JZPKIL

来源：互联网发布：鬼武者系列知乎编辑：程序博客网时间：2024/06/05 10:50
终于A掉JZPKIL了。。感觉以后可以当数论题模板了(Bernoulli数、Pollard's rho、Miller Rabin、*快速乘、快速幂、逆元、组合数)
推导见 ydc blog
最神奇的是gyz的O(1)快速乘
证明见:http://tieba.baidu.com/p/2175344676
刷到R1了好开心
#include <iostream>#include <cstdio>#include <algorithm>#include <set>#include <cmath>#include <cstdlib>#define mod 1000000007#define abs(x) ((x)>0?(x):(-(x)))using namespace std;typedef long long LL;LL p[20000],prime[20000];int num[20000],tot;LL P[20000],Q[20000],B[20000];LL b[3010];int t[3010][3010],c[3010][3010];int cnt;LL mul(LL x,LL y,LL z){return (x*y-(LL)(x/(long double)z*y+1e-3)*z+z)%z;}inline LL power(LL p,LL n,LL mo){LL ans=1;for(;n;n>>=1,p=mul(p,p,mo))if(n&1)ans=mul(ans,p,mo);return ans;}LL power1(LL a,LL b,LL mo){    if (b==0) return 1;    if (b==1) return a%mo;    LL m=power1(a,b/2,mo);    m=(m*m)%mo;    if (b&1) m=(m*a)%mo;    return m;}inline bool Miller_Rabin(int p,LL n){int s=0;LL d=n-1;while(!(d&1))++s,d>>=1;LL w=power(p,d,n);if(w==1)return true;for(int i=0;i<=s-1;++i){if(w==n-1)return true;w=mul(w,w,n);}return false;}inline bool isPrime(LL n){static int nd[]={2,3,5,7,11,13,17,19,23};if(n==1)return false;if(find(nd,nd+9,n)!=nd+9)return true;if(!(n&1))return false;for(int i=0;i<9;++i)if(!Miller_Rabin(nd[i],n))return false;return true;}LL gcd(LL a,LL b){    if (b==0) return a;    else return gcd(b,a%b);}void frac(LL n){if(isPrime(n)){p[++cnt]=n;return ;}int c=16831;while(1){LL x1,x2;int i=1,k=2;x1=x2=1;while(1){x1=mul(x1,x1,n)+c;x1%=n;LL d=gcd(abs(x1-x2),n);if(d!=1&&d!=n){frac(d),frac(n/d);return ;}if(x1==x2)break;if(++i==k)k<<=1,x2=x1;}c--;}}LL exgcd(LL a,LL b,LL &x,LL &y){    if (b==0)    {        x=1,y=0;        return a;    }    LL r=exgcd(b,a%b,x,y);    LL t=x;    x=y;    y=t-a/b*y;    return r;}LL anti(LL n){    LL x,y;    exgcd(n,mod,x,y);    x=x%mod;    if (x<0) x+=mod;    return x;}void init(){    int n=3000;    c[0][0]=1;    for (int i=1;i<=n+1;++i)    {        c[i][0]=1;        for (int j=1;j<=i;++j)        {            c[i][j]=c[i-1][j]+c[i-1][j-1];            if (c[i][j]>=mod) c[i][j]-=mod;        }    }    for (int i=0;i<=n;++i)    {        b[i]=i+1;        for (int j=0;j<i;++j)        {            b[i]=(b[i]-b[j]*c[i+1][j])%mod;            if (b[i]<0) b[i]+=mod;        }        b[i]=(b[i]*anti(c[i+1][i]))%mod;    }    for (int i=0;i<=n;++i)    {        LL a=anti(i+1);        for (int j=0;j<=i;++j)            t[i][j]=(a*c[i+1][j]%mod*b[j])%mod;    }}void add(LL n,LL m){    prime[++tot]=m;    num[tot]=0;    while (n%m==0)    {        n/=m;        num[tot]++;    }}LL kil[10000],jzp[10000],m[10000],super[100][100];LL calc(LL n,LL k,int x,int y){    LL res=1LL;    for (int i=1;i<=tot;++i)    {        LL d=1LL,tmp=0;        LL kk=1;        for (int j=0;j<=num[i];++j)        {            LL a=super[i][num[i]-j];            if (num[i]-j>=1) a-=(kil[i]*super[i][num[i]-j-1])%mod;            a%=mod;            if (a<0) a+=mod;            tmp+=kk*a;            kk*=jzp[i];            kk%=mod;            tmp%=mod;            d*=prime[i];        }        res=res*tmp;        res%=mod;    }    return res;}LL work(LL n,int x,int y){    cnt=0,tot=0;    frac(n);    sort(p+1,p+cnt+1);    add(n,p[1]);    for (int i=2;i<=cnt;++i)        if (p[i]!=p[i-1]) add(n,p[i]);    LL ans=0;for (int i=1;i<=tot;++i){        kil[i]=power1(prime[i]%mod,y,mod);        jzp[i]=power1(prime[i]%mod,x,mod);        m[i]=power1(prime[i],num[i],n+1);LL d=1;for (int j=0;j<=num[i];++j,d*=prime[i]) super[i][j]=d%mod;}    for (int i=y;i>=0;--i)    {        ans+=(t[y][i]*calc(n,i,x,y))%mod;        ans%=mod;        if (ans<0) ans+=mod;for (int i=1;i<=tot;++i){LL d=1;for (int j=0;j<=num[i];++j,d*=prime[i]) super[i][j]=(super[i][j]*(d%mod))%mod;}    }n%=mod;    ans=(ans*power1(n,y,mod))%mod;    return ans;}LL kill(LL n,int x,int y){    LL ans=0,p=n%mod,d=p;    for(int i=y;i>=0;--i)    {        ans=ans+t[y][i]*d;        if(ans>=mod)            ans%=mod;d=d*p;if(d>=mod)d%=mod;    }    return ans*power1(p,y,mod)%mod;}int main(){    int T;LL n;    int x,y;    init();    cin>>T;    while (T--)    {        cin>>n>>x>>y;//frac(n);        if (x==y) cout<<kill(n,x,y)<<endl;        else cout<<work(n,x,y)<<endl;    }    return 0;}
0 0