PWJ的数论模板整理

Posted 2022-11-11 lmcc1108

一些还没学到，但已经听说的就先copy其他博客的

数论

欧拉降幂

求a1^a2^a3^a4^a5^a6 mod m

#include<cstdio>
#include<cstring>
const int N=1e4+11;
typedef long long ll;
char s[10]; 
int n,lens,phi[N];
ll md,a[15];
void init()
    for(int i=1;i<N;i++) phi[i]=i;
    for(int i=2;i<N;i+=2) phi[i]/=2;
    for(int i=3;i<N;i+=2) if(phi[i]==i)
        for(int j=i;j<N;j+=i) phi[j]=phi[j]/i*(i-1);

ll poww(ll x,ll y,ll z)
    ll ans=1;
    if(x>=z) x=x%z+z;
    while(y)
        if(y&1)
            ans*=x;
            if(ans>=z) ans=ans%z+z;
        
        x*=x;
        if(x>=z) x=x%z+z;
        y>>=1;
    
    return ans;

ll dfs(int x,ll y)
    if(x==n-1) return a[x]>=y ? a[x]%y+y : a[x];
    if(a[x]%y==0) return 0;
    return poww(a[x],dfs(x+1,phi[y]),y);

int main()
    init();
    int t=1;
    while(~scanf("%s",s)&&s[0]!=‘#‘)
        md=0;lens=strlen(s);
        for(int i=0;i<lens;i++) md=md*10+(s[i]-‘0‘); 
        scanf("%d",&n);
        for(int i=0;i<n;i++) scanf("%lld",&a[i]);
        printf("Case #%d: %lld\n",t++,dfs(0,md)%md); 
    
    return 0;

广义斐波那契循环节

int pn,pri[100000+10];
void init()
    ll lim;
    for(int i=2;i<=100000;i++)
        if(!pri[i]) pri[pn++]=i;
        for(int j=0;j<pn;j++)
            lim=1ll*i*pri[j];
            if(lim>100000) break;
            pri[lim]=1;
            if(i%pri[j]==0) break; 
        
    

ll solve(ll x)
    ll ans1=1,ans2=1,xx=x;
    for(int i=0;i<pn;i++)
        if(1ll*pri[i]*pri[i]>x) break;
        if(x%pri[i]==0)
            ans1*=(pri[i]-1)*(pri[i]+1);
            ans2*=pri[i];
            while(x%pri[i]==0) x/=pri[i];
        
    
    if(x>1)
        ans1*=(x-1)*(x+1);
        ans2*=x;
    
    return xx/ans2*ans1;

二次剩余

求x²Ξa(mod m)的解x

#include<cstdio>
#include<cstdlib>
#include<ctime>
struct Ima
    int x,y;
;
int p,w;
Ima muli(const Ima &i1,const Ima &i2)
    Ima ans;
    ans.x=(i1.x*i2.x%p+i1.y*i2.y%p*w%p)%p;
    ans.y=(i1.x*i2.y%p+i1.y*i2.x%p)%p;
    return ans;

Ima powi(Ima a,int b)
    Ima ans;
    ans.x=1,ans.y=0;
    while(b)
        if(b&1) ans=muli(ans,a);
        a=muli(a,a);
        b>>=1;
    
    return ans;

int poww(int a,int b)
    int ans=1;
    a%=p;
    while(b)
        if(b&1) ans=ans*a%p;
        a=a*a%p;
        b>>=1;
    
    return ans;

int Cipolla(int n)
    if(p==2) return 1;
    if(poww(n,(p-1)>>1)+1==p) return -1;
    int a;
    while(true)
        a=rand()%p;
        w=((a*a%p-n)%p+p)%p;
        if(poww(w,(p-1)>>1)+1==p) break;
    
    Ima ans;
    ans.x=a,ans.y=1;
    ans=powi(ans,(p+1)>>1);
    return ans.x;

int main()
    int t,n,ans1,ans2;
    srand(time(NULL));
    scanf("%d",&t);
    while(t--)
        scanf("%d%d",&n,&p);
        n%=p;
        ans1=Cipolla(n),ans2=p-ans1;
        if(ans1==-1) printf("No root\n");
        else if(ans1==ans2) printf("%d\n",ans1);
        else if(ans1<ans2) printf("%d %d\n",ans1,ans2);
        else printf("%d %d\n",ans2,ans1); 
    
    return 0;

大素数判断

typedef unsigned long long  ll;
const int S=20;
ll mult_mod(ll a, ll b, ll c)
    a%=c;
    b%=c;
    ll ret=0,tmp=a;
    while(b)
        if(b&1)
            ret+=tmp;
            if(ret>c) ret-=c;
        
        tmp<<=1;
        if(tmp>c) tmp-=c;
        b>>=1;
    
    return ret;

ll pow_mod(ll a, ll n, ll mod)
    ll ret=1,tmp=a%mod;
    while(n)
        if(n&1)    ret=mult_mod(ret,tmp,mod);
        tmp=mult_mod(tmp,tmp,mod);
        n >>=1;
    
    return ret;

bool check(ll a, ll n, ll x, ll t)
    ll ret=pow_mod(a,x,n);
    ll last=ret;
    for(int i=1;i<=t;i++)
        ret=mult_mod(ret,ret,n);
        if(ret==1&&last!=1&&last!=n-1)    return true;
        last=ret;
    
    if(ret!=1) return true;
    else return false;

bool Miller_Rabin(ll n)
    if(n<2)    return false;
    if(n==2) return true;
    if((n&1)==0) return false;
    ll x=n-1,t=0;
    while((x&1)==0)
        x>>=1;
        t++;
    
//    srand(time(NULL));
    for(int i=0;i<S;i++)
        ll a=rand()%(n-1)+1;
        if(check(a,n,x,t))    return false;
    
    return true;

质因子分解

long long factor[100];//质因数分解结果（刚返回时是无序的）
int tol;//质因数的个数。数组小标从0开始

long long gcd(long long a,long long b)

    if(a==0)return 1;//???????
    if(a<0) return gcd(-a,b);
    while(b)
    
        long long t=a%b;
        a=b;
        b=t;
    
    return a;


long long Pollard_rho(long long x,long long c)

    long long i=1,k=2;
    long long x0=rand()%x;
    long long y=x0;
    while(1)
    
        i++;
        x0=(mult_mod(x0,x0,x)+c)%x;
        long long d=gcd(y-x0,x);
        if(d!=1&&d!=x) return d;
        if(y==x0) return x;
        if(i==k)y=x0;k+=k;
    

//对n进行素因子分解
void findfac(long long n)

    if(Miller_Rabin(n))//素数
    
        factor[tol++]=n;
        return;
    
    long long p=n;
    while(p>=n)p=Pollard_rho(p,rand()%(n-1)+1);
    findfac(p);
    findfac(n/p);

中国剩余定理

//n个方程：x=a[i](mod m[i]) (0<=i<n)
LL china(int n, LL *a, LL *m)
    LL M = 1, ret = 0;
    for(int i = 0; i < n; i ++) M *= m[i];
    for(int i = 0; i < n; i ++)
        LL w = M / m[i];
        ret = (ret + w * inv(w, m[i]) * a[i]) % M;
    
    return (ret + M) % M;

扩展中国剩余定理

typedef long long LL;
typedef pair<LL, LL> PLL;
PLL linear(LL A[], LL B[], LL M[], int n) //求解A[i]x = B[i] (mod M[i]),总共n个线性方程组 
    LL x = 0, m = 1;
    for(int i = 0; i < n; i ++) 
        LL a = A[i] * m, b = B[i] - A[i]*x, d = gcd(M[i], a);
        if(b % d != 0)  return PLL(0, -1);//答案不存在，返回-1 
        LL t = b/d * inv(a/d, M[i]/d)%(M[i]/d);
        x = x + m*t;
        m *= M[i]/d;
    
    x = (x % m + m ) % m;
    return PLL(x, m);//返回的x就是答案，m是最后的lcm值 

java实现

import java.util.*;
import java.math.*;
public class Main 
    static BigInteger x,y;
    public static void main(String[] args)
        Scanner cin=new Scanner(System.in);
        int n=cin.nextInt();
        BigInteger m=cin.nextBigInteger();
        BigInteger x=BigInteger.ZERO,mo=BigInteger.ONE;
        BigInteger B,M,b,d,t;
        int i;
        bool flag=true;
        for(i=0;i<n;++i)
            M=cin.nextBigInteger();
            B=cin.nextBigInteger();
            b=B.subtract(x);
            d=gcd(M,mo);
            if(b.mod(d).compareTo(BigInteger.ZERO)!=0) flag=false;
            if(!flag) continue;
            t=((b.divide(d)).multiply((mo.divide(d)).modInverse(M.divide(d)))).mod(M.divide(d));
            x=x.add(mo.multiply(t));
            mo=mo.multiply(M.divide(d));
        
        if(!flag) System.out.println("无解");
        else
            x=((x.mod(mo)).add(mo)).mod(mo);
            System.out.println(x);
        
    
    static BigInteger gcd(BigInteger a,BigInteger b)
        if(b.compareTo(BigInteger.ZERO)==0)return a;
        return gcd(b,a.mod(b));
    

一阶线性同余方程

通解为r+a*k r为最小非负整数解 a为lcm(a1,a2,a3...,an)

#include<cstdio>
typedef long long ll;
const int N=1e4+11; 
int n;
ll aa[N],rr[N];
ll exgcd(ll a,ll b,ll &x,ll &y)
    if(!b)
        x=1;
        y=0;
        return a;
    
    ll g=exgcd(b,a%b,y,x);
    y-=a/b*x;
    return g;

ll linemod()
    ll a,b,c,r,x,y,g,k;
    a=aa[0],r=rr[0];
    for(int i=1;i<n;i++)
        b=aa[i];
        c=rr[i]-r;
        g=exgcd(a,b,x,y);
        if(c%g!=0) return -1;
        k=b/g;
        x=(x*(c/g)%k+k)%k;
        r=x*a+r;
        a=a/g*b;
    
    return r;

int main()
    while(~scanf("%d",&n))
        for(int i=0;i<n;i++) scanf("%lld%lld",&aa[i],&rr[i]); 
        printf("%lld\n",linemod());
    
    return 0;

拔山盖世算法

求a^x=b(mod n)的x

#include <cstdio>
#include <cstring>
#include <cmath>
#include <map>
#include <iostream>
#include <algorithm>
using namespace std;
#define LL __int64
LL gcd(LL a,LL b)

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

//拓展欧几里得定理，求ax+by=gcd(a,b)的一组解(x,y),d=gcd(a,b)
void gcd_mod(LL a,LL b,LL &d,LL &x,LL &y)

    if(!b)d=a;x=1;y=0;
    elsegcd_mod(b,a%b,d,y,x);y-=x*(a/b);

//求解模方程d*a^(x-c)=b(mod n)。d,a和n互质，无解返回-1
LL log_mod (LL a,LL b,LL n,LL c,LL d)

    LL m,v,e=1,i,x,y,dd;
    m=ceil(sqrt(n+0.5));     //x=i*m+j
    map<LL,LL>f;
    f[1]=m;
    for(i=1;i<m;i++)  //建哈希表，保存a^0,a^1,...,a^m-1
    
        e=(e*a)%n;
        if(!f[e])f[e]=i;
    
    e=(e*a)%n;//e=a^m
    for(i=0;i<m;i++)//每次增加m次方，遍历所有1<=f<=n
    
        gcd_mod(d,n,dd,x,y);//d*x+n*y=1-->(d*x)%n=1-->d*(x*b)%n==b
        x=(x*b%n+n)%n;
        if(f[x])
        
            LL num=f[x];
            f.clear();//需要清空，不然会爆内存
            return c+i*m+(num==m?0:num);
        
        d=(d*e)%n;
    
 
    return -1;

int main()

    LL a,b,n;
    while(scanf("%I64d%I64d%I64d",&a,&n,&b)!=EOF)
    
        if(b>=n)
        
            printf("Orz,I can’t find D!\n");
            continue;
        
        if(b==0)
        
            printf("0\n");
            continue;
        
        LL ans=0,c=0,d=1,t;
        while((t=gcd(a,n))!=1)
        
            if(b%t)ans=-1;break;
            c++;
            n=n/t;
            b=b/t;
            d=d*a/t%n;
            if(d==b)ans=c;break;//特判下是否成立。
        
        if(ans!=0)
        
            if(ans==-1)printf("Orz,I can’t find D!\n");
            else printf("%I64d\n",ans);
        
        else
        
            ans=log_mod(a,b,n,c,d);
            if(ans==-1)printf("Orz,I can’t find D!\n");
            else printf("%I64d\n",ans);
        
    
    return 0;

/*
    求解模方程a^x=b(mod n)，n不为素数。拓展Baby Step Giant Step
    模板题。
    
    方法：
    初始d=1,c=0,i=0;
    1.令g=gcd(a,n),若g==1则执行下一步。否则由于a^x=k*n+b;(k为某一整数),则(a/g)*a^k=k*(n/g)+b/g,(b/g为整除，若不成立则无解)
令n=n/g，d=d*a/g，b=b/g,c++则d*a^(x-c)=b(mod n),接着重复1步骤。
    2.通过1步骤后，保证了a和d都与n互质，方程转换为d*a^(x-c)=b(mod n)。由于a和n互质，所以由欧拉定理a^phi(n)==1(mod n),(a,n互质)
可知，phi(n)<=n,a^0==1(mod n),所以构成循环，且循环节不大于n。从而推出如果存在解，则必定1<=x<n。(在此基础上我们就可以用
Baby Step Giant Step方法了)
    3.令m=ceil(sqrt(n)),则m*m>=n。用哈希表存储a^0,a^1,..,a^(m-1)，接着判断1~m*m-1中是否存在解。
    4.为了减少时间，所以用哈希表缩减复杂度。分成m次遍历，每次遍历a^m长度。由于a和d都与n互质，所以gcd(d,n)=1，
所以用拓展的欧几里德定理求得d*x+n*y=gcd(d,n)=1,的一组整数解(x,y)。则d*x+n*y=1-->d*x%n=(d*x+n*y)%n=1-->d*(x*b)%n=((d*x)%n*b%n)%n=b。
所以若x*b在哈希表中存在，值为k，则a^k*d=b(mod n),答案就是ans=k+c+i*m。如果不存在，则令d=d*a^m,i++后遍历下一个a^m，直到遍历a^0到a^(m-1)还未找到，则说明不解并退出。

哈希表实现

struct HashMap//哈希表
    static const int Hash=999917,maxn=46340;
    int num,link[Hash],son[maxn+5],next[maxn+5],w[maxn+5];
    int top,Stack[maxn+5];
    void clear()//清空表
        num=0;
        while(top)
            link[Stack[top--]]=0;
    
    void add(int x,int y)//添加键值元素
        son[++num]=y;
        next[num]=link[x];
        w[num]=INF;
        link[x]=num;
    
    bool count(int y)//判断表中是否有对应值
        int x=y%Hash;
        for(int j=link[x];j;j=next[j])
            if(y==son[j])
                return true;
        return false;
    
    int &operator [](int y)//获取键的对应值
        int x=y%Hash;
        for(int j=link[x];j;j=next[j])
            if(y==son[j])
                return w[j];
        add(x,y);
        Stack[++top]=x;
        return w[num];
    
hashMap;
int GCD(int a,int b)
    if(b==0)
        return a;
    return GCD(b,a%b);

int extendedGCD(int x,int y,int &a,int &b)
    if(y==0)
        a=1;
        b=0;
        return x;
    
 
    int temp;
    int gcd=extendedGCD(y,x%y,a,b);
    temp=a;
    a=b;
    b=temp-x/y*b;
    return gcd;

 
int extendBSGS(int A,int B,int C)
    //三种特判
    if(C==1)
        if(!B)
            return A!=1;
        return -1;
    
    if(B==1)
        if(A)
            return 0;
        return -1;
    
    if(A%C==0)
        if(!B)
            return 1;
        return -1;
    
 
    int gcd=GCD(A,C);
    int D=1,num=0;
    while(gcd!=1)//把A,C变成(A,C)=1为止
        if(B%gcd)
            return -1;
 
        B/=gcd;//从B中约去因子
        C/=gcd;//约C中约去因子
        D=((LL)D*A/gcd)%C;//将多出的乘给D
        num++;//统计约去次数
        gcd=GCD(A,C);
    
 
    int now=1;
    for(int i=0;i<=num-1;i++)//枚举0~num-1
        if(now==B)
            return i;
        now=((LL)now*A)%C;
    
 
    hashMap.clear();
    int Size=ceil(sqrt(C)),Base=1;
    for(int i=0;i<=Size-1;i++)//将A^j存进哈希表
        hashMap[Base]=min(hashMap[Base],i);//存储答案最小的
        Base=((LL)Base*A)%C;
    
    for(int i=0;i<=Size-1;i++)//扩展欧几里得求A^j
        int x,y;
        int gcd=extendedGCD(D,C,x,y);//求出x、y
        x=((LL)x*B%C+C)%C;
        if(hashMap.count(x))
            return i*Size+hashMap[x]+num;//加回num
        D=((LL)D*Base)%C;
    
    return -1;//无解

 
int main()
    int A,B,C;
    while(scanf("%d%d%d",&A,&B,&C)!=EOF&&(A+B+C))
        int res=extendBSGS(A,B,C);
        if(res==-1)
            printf("No Solution\n");
        else
            printf("%d\n",res);
    
    return 0;