Pollard_Rho质因数分解 似乎有点卡常？
#include"bits/stdc++.h" using namespace std; const int N=3005,P=1000000007; typedef long long LL; int inv(LL a,LL t=P-2){ int r=1;t%=P-1,a%=P; if(t<0)t+=P-1; while(t){ if(t&1)r=1LL*r*a%P; a=1LL*a*a%P;t>>=1; } return r; } int b[N],f[N][N],val[N][N],C[N][N]; void pre(){ b[0]=C[0][0]=1; for(int i=1;i<=3001;i++)for(int j=0;j<=i;j++)C[i][j]=(C[i-1][j-1]+C[i-1][j])%P; for(int d=1;d<=3000;d++){ for(int k=0;k<d;k++)b[d]=(b[d]-1LL*b[k]*C[d+1][k]%P+P)%P; b[d]=1LL*b[d]*inv(d+1)%P; } for(int d=0;d<=3000;d++){ int inv_d=inv(d+1); for(int k=0;k<=d;k++)val[d][d+1-k]=1LL*b[k]*C[d+1][k]%P*inv_d%P; (val[d][d]+=1)%=P; } } LL mult(LL x,LL y,LL z){ //__int128 r=x;r=r*y; //return (LL)(r%z); return (x*y-(LL)(((long double)x*y+0.5)/(long double)z)*z+z)%z; } //LL mult(LL a,LL b,LL c){return (LL)((long double)a*b-(LL)(((long double)a*b+0.5)/c)*(long double)c);} LL ksm(LL a,LL t,LL P){ LL r=1; while(t){ if(t&1)r=mult(r,a,P); a=mult(a,a,P);t>>=1; } return r; } const int JudgeP[] = {2,3,5,7,11,13,17,19,23}; inline int irand(){return (rand()<<16)^rand();} bool Miller_Rabin(LL n){ if(n==2||n==3||n==5||n==7||n==11||n==13||n==17||n==19||n==23)return true; if((n&1)==0)return false; LL r=n-1;int k=0; while((r&1)==0)r>>=1,k++; for(int i=0;i<9;i++){ LL x=ksm(JudgeP[i],r,n),y; for(int i=0;i<k;i++){ y=mult(x,x,n); if(y==1&&x!=1&&x!=n-1)return false; x=y; } if(y!=1)return false; } return true; } int T,numc; LL num[70]; LL gcd(LL x,LL y){return y?gcd(y,x%y):x;} #define Func(x) (mult(x,x,n)+1) LL getnum(LL n){ LL x=irand()%n+1,y=Func(x); while(x!=y){ LL k=gcd(n,(y-x+n)%n); if(k!=1&&k!=n)return k; x=Func(x),y=Func(Func(y)); } return -1; } void Pollard_Rho(LL n){ while(n%2==0){ num[++numc]=2; n>>=1; } if(n==1)return; if(Miller_Rabin(n)){ num[++numc]=n; return; } LL x; do x=getnum(n);while(x==-1); Pollard_Rho(x); Pollard_Rho(n/x); } LL n; void work(){ srand(19260817); scanf("%d",&T); int x,y; while(T--){ scanf("%lld%d%d",&n,&x,&y); numc=0; if(n==1){puts("1");continue;} Pollard_Rho(n); sort(num+1,num+numc+1); int l=unique(num+1,num+numc+1)-num-1; LL r=0,w,m,s,k,g; for(int i=1;i<=y+1;i++){ w=1,m=n; for(int j=1;j<=l;j++){ s=1,k=0; while(m%num[j]==0)m/=num[j],s*=num[j],k++; g=0; for(int t=0;t<=k;t++)g=(g+inv(num[j],1LL*t*(x-i)+1LL*k*i))%P; for(int t=0;t<k;t++)g=(g-inv(num[j],1LL*t*(x-i)+1LL*(k-1)*i+y))%P; w=1LL*w*(g+P)%P; } r=((r+1LL*val[y][i]*w)%P+P)%P; } r=1LL*r*inv(n,y)%P; printf("%lld\n",r); } } int main(){ pre(); work(); return 0; }
