数论总结与代码

最新推荐文章于 2023-03-17 11:37:53 发布

潇洒走一回LW

最新推荐文章于 2023-03-17 11:37:53 发布

阅读量1.2k

点赞数

CC 4.0 BY-SA版权

本文链接：https://blog.youkuaiyun.com/u011466175/article/details/20218021

本文深入讲解了多种数论算法，包括素数打表、欧几里得算法、快速幂取模、欧拉phi函数、乘法逆元、模方程求解、Lucas定理、指数循环节、高斯消元法以及大素数判断算法。通过这些算法的学习，读者将能够解决复杂的数学问题。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

一.素数打表：

1.平常打表：

for(i=2;i<=n;i++)    
        if(!s[i])    
        {    
            for(j=2*i;j<=n;j+=i)    
                s[j]=1;    
        }

2.线性打表：

void get_prime()  
{  
    int cnt = 0;  
    for (int i = 2; i < N; i++)  
    {  
        if (!tag[i])  
            p[cnt++] = i;  
        for (int j = 0; j < cnt && p[j] * i < N; j++)  
        {  
            tag[i*p[j]] = 1;  
            if (i % p[j] == 0)  
                break;  
        }  
    }  
}

二.欧几里得与扩展欧几里得算法：

1.欧几里得：

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

2.扩展欧几里得：

//求整数x和y，使得ax+by=d，且|x|+|y|最小。其中d=gcd(a,b)
void extend_gcd(int a,int b,int &d,ll &x,ll &y)
{
    if(!b) d=a,x=1,y=0;
    else {extend_gcd(b,a%b,d,y,x);y-=x*(a/b);}
}

三.快速幂取模：

//求a^b%n
ll mod(ll a,ll b,int n)
{
    if(b==0) return 1;
    ll ans=mod(a,b/2,n);
    ans=ans*ans%n;
    if(b&1) ans=ans*a%n;
    return ans;
}

四.欧拉phi函数：

欧拉函数phi(n)等于不超过n且和n互素的整数个数。互素：两个数的最大公约数为1的数称为互素数。

1.单个欧拉函数求法代码：

int euler_phi(int n)
{
    int m=sqrt(n+0.5),ans=n;
    for(int i=2;i<=m;i++)
        if(n%i==0)
        {
            ans=ans/i*(i-1);
            while(n%i==0) n/=i;
        }
    if(n>1) ans=ans/n*(n-1);
    return ans;
}

2.函数表：多个欧拉函数一起求，类似筛选法计算phi(1),phi(2),……phi(n).

int phi[MM]
void phi_table(int n)
{
    mem(phi,0); 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);
            }
}

五.乘法逆：

在n的同余类中的两个数a和b满足ab=1,比如在15的同余类中7*13=1，在这种情况下，说a和b互为乘法的逆。即7的逆的13，13的逆为7；如果知道a就可以求b，知道b也就可以求a。代码如下：

1.用扩展欧几里得求：

//计算模n下a的逆。如果不存在逆，返回-1
int inv(int a,int n)
{
    int d,x,y;
    extend_gcd(a,n,d,x,y); //扩展欧几里得函数
    return d==1?(x+n)%n:-1;
}

2.利用欧拉定理求：

//a的逆就是mod(a,n-2,n)
ll mod(ll a,ll n-2,int n)
{
    if(b==0) return 1;
    ll ans=mod(a,b/2,n);
    ans=ans*ans%n;
    if(b&1) ans=ans*a%n;
    return ans;
}

六.模方程：

1.线性模方程：axºb(mod n)

给a线性模方程：ax ºb(mod n)。求出x,因为x有多个解，所以代码也有多个解。把其化成ax-ny=b，d=gcd(a,n)不是b的约数时无解，否则两边同时除以d……即利用扩展欧几里得求解。

//返回0时不存在解，为1有解
int mod_gcd(int a,int b,int n)
{
    int x,y,d,x0,i;
    extend_gcd(a,n,d,x,y);
    if(b%d) return 0;
    x0=x*(b/d)%n;
    for(i=1;i<d;i++)
        printf("%d\n",(x0+i*(n/d))%n);
    return 1;
}

2.费马小定理：a^(p-1) ≡1（mod p）

费马小定理是数论中的一个重要定理，其内容为：假如p是质数，且(a,p)=1，那么 a^(p-1) ≡1（mod p）。

即：假如p是质数，且a,p互质，那么a的(p-1)次方除以p的余数恒等于1。

3.中国剩余定理：xºai(mod mi)

如果线性模方程有多个，xºai(mod mi)该怎么做呢？这里就可以用中国剩余定理了。令M为所有mi的乘积，则在模M的剩余系下原方程组有唯一解。

但是用中国剩余定理的条件是mi之间是两两互质的数：

//n个方程：x=a[i](mod m[i]) (0<=i<n)
ll china(int n,int *a,int *m)
{
    ll M=1,d,y,x=0;
    for(int i=0;i<n;i++)
        M*=m[i];
    for(int i=0;i<n;i++)
    {
        ll w=M/m[i];
        gcd(m[i],w,d,d,y);
        x=(x+y*w*a[i])%M;
    }
    return (x+M)%M;
}

如果mi不是两两互质的数，那该怎么办呢？那就得要扩展欧几里得算法把两个等式合并成一个了，然后继续下去求了……（POJ 2891)

ll gcd(ll a,ll b)  
{  
    ll t=(!b?a:gcd(b,a%b));   //记忆递归省点时间  
    return t;  
}  
void exgcd(ll a,ll b,ll &x,ll &y)  
{  
    if(!b) {x=1;y=0;return;}  
    exgcd(b,a%b,y,x);  //这里用记忆化递归时间还是一样的，所以不用  
    y-=a/b*x;  
}  
int main()  
{  
    ll i,k,a1,r1,a2,r2,c,x,y,l,t;  
    while(cin>>k)  
    {  
        int flag=0;  
        cin>>a1>>r1;  
        for(i=1;i<k;i++)  
        {  
            cin>>a2>>r2;  
            if(flag) continue;  
            c=r2-r1;  
            l=gcd(a1,a2);  
            if(c%l!=0) {flag=1;continue;}  
            exgcd(a1,a2,x,y);  
            t=a2/l;  
            x=((c/l*x)%t+t)%t;  
            r1+=a1*x;  
            a1*=t;  
        }  
        if(flag) cout<<-1<<endl;  
        else cout<<r1<<endl;  
    }  
    return 0;  
}

4.高次同余模方程求解：a^x ºb(mod n)

利用欧拉定理求解模方程，当n为素数时，解模方程a^xºb(mod n)。根据欧拉定理，只需检查x=0,1,2,……,n-1是不是解即可因为a^(n-1)º1(mod n),当x超过n-1时a^x就开始循环了。

//求解模方程a^xb(mod n)。n为素数。无解返回0
int log_mod(int a,int b,int n)
{
    int m,v,e=1,i;
    m=sqrt(n+0.5);
    v=inv(mod(a,m,n),n);
    map<int,int>x;
    x[1]=0;
    for(i=1;i<m;i++)
    {
        e=e*a%n;
        if(!x.count(e)) x[e]=i;
    }
    for(i=0;i<m;i++)
    {
        if(x.count(b)) return i*m+x[b];
        b=b*v%n;
    }
    return 0;
}

哈希法的高次同余方程模板：优点快速，代码少。

#define MOD 76543
int hs[MOD],head[MOD],next[MOD],id[MOD],top;
void insert(int x,int y)
{
    int k = x%MOD;
    hs[top] = x, id[top] = y, next[top] = head[k], head[k] = top++;
}
int find(int x)
{
    int k = x%MOD;
    for(int i = head[k]; i != -1; i = next[i])
        if(hs[i] == x)
            return id[i];
    return -1;
}
int BSGS(int a,int b,int n) //a^L=b(mod n) 求L
{
    memset(head,-1,sizeof(head));
    top = 1;
    if(b == 1)return 0;
    int m = sqrt(n*1.0), j;
    long long x = 1, p = 1;
    for(int i = 0; i < m; ++i, p = p*a%n)insert(p*b%n,i);
    for(long long i = m; ; i += m)
    {
        if( (j = find(x = x*p%n)) != -1 )return i-j;
        if(i > n)break;
    }
    return -1;  //L不存在的情况
}

还有就是算法经典入门里的模板，用的是map<int,int>x，缺点：太慢，而且代码又多！在POJ 2417上TLE了。

void extend_gcd(ll a,ll b,ll &d,ll &x,ll &y)
{
    if(!b) d=a,x=1,y=0;
    else {extend_gcd(b,a%b,d,y,x);y-=x*(a/b);}
}
ll mod(ll a,ll b,ll n)
{
    if(b==0) return 1;
    ll ans=mod(a,b/2,n);
    ans=ans*ans%n;
    if(b&1) ans=ans*a%n;
    return ans;
}
ll inv(ll a,ll n)
{
    ll d,x,y;
    extend_gcd(a,n,d,x,y); //扩展欧几里得函数
    return d==1?(x+n)%n:-1;
}
ll log_mod(ll a,ll b,ll n)
{
    ll m,v,e=1,i;
    m=sqrt(n+0.5);
    v=inv(mod(a,m,n),n);
    map<int,int>x;
    x[1]=0;
    for(i=1;i<m;i++)
    {
        e=e*a%n;
        if(!x.count(e)) x[e]=i;
    }
    for(i=0;i<m;i++)
    {
        if(x.count(b)) return i*m+x[b];
        b=b*v%n;
    }
    return -1;
}

七.Lucas定理：

A、B是非负整数，p是质数。A B写成p进制：A=a[n]a[n-1]...a[0]，B=b[n]b[n-1]...b[0]。
则组合数C(A,B)与C(a[n],b[n])*C(a[n-1],b[n-1])*...*C(a[0],b[0]) mod p同余

即：Lucas(n,m,p)=C(n%p,m%p)*Lucas(n/p,m/p,p)

求大组合数取模：C(n,m)%p，大组合数取模。

ll fac[100003];
void init(ll p)
{
    fac[0] = 1;
    for (int i=1; i<=p; i++) fac[i] = fac[i-1]*i%p;
}
ll PowerMod(ll a, ll b, ll k)
{
    ll tmp = a, ret = 1;
    while (b)
    {
        if (b & 1) ret = ret * tmp % k;
        tmp = tmp * tmp % k;
        b >>= 1;
    }
    return ret;
}
ll Lucas(ll n, ll m, ll p)
{
    ll ret = 1;
    while (n && m)
    {
        ll nn = n%p, mm = m%p;
        if (nn < mm) return 0;
        ret = ret*fac[nn]*PowerMod(fac[mm]*fac[nn-mm]%p, p-2, p)%p;
        n /= p;
        m /= p;
    }
    return ret;
}

八.指数循环节：

即求：a[0]^a[1]^……^a[n]的值，比如：2^4^3^5=2。

//求a[0]^a[1]^……^a[n]%m的值（指数循环节）
LL a[20],n,m;
LL euler(LL n)
{
    LL m=(LL)sqrt(n+0.5);
    LL ans=n;
    for(LL i=2; i<=m; i++)
        if(n%i==0)
        {
            ans=ans/i*(i-1);
            while(n%i==0)n/=i;
        }
    if(n>1)ans=ans/n*(n-1);
    return ans;
}
LL pmod(LL a,LL n,LL m)
{
    LL now = 1;
    for(LL i=0; i<n; i++)
    {
        now *= a;
        if(now>=m)  break;
    }
    if(now >= m)   now = m;
    else now=0;

    if(n==0)return 1;
    LL x=pmod(a,n/2,m);
    LL ans=(LL)x*x%m;
    if(n%2==1) ans=ans*a%m;
    return ans+now;
}
LL solve(LL cur,LL n,LL m)
{
    if(cur==n-1)
    {
        if(a[cur]>=m)
            return a[cur]%m+m;
        else
            return a[cur];
    }
    LL M=euler(m);
    LL p=solve(cur+1,n,M);
    LL ans=pmod(a[cur],p,m);
    return ans;
}
//输出的是solve(0,n,m)%m

九.高斯消元法：

//要求系数矩阵可逆
//这里的A是增广矩阵，即A[i][n]是第i个方程右边的常数bi。
//运行结束后A[i][n]是第i个未知数的值
typedef double Matrix[MM][MM];
void gauss(Matrix A,int n)
{
    int i,j,k,r;
    //消元过程
    for(i=0;i<n;i++)
    {
        //选一行r并与第i交换
        r=i;
        for(j=i+1;j<n;j++)
            if(fabs(A[i][j])>fabs(A[r][i])) r=j;
        if(r!=i)
        {
            for(j=0;j<=n;j++)
                swap(A[r][j],A[i][j]);
        }
        //与第i+1~n行进行消元
        for(j=n;j>=i;j--)
            for(k=i+1;k<n;k++)
                A[k][j]-=A[k][i]/A[i][i]*A[i][j];
    }
    for(i=n-1;i>=0;i--)
    {
        for(j=i+1;j<n;j++)
            A[i][n]-=A[j][n]*A[i][j];
        A[i][n]/=A[i][i];
    }
}

十.大素数判断Miller_Rabin 算法：

判断的句子是：if(Is_Prime(n)) cout<<"YES"<<endl;

const int TIMES=200;  //判断次数可以减少，但是越多，判断数组的准确率越高
long long Gcd(long long a,long long b)
{
    if(b==0)
    {
        return a;
    }
    else
    {
        return Gcd(b,a%b);
    }
}
long long Cheng_Mod(long long a,long long x,long long p)
{
    long long Ret=0,s=a;
    while(x>0)
    {
        if(x%2==1)
        {
            Ret=(Ret+s)%p;
        }
        s=(s+s)%p;
        x/=2;
    }
    return Ret;
}
long long Mult_Mod(long long a,long long x,long long p)
{
    long long Ret=1,s=a;
    while(x>0)
    {
        if(x%2==1)
        {
            Ret=Cheng_Mod(Ret,s,p);
        }
        s=Cheng_Mod(s,s,p);
        x/=2;
    }
    return Ret;
}
bool Is_Prime(long long p)
{
    if(p==1)
    {
        return 0;
    }
    if(p<=3)
    {
        return 1;
    }
    int i,j;
    long long b=0,m=p-1,d=p-1,x;
    while(m%2==0)
    {
        b++;
        m/=2;
    }
    for(i=1; i<=TIMES; i++)
    {
        x=rand()%(p-3)+2;
        d=Mult_Mod(x,m,p);
        if(d==1||d==p-1)
        {
            continue;
        }
        for(j=1; j<=b; j++)
        {
            d=Cheng_Mod(d,d,p);
            if(d==1)
            {
                return 0;
            }
            if(d==p-1)
            {
                break;
            }
        }
        if(d!=p-1)
        {
            return 0;
        }
    }
    return 1;
}