数论总结——为搞事而生

目录

素数的判定

一个数 n 如果是合数，那么它的所有的因子不超过sqrt(n)，复杂度是o(n*sqrt(n))：

bool prime(int x){//判断x是不是质数，是返回true，不是返回false     if(x <= 1) return false;     for(int i = 2; i <= sqrt(x + 0.5); i ++){//0.5是防止根号的精度误差         if(x % i == 0) return false;    }    return true;}

或者是这样写也行：

bool prime(int x){//判断x是不是质数，是返回true，不是返回false     if(x <= 1) return false;     for(int i = 2; i * i <= x; i ++){//用乘法避免根号的精度误差         if(x % i == 0) return false;    }    return true;}//根据题目不同，如果i*i会爆int，记得开longlong

引入素数表来进行维护：

bool IsPrime(unsigned n){    if ( n < 2 )    { // 小于2的数即不是合数也不是素数        throw 0;    }    static unsigned aPrimeList[] = { // 素数表        1, 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41,        43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 113,         193, 241, 257, 337, 353, 401, 433, 449, 577, 593, 641,         673, 769, 881, 929, 977, 1009, 1153, 1201, 1217, 1249,         1297,1361, 1409, 1489, 1553, 1601, 1697, 1777, 1873,         1889, 2017, 2081, 2113, 2129, 2161, 2273, 2417, 2593,         2609, 2657, 2689, 2753, 2801, 2833, 2897, 3041, 3089,         3121, 3137, 3169, 3217, 3313, 3329, 3361, 3457, 3617,         3697, 3761, 3793, 3889, 4001, 4049, 4129, 4177, 4241,         4273, 4289, 4337, 4481, 4513, 4561, 4657, 4673, 4721,         4801, 4817, 4993, 5009, 5153, 5233, 5281, 5297, 5393,         5441, 5521, 5569, 5857, 5953, 6113, 6257, 6337, 6353,         6449, 6481, 6529, 6577, 6673, 6689, 6737, 6833, 6961,         6977, 7057, 7121, 7297, 7393, 7457, 7489, 7537, 7649,         7681, 7793, 7841, 7873, 7937, 8017, 8081, 8161, 8209,         8273, 8353, 8369, 8513, 8609, 8641, 8689, 8737, 8753,         8849, 8929, 9041, 9137, 9281, 9377, 9473, 9521, 9601,         9649, 9697, 9857     };    const int nListNum = sizeof(aPrimeList)/sizeof(unsigned);//计算素数表里元素的个数    for (unsigned i=2;i<nListNum;++i )    {         if(n/2+1<aPrimeList[i])        {            return true;        }        if(0==n%aPrimeList[i])        {            return false;        }    }    for (unsigned i=aPrimeList[nListNum-1];i<n/2+1;i++ )    {         if (0==n%i)        { // 除尽了，合数             return false;        }    }    return true; } 

二次探测定理：如果p是素数，且0<x<p，则方程x^2=1(mod p)的解为1或p-1。
最后一种情况，面对大数时采用Miller-Rabbin素数测试法：
Miller-Rabbin素数测试是一个不确定的算法，只能从概率意义上判定一个数可能是素数，但并不能确保。算法流程如下:
1.选择T个随机数A，并且有A

int Miller_Rabbin(long long n)    {        long long i,s,a;        s=10;    //s的值可以根据需要变大     // randomize();        for(i=0;i <s;i++)        {            a=long long(rand()%(n-1)+2); //自动生成受限            if(modular_exp(a,n-1,n)>1)                return 0;        }        return 1;    }    long long modular_exp(long long a,long long b,long long c)//求a^b%c该函数受限    {        if(a==0)            return 0;        if(b==0)            return 1;        if(b==1)            return a%c;        return (a*modular_exp(a,b-1,c))%c;    } 

完整版代码：

#include <stdio.h>#include <stdlib.h>#include <string.h>#include <algorithm>#include <iostream>#include <math.h>using namespace std;const int Times = 10;typedef long long LL;LL multi(LL a, LL b, LL m){    LL ans = 0;    a %= m;    while(b)    {        if(b & 1)        {            ans = (ans + a) % m;            b--;        }        b >>= 1;        a = (a + a) % m;    }    return ans;}LL quick_mod(LL a, LL b, LL m){    LL ans = 1;    a %= m;    while(b)    {        if(b & 1)        {            ans = multi(ans, a, m);            b--;        }        b >>= 1;        a = multi(a, a, m);    }    return ans;}bool Miller_Rabin(LL n){    if(n == 2) return true;    if(n < 2 || !(n & 1)) return false;    LL m = n - 1;    int k = 0;    while((m & 1) == 0)    {        k++;        m >>= 1;    }    for(int i=0; i<Times; i++)    {        LL a = rand() % (n - 1) + 1;        LL x = quick_mod(a, m, n);        LL y = 0;        for(int j=0; j<k; j++)        {            y = multi(x, x, n);            if(y == 1 && x != 1 && x != n - 1) return false;            x = y;        }        if(y != 1) return false;    }    return true;}int main(){    int T;    scanf("%d",&T);    while(T--)    {        LL n;        scanf("%I64d",&n);        if(Miller_Rabin(n)) puts("Yes");        else puts("No");    }    return 0;}

质数检测Java大数版：

import java.io.*;import java.util.*;import java.math.BigInteger;public class Main{        public static final int Times = 10;        public static BigInteger quick_mod(BigInteger a,BigInteger b,BigInteger m){                BigInteger ans = BigInteger.ONE;                a = a.mod(m);                while(!(b.equals(BigInteger.ZERO))){                        if((b.mod(BigInteger.valueOf(2))).equals(BigInteger.ONE)){                                ans = (ans.multiply(a)).mod(m);                                b = b.subtract(BigInteger.ONE);                        }                        b = b.divide(BigInteger.valueOf(2));                        a = (a.multiply(a)).mod(m);                }                return ans;        }        public static boolean Miller_Rabin(BigInteger n){                if(n.equals(BigInteger.valueOf(2))) return true;                if(n.equals(BigInteger.ONE)) return false;                if((n.mod(BigInteger.valueOf(2))).equals(BigInteger.ZERO)) return false;                BigInteger m = n.subtract(BigInteger.ONE);                BigInteger y = BigInteger.ZERO;                int k = 0;                while((m.mod(BigInteger.valueOf(2))).equals(BigInteger.ZERO)){                        k++;                        m = m.divide(BigInteger.valueOf(2));                }                Random d = new Random();                for(int i=0;i<Times;i++){                        int t = 0;                        if(n.compareTo(BigInteger.valueOf(10000)) == 1){                                t = 10000;                        }else{                                t = n.intValue() - 1;                        }                        int a = d.nextInt(t) + 1;                        BigInteger x = quick_mod(BigInteger.valueOf(a),m,n);                        for(int j=0;j<k;j++){                                y = (x.multiply(x)).mod(n);                                if(y.equals(BigInteger.ONE) && !(x.equals(BigInteger.ONE)) && !(x.equals(n.subtract(BigInteger.ONE)))) return false;                                x = y;                        }                        if(!(y.equals(BigInteger.ONE))) return false;                }                return true;        }        public static void main(String[] args){                Scanner cin = new Scanner(System.in);                while(cin.hasNextBigInteger()){                        BigInteger n = cin.nextBigInteger();                        if(Miller_Rabin(n)) System.out.println("Yes");                        else System.out.println("No");                }        }}

素数的筛选（构造素数表）

不采用筛素：假设我们已经我素数序列: p1, p2, .. pn
现在要判断pn+1是否是素数, 则需要(1, sqrt(pn+1)]范围内的所有素数序列,
而这个素数序列显然已经作为p1, p2, .. pn的一个子集被包含了

void makePrimes(int primes[], int num){    int i, j, cnt;    primes[0] = 2;    primes[1] = 3;    for(i = 5, cnt = 2; cnt < num; i += 2)    {        int flag = true;        for(j = 1; primes[j]*primes[j] <= i; ++j)        {            if(i%primes[j] == 0)                {                    flag = false; break;                }        }        if(flag) primes[cnt++] = i;    }}

可以使用筛素法，复杂度O(nloglogn)：

bool prime[N];void init(){    for(int i = 2; i < N; i ++) prime[i] = true;    for(int i = 2; i*i < N; i ++){//判断改成i*i<N         if(prime[i]){            for(int j = i*i; j < N; j += i){//从i*i开始就可以了                 prime[j] = false;              }        }    }}

或者采用线性筛素法O(n)：

#include<cstdio>const int N = 100000 + 5;bool prime[N];//prime[i]表示i是不是质数 int p[N], tot;//p[N]用来存质数 void init(){    for(int i = 2; i < N; i ++) prime[i] = true;//初始化为质数     for(int i = 2; i < N; i++){        if(prime[i]) p[tot ++] = i;//把质数存起来         for(int j = 0; j < tot && i * p[j] < N; j++){            prime[i * p[j]] = false;            if(i % p[j] == 0) break;//保证每个合数被它最小的质因数筛去         }    }    }int main(){    init();}

Meisell-Lehmer 模板，求1-n范围内的素数个数：

#include<bits/stdc++.h>using namespace std;typedef long long LL;const int N = 5e6 + 2;bool np[N];int prime[N], pi[N];int getprime() {    int cnt = 0;    np[0] = np[1] = true;    pi[0] = pi[1] = 0;    for(int i = 2; i < N; ++i) {        if(!np[i]) prime[++cnt] = i;        pi[i] = cnt;        for(int j = 1; j <= cnt && i * prime[j] < N; ++j) {            np[i * prime[j]] = true;            if(i % prime[j] == 0)   break;        }    }    return cnt;}const int M = 7;const int PM = 2 * 3 * 5 * 7 * 11 * 13 * 17;int phi[PM + 1][M + 1], sz[M + 1];void init() {    getprime();    sz[0] = 1;    for(int i = 0; i <= PM; ++i)  phi[i][0] = i;    for(int i = 1; i <= M; ++i) {        sz[i] = prime[i] * sz[i - 1];        for(int j = 1; j <= PM; ++j) {            phi[j][i] = phi[j][i - 1] - phi[j / prime[i]][i - 1];        }    }}int sqrt2(LL x) {    LL r = (LL)sqrt(x - 0.1);    while(r * r <= x)   ++r;    return int(r - 1);}int sqrt3(LL x) {    LL r = (LL)cbrt(x - 0.1);    while(r * r * r <= x)   ++r;    return int(r - 1);}LL getphi(LL x, int s) {    if(s == 0)  return x;    if(s <= M)  return phi[x % sz[s]][s] + (x / sz[s]) * phi[sz[s]][s];    if(x <= prime[s]*prime[s])   return pi[x] - s + 1;    if(x <= prime[s]*prime[s]*prime[s] && x < N) {        int s2x = pi[sqrt2(x)];        LL ans = pi[x] - (s2x + s - 2) * (s2x - s + 1) / 2;        for(int i = s + 1; i <= s2x; ++i) {            ans += pi[x / prime[i]];        }        return ans;    }    return getphi(x, s - 1) - getphi(x / prime[s], s - 1);}LL getpi(LL x) {    if(x < N)   return pi[x];    LL ans = getphi(x, pi[sqrt3(x)]) + pi[sqrt3(x)] - 1;    for(int i = pi[sqrt3(x)] + 1, ed = pi[sqrt2(x)]; i <= ed; ++i) {        ans -= getpi(x / prime[i]) - i + 1;    }    return ans;}LL lehmer_pi(LL x) {    if(x < N)   return pi[x];    int a = (int)lehmer_pi(sqrt2(sqrt2(x)));    int b = (int)lehmer_pi(sqrt2(x));    int c = (int)lehmer_pi(sqrt3(x));    LL sum = getphi(x, a) + LL(b + a - 2) * (b - a + 1) / 2;    for (int i = a + 1; i <= b; i++) {        LL w = x / prime[i];        sum -= lehmer_pi(w);        if (i > c) continue;        LL lim = lehmer_pi(sqrt2(w));        for (int j = i; j <= lim; j++) {            sum -= lehmer_pi(w / prime[j]) - (j - 1);        }    }    return sum;}int main() {    init();    LL n;    while(cin >> n) {        cout << lehmer_pi(n) << endl;    }    return 0;}

PollardRho大整数分解

Pollard rho算法的原理就是通过某种方法得到两个整数a和b，而待分解的大整数为n，计算p=gcd(a-b,n)，直到p不为1，或者a，b出现循环为止。然后再判断p是否为n，如果p=n成立，那么返回n是一个质数，否则返回p是n的一个因子，那么我们又可以递归的计算Pollard(p)和Pollard(n/p)，这样，我们就可以求出n的所有质因子。
模板如下：

#include <iostream>#include <stdlib.h>#include <string.h>#include <algorithm>#include <stdio.h>const int Times = 10;const int N = 5500;using namespace std;typedef long long LL;LL ct, cnt;LL fac[N], num[N];LL gcd(LL a, LL b){    return b? gcd(b, a % b) : a;}LL multi(LL a, LL b, LL m){    LL ans = 0;    a %= m;    while(b)    {        if(b & 1)        {            ans = (ans + a) % m;            b--;        }        b >>= 1;        a = (a + a) % m;    }    return ans;}LL quick_mod(LL a, LL b, LL m){    LL ans = 1;    a %= m;    while(b)    {        if(b & 1)        {            ans = multi(ans, a, m);            b--;        }        b >>= 1;        a = multi(a, a, m);    }    return ans;}bool Miller_Rabin(LL n){    if(n == 2) return true;    if(n < 2 || !(n & 1)) return false;    LL m = n - 1;    int k = 0;    while((m & 1) == 0)    {        k++;        m >>= 1;    }    for(int i=0; i<Times; i++)    {        LL a = rand() % (n - 1) + 1;        LL x = quick_mod(a, m, n);        LL y = 0;        for(int j=0; j<k; j++)        {            y = multi(x, x, n);            if(y == 1 && x != 1 && x != n - 1) return false;            x = y;        }        if(y != 1) return false;    }    return true;}LL pollard_rho(LL n, LL c){    LL i = 1, k = 2;    LL x = rand() % (n - 1) + 1;    LL y = x;    while(true)    {        i++;        x = (multi(x, x, n) + c) % n;        LL d = gcd((y - x + n) % n, n);        if(1 < d && d < n) return d;        if(y == x) return n;        if(i == k)        {            y = x;            k <<= 1;        }    }}void find(LL n, int c){    if(n == 1) return;    if(Miller_Rabin(n))    {        fac[ct++] = n;        return ;    }    LL p = n;    LL k = c;    while(p >= n) p = pollard_rho(p, c--);    find(p, k);    find(n / p, k);}int main(){    LL n;    while(cin>>n)    {        ct = 0;        find(n, 120);        sort(fac, fac + ct);        num[0] = 1;        int k = 1;        for(int i=1; i<ct; i++)        {            if(fac[i] == fac[i-1])                ++num[k-1];            else            {                num[k] = 1;                fac[k++] = fac[i];            }        }        cnt = k;        for(int i=0; i<cnt; i++)            cout<<fac[i]<<"^"<<num[i]<<" ";        cout<<endl;    }    return 0;}

幂乘快速计算（二分思想）

逆向思维，不断地二分，累积底数的平方即可：

typedef long long LL;LL pow_mod(LL a, LL b, LL p){//a的b次方求余p     LL ret = 1;    while(b){        if(b & 1) ret = (ret * a) % p;        a = (a * a) % p;        b >>= 1;    }    return ret;}

衍生出ll的快速乘：

LL mul(LL a, LL b, LL p){//快速乘，计算a*b%p     LL ret = 0;    while(b){        if(b & 1) ret = (ret + a) % p;        a = (a + a) % p;        b >>= 1;    }    return ret;}

关于计算结果过除以MOD取余的问题：
1、计算加法每相加一次执行一次取余；
2、计算减法给被减数加上MOD之后先减法后取余；
3、计算乘法每相乘一次取余一次。

矩阵快速幂

算一个n*n矩阵的的x次方对1e9+7求余：

#include <bits/stdc++.h>using namespace std;typedef long long ll;const int maxn = 110;const int MOD = 1e9+7;#define mod(x) ((x)%MOD)int n;struct mat{    int m[maxn][maxn];}unit;mat operator * (mat a,mat b){    mat ret;    ll x;    for(int i = 0;i < n;i++){        for(int j = 0;j < n;j++){            x = 0;            for(int k = 0;k < n;k++)                x += mod((ll)a.m[i][k]*b.m[k][j]);            ret.m[i][j] = mod(x);        }    }    return ret;}void init_unit(){    for(int i = 0;i < maxn;i++)        unit.m[i][i] = 1;    return ;}mat pow_mat(mat a,ll n){    mat ret = unit;    while(n){        if(n&1) ret = ret*a;        a = a*a;        n >>= 1;    }    return ret;}int main(){    ll x;    init_unit();    while(cin>>n>>x){        mat a;        for(int i = 0;i < n;i++)            for(int j = 0;j < n;j++)                cin>>a.m[i][j];        a = pow_mat(a,x);        for(int i = 0;i < n;i++){            for(int j = 0;j < n;j++){                if(j + 1 == n)  cout << a.m[i][j] << endl;                else cout<<a.m[i][j]<<" ";            }        }    }    return 0;}

GCD&LCM的相关

计算gcd循环写法和递归写法：

LL gcd(LL a, LL b){    if(b == 0) return a;    else return gcd(b, a%b);}LL gcd(LL a, LL b){    return b ? gcd(b, a%b) : a;}LL gcd(LL a, LL b){    LL t;    while(b){        t = b;        b = a % b;        a = t;    }    return a;}

相关的公式：

a*b = gcd(a,b) * lcm(a,b)gcd(ka, kb) = k * gcd(a, b)lcm(ka, kb) = k * lcm(a, b)lcm(S/a, S/b) = S/gcd(a, b)

扩展欧几里德算法

那么已知 a，b 求 一组解 x，y 满足 ax+by = gcd(a, b) 这个公式:

代码模板：

#include<cstdio>typedef long long LL;void ex_Eulid(LL a, LL b, LL &x, LL &y, LL &d){    if (!b) {d = a, x = 1, y = 0;}    else{        extend_Eulid(b, a % b, y, x, d);        y -= x * (a / b);    }}int main(){    LL a, b, d, x, y;    while(~scanf("%lld%lld", &a, &b)){        ex_Eulid(a, b, x, y, d);        printf("%lld*a + %lld*b = %lld\n", x, y, d);    }}

威尔逊定理

当且仅当p为素数时：( p -1 )! ≡ -1 ( mod p )
或者这么写( p -1 )! ≡ p-1 ( mod p )
若p为质数，则p能被(p-1)!+1整除
在初等数论中，这是威尔逊给出了判定一个自然数是否为素数的充分必要条件。

欧拉定理

欧拉定理，也称费马-欧拉定理
若n,a为正整数，且n,a互质，即gcd(a,n) = 1，则a^φ(n) ≡ 1 (mod n)

φ(n) 是欧拉函数，欧拉函数是求小于等于n的数中与n互质的数的数目
欧拉函数是求 （小于n的数 ）中 （与n互质的数 ）的数目或者说欧拉函数是求 1到n-1 中 与n互质的数的数目。

如果n是质数那么1到n-1所有数都是与n互质的，所以φ(n) = n-1
欧拉函数模板：

#include<cstdio>const int N = 100000 + 5;int phi[N];void Euler(){    phi[1] = 1;    for(int i = 2; i < N; i ++){        if(!phi[i]){            for(int j = i; j < N; j += i){                if(!phi[j]) phi[j] = j;                phi[j] = phi[j] / i * (i-1);            }        }    }}int main(){    Euler();}

另一种，比上面更快的方法需要用到如下性质

p为质数时有：
phi(p)=p-1 因为质数p除了1以外的因数只有p，故1至p的整数只有p与p不互质 如果i mod p = 0, 那么 phi(i * p)=phi(i) * p
若i mod p ≠0, 那么 phi( i * p )=phi(i) * ( p-1 )

#include<cstdio>using namespace std;const int N = 1e6+10 ;int phi[N], prime[N];int tot;//tot计数，表示prime[N]中有多少质数 void Euler(){    phi[1] = 1;    for(int i = 2; i < N; i ++){        if(!phi[i]){            phi[i] = i-1;            prime[tot ++] = i;        }        for(int j = 0; j < tot && 1ll*i*prime[j] < N; j ++){            if(i % prime[j]) phi[i * prime[j]] = phi[i] * (prime[j]-1);            else{                phi[i * prime[j] ] = phi[i] * prime[j];                break;            }        }    }}int main(){    Euler();}

中国剩余定理

CRT问题代码描述如下：

void extend_Euclid(LL a, LL b, LL &x, LL &y){    if(b == 0)    {        x = 1;        y = 0;        return;    }    extend_Euclid(b, a % b, x, y);    LL tmp = x;    x = y;    y = tmp - (a / b) * y;}int CRT(int a[],int m[],int n){    int M = 1;    int ans = 0;    for(int i=1; i<=n; i++)        M *= m[i];    for(int i=1; i<=n; i++)    {        int x, y;        int Mi = M / m[i];        extend_Euclid(Mi, m[i], x, y);        ans = (ans + Mi * x * a[i]) % M;    }    if(ans < 0) ans += M;    return ans;}

费马小定理

假如p是质数：
若p不能整除a，则 a^(p-1) ≡1（mod p）；
若p能整除a，则a^(p-1) ≡0（mod p）。
或者说，若p是质数，且a,p互质，那么 a的(p-1)次方除以p的余数恒等于1。