红书数论模板

long long pow_mod(long long a,long long i,long long n)
{
    if(i==0)
        return 1%n;
    long long temp=pow_mod(a,i>>1,n)%n;
    temp=temp*temp%n;
    if(i&1)
        temp=(long long)(temp%n)*(a%n)%n;
    return temp;
}

long long extend_gcd(long long a,long long b,long long &x,long long &y)
{
    if(b==0)
    {
        x=1;
        y=0;
        return a;
    }
    else 
    {
        long long r=extend_gcd(b,a%b,y,x);
        y-=x*(a/b);
        return r;
    }
}

vector<long long> a;
bool g_test(long long g,long long p)
{
    for(long long i=0;i<a.size();i++)
        if(pow_mod(g,(p-1)/a[i],p)==1)
            return 0;
    return 1;
}

long long primitive_root(long long p) //求mod p的原根
{
    long long tmp=p-1;
    for(long long i=2;i<=tmp/i;i++)
        if(tmp%i==0)
        {
            a.push_back(i);
            while(tmp%i==0)
                tmp/=i;
        }
    if(tmp!=1)
        a.push_back(tmp);
    long long g=1;
    while(true)
    {
        if(g_test(g,p))
            return g;
        ++g;
    }
}

long long discrete_log(long long x,long long n,long long m) //求x^y=n(mod m)的解 m是素数
{
    map<long long,int> rec;
    long long s=(long long)(sqrt((double)m));
    for(;(long long)s*s<=m;)
        s++;
    long long cur=1;
    for(long long i=0;i<s;i++)
    {
        rec[cur]=i;
        cur=cur*x%m;
    }
    long long mul=cur;
    cur=1;
    for(long long i=0;i<s;i++)
    {
        long long more=(long long)n*pow_mod(cur,m-2,m)%m;
        if(rec.count(more))
        {
            return i*s+rec[more];
        }
        cur=cur*mul%m;
    }
    return -1;
}

vector <long long> residue(long long p,long long N,long long a) //求x^n=a(mod p)所有解 p是素数
{
    long long g=primitive_root(p);
    long long m=discrete_log(g,a,p);
    vector<long long> ret;
    if(a==0)
    {
        ret.push_back(0);
        return ret;
    }
    if(m==-1)
    {
        return ret;
    }
    long long A=N,B=p-1,C=m,x,y;
    long long d=extend_gcd(A,B,x,y);
    if(C%d!=0)
        return ret;
    x=x*(C/d)%B;
    long long delta=B/d;
    for(long long i=0;i<d;++i)
    {
        x=((x+delta)%B+B)%B;
        //    printf("%I64d %I64d %I64d\n",g,x,p);
        //    printf("---%I64d\n",pow_mod(g,x,p));
        ret.push_back(pow_mod(g,x,p));
    }
    sort(ret.begin(),ret.end());
    ret.erase(unique(ret.begin(),ret.end()),ret.end());
    return  ret;
}

long long CRT(long long a[],long long m[],int n) //求方程组x=ai(mod mi) mi两两互质
{
    long long M=1;
    for(int i=0;i<n;i++)
        M*=m[i];
    long long ret=0;
    for(int i=0;i<n;i++)
    {
        long long x,y;
        long long tm=M/m[i];
        extend_gcd(tm,m[i],x,y);
        ret=(ret+tm*x*a[i])%M;    
    }
    return (ret+M)%M;
}

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

void line_mod(long long a,long long b,long long n) //求解ax=b(mod n)
{
    long long x,y;
    long long d=extend_gcd(a,n,x,y);
    vector<long long> ans;
    ans.clear();
    if(b%d==0)
    {
        x%=n;
        x+=n;
        x%=n;
        ans.push_back(x*(b/d)%(n/d));
        for(long long i=1;i<d;i++)
        {
            
            ans.push_back((ans[0]+i*n/d)%n);
        }
    }
}