Acdreamer博客数论学习(1)

Acdreamer博客数论学习

Day One

等比数列求和．

求解$S_n = (a+a^2+…+a^n)modM$

采用二分．n%2==0时 $ S_n = (1+a^{(\frac{n}{2})} )S_{\frac{n}{2}} $ ,n%2==1时，$S_n = (1+a^{\frac{n+1}{2}})S_{\frac{n-1}{2}} + a^{\frac{n+1}{2}}$

#include<cstdio>
#include<iostream>
#define ll long long
using 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;
}

素数定理

$$
\pi(n)为小于等于n的素数的个数，有lim_{n→∞}\frac{\pi(n)}{n/lnn}=1
$$

定理$gcd(a^n-1,a^m-1)=a^{gcd(n,m)}-1$

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

#include<cstdio>
#include<iostream>
#define ll long long
using 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(a^m-b^m,a^n-b^n)=a^{gcd(m,n)}-b^{gcd(m,n)}$

定理$G=gcd(C_n^1,C_n^2,C_n^3,…,Cn^{n-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 1000100
using 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;
}

定理$F_n 为斐波那契额数列 ，gcd(F_m,F_n)=F_{gcd(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=$\frac{(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;
}

定理$\sum{gcd(i,N)=\sum{d\phi(\frac{N}{d})}}$

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

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

#include<cstdio>
#include<cstring>
#include<iostream>
#define ll long long
#define maxn 50000
using 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(C_n^0,C_n^1,C_n^2,…,C_N^{n-1},C_n^N)=lcm(1,2,…,n+1)
$$

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

结论：
$$
如果p是素数，那么C_p^1,C_p^2,C_p^3,…C_p^{p-1},均能被p整除\
如果p是素数，那么有(x+y+..+w)^p≡x^p+y^p+…+w^p(mod p)\
$$

Day Two

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

$$
梅森素数M_p =2^p-1，p为素数．，定义序列{s_i}，所有的i ≥ 0, \
si= \begin{cases} 4, i=0 \ s_{i-1}^{2}-2,i>0 \end{cases}
M_p是素数当且仅当s_{p-2} \equiv 0(modM_p)
$$

#include<cstdio>
#include<iostream>
#define ll long long
using 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，有$$a^p ≡ a(mod p)$$（同余）。反过来，满足$$a^p ≡ a(mod p)$$，p也几乎一定是素数。

　　伪素数：如果n是一个正整数，如果存在和n互素的正整数a满足 $a^{n-1} ≡ 1(mod n)$，我们说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次)，计算是否每次都有$$b^{n-1} ≡ 1(mod n)$$，若每次都成立则n是素数，否则为合数。　

​ 二次探测定理 ：p为素数，０<x<p,则方程$x^2 ≡ 1(mod p)的解为１或者p-1$

#include<cstdio>
#include<cmath>
#include<cstdlib>
#include<iostream>
#define ll long long
using 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");
        }
    }
}

扩展欧几里得

二元一次方程$$ax_1+by_1=gcd(a,b)$$;

a>b时，b==0时，x=1,y=0,　否则设$bx_2+(a\ mod\ b)y_2 = gcd(b, a\ mod \ b)=gcd(a,b)$

又有$ax_1 +by_1=bx_2+(a-[\frac{a}{b}]b)y_2$

则有$x_1=y_2,y_1=x_2-[\frac{a}{b}]y_2$,　用递推实现．

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 long
using 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;
}

毕达哥拉斯三元组的解

$x^2+y^2=z^2$，只考虑本原解．不妨假设x为偶数,y,z为奇数.这样就必有$x=2mn,y=m^2-n^2,z=m^2+n^2$．

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

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

#include<cstdio>
#include<cstring>
#include<iostream>
#define maxn 1000111
using 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 1001000
using 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互素，\phi(mn)=\phi(m)\phi(n)$$, 所以n为奇数时，我们有$\phi(2n)=\phi(n)$

定理二: p是素数，有$\phi(p^a)=p^a-p^{a-1}$,容斥原理．

定理三: $n=p_1^{\alpha_1}p_2^{\alpha_2}…p_k^{\alpha_k}时，\phi(n)=n(1-\frac{1}{p_1})(1-\frac{1}{p_2})…(1-\frac{1}{p_k})$　定理１＋３

定理四: $\sum_{d|n}{\phi(d)}= n$ 莫比乌斯反演．

定理五: (a,m)=1,则$x \equiv ba^{\phi(m)-1} 是同于方程 ax \equiv b(mod \ m)$

定理六: n>2时，$\phi(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 100000
int 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 50000
int 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 \equiv 1(mod \ m)$,这个同于方程中x的最小正整数解叫做a模m的逆元．

逆元一般用扩展欧几里得算法来求．m为素数时，由费马小定理得到逆元为$a^{m-2}modm$

逆元常见问题

已知b|a，求解$ans = \frac{a}{b}mod \ m$

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

除法取模求逆元的通用方法： $ans =\frac{ a \ mod(mb)}{b}$
$$
\frac{a}{b}=km+x(x<m),a=kbm+bx,a\ mod(bm)=bx,x = \frac{a\ mod(bm)}{b}
$$
例题：https://vjudge.net/problem/POJ-1845

求$A^B$的所有因子和对9901取余的结果．

考虑$$A=p_1^{\alpha_1}p_2^{\alpha_2}…p_k^{\alpha_k} ,A^B=p_1^{\alpha_1 B}p_2^{\alpha_2B}…p_k^{\alpha_kB}$$

所有的因子和：$$sum=\prod_{1\le i \le k} \sum_{0 \le j \le \alpha_iB} p_i^{j}$$

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

#include<cstdio>
#include<iostream>
#include<cstring>
#include<cmath>
#define ll long long
#define maxn 10000
using 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,有ti+k \equiv 0 (modM), -ti \equiv k(modM)\ 两边同时除ik记作:-tinv[k] \equiv invi\ 即　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互素的个数为$\frac{n}{m}\phi(m)$

题目中　$\frac{N!}{M!}\phi(M!) = \frac{N!}{M!}M!(1-\frac{1}{p_1})(1-\frac{1}{p_2})…(1-\frac{1}{p_k})$.　　

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

#include<cstdio>
#include<bitset>
#define ll long long
#define maxn 10000005
using 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]]，利用了欧拉函数的积性    
       }    
   }    
}