博客

  • 快速幂与矩阵快速幂

    首先介绍快速幂:

    快速幂可以在常数时间内计算$$a^b$$,而朴素的算法则需要$$O(b)$$

    思想:将b转换为二进制,则

    $$b=1+2+4+8+……$$

    $$a^b=a^{1}*a^{2}*a^{4}*….=a^{1}*(a^{1})^{2}*((a^{1})^{2})^{2}*…$$

    例如:

    求$$3^5$$

    $$\because 5=1+4$$

    $$\therefore 3^5=3^{1+4}=3^1*3^4$$

    先计算$$3^1=3$$,结果乘上3 然后计算$$(3^1)^2=3^2=9$$,然后计算$$(3^2)^2=3^4=81$$,结果乘上81,最终结果为$$3*81=243$$

    代码:

     

    #include 
    #include 
    #include 
    #include 
    #include 
    #include 
    #include 
    using namespace std;
    using int_t = long long int;
    
    int_t fastPower(int_t base, int_t index) {
        //结果
        int_t result = 1;
        //只要指数不为0就一直循环
        while (index) {
            //index的最低位为1,结果中应该乘上base的对应次幂
            if (index & 1) {
                result = result*base;
            }
            base = base*base;
            //index右移一位
            index >>= 1;
        }
        return result;
    }
    
    int main() {
        int_t base, index;
        while (true) {
            cin >> base>>index;
            if (base == 0 && index == 0) break;
            cout << fastPower(base, index) << endl;
    
        }
        return 0;
    }
    

    矩阵快速幂:

    使用快速幂的思想来计算矩阵的幂

     

    #include 
    #include 
    #include 
    #include 
    #include 
    #include 
    #include 
    #include 
    #include 
    using namespace std;
    using int_t = unsigned long long int;
    const int_t mod = (int_t) 1e9 + 7;
    //Matrix的前置声明,用于下面两个函数
    template 
    class Matrix;
    //友元函数的声明
    //模板类的友元函数必须显式的写明模板
    template 
    ostream& operator<<(ostream & os, const Matrix & mat);
    template 
    istream& operator>>(istream & os, const Matrix & mat);
    
    template 
    class Matrix {
    private:
        ValType* data;
        int_t row;
        int_t col;
        int_t size;
        int_t num;
    public:
        //析构函数,释放内存
        ~Matrix() {
            delete[] data;
        }
        //拷贝构造函数
        Matrix(const Matrix& other) :
        row(other.row), col(other.col), size(other.size), num(other.num) {
            data = new ValType[other.num];
            memcpy(data, other.data, other.size);
        }
        //默认构造函数,动态分配内存
        Matrix(int_t r, int_t c) :
        row(r), col(c) {
            this->row = r;
            this->col = c;
            num = (row + 1)*(col + 1);
            data = new ValType[num];
            size = sizeof (ValType) * num;
            memset(data, 0, size);
        }
        //访问矩阵的某个元素
        ValType& at(int_t r, int_t c) {
            return data[c + r * col];
        }
        //重载的*运算符,矩阵无法相乘时会抛出异常
        Matrix operator*(Matrix& mat) throw (const char*) {
            if (this->col != mat.row) throw "Column number of matrix1 doesn't equals that in matrix2";
            Matrix result(this->row, mat.col);
            for (int_t i = 1; i <= this->row; i++) {
                for (int_t j = 1; j <= mat.col; j++) {
                    int_t sum = 0;
                    for (int_t x = 1; x <= this->col; x++) {
                        sum = (sum % mod + at(i, x) % mod * mat.at(x, j) % mod) % mod;
                    }
                    result.at(i, j) += sum % mod;
                }
            }
            return result;
        }
    
        int_t getCol() const {
            return col;
        }
    
        int_t getRow() const {
            return row;
        }
        //赋值构造函数
        Matrix & operator=(const Matrix & other) {
            this->col = other.col;
            this->num = other.num;
            this->row = other.row;
            this->size = other.size;
            this->data = new ValType[num];
            memcpy(data, other.data, size);
        }
        //友元,函数名后必须写<>以告诉编译器这个友元使用了模板,如果不能从参数中推断模板类型那么还必须在<>内写上模板类型
        friend ostream& operator<< <>(ostream & os, const Matrix & mat);
        friend istream& operator>> <>(istream & os, const Matrix & mat);
    };
    //友元函数的实现
    template 
    ostream & operator<<(ostream & os, Matrix & mat) {
        for (int_t i = 1; i <= mat.getRow(); i++) {
            for (int_t j = 1; j <= mat.getCol(); j++) {
                os << mat.at(i, j) % mod << " ";
            }
            os << endl;
        }
        return os;
    }
    
    template 
    istream & operator>>(istream & is, Matrix & mat) {
        for (int_t i = 1; i <= mat.getRow(); i++) {
            for (int_t j = 1; j <= mat.getCol(); j++) {
                is >> mat.at(i, j);
                mat.at(i, j) %= mod;
            }
        }
        return is;
    }
    using MatType = Matrix;
    
    int main() {
    
    
        int_t n, index;
        cin >> n>>index;
        MatType mat(n, n);
        cin>>mat;
        MatType result = mat;
        //普通的快速幂
        for (index--; index; index>>=1,mat = mat * mat) {
            if (index & 1) result = result * mat;
        }
        cout << result;
        return 0;
    }
    

    此外,矩阵$$A=\begin{bmatrix}
    0&1 \\
    1&1
    \end{bmatrix}$$的n次方所得的矩阵的对角线上的元素即为斐波那契数列的第n项。

  • 乘法逆元

    如果$$ax\equiv 1(mod n)$$则称x为a模n意义下的逆

    则有$$ax-ny=1$$,然后用扩展欧几里得算法求解即可。

  • 唯一分解定理

    任何一个正整数都可以写成几个质数的幂次方之积,即
    $$n={p_1}^{k_1} {p_2}^{k_2} {p_3}^{k_3} {p_4}^{k_4} …..{p_m}^{k_m} $$ 其中$$p_m$$为素数,$$k_m$$为非负整数

    例如

    $$30=2^1*3^1*5^1$$

    $$40=5*8=5^{1}*2^{3}$$

    $$72=2*36=2^1*6^2$$

    代码实现:

     

    #include 
    #include 
    #include 
    #include 
    #include 
    #include 
    #include 
    #include 
    using namespace std;
    using int_t = long long int;
    
    map func(int_t x) {
        //map
        //x=p1^k1*p2^k2...........pn^kn
        map result;
        for (int_t i = 2; i <= x; i++) {
            //如果x的可以整除i,则结果中i的指数+1
            while (x % i == 0) {
                result[i]++;
                x /= i;
            }
        }
        return result;
    }
    
    int main() {
        int_t x;
        while (cin >> x) {
            auto result = func(x);
            cout << x << " = ";
            //请手动忽略最后多余的乘号
            for (auto p : result) {
                cout << p.first << "^" << p.second << " * ";
            }
            cout << endl;
        }
    }
    
  • 欧拉函数

    设函数$$\phi (x)$$为不超过x且与x互素的正整数的个数(两个数互素,则两个数的最大公因数为1)

    则$$\phi(x)=x(1-\frac {1} {p_1})(1-\frac {1} {p_2})(1-\frac {1} {p_3})(1-\frac {1} {p_4})·····(1-\frac {1} {p_n})$$

    其中$$p_n$$为将x通过唯一分解定理所得的式子的底数,即$$n={p_1}^{k_1} {p_2}^{k_2} {p_3}^{k_3} {p_4}^{k_4} …..{p_m}^{k_m} $$ 其中$$p_n$$为素数

     

    #include 
    #include 
    #include 
    #include 
    #include 
    #include 
    #include 
    #include 
    using namespace std;
    using int_t = long long int;
    
    int_t func(int_t x) {
        int_t result = x;
        for (int_t i = 2; i <= x; i++) {
            //如果x的可以整除i,则结果中i的指数+1
            if (x % i == 0) {
                //    result = result * (i - 1) / i;
                //先乘后除防止溢出
                result = result / i * (i - 1);
                while (x % i == 0) x /= i;
            }
        }
        if (x > 1) result = result / x * (x - 1);
        return result;
    }
    
    int main() {
        int_t x;
        while (cin >> x) {
            int_t result = func(x);
            cout<
    

    或者可以通过类似于线性筛素数的方式求出一定范围内所有的数的欧拉函数值

    #include 
    #include 
    #include 
    #include 
    #include 
    #include 
    #include 
    #include 
    
    using namespace std;
    using int_t = long long int;
    const int_t MAX = 100000;
    int_t phi[MAX + 1];
    int main() {
        int_t x;
        cin>>x;
        memset(phi, 0, sizeof (phi));
        phi[1] = 1;
        for (int_t i = 2; i <= x; i++) {
            if (phi[i] == 0) {
                for (int_t j = i; j <= x; j = j + i) {
                    if (phi[j] == 0) phi[j] = j;
                    phi[j] = phi[j] / i * (i - 1);
                }
            }
    
        }
        for(int_t i=1;i<=x;i++) cout<<"phi("<<i<<") = "<<phi[i]<<endl;
    
    }
    
  • 线性筛素数

    使用欧拉筛法可以接近线性的时间内筛出一定范围内所有素数。

    基本思想:如果$$a$$是素数,那么$$a*1,a*2,a*3,a*4,a*5…..$$都不是素数,这样就可以对一定范围内的非素数进行标记。

    代码:

     

    #include 
    #include 
    #include 
    #include 
    #include 
    #include 
    #include 
    using namespace std;
    using int_t = long long int;
    const int_t MAX = 10000;
    char isPrime[MAX + 1];
    vector primeTable;
    
    int main() {
        //初始时,假设所有数都是素数
        memset(isPrime, 1, sizeof (isPrime));
        int_t n;
        cin>>n;
        //1当然不是素数
        isPrime[1] = 0;
        //i只需要枚举到sqrt(n),因为i>sqrt(n)的情况已经被考虑过了
        for (int_t i = 2; i <= sqrt(n + 0.5); i++) {
            //当外层循环执行到i时,如果i仍未被标记非素数,就说明i真的是素数
            //如果i不是素数,那样就不需要继续标记i*2,i*3,i*4....了 因为这些情况已经被标记了
            if (isPrime[i]) {
                for (int_t j = i * i; j <= n; j =j+i) {
                    isPrime[j] = false;
                }
            }
        }
        for (int_t i = 1; i <= n; i++) {
            if (isPrime[i]) {
                primeTable.push_back(i);
            }
        }
        for (int_t prime : primeTable) {
            cout << prime << endl;
        }
    
        return 0;
    }
    
  • 求最长上升子序列的O(n log n)算法

    上升子序列:从一个序列中按照顺序挑出一些数,使得后一个数大于前一个数

    例如对于序列$$1,5,3,7,7,2$$

    $$1,5,7$$ $$5,7$$ $$1,3$$都是其上升子序列

    但是$$1,5,7,7$$不是上升子序列(不满足后一个数大于前一个数)

    最长上升子序列:一个序列中最长的上升子序列

    对于序列$$1,5,3,7,7,2$$ 其最长上升子序列是$$1,3,7$$或者$$1,5,7$$

     

    求最长上升子序列的算法

    1.朴素的O(n²)的动态规划

    设dp[tail]为序列中以第tail个元素结尾的最长上升子序列的长度

    则有dp[1]=1 (第一个数自成一个最长上升子序列)

    tail>1时$$dp[tail]=max\{0,dp[1],dp[2],dp[3]…dp[i]…dp[tail-1]|满足seq[i]<seq[tail]\}+1$$

    结果为$$max{dp[i]|1<=i<=序列长度}$$

    解释:求dp[tail]时,从序列中tail前的所有数中找到比seq[tail]小的数,在满足这个数比seq[tail]小的情况下使dp[tail]最大

    举例:对于序列

    $$1,5,3,2,5,6$$

    使用seq[1]表示序列中的第一项,seq[2]表示序列中的第二项,以此类推

    初始时dp[1]=1;

    然后开始求dp[2] 发现seq[1]<seq[2] 所以dp[2]=dp[1]+1=2

    开始求dp[3] 发现seq[1]<seq[3] 更新dp[3]=dp[1]+1=2

    开始求dp[4] 发现seq[1]<seq[4] 更新dp[4]=dp[1]+1=2;

    开始求dp[5] 发现seq[1]<seq[5] 更新dp[5]=dp[1]+1=2

    又发现seq[3]<seq[5] 再次更新dp[5]=dp[3]+1=3

    开始求dp[6] 发现seq[1]<seq[6] 更新dp[6]=dp[1]+1=2

    又发现seq[2]<seq[6] 更新dp[6]=dp[2]+1=3

    又发现seq[3]<seq[6] 更新dp[6]=dp[3]+1=3

    以此类推

     

    最终 dp[1]=1 dp[2]=2 dp[3]=2 dp[4]=2 dp[5]=3 dp[6]=3

     

    最终结果为3

     

    程序:

    #include 
    
    using namespace std;
    typedef long long int_t;
    int_t seq[1001];
    int_t dp[1001];
    
    int main() {
        dp[1] = 1;
        int_t n;
        cin>>n;
        for (int_t i = 1; i <= n; i++) cin >> seq[i];
        //一次性计算出LIS
        for (int_t i = 1; i <= n; i++) {
            int_t length = 0;
            for (int_t j = 1; j <= i; j++) { if (seq[i] > seq[j]) {
                    if (dp[j] > length) length = dp[j];
                }
            }
            dp[i] = length + 1;
    
        }
        int_t result = 0;
        for (int_t i = 1; i <= n; i++) {
            result = max(result, dp[i]);
        }
        cout << result;
        return 0;
    }
    

    2.复杂度为O(n log n)的算法

    设dp[i]为以序列中第i个元素结尾的最长上升子序列的长度,那么如果对于序列中任意两个数seq[i]和seq[j] 如果满足dp[i]==dp[j]且seq[i]<seq[j] 那么仅保存i一定不会丢失最优解(因为seq[i]<seq[j],所以可以接到seq[j]后面的数一定可以接到seq[i]后面,然而可以接到seq[i]后面的数却不一定能接到seq[j]后面)

    设g[x]为满足dp[i]==x的最小的seq[i]

    下一句的最长上升子序列长度是指以该数结尾的最长上升子序列的长度

    例:g[3]为满足最长上升子序列长度为3的序列中最小的数

    对于以下序列$$1,5,3,6,7$$

    $$

    g[1]=1

    g[2]=3 (不能是5,因为3比5小)

    g[3]=6

    g[4]=7

    g[5]=\infty (因为不存在长度的5的最长上升子序列)

    $$

    同时一定满足

    $$g[1]<=g[2]<=g[3]<=g[4]……..<=g[n]$$

    所以为了确定序列中的一个数seq[i]的最长上升子序列长度,只需要在g数组中进行二分查找即可。查找后,查找结果即为以seq[i]结尾的最长上升子序列的长度

    但要注意:

    假设在g中二分查找seq[i]所得到的结果是index,那么还需要把g[index]改为seq[i],以此维护g数组的单调性(因为seq[i]一定比g[index]小,如果seq[i]比g[index]大,那么二分查找的结果就不是index了)

    代码:注意,初始时g数组需要全部设置为$$\infty$$

     

     

    #include 
    #include 
    #include 
    using namespace std;
    using int_t = long long int;
    int_t seq[100000 + 1];
    int_t g[100000 + 1];
    int_t result[100000 + 1];
    
    int main() {
        int_t n;
        cin>>n;
        for (int_t i = 1; i <= n; i++) cin >> seq[i];
        memset(g, 0x7f, sizeof (g));
    
        for (int_t i = 1; i <= n; i++) {
            int_t index = lower_bound(g + 1, g + 1 + n, seq[i]) - g;
            result[i] = index;
            g[index] = seq[i];
        }
        int_t r = 0;
        for (int_t i = 1; i <= n; i++) r = max(r, result[i]);
        cout << r;
        return 0;
    }
    

    因为二分查找的复杂度为$$O(log n)$$,所以该算法的复杂度为$$O(n log n)$$

     

    两种算法的效率对比:

    数据量为10000时:

     
    每组数据第一行为算法2,第二行为算法1

    1:
            0
            0.672
    2:
            0
            0.656
    3:
            0.016
            0.641
    4:
            0.015
            0.656
    5:
            0.016
            0.641
    6:
            0.015
            0.719
    7:
            0
            0.656
    8:
            0.016
            0.609
    9:
            0
            0.687
    10:
            0.016
            0.672
    
    

    数据量为100000时

    1:
            0.094
            65.765
    2:
            0.093
            63.454
    3:
            0.078
            65.172
    4:
            0.093
            65.391
    5:
            0.079
            61.999
    6:
            0.094
            63.703
    7:
            0.094
            65.047
    8:
            0.094
            65.124
    9:
            0.109
            64.954
    10:
            0.062
            63.766
    
    可以看到,n log n的算法明显要比n^2的算法快
    
  • 欧几里得算法与扩展欧几里得算法

    欧几里得算法,即人教版高中数学必修三上的辗转相除法

    用于求两个整数的最大公因数(即gcd(a,b))
    注:a b的最小公倍数等于$$\frac {ab}{\gcd \left(a,b\right)}$$

    int gcd(int a, int b) {
        //边界b==0 此时a就是结果
        if (b == 0) {
            return a;
        } else {
            return gcd(b, a % b);
        }
    }
    

    扩展欧几里得算法:在求最大公约数的同时求出方程$$ax+by=gcd(a,b)$$的一组整数解,使|x|+|y|最小

    //a b 方程中的a b 
    //gcd 用来存放gcd(a,b)
    //x y 用来存放解
    void exgcd(int a, int b, int & gcd, int& x, int& y) {
        //边界b==0 此时a就是结果
        if (b == 0) {
            x = 1;
            y = 0;
            gcd = a;
        } else {
            exgcd(b, a % b, gcd, y, x);
            y -= (a / b) * x;
    
        }
    }
    

    求出$$ax+by=gcd(a,b)$$的解之后,如果$$c$$可以整除$$gcd(a,b)$$ 就说明方程$$ax+by=c$$存在整数解$$x_0,y_0$$,而且这组解为$$x_1=x_0*c/gcd(a,b),y_1=y_0*c/gcd(a,b)$$

    然后,如果$$x_1,y_1$$是方程$$ax+by=c$$的解,则$$x_1+k*a/gcd(a,b),y_1-k*b/gcd(a,b)$$都是方程$$ax+by=c$$的解

  • 省选准备计划

    1.网络流 Dinic ISAP 二分图等

    2.各种平衡树 Splay Treap 等

    3.博弈论

    4.计算几何

    5.字符串

    6.可持久化数据结构

    7.数论