数论算法模板（不定期更新）

/**********/

#include<iostream>
#include<cstdio>
#include<cmath>
#include<algorithm>
#include<cstring>
#include<string>
#include<cstdlib>
#include<vector>
#include<stack>
#include<map>
using namespace std;
typedef long long ll;

/***********GCD**********************/

ll gcd(ll a,ll b)
{
    if(b==0)
        return a;
    return gcd(b,a%b);
}

/*************快速乘***************/

ll mult_mod(ll a,ll b,ll mod)
{
    a%=mod; b%=mod;
    ll res=0;
    while(b)
    {
        if(b&1)
        {
            res+=a;
            res%=mod;
        }
        a<<=1;
        if(a>=mod) a%=mod;
        b>>=1;
    }
    return res;
}

/*************快速幂**************/

ll pow_mod(ll x,ll n,ll mod)
{
    if(n==1) return x%mod;
    x%=mod;
    ll t=x,res=1;
    while(n)
    {
        if(n&1) res=mult_mod(res,t,mod);
        t=mult_mod(t,t,mod);
        n>>=1;
    }
    return res;
}

/************扩展欧几里德****************/

void extend_gcd(ll a,ll b,ll &x,ll &y)
{
    ll d; //d=gcd(a,b)
    if(b==0)
    {
        x=1,y=0;
        d=a;
    }
    else
    {
        d=extend_gcd(b,a%b,y,x);
        ll xx=x,yy=y;
        x=yy;
        y=xx-(a/b)*yy;
    }
}

/************素数筛法******************/

ll l,u,prime[N];
int tot;
int vis[N],ans[10000005];

void isPrime()
{
    tot=0;
    memset(vis,0,sizeof(vis));
    memset(prime,0,sizeof(prime));
    for(ll i=2;i<N;i++)
    {
        if(!vis[i])
        {
            prime[tot++]=i;
            for(ll j=i*i;j<N;j+=i)
                vis[j]=1;
        }
    }
}

/*********素数判定Miller_Rabin***********/

bool test(ll a,ll n) //Miller_Rabin算法的核心
{
    ll x=n-1,t=0,res,last;
    while((x&1)==0)
    {
        x>>=1;
        t++;
    }
    last=pow_mod(a,x,n);

    for(int i=0;i<t;i++)
    {
        res=pow_mod(last,2,n);
        if(res==1&&last!=1&&last!=n-1)
            return true;
        last=res;
    }
    if(res!=1) return true;
    return false;
}

bool millier_rabin(ll n)
{
    if(n==2) return true;
    if(n==1||(n&1)==0) return false;
    for(int i=0;i<times;i++)
    {
        ll a=rand()%(n-1)+1;
        if(test(a,n))
            return false;
    }
    return true;
}

/*********大数因数分解Pollard_rho*************/

ll pollard_rho(ll x,ll c)
{
    ll x0,y,i=1,k=2;
    x0=rand()%x;
    y=x0;
    while(1)
    {
        i++;
        x0=(mult_mod(x0,x0,x)+c)%x;
        ll d=gcd(y-x0,x);
        if(d>1&&d<x) return d;
        if(y==x0) break;
        if(i==k)
        {
            y=x0;
            k+=k;
        }
    }
    return x;
}

void find_fac(ll n,int c)
{
    if(n==1) return;
    if(millier_rabin(n))
    {
        factor[tot++]=n;
        return;
    }
    ll p=n;
    while(p>=n)
        p=pollard_rho(p,c--);
    find_fac(p,c);
    find_fac(n/p,c);
}

/**********欧拉函数**********************/

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

/**********欧拉函数**************/

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

/******素数筛法+欧拉打表**********/

ll e[N+5],p[N+5];
bool vis[N+5];

void init()
{
    memset(e,0,sizeof(e));
    memset(p,0,sizeof(p));
    memset(vis,false,sizeof(vis));
    ll i,j;
    p[0]=1;//记录素数总个数
    p[1]=2;

    for(i=3;i<N;i+=2)
    {
        if(!vis[i])
        {
            p[++p[0]]=i;
            for(j=i*i;j<N;j+=i)
                vis[j]=true;
        }
    }

    e[1]=1;
    for(i=1;i<=p[0];i++)
        e[p[i]]=p[i]-1; //为什么要-1？

    for(i=2;i<N;i++)
    {
        if(!e[i])
        {
            for(j=1;j<=p[0];j++)
            {
                if(i%p[j]==0)
                {
                    if(i/p[j]%p[j])
                        e[i]=e[i/p[j]]*e[p[j]];
                    else
                        e[i]=e[i/p[j]]*p[j];
                    break;
                }
            }
        }
    }
}

/*************欧拉打表*更快版***************/

#define N 1000000

int e[N+5];
void init()
{
    int i,j;
    memset(e,0,sizeof(e));
    for(i=2;i<=N;i++)
    {
        if(!e[i])
        {
            for(j=i;j<=N;j+=i)
            {
                if(!e[j])
                    e[j]=j;
                e[j]=e[j]/i*(i-1);
            }
        }
    }
}