Acdreamer博客数论学习(1)

Acdreamer博客数论学习

Day One

等比数列求和．

求解Sn=(a+a2+...+an)modM

采用二分．n%2==0时 Sn=(1+a(n2))Sn2 ,n%2==1时，Sn=(1+an+12)Sn−12+an+12

#include<cstdio>#include<iostream>#define ll long longusing namespace std;const int mod = 7;ll pow(int a, int k){    ll ans = 1;    while(k){        if(k&1)ans=ans*a%mod;        a=a*a%mod;        k>>=1;    }    return ans;}ll S(int a, int n){    if(n==1)return a%mod;    if(n&1){        return ((1+pow(a,(n+1)/2))%mod*S(a,(n-1)/2)%mod + pow(a,(n+1)/2)%mod)%mod;    }else{        return (1+pow(a,n/2))%mod*S(a,n/2)%mod;    }}int main(){    int a,n;    while(~scanf("%d%d",&a,&n)){        printf("%lld\n",S(a,n));    }    return 0;}

题目：https://vjudge.net/problem/POJ-3233．

计算矩阵１－n次方和

int n,k,m;int mod = 1e9+7;#include<cstdio>#include<iostream>#include<cstring>using namespace std;struct Matrix{    int mp[35][35];    Matrix(){        memset(mp,0,sizeof mp);    }    Matrix(int v){        memset(mp,0,sizeof mp);        for(int i = 1; i <= n; i++){            mp[i][i] = v;        }    }};Matrix multi(Matrix a, Matrix b){    Matrix c;    for(int i = 1; i <= n; i++){        for(int k = 1; k <= n; k++){if(a.mp[i][k])            for(int j = 1; j <= n; j++){if(b.mp[k][j])                c.mp[i][j] = c.mp[i][j] + a.mp[i][k]*b.mp[k][j];                c.mp[i][j] %= mod;            }        }    }    return c;}Matrix pls(Matrix a, Matrix b){    Matrix c;    for(int i = 1; i <= n; i++){        for(int j = 1; j <= n; j++){            c.mp[i][j] = a.mp[i][j]+b.mp[i][j];            c.mp[i][j] %= mod;        }    }    return c;}Matrix powmul(Matrix a, int k){    Matrix c = Matrix(1);    for(int i = 1; i <= n; i++)        c.mp[i][i] = 1;    while(k){        if(k&1)c = multi(c,a);        a = multi(a,a);        k>>=1;    }    return c;}Matrix S(Matrix a, int k){    Matrix temp = Matrix(1);    if(k==1)return a;    if(k&1){        Matrix temp2 = powmul(a,(k+1)/2);        return pls(multi(pls(temp, temp2), S(a,(k-1)/2)), temp2);    }else{        Matrix temp2 = powmul(a,k/2);        return multi(pls(temp, temp2), S(a,k/2));    }}int main(){    while(~scanf("%d%d%d",&n,&k,&m)){        Matrix a;        for(int i = 1; i <= n; i++){            for(int j = 1; j <= n; j++){                scanf("%d",&a.mp[i][j]);            }        }        mod = m;        a = S(a,k);        for(int i = 1 ; i <= n; i++){            for(int j = 1; j <= n; j++){                printf("%d",a.mp[i][j]);                if(j!=n)printf(" ");            }            puts("");        }    }    return 0;}

素数定理

π(n)为小于等于n的素数的个数，有limn→∞π(n)n/lnn=1

定理gcd(an−1,am−1)=agcd(n,m)−1

在a>1,m,n>0时成立．

#include<cstdio>#include<iostream>#define ll long longusing namespace std;int mod = 1e9+7;ll gcd(ll a, ll b){    return b?gcd(b,a%b):a;}ll pow(ll a, ll k){    ll ans = 1;    while(k){        if(k&1)ans=ans*a%mod;        a = a*a%mod;        k>>=1;    }    return ans;}int main(){    int a,m,n,k1,t;scanf("%d",&t);    while(t--){        scanf("%d%d%d%d",&a,&m,&n,&k1);        mod = k1;        ll ans = pow(a,gcd(n,m))-1+mod;        printf("%lld\n",ans%mod);    }    return 0;}

定理扩展a>b,gcd(a,b)=1,gcd(am−bm,an−bn)=agcd(m,n)−bgcd(m,n)

定理G=gcd(C1n,C2n,C3n,...,Cnn−1)

1. n为素数，答案为n；
2. n有多个素因子，答案为１；
3. n只有一个素因子,答案就是该素因子(其实包括了第一种情况

题目：https://vjudge.net/problem/HDU-2582

#include<cstdio>#include<iostream>#include<cstring>#define ll long long#define maxn 1000100using namespace std;bool judge[maxn];int prime[maxn/10];ll biao[maxn];int tot = 0;void init(){    memset(prime,0,sizeof prime);    memset(judge,0,sizeof judge);    for(int i = 2; i < maxn; i++){        if(!judge[i]){            for(int j = i+i; j < maxn; j+=i){                judge[j]=1;            }            prime[tot++]=i;        }    }    fill(biao,biao+maxn,1);    for(int i = 0; i < tot; i++){        for(ll j = prime[i]; j < maxn; j*=prime[i]){            biao[j]=prime[i];        }    }    for(int i = 4; i < maxn; i++){        biao[i]+=biao[i-1];    }}int main(){    int n;    init();    while(~scanf("%d",&n)){        printf("%lld\n",biao[n]);    }    return 0;}

定理Fn为斐波那契额数列，gcd(Fm,Fn)=Fgcd(m,n)

http://acm.nyist.net/JudgeOnline/problem.php?pid=468

此处需要学会一种新的判断素数的方式叫做Miller_Rabin素数测试．在DayTwo中有．

#include<stdio.h>#include<algorithm>#include<iostream>using namespace std;long long multi(long long a, long long b, long long mod){    long long temp = a, sum = 0;    while(b){        if(b & 1) sum = (sum + temp) % mod;        temp = (temp + temp) % mod;        b = b >> 1;    }    return sum;}long long Modular_exponent(long long a, long long x, long long mod){    long long t = a % mod, r = 1;    while(x){        if(x & 1) r = multi(r, t, mod);        t = multi(t, t, mod); x = x >> 1;    }    return r;}bool Miller_Rabin(long long n){    long long time = 20;    if(n < 2) return false;    if(n == 2) return true;    bool flag = false;    for(long long k = 0; k <= time; ++k){        flag = false;        long long d = n - 1, r = 0, t, a = rand() % (n - 2) + 2;        while((d & 1) == 0){            d = d>>1;            r++;        }        t = Modular_exponent(a, d, n);        if(t == 1 || t == n-1) {flag = true; continue;}        for(long long i = 1; i < r; i++){            t = multi(t, t, n);            if(t == 1) {flag = false; return flag;}            if(t == n-1) {flag = true; break;}        }        if(!flag) break;    }    return flag;}int main(){    long long n;    while(scanf("%lld",&n)!=EOF)    puts((n == 4 || (n != 2 && Miller_Rabin(n))) ? "Yes" : "No");    return 0;} 

定理：给定两个互素的正整数Ａ和Ｂ，那么他们最大不能组合的数为AB-A-B,不能组合的个数为num=(A−1)(B−1)2

ax+by=c的存在非负解的条件就是c>AB-A-B.

HDU - 1792

#include<cstdio>#include<iostream>using namespace std;int main(){    int a,b;    while(~scanf("%d%d",&a,&b)){        printf("%d %d\n",a*b-a-b,(a-1)*(b-1)/2);    }    return 0;}

定理∑gcd(i,N)=∑dϕ(Nd)

数论专题里也有一题如是说．

https://vjudge.net/problem/POJ-2480

#include<cstdio>#include<cstring>#include<iostream>#define ll long long#define maxn 50000using namespace std;bool judge[maxn];int prime[maxn];int tot;void init(){    tot = 0;    memset(judge,0,sizeof judge);    memset(prime, 0, sizeof prime);    for(int i = 2; i < maxn; i++){        if(!judge[i]){            for(int j = i+i;  j < maxn; j+=i){                judge[j] = 1;            }            prime[tot++]=i;        }    }}int main(){    init();    ll n;    while(~scanf("%lld",&n)){        ll ans = n;        for(int i = 0; i < tot && prime[i]*prime[i] <= n; i++){            if(n%prime[i]==0){                int cot = 0;                while(n%prime[i]==0){cot++;n/=prime[i];}                ans += ans * cot * (prime[i]-1)/prime[i];            }        }        if(n!=1)ans += ans*(n-1)/n;        printf("%lld\n",ans);    }}

剩下没有题目的定理

(n+1)lcm(C0n,C1n,C2n,...,Cn−1N,CNn)=lcm(1,2,...,n+1)

定理：任何n个连续的正整数的乘积均可被n!整除

结论：

如果p是素数，那么C1p,C2p,C3p,...Cp−1p,均能被p整除如果p是素数，那么有(x+y+..+w)p≡xp+yp+...+wp(modp)

Day Two

卢卡斯－莱默检验法 检验梅森素数

梅森素数Mp=2p−1，p为素数．，定义序列{si}，所有的i≥0,si={4,i=0s2i−1−2,i>0Mp是素数当且仅当sp−2≡0(modMp)

#include<cstdio>#include<iostream>#define ll long longusing namespace std;ll T,p;ll multi(ll a, ll b, ll m){    ll ans = 0;    while(b){        if(b&1)ans = (ans+a)%m;        a = (a+a)%m;        b>>=1;    }    return ans;}int main(){    scanf("%lld",&T);    while(T--){        scanf("%lld",&p);        ll n = ((ll) 1<<p)-1;        ll r = 4;        if(p==2)puts("yes");        else{            while((--p)!=1)r=(multi(r,r,n)-2+n)%n;            if(r%n==0)puts("yes");            else puts("no");        }    }}

Miller-Rabin素数测试

​ 费马小定理：对于素数p和任意整数a，有

ap≡a(modp)

（同余）。反过来，满足

ap≡a(modp)

，p也几乎一定是素数。

　　伪素数：如果n是一个正整数，如果存在和n互素的正整数a满足 an−1≡1(modn)，我们说n是基于a的伪素数。如果一个数是伪素数，那么它几乎肯定是素数。

​ 序列:
​ a^m%n; a^(2m)%n; a^(4m)%n; …… ;a^(m*2^q)%n把上述测试序列叫做Miller测试,关于Miller测试,有下面的定理:

定理:若n是素数,a是小于n的正整数,则n对以a为基的Miller测试,结果为真.Miller测试进行k次,将合数当成素数处理的错误概率最多不会超过4^(-k).

　　Miller-Rabin测试：不断选取不超过n-1的基b(s次)，计算是否每次都有

bn−1≡1(modn)

，若每次都成立则n是素数，否则为合数。　

​ 二次探测定理 ：p为素数，０

#include<cstdio>#include<cmath>#include<cstdlib>#include<iostream>#define ll long longusing namespace std;const int Times = 10;ll multi(ll a, ll b ,ll m){    ll ans = 0;    a%=m;    while(b){        if(b&1)ans = (ans+a)%m;        a = (a+a)%m;        b >>= 1;    }    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);        a = multi(a,a,m);        b>>=1;    }    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("%lld",&n);        if(Miller_Rabin(n))puts("Yes");        else puts("No");    }    return 0;}

例题：[http://www.51nod.com/onlineJudge/questionCode.html#!problemId=1186

java写大数：

import java.math.BigInteger;import java.util.Random;import java.util.Scanner;/** * Created by huchi on 17-7-18. */public class test {    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 < k; 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 in = new Scanner(System.in);        while(in.hasNextBigInteger()){            BigInteger n = in.nextBigInteger();            if(Miller_Rabin(n))System.out.println("Yes");            else System.out.println("No");        }    }}

扩展欧几里得

二元一次方程

ax1+by1=gcd(a,b)

;

a>b时，b==0时，x=1,y=0,　否则设bx2+(a mod b)y2=gcd(b,a mod b)=gcd(a,b)

又有ax1+by1=bx2+(a−[ab]b)y2

则有x1=y2,y1=x2−[ab]y2,　用递推实现．

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

代表的是ax+by=gcd(a,b)的解x,y,和d为gcd(a,b);，得到的x,y中a>b则x为正．在gcd(a,b)==1的时候，可以领y=(y%a+a)％a把y变成正的，再根据等式把x变成负的就行．

例题：https://vjudge.net/problem/POJ-2142

#include<cstdio>#include<iostream>#include<cmath>using namespace std;int gcd(int a, int b){    return b?gcd(b,a%b):a;}void e_gcd(int a, int b, int &d, int &x, int &y){    if(!b){        d=a,x=1,y=0;        return;    }else{        e_gcd(b,a%b,d,y,x);        y -= (a/b)*x;    }}int main(){    int a,b,c;    while(~scanf("%d%d%d",&a,&b,&c) && (a+b+c)){        int x, y, d, tx, ty;        int gd = gcd(a,b);        a /= gd; b/= gd; c/= gd;        e_gcd(a,b,d,x,y);        //获取x非负的时候的情况        tx = x*c;        tx = (tx%b + b)%b;        ty = (c-tx*a)/b;        if(ty < 0)ty = -ty;        //获取y非负的情况        y = y*c;        y = (y%a + a)%a;        x = (c-b*y)/a;        if(x<0)x =-x;        if(x+y>tx+ty){            x=tx;            y=ty;        }        printf("%d %d\n",x,y);    }    return 0;}

https://vjudge.net/problem/HIT-2815

题目没读清，特判就写错了，扩展欧几里得和上面一道题一样的用．

#include<cstdio>#include<iostream>#define ll long longusing namespace std;ll gcd(ll a, ll b){    return b?gcd(b, a%b):a;}void e_gcd(ll a, ll b, ll &d, ll &x, ll &y){    if(!b){        d = a; x = 1; y = 0;        return ;    }else{        e_gcd(b,a%b,d,y,x);        y -= (a/b)*x;    }}int main(){    ll t,a,b,d,x,y;scanf("%lld",&t);    while(t--){        scanf("%lld%lld",&a,&b);        if(a*b==0){            if(a==1 || b==1){                puts("1");            }else{                puts("-1");            }        }else if(a==1 || b==1){            if(max(a,b)==2)                puts("1");            else puts("2");        }        else if(gcd(a,b)!=1){            puts("-1");        }else{            e_gcd(a,b,d,x,y);            ll tx,ty;            tx = (x%b+b)%b;            ty = (1-a*x)/b;            if(ty<0)ty = -ty;            y = (y%a+a)%a;            x = (1-b*y)/a;            if(x<0)x=-x;            if(x+y>tx+ty){                x = tx;                y = ty;            }            printf("%lld\n",x+y-1);        }    }    return 0;}

毕达哥拉斯三元组的解

x2+y2=z2，只考虑本原解．不妨假设x为偶数,y,z为奇数.这样就必有x=2mn,y=m2−n2,z=m2+n2．

例题：https://vjudge.net/problem/POJ-1305

输出两个数据，一个是z不大于n的本原解的个数，和本原解中没有用到的数的个数．

#include<cstdio>#include<cstring>#include<iostream>#define maxn 1000111using namespace std;int biao[maxn];bool flag[maxn];int gcd(int a, int b){    return b?gcd(b,a%b):a;}int main(){    int n;    memset(flag,0,sizeof flag);    for(int i = 2; i*i < maxn; i++){        for(int j = 1; j < i; j++){            if(gcd(i,j)!=1)continue;            if((i&1)&&(j&1))continue;            if(i*i+j*j<maxn)biao[i*i+j*j]++;        }    }    for(int i = 1; i < maxn; i++){        biao[i]+=biao[i-1];    }    while(~scanf("%d",&n)){    memset(flag,0,sizeof flag);        int ans = 0;        for(int i = 2; i*i < maxn; i++){            for(int j = 1; j < i; j++){                if(gcd(i,j)!=1)continue;                if((i&1)&&(j&1))continue;                for(int k = 1; k*(i*i+j*j)<=n; k++){                    flag[2*i*j*k]=flag[k*(i*i-j*j)]=flag[k*(i*i+j*j)]=1;                }            }        }        for(int i = 1; i <= n; i++)            if(flag[i])ans++;        printf("%d %d\n",biao[n],n-ans);    }    return 0;}

例题二:

利用毕达哥拉斯三元组＋欧拉函数＋容斥原理

https://vjudge.net/problem/HDU-3939???

寻找满足的i,j满足gcd(i,j)=1,dfs()容斥原理来找有多少互斥的数．

这里需要判断i<＝j和i>j的情况．其中,i<=j时，i为偶数,则直接加phi[i]，i为奇数时,由于j必须为偶数,不能用phi[i]，所以用容斥原理，一共i>>1个数(晒去了２的倍数)．大于时同理．

#include<cstdio>#include<iostream>#include<cstring>#include<cmath>#define ll long long#define maxn 1001000using namespace std;bool judge[maxn];int prime[maxn/10];int phi[maxn];int check[40];int tot,num;ll ans;void init(){    tot = 0;    memset(judge, 0, sizeof judge);    for(int i = 2; i < maxn; i++){        if(!judge[i]){            for(int j = i+i; j < maxn; j+=i){                judge[j]=1;            }            prime[tot++]=i;        }    }    for(int i = 1; i < maxn; i++)phi[i]=i;    for(int i = 2; i < maxn; i+=2)phi[i]>>=1;    for(int i = 3; i < maxn; i+=2)if(phi[i]==i){        for(int j = i; j < maxn; j+=i){            phi[j] = phi[j] - phi[j]/i;        }    }}void prime_check(int n){    num = 0;    if(!judge[n]){        check[num++]=n;        return ;    }    for(int i = 0; i < tot && n > 1; i++){        if(n%prime[i]==0){            check[num++]=prime[i];            while(n%prime[i]==0) n/=prime[i];            if(n>1 && !judge[n]){                check[num++]=n;                return ;            }        }    }}//容斥原理计算无关的数fun(j)=j-(j/c1+j/c2+j/c3)+(j/(c1*c2)+j/(c1*c3)+j/(c2*c3))-(j/(c1*c2*c3))。void dfs(int k, int r, int s, int n){    if(k==num){        if(r&1) ans -= n/s;        else ans+= n/s;        return;    }    dfs(k+1, r, s, n);    dfs(k+1, r+1, s*check[k], n);}ll L;int main(){    init();    int t;scanf("%d",&t);    while(t--){        ans = 0;        scanf("%lld",&L);        int m = (int)sqrt(1.0*L);        for(int i = m; i > 0; i--){            int p = (int)sqrt(L - (ll)i*i);            //需要寻找gcd(i,j)=1的个数            if(i<=p){                if(i&1){                    prime_check(i);                    dfs(0,0,1,i>>1);                }else ans+=phi[i];            }else{                // i>p时，j所取的范围在[1,p]                prime_check(i);                if(i&1)dfs(0,0,1,p>>1);                else dfs(0,0,1,p);            }        }        printf("%lld\n",ans);    }}

欧拉函数与欧拉定理

定理一：

m与n互素，ϕ(mn)=ϕ(m)ϕ(n)

, 所以n为奇数时，我们有ϕ(2n)=ϕ(n)

定理二: p是素数，有ϕ(pa)=pa−pa−1,容斥原理．

定理三: n=pα11pα22...pαkk时，ϕ(n)=n(1−1p1)(1−1p2)...(1−1pk)　定理１＋３

定理四: ∑d|nϕ(d)=n 莫比乌斯反演．

定理五: (a,m)=1,则x≡baϕ(m)−1是同于方程ax≡b(mod m)

定理六: n>2时，ϕ(n)为偶数

定理三求解欧拉函数：

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

容斥原理，筛法去减不满足的．

#define maxn 100000int p[maxn];void phi_table(){    memset(p,0,sizeof p);    for(int i = 1; i < maxn; i++)p[i]=i;    for(int i = 2; i < maxn; i+=2)p[i]>>=1;    for(int i = 3; i < maxn; i+=2){        if(p[i]==i){            for(int j = i; j < maxn; j+=i)                p[j]=p[j]-p[j]/i;        }    }}

二次筛法

这道题有类似的在kuangbin的数论基础专题里有出现过．

例题：https://vjudge.net/problem/POJ-2689

注意1　2这个样例．．．还有晒的时候不能把质数也晒了．．．．

#include<cstdio>#include<iostream>#include<cstring>#define maxn 50000int prime[maxn];bool judge2[21*maxn];bool judge1[maxn];using namespace std;int tot=0;void init(){    memset(judge1,0,sizeof judge1);    for(int i = 2; i<maxn; i++){        if(!judge1[i]){            prime[tot++] = i;            for(int j = i+i; j < maxn; j+=i){                judge1[j]=true;            }        }    }}int main(){    init();    unsigned int u,v;    while(~scanf("%d%d",&u,&v)){        if(u==1)u=2;        memset(judge2,0,sizeof judge2);        for(int i = 0; i < tot; i++){            if(prime[i]>v)break;            unsigned int j=u/prime[i]*prime[i];            if(u%prime[i]==0 && u!=prime[i])judge2[0]=1;            while(j<=u)j+=prime[i];            if(j==prime[i])j+=j;            for(;j<=v;j+=prime[i]){                judge2[j-u]=1;            }        }        int st,stlen=maxn,ed,edlen=0;        int cot=0,len=0;        for(int i = 0; i <= v-u; i++){            if(!judge2[i]){                if(len &&len<stlen){                    st = i-len;                    stlen = len;                }                if(len && len>edlen){                    ed = i-len;                    edlen = len;                }                cot++;                len = 0;            }            len ++ ;        }        if(cot!=2 && st==ed){            ed = st+stlen;        }        if(cot>1)printf("%d,%d are closest, %d,%d are most distant.\n",u+st,u+st+stlen,u+ed,u+ed+edlen);        else puts("There are no adjacent primes.");    }    return 0;}

逆元

正整数a,m，如果有ax≡1(mod m),这个同于方程中x的最小正整数解叫做a模m的逆元．

逆元一般用扩展欧几里得算法来求．m为素数时，由费马小定理得到逆元为am−2modm

逆元常见问题

已知b|a，求解ans=abmod m

常见解法:扩展欧几里得．m是素数时，直接用费马小定理．但他们都需要a与m互素．

除法取模求逆元的通用方法： ans=a mod(mb)b

ab=km+x(x<m),a=kbm+bx,a mod(bm)=bx,x=a mod(bm)b

例题：https://vjudge.net/problem/POJ-1845

求AB的所有因子和对9901取余的结果．

考虑

A=pα11pα22...pαkk,AB=pα1B1pα2B2...pαkBk

所有的因子和：

sum=∏1≤i≤k∑0≤j≤αiBpji

可以用二分求等比数列和，也可以用等比数列公式(需要用逆元)直接求.分别写在了Method_one和Method_two中．

#include<cstdio>#include<iostream>#include<cstring>#include<cmath>#define ll long long#define maxn 10000using namespace std;ll mod =9901;bool judge[maxn];ll check[maxn][2];int prime[maxn/2];int tot,num;ll a,b;void init(){    tot=0;    memset(judge, 0, sizeof judge);    for(int i = 2; i < maxn; i++){        if(!judge[i]){            for(int j = i+i; j < maxn; j+=i){                judge[j]=1;            }            prime[tot++]=i;        }    }}ll multi(ll a, ll b, ll m){    ll ans = 0;    a %= m;    while(b){        if(b&1)ans = (ans+a)%m;        a = (a+a)%m;        b >>= 1;    }    return ans;}ll largequickmod(ll p, ll k, ll m){    ll ans = 1;    while(k){        if(k&1)ans = multi(ans,p,m);        p = multi(p,p,m);        k>>=1;    }    return ans;}ll quickmod(ll p, ll k){    ll ans=1;    p%=mod;    while(k){        if(k&1)ans = ans*p%mod;        p = (p*p)%mod;        k>>=1;    }    return ans;}ll S(ll p, ll k){    if(p==0)return 0;    if(k==0)return 1;    ll temp = quickmod(p,(k>>1))%mod;    if(k&1){        return S(p,k>>1)%mod*(1+temp*p%mod)%mod;    }else{        return (S(p,(k>>1)-1)%mod*(1+temp*p%mod)%mod+temp)%mod;    }}void Method_one(){    ll ans = 1;    for(int i = 0; i < num; i++){        ans = (ans * S(check[i][0], b*check[i][1]))%mod;    }    if(a==0)ans=0;    printf("%lld\n",ans);}void Method_two(){    ll ans = 1;    for(int i = 0; i < num; i++){        ll M = mod*(check[i][0]-1);        ans *= (largequickmod(check[i][0], check[i][1]*b+1, M)+M-1)/(check[i][0]-1);        ans %= mod;    }    printf("%lld\n",ans);}int main(){    init();    while(~scanf("%lld%lld",&a,&b)){        ll temp = a;        num=0;        memset(check, 0, sizeof check);        for(int i = 0; i < tot; i++){            if(temp<=1)break;            if(temp%prime[i]==0){check[num++][0]=prime[i];}            while(temp%prime[i]==0){                temp/=prime[i];                check[num-1][1]++;            }        }        if(temp>1){            check[num][0]=temp;            check[num++][1]=1;        }        //Method_one();        Method_two();    }    return 0;}

1→p模p的所有逆元，这里p为奇质数．用快速幂的复杂度Ｏ(plog(p)),用递推式可以优化到Ｏ(p)：

inv[i]=(M−M/i)∗inv[M%i]%M推导：t=M/i,k=M%i,有t∗i+k≡0(modM),−t∗i≡k(modM)两边同时除i∗k记作:−t∗inv[k]≡inv[i](modM)即　inv[i]=(M−M/i)∗inv[M%i]%M

初始化inv[1]=1,就可以通过递推法求出1→p模奇素数的所有逆元．

http://www.lydsy.com/JudgeOnline/problem.php?id=2186

例题：1-N!中与M!互质的数的个数，M<=N;M!|N!, ．当n是m的倍数时，1-n中与m互素的个数为nmϕ(m)

题目中　N!M!ϕ(M!)=N!M!M!(1−1p1)(1−1p2)...(1−1pk).　　

时间卡的特别紧．bool开的就超时了．inv[i]算逆元．multi[maxn]算阶乘取模，要求p与m互质，否则需要用除法取模．

#include<cstdio>#include<bitset>#define ll long long#define maxn 10000005using namespace std;ll mod;bitset<maxn> judge;void init(){    judge.set();    for(int i = 4; i < maxn; i+=2)judge[i]=0;    for(int i = 3; i < maxn; i+=2){        if(judge[i]){            for(int j = i+i; j < maxn; j+=i){                judge[j]=0;            }        }    }}ll multi[maxn];int inv[maxn];ll ans[maxn];int main(){    init();    int T,R;    scanf("%d%d",&T,&R);    mod = R;    multi[0]=1;    for(int i = 1; i < maxn; i++){        multi[i]=multi[i-1]*i%mod;    }    inv[1]=1;    for(int i = 2; i < mod && i < maxn; i++){        inv[i]=(mod-mod/i)*inv[mod%i]%mod;    }    ans[1]=1;    for(int i = 2; i < maxn; i++){        if(judge[i]){            ans[i] = ans[i-1]*(i-1)%mod;            ans[i] = ans[i]*inv[i%mod]%mod;        }else{            ans[i] = ans[i-1];        }    }    int n,m;    while(T--){        scanf("%d%d",&n,&m);        ll anss = 1LL*multi[n]*ans[m]%mod;        printf("%lld\n",anss);    }    return 0;}

欧拉筛

void getphi()    {       int i,j;       phi[1]=1;       for(i=2;i<=N;i++)//相当于分解质因式的逆过程       {           if(!mark[i])               {                 prime[++tot]=i;//筛素数的时候首先会判断i是否是素数。                 phi[i]=i-1;//当 i 是素数时 phi[i]=i-1                 }           for(j=1;j<=tot;j++)           {              if(i*prime[j]>N)  break;              mark[i*prime[j]]=1;//确定i*prime[j]不是素数              if(i%prime[j]==0)//接着我们会看prime[j]是否是i的约数              {                 phi[i*prime[j]]=phi[i]*prime[j];break;              }              else  phi[i*prime[j]]=phi[i]*(prime[j]-1);//其实这里prime[j]-1就是phi[prime[j]]，利用了欧拉函数的积性           }       }    }