基础数学算法

这里介绍的是，算法题目中的一些，基础的考察数学问题的算法，主要包括gcd，素数问题，快速幂这三部分。

一、gcd和extgcd：
首先是gcd问题，求解的思路是辗转相除法。下面看书上的例题：给定平面上的两个格点P1=(x1,y1)和P2=(x2,y2)，求线段P1P2上，除P1和P2以外一共有几个格点？
$-10^9 <= x1,x2,y1,y2 <= 10^9$
数据范围很大，遍历横坐标，计算纵坐标的暴力方法肯定是不能用的。这里先说一下gcd的做法。比如给出的例子是(9,1)和(1,13)这两个点，明显斜率是 $-\frac{12}{8}$ ，也就是 $-\frac{3}{2}$ 。此时有一个很重要的信息，x轴变化两个单位，y轴会随着变化3个单位。所以这个时候有一个很明显的思路，计算x轴从1出发，经过了几个2会到9。答案是（9-1）/ 2 = 4个，去掉端点的一个，中间有3个整数点（一定是整数点才是格点）。
回过头来想一想，这个题究竟是怎么做的。1. 求得了斜率，并且化简。2. 得出直线在x轴变化时，y轴点的变化规律，核心是确定整数变化规律。这个思路是没问题的，所以核心就是找到这个化简后的分母或者分子，由此将问题变成了找最大公约数(gcd)。这道题的意义在于，让我们发现了gcd的又一个作用——化简到最小整数比。然后就可以找到点在直线上运动时，x轴和y轴都为整数的点。下面贴出gcd的代码模板：

int gcd(int a,int b){
    if(b == 0)return a;
    return gcd(b, a%b);
}

接下来是扩展欧几里得算法（extgcd)，也就是扩展的辗转相除法。这个算法解决的主要是线性同余方程那里的问题，即： $ax+by=gcd(a,b)$ 先说明一下extgcd这个函数，它的输入是(a,b,x,y)，x和y只有是一个开了空间的变量就好，值不用初始化，返回值是gcd(a,b)，并且x和y已经更新为所求的值(x和y有可能取负）。这里的求解方法不再赘述，原理的话我也没有看太多，我觉得这里关键还是要了解怎么用，之后在数学问题中的取模部分还会对这里进行分析，现在先直接看代码：

int extgcd(int a,int b,int& x,int& y){
    int d = a;
    if(b != 0){
        d = extgcd(b,a%b,y,x);
        y -= (a/b) * x;
    }
    else{
        x=1;y=0;
    }
    return d;
}

关于辗转相除法的时间复杂度，没有一个准确的估计。不过可以粗略确定的是，这个算法的时间复杂度在 $O(logmax(a,b))$ 以内。

二、素数问题
素数问题包含两类，也可以说是做题时经常会有的两类需求。第一类是判断一个数是否是素数，当判断的次数较少时，会这样判断；第二类则是直接筛出一个区间内所有的素数，需要判断的次数较多时，会用筛法来判断。这两类问题其实很常见，前者往往速度快一些，但是只能处理一个数；后者往往慢一点，但是可以批量处理出一个区间的数。
首先是判断一个数n是否为素数，这个其实很简单，最直接的思路就是遍历一下[2,n-1]，看是否有一个数m满足m|n(n % m == 0)。但是其实[2,n-1]的区间有点大，下一步最直接的想法是[2,n/2]，但是其实还，可以再进一步优化到[2, $\sqrt{n}$ ]。比如说8=2*4，如果判断到4，得到8%4 = 0，那么在此之前，一定已经根据8 % 2 = 0得到了结果。所以也就是只需要判断约数中小的那一部分即可，判断到 $\sqrt{n}$ 刚好到临界。下面是素数判定的代码：

bool is_prime(int n){
    for(int i = 2;i*i <= n;i++){
        if(n % i == 0)return false;
    }
    return true;
}

接下来是素数筛法，一般的题目用埃式筛法已经够了，这个筛法的时间复杂度是O(n*lg(lgn))，还有一种欧拉筛法，复杂度是O(n)，但是其实lg(lgn)的复杂度，也和常数的复杂度没太大区别了。n就算有 $10^9$ 这么大（再大O(n)都会卡）,lg(lgn)也就是lg(30)大约是5，代价并不大。下面是埃式筛法的代码，我相信只看代码足以看懂思路：

int prime[max_n];
bool is_prime[max_n+1];

int sieve(int n){
    int p = 0;
    for(int i = 0;i <= n;i++)is_prime = true;
    is_prime[0] = is_prime[1] = false;
    for(int i = 2;i <= n;i++){
        if(is_prime[i]){
            prime[p] = i;
            p++;                //这里我习惯把赋值和更新分开
            for(int j = 2*i;j <= n;j += i)is_prime[j] = false;
        }
    }
    return p;                  //返回素数的个数
}

这里还有几个扩展，首先还是筛选素数的问题。如果要询问一段区间[a,b)内素数的个数，a<b<=10¹²，b-a <= $10^6$ 。此时就不能先筛一遍a，再筛一遍b这样做了，因为筛10¹²一定会超时。
这个时候用的方法其实也简单，首先可以猜到，筛的区间一定不再是到a或b，而是直接筛[a,b)这个区间。这个时候其实最关键的问题是，拿谁来筛？这里的方法其实是借助了素数判定时的思路，判断b是不是约数，在[2, $\sqrt{b}$ ]来找约数就好了。这里也是一样，用[2, $\sqrt{b}$ ]来筛，超过这个区间的数能筛的情况，这个区间的数也能筛。有了这个思路，看下面的代码：

typedef long long ll;

bool is_prime[max_l];
bool is_prime_small[max_sqrt_b];

void segment_sieve(ll a,ll b){
    for(int i = 0;(ll)i*i < b;i++)is_prime_small[i] = true; //一定要转换类型，不然会wa的很惨
    for(int i = 0;i < b-a;i++)is_prime[i] = true;
    
    for(int i = 2;(ll)i*i < b;i++){
        if(is_prime_small[i]){
            for(int j = 2*i;(ll)j*j < b;j+=i)is_prime_small[j] = false; // 筛出小区间的素数
            /* 这里其实难的是看j的起点，(a+i-1)/i*i可以保证，若a % i == 0,取a,
            否则是[a,b)的首个能被i整除的数。max是为了可以确保不是从i开始的，因为0*i或者
            1*i开始都是没意义的，最少要2*i开始*/
            for(ll j = max(2LL, (a+i-1)/i)*i,j < b;j += i)is_prime[j-a] = false;
        }
    }
}

这里的代码是真的很厉害，首先2LL给出了常数用max和long long类型的数做比较时的方法，用LL转化为long long类型；然后就是后面(a+i-1)*i/I的设计，得到[a,b)的首个可以被i整除的数，真的很厉害。
另外补充两份代码，一个是枚举n所有的约数，用前面讲的[2, $\sqrt{n}$ ]的想法可以在 $\sqrt{n}$ 的时间完成，思路也不复杂；另一个是将整数分解为质因数的乘积，这里有一步很出彩，值得学习，两份代码如下：

vector<int> divisor(int n){
    vector<int> res;
    for(int i = 1;i*i <= n;i++){
        if(n%i == 0){
            res.push_back(i);
            if(i != n/i)res.push_back(n/i);
        }
    }
    return res;
}

map<int,int> prime_factor(int n){
    map<int,int> res;
    for(int i = 2;i*i <= n;i++){
        while(n % i == 0){
            res[i]++;
            /*这一步很厉害，试想如果是找出所有质因数，那么找到一个质因数之后，
            直接通过这一步将该质因数除干净了，再去找下一个质因数即可，比如说
            处理12，i = 2时，直接12/2=6，6/2=3，然后再拿3去找下一个质因数即可，
            不需要再去管质因数有几个这样的问题。*/
            n /= i;
        }
    }
    return res;
}

三、快速幂
快速幂的功能其实也比较简单，就是输入(n,m)求出 $n^m$ 的值。如果是普通的算法，连乘m次n，需要O(m)的时间复杂度。还有一种想法是倍增，从res=n开始，res=res * res这样去处理，可以降到O(lgn)，但是要解决找过了头的问题。比如m=8，这种情况好说，三次res = res * res就好了。但是如果m=9，三次res = res * res之后，此时是8次方，还差一次方，要么是再来一次res = res * res，然后再除回来；要么是直接乘上去，但是其实中间的过程都挺复杂的。快速幂的出现，解决了这个问题。
快速幂降低时间复杂度的想法其实和倍增类似，都是用现有的res来平方，从而达到倍增的目的。但是快速幂有一点另辟蹊径的地方在于，运用位运算。比如说要求出n¹⁰，10用二进制表示是1010，一开始的res=n，代表的就是10变成二进制以后的第一位，也就是0；然后res = res * res，此时是 $n^2$ ，其实就对应第二位，是1，然后再来一个res = res * res，是第三位0，以此类推下去，res=res*res其实可以表示幂次m转变为二进制的每一位。接下来，如果该位取1，那么就代表这一位需要留下，为0说明结果不需要这一位。理解到这里，就可以看下面的代码了，代码包含了取余的部分，因为快速幂的结果一般很大，留的是取余后的结果：

typedef long long ll;

ll mod_pow(ll x,ll n,ll mod){
    ll res = 1;
    while(n > 0){
        if(n & 1)res = res*x%mod;
        x = x*x % mod;
        n = n >> 1;
    }
    return res;
}