BZOJ2627 JZPKIL
一句话题意：求$\sum_{i=1}^ngcd(i,n)^xlcm(i,n)^y$，多组数据，模$1e9+7$。
30%的数据，$x=y$
另30%的数据，$n<=10^9$, $x, y<=100$
100%的数据，$T<=100, 1<=n<=10^{18}, 0<=x, y<=3000$
orz出题人顾宇宙！！！orz全网题解的源头ydc！！！感觉ydc的题解已经写得很清楚了，但是为了内部交流还是写一写……下面的公式中$[x]$表示$x==1$。
推导部分
Part1 原式整理$$\sum_{i=1}^ngcd(i,n)^xlcm(i,n)^y$$$$=\sum_{i=1}^n(in)^ygcd(i,n)^{x-y}$$$$=n^y\sum_{i=1}^ni^ygcd(i,n)^{x-y}$$$$=n^y\sum_{d|n}d^{x-y}\sum_{k=1,[gcd(k,\frac{n}{d}))]}^{\frac{n}{d}}(kd)^y$$$$=n^y\sum_{d|n}d^x\sum_{k=1,[gcd(k,\frac{n}{d})]}^{\frac{n}{d}}k^y$$
Part2 套路反演
这一部分我们来整理上式最后的$\sum_{k=1,[gcd(k,\frac{n}{d})]}^{\frac{n}{d}}k^y$
$$\sum_{i=1,[gcd(i,n)]}^ni^y$$$$=\sum_{i=1}^ni^y\sum_{d|gcd(i,n)}\mu(d)$$$$=\sum_{d=1}^n\mu(d)\sum_{i=1}^{\frac{n}{d}}(id)^y$$$$=\sum_{d=1}^n\mu(d)d^y\sum_{i=1}^{\frac{n}{d}}i^y$$
Part3 伯努利数
可以发现Part2的最后我们把Part1的幂之和化为了从1开始的连续一段……这有什么用呢？我们可以套结论了！有一个叫伯努利数的东西（下文推导中用$B_k$表示第$k$个伯努利数），在网上搜一搜的话资料不少，百度百科里基本的公式也给了一些。刚开始看证明说什么从生成函数定义的角度可以证明它和自然数幂和的数学关系，但是我后来看到它就是在求自然数幂和的过程中被发明的就感觉证这个没必要了……求伯努利数的方法放在实现部分里，现在我们继续关注推导过程。
$$\sum_{i=1}^ni^m=\frac{1}{m+1}\sum_{k=0}^mC_{m+1,k}B_kn^{m-k+1}$$
大多数情况下伯努利数都是这么用的……可以看成公式
你可能会发现伯努利数的公式很乱，有人用$n$的幂有人用$n+1$的幂还有人带-1的幂……然而因为不会证我也不能确切地说哪一种是对的。维基百科上大概是说伯努利数有两类，两类伯努利数的区别在于B1。第一类伯努利数，由美国国家标准技术研究所 (NIST)制定，在这标准下B1=-1/2;第二类伯努利数，又被称为是“原始的伯努利数”，在这标准下B1=1/2。$n+1$的式子是关于第一类伯努利数的，而$n$的式子是关于第二类伯努利数的;第一类伯努利数乘上-1的幂即为第二类伯努利数……假装他们都是对的好了编套理论真难啊
设$t_k=\frac{1}{m+1}C_{m+1,k}B_k$，则$\sum_{i=1}^ni^m=\sum_{k=0}^mt_kn^{m-k+1}$
Part4 积性函数
即使化成这步代入原式它也不好处理……不过看一眼$10^{18}$的数据范围，分解质因数来做好像是一种历史必然？so现在我们试一试把积性函数的部分整理到一起会发生什么……
$$n^y\sum_{d|n}d^x\sum_{d_1|\frac{n}{d}}\mu(d_1)d_1^y\sum_{k=0}^yt_k(\frac{n}{dd_1})^{y-k+1}$$
$$=\sum_{k=0}^yt_kn^y\sum_{d|n}d^x\sum_{d_1|\frac{n}{d}}\mu(d_1)d_1^y(\frac{n}{dd_1})^{y-k+1}$$
这时候除去$t_k$的部分就是一个积性函数了，$y$很小只有$3000$，我们可以枚举$k$然后分解质因数求后面的部分。那么就变成了对每一个质因子$p_i$的$k_i$幂$p_k$求$p_k^y\sum_{d|p_k}d^x\sum_{d_1|\frac{p_k}{d}}\mu(d_1)d_1^y(\frac{p_k}{dd_1})^{y-k+1}$，又因为$\mu(d_1)$只有在$d_1=1$时为1、$d_1=p_i$时为-1而其他情况下都为0，式子可以直接变成$p_k^y\sum_{d|p_k}d^x((\frac{p_k}{d})^{y-k+1}-p_i^y\frac{p_k}{dp_i}^{y-k+1})$，这就是推式子的终点了。
实现部分
Part1 伯努利数
你问我伯努利数怎么求啊……
$\sum_{k=0}^nC_{n+1,k}B_k=n+1$
所以杨辉三角求一下组合数就可以$O(k^2)$求出所有伯努利数了。
Part2 分解质因数
给$10^{18}$分解质因数？需要$n^{\frac{1}{4}}$的算法！Miller-Rabin和Pollard-Rho可以解决，关于这两个算法我们的朋友Doggu的blog堪称全网最良心。
Part3 快速乘
$10^{18}*10^{18}$应该怎么做？快速乘显然是必要的。这是ydc改良过的快速乘板子。
1 inline LL mul(LL x,LL y,LL z) 2 { 3 return (x*y-(LL)(((long double)x*y+0.5)/(long double)z)*z+z)%z; 4 }
好像少了点什么东西
没有错，这题是有部分分的吧？对于$x=y$的$30$分，式子就是$\sum_{i=1}^n(in)^y$，如果你会用伯努利数就可以拿到这一部分的分。对于$n<=10^9$的$30$分，伯努利数仍然是必要的，不过分解质因数可以$\sqrt{n}$来做。看来出题人只打算给会伯努利数的人分数（显然我不在此之列，所以部分分和正解一样困难……）。想打部分分的同学可以去cogs上提交。
最后放代码，完结撒花。
1 #include<iostream> 2 #include<cstdio> 3 #include<cstring> 4 #include<cstdlib> 5 #include<algorithm> 6 #define ll long long 7 using namespace std; 8 const int mod=1e9+7; 9 const int sj=3010; 10 const int ss=110; 11 inline ll read() 12 { 13 ll jg=0;int jk=getchar()-'0'; 14 while(jk<0||jk>9) jk=getchar()-'0'; 15 while(jk>=0&&jk<=9) jg*=10,jg+=jk,jk=getchar()-'0'; 16 return jg; 17 } 18 int ca,ai,bi,cnt,tot,s[ss],sp,t[sj][sj],c[sj][sj]; 19 ll n,a1,a2,b[sj],pi[ss],pk[ss]; 20 ll p[sj],prime[10]={2,3,5,7,11,13,17,19,23,29},fac[ss]; 21 ll pw[ss][ss],pw1[ss]; 22 bool vi[sj]; 23 ll gcd(ll x,ll y) 24 { 25 if(!y) return x; 26 return gcd(y,x%y); 27 } 28 inline ll mul(ll x,ll y,ll z) 29 { 30 return (x*y-(ll)(((long double)x*y+0.5)/(long double)z)*z+z)%z; 31 } 32 ll qpow(ll x,ll y,ll z) 33 { 34 ll ret=1;x%=z; 35 while(y) 36 { 37 if(y&1) ret=mul(ret,x,z); 38 y>>=1,x=mul(x,x,z); 39 } 40 return ret; 41 } 42 ll ksm(ll x,int y) 43 { 44 ll ret=1;x%=mod; 45 while(y) 46 { 47 if(y&1) ret=ret*x%mod; 48 x=x*x%mod,y>>=1; 49 } 50 return ret; 51 } 52 void pre() 53 { 54 for(int i=0;i<=3005;i++) 55 { 56 c[i][0]=1; 57 for(int j=1;j<=i;j++) 58 { 59 c[i][j]=c[i-1][j]+c[i-1][j-1]; 60 if(c[i][j]>=mod) c[i][j]-=mod; 61 } 62 } 63 for(int i=0;i<=3005;i++) 64 { 65 b[i]=i+1; 66 for(int j=0;j<i;j++) 67 { 68 b[i]=(b[i]-c[i+1][j]*b[j])%mod; 69 if(b[i]<0) b[i]+=mod; 70 } 71 b[i]=b[i]*ksm(c[i+1][i],mod-2)%mod; 72 } 73 for(int i=0;i<=3005;i++) 74 { 75 a1=ksm(i+1,mod-2); 76 for(int j=0;j<=i;j++) 77 t[i][j]=c[i+1][j]*b[j]%mod*a1%mod; 78 } 79 for(int i=2;i<=1000;i++) 80 for(int j=2;j<=1000;j++) 81 if(i*j<=1000) 82 vi[i*j]=1; 83 for(int i=2;i<=1000;i++) 84 if(!vi[i]) p[++cnt]=i; 85 } 86 inline ll calc(ll x) 87 { 88 ll ret=0,q=x%mod; 89 ll nt=q; 90 for(int i=bi;i>=0;i--) 91 { 92 ret=(ret+t[bi][i]*nt)%mod; 93 nt=nt*q%mod; 94 } 95 return ret*ksm(q,bi)%mod; 96 } 97 bool Miller_Rabin(ll x) 98 { 99 if(x==1) return 0; 100 for(int j=0;j<=9;j++) 101 if(x==prime[j]) 102 return 1; 103 ll y=x-1; 104 int k=0; 105 while((y&1)==0) k++,y>>=1; 106 for(int i=0;i<10;i++) 107 { 108 ll z=rand()%(x-1)+1; 109 ll c=qpow(z,y,x),d; 110 for(int j=0;j<k;j++) 111 { 112 d=mul(c,c,x); 113 if(d==1&&c!=1&&c!=x-1) return 0; 114 c=d; 115 } 116 if(d!=1) return 0; 117 } 118 return 1; 119 } 120 ll pollard_rho(ll x,ll y) 121 { 122 ll i=1,k=2,c=rand()%(x-1)+1; 123 ll d=c; 124 while(1) 125 { 126 i++; 127 c=(mul(c,c,x)+y)%x; 128 ll g=gcd((d-c+x)%x,x); 129 if(g!=1&&g!=x) return g; 130 if(c==d) return x; 131 if(i==k) d=c,k<<=1; 132 } 133 } 134 void find(ll x,int tim) 135 { 136 if(x==1) return; 137 if(Miller_Rabin(x)) 138 { 139 fac[++tot]=x; 140 return; 141 } 142 ll yz=x,mk=tim; 143 while(yz>=x) yz=pollard_rho(yz,tim--); 144 find(yz,mk),find(x/yz,mk); 145 } 146 ll solve(ll x,ll y,ll z) 147 { 148 ll ret=0,qwq,s1,s2; 149 tot=sp=0; 150 for(int j=1;j<=cnt;j++) 151 while(x%p[j]==0) 152 fac[++tot]=p[j],x/=p[j]; 153 find(x,120); 154 sort(fac+1,fac+tot+1); 155 for(int i=1;i<=tot;i++) 156 { 157 if(fac[i]!=fac[i-1]) sp++,pi[sp]=pk[sp]=fac[i],s[sp]=1; 158 else pk[sp]*=fac[i],s[sp]++; 159 } 160 for(int i=1;i<=sp;i++) 161 { 162 pw[i][0]=1,a1=ksm(pi[i],y); 163 for(int j=1;j<=s[i];j++) 164 pw[i][j]=pw[i][j-1]*a1%mod; 165 } 166 for(int i=0;i<=z;i++) 167 { 168 qwq=1; 169 for(int j=1;j<=sp;j++) 170 { 171 s1=s2=0,pw1[0]=1,a1=ksm(pi[j],z),a2=ksm(pi[j],z-i+1); 172 for(int k=1;k<=s[j];k++) pw1[k]=pw1[k-1]*a2%mod; 173 for(int k=0;k<=s[j];k++) s1=(s1+pw[j][k]*pw1[s[j]-k])%mod; 174 for(int k=0;k<s[j];k++) s2=(s2+pw[j][k]*pw1[s[j]-k-1]%mod*a1)%mod; 175 qwq=qwq*(s1-s2+mod)%mod*ksm(pk[j],z)%mod; 176 } 177 ret=(ret+qwq*t[z][i])%mod; 178 } 179 return ret; 180 } 181 int main() 182 { 183 pre(); 184 ca=read(); 185 for(int i=1;i<=ca;i++) 186 { 187 n=read(),ai=read(),bi=read(); 188 if(ai==bi) printf("%lld\n",calc(n)); 189 else printf("%lld\n",solve(n,ai,bi)); 190 } 191 return 0; 192 }
