码迷，mamicode.com
首页 > 其他好文 > 详细

PWJ的数论模板整理

时间：2019-09-19 12:16:40      阅读：86      评论：0      收藏：0      [点我收藏+]

标签：unsigned   模板   sub   dfa   ctime   方程组   因子   复杂   color   

一些还没学到，但已经听说的就先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;
}

二次剩余

求x2Ξ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;}
    else{gcd_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;
}

 

PWJ的数论模板整理

标签：unsigned   模板   sub   dfa   ctime   方程组   因子   复杂   color   

原文地址：https://www.cnblogs.com/LMCC1108/p/11548240.html
踩
(0)
赞
(0)
   
举报
评论 一句话评论（0）
登录后才能评论！
分享档案
更多>
2021年07月29日 (22)
2021年07月28日 (40)
2021年07月27日 (32)
2021年07月26日 (79)
2021年07月23日 (29)
2021年07月22日 (30)
2021年07月21日 (42)
2021年07月20日 (16)
2021年07月19日 (90)
2021年07月16日 (35)
周排行
mamicode.com排行更多图片
更多
友情链接
兰亭集智   国之画   百度统计   站长统计   阿里云   chrome插件   新版天听网
关于我们 - 联系我们 - 留言反馈
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码！