题面:
2627: JZPKIL
Time Limit: 50 Sec Memory Limit: 128 MBSubmit: 347 Solved: 104
[Submit][Status][Discuss]
Description
Input
第一行,询问个数T。
下面T行,每行三个整数,n, x, y。
Output
Sample Input
6 0 0
6 0 1
6 1 0
6 1 1
1000000000 50 50
Sample Output
66
15
126
393442025
HINT
数据规模和约定
30%的数据,x=y
另30%的数据,n<=10^9, x, y<=100
100%的数据,T<=100, 1<=n<=10^18, 0<=x, y<=3000
这题是个50sec卡评测反演好题。
$\quad \sum_{i=1}^{n}(i,n)^{x}[i,n]^{y}$
$=\sum_{i=1}^{n}(i,n)^{x-y}(i\cdot n)^{y}$
$=n^{y}\sum_{i=1}^{n}(i,n)^{x-y}i^{y}$
$=n^{y}\sum_{d|n}d^{x}\sum_{i=1}^{\frac{n}{d}}\sum_{k|i,k|\frac{n}{d}}\mu(k)i^{y}$
$=n^{y}\sum_{d|n}d^{x}\sum_{k|\frac{n}{d}}\mu(k)k^{y}\sum_{i=1}^{\frac{n}{k\cdot d}i^{y}}$
$=n^{y}\sum_{d|n}d^{x}\sum_{k|\frac{n}{d}}\mu(k)k^{y}\sum_{i=0}^{y}\frac{1}{y+1}\dbinom{y+1}{i}B_{i}(\frac{n}{k\cdot d})^{y+1-i}$
$=\sum_{i=0}^{y}t_{i}n^{y}\sum_{d|n}d^{x}\sum_{k|\frac{n}{d}}\mu(k)k^{y}(\frac{n}{k\cdot d})^{y+1-i}$
$t_{i}=\frac{1}{y+1}\dbinom{y+1}{i}B_{i}$ ($B_{i}$为伯努利数)
将$n$分解质因数$n=p_{1}^{k_1}\cdot p_{2}^{k_2}\cdot p_{3}^{k_3}\cdots $
令$f(p_{j}^{k_j},i)=\sum_{d|p_{j}^{k_j}}d^{p_{j}^{k_j}}\sum_{k|\frac{p_{j}^{k_j}}{d}}\mu(k)k^{y}(\frac{p_{j}^{k_j}}{k\cdot d})^{y+1-i}$
$f$为积性函数,可以用每个质因子的答案相乘得到$n$的答案。
但是由于$n$很大,我们只能用$Pollard_rho$分解大数因子。
时间复杂度为$O(y^{2}+T\cdot n^{\frac{1}{4}}+T\cdot y\cdot log\;n)$
![](https://i-blog.csdnimg.cn/blog_migrate/8f900a89c6347c561fdf2122f13be562.gif)
![](https://i-blog.csdnimg.cn/blog_migrate/961ddebeb323a10fe0623af514929fc1.gif)
1 #include <iostream> 2 #include <stdio.h> 3 #include <algorithm> 4 using namespace std; 5 #define LL long long 6 template<typename __> 7 inline void read(__ &s) 8 { 9 char ch=getchar(); 10 for(s=0;ch<'0'||ch>'9';ch=getchar()); 11 for(;ch>='0'&&ch<='9';s=s*10+ch-'0',ch=getchar()); 12 } 13 namespace Pollard_rho 14 { 15 inline LL mul(LL x,LL y,LL mod) 16 { 17 return (x*y-(LL)(x/(long double)mod*y+0.5)*mod+mod)%mod; 18 } 19 inline LL qpow(LL x,LL y,LL mod) 20 { 21 LL ans=1; 22 for(;y;y>>=1,x=mul(x,x,mod)) 23 if(y&1) 24 ans=mul(ans,x,mod); 25 return ans; 26 } 27 inline bool judge(LL n,LL p) 28 { 29 if(n<=p) 30 return 1; 31 if(n%p==0) 32 return 0; 33 LL x=n-1; 34 int s=-1; 35 while(!(x&1)) 36 x>>=1,++s; 37 x=qpow(p,x,n); 38 if(x==1||x==n-1) 39 return 1; 40 while(s--) 41 { 42 x=mul(x,x,n); 43 if(x==1) 44 return 0; 45 if(x==n-1) 46 return 1; 47 } 48 return 0; 49 } 50 inline bool miller_rabin(LL n) 51 { 52 return judge(n,2)&&judge(n,3)&&judge(n,5)&&judge(n,7)&&judge(n,11)&&judge(n,13)&&judge(n,17)&&judge(n,19)&&judge(n,23); 53 } 54 LL *P; 55 int Pcnt; 56 inline LL pollard_rho(LL n) 57 { 58 LL x=rand()%(n-1)+1,xx=mul(x,x,n)+1,t=1; 59 while(t==1) 60 { 61 x=mul(x,x,n)+1; 62 xx=mul(xx,xx,n)+1; 63 xx=mul(xx,xx,n)+1; 64 t=__gcd(x-xx+n,n); 65 } 66 return t; 67 } 68 inline void divide(LL n) 69 { 70 if(n==1) 71 return ; 72 else 73 if(n%2==0) 74 P[++Pcnt]=2,divide(n/2); 75 else 76 if(miller_rabin(n)) 77 P[++Pcnt]=n; 78 else 79 { 80 LL x=pollard_rho(n); 81 divide(x); 82 divide(n/x); 83 } 84 } 85 inline void init(LL n,LL *p,int &cnt) 86 { 87 P=p; 88 Pcnt=0; 89 divide(n); 90 cnt=Pcnt; 91 } 92 } 93 #define mod 1000000007 94 inline int mul(int x,int y) 95 { 96 return ((LL)x*y)%mod; 97 } 98 inline int add(int x,int y) 99 { 100 x+=y; 101 if(x>=mod) 102 x-=mod; 103 return x; 104 } 105 inline int dec(int x,int y) 106 { 107 x-=y; 108 if(x<0) 109 x+=mod; 110 return x; 111 } 112 inline int qpow(int x,int y) 113 { 114 x%=mod; 115 int ans=1; 116 for(;y;y>>=1,x=mul(x,x)) 117 if(y&1) 118 ans=mul(ans,x); 119 return ans; 120 } 121 #define maxn 3002 122 int B[maxn],C[maxn][maxn],a[maxn],pw[20][maxn],iv[20][maxn]; 123 int p[20],c[20],x,y,N; 124 LL n; 125 int tot; 126 inline void init() 127 { 128 int i,j; 129 C[0][0]=1; 130 for(i=1;i<=3001;i++) 131 { 132 C[i][0]=1; 133 for(j=1;j<=i;j++) 134 C[i][j]=add(C[i-1][j],C[i-1][j-1]); 135 } 136 B[0]=1; 137 for(i=1;i<=3001;i++) 138 { 139 for(j=0;j<i;++j) 140 B[i]=dec(B[i],mul(C[i+1][j],B[j])); 141 B[i]=mul(B[i],qpow(i+1,mod-2)); 142 } 143 } 144 inline void pre() 145 { 146 for(int i=0;i<=y+1;i++) 147 a[i]=0; 148 for(int i=0;i<=y;i++) 149 a[y+1-i]=mul(C[y+1][i],B[i]); 150 int tmp=qpow(y+1,mod-2); 151 for(int i=0;i<=y+1;i++) 152 a[i]=mul(a[i],tmp); 153 a[y]=add(a[y],1); 154 if(!y) 155 a[0]=dec(a[0],1); 156 } 157 inline void get_frac() 158 { 159 static LL frac[97]; 160 int cnt=0; 161 tot=0; 162 Pollard_rho::init(n,frac,cnt); 163 sort(frac+1,frac+cnt+1); 164 for(int i=1;i<=cnt;i++) 165 if(frac[i]!=frac[i-1]) 166 p[++tot]=frac[i]%mod,c[tot]=1; 167 else 168 ++c[tot]; 169 for(int i=1;i<=tot;i++) 170 { 171 pw[i][0]=iv[i][0]=1; 172 int t=max(max(x,y),c[i])+1; 173 for(int j=1;j<=t;j++) 174 pw[i][j]=mul(pw[i][j-1],p[i]); 175 int v=qpow(p[i],mod-2); 176 for(int j=1;j<=t;j++) 177 iv[i][j]=mul(iv[i][j-1],v); 178 } 179 } 180 inline int geo(int i,int x) 181 { 182 int c=::c[i]; 183 int p; 184 if(x<0) 185 p=iv[i][-x]; 186 else 187 p=pw[i][x]; 188 if(p==1) 189 return c; 190 int ans=0,t=p; 191 for(int j=1;j<=c;j++) 192 { 193 ans=add(ans,t); 194 t=mul(t,p); 195 } 196 return ans; 197 } 198 inline void input() 199 { 200 read(n); 201 read(x); 202 read(y); 203 unsigned int rand_seed=18230742; 204 rand_seed+=(n+(x<<16)+y)+(rand_seed<<2); 205 srand(rand_seed); 206 } 207 inline int calc(int x,int y) 208 { 209 int ans=x<0?qpow(qpow(N,mod-2),-x):qpow(N,x); 210 for(int i=1;i<=tot;i++) 211 { 212 int t1=geo(i,-x); 213 int t2=dec(1,y>=0?pw[i][y]:qpow(p[i],mod-2)); 214 ans=mul(ans,add(1,mul(t1,t2))); 215 } 216 return ans; 217 } 218 void solve() 219 { 220 pre(); 221 get_frac(); 222 N=n%mod; 223 if(!N) 224 return (void)(puts("0")); 225 int ans=0; 226 int t=1,tmp; 227 for(int i=0;i<=y+1;i++) 228 { 229 tmp=calc(x-i,y-i); 230 ans=add(ans,mul(mul(a[i],t),tmp)); 231 t=mul(t,N); 232 } 233 ans=mul(ans,qpow(N,y)); 234 printf("%d\n",ans); 235 } 236 int main() 237 { 238 init(); 239 int T; 240 read(T); 241 while(T--) 242 { 243 input(); 244 solve(); 245 } 246 }
PS:$Pollard_rho$写丑了会TLE,需要大力卡一波常数才能过。