【BZOJ2627】JZPKIL（数论大杂烩）
大致题意： 求\(\sum_{i=1}^ngcd(i,n)^xlcm(i,n)^y\)。
前言
该怎么说呢？不管怎样，至少最近在数论方面还是有点长进的吧，除了伯努利数那个地方实在是不知道，其他式子都是自己一步一步推出来的。
尽管这可能很大程度是归功于做过一道这题的简单版：【BZOJ3601】一个人的数论。。。（那道题十分良心，题目中给出的是\(n\)的质因数分解，不需要\(Pollard-Rho\)）
这是一道数论大杂烩，涉及到莫比乌斯反演、伯努利数、积性函数、\(Pollard-Rho\)（说到它，自然也包含\(MillerRabin\)）等众多知识，请做好心理准备。
第一次（大概吧）能把数论题的代码写到\(100\)行。。。
推式子+各种准备工作共计半小时，写代码一个小时，调试半小时，写博客又快一个小时，共计三小时。。。
莫比乌斯反演
众所周知，\(gcd(i,n)\cdot lcm(i,n)=i\cdot n\)，因此我们完全可以去除\(lcm(i,n)\)，得到：
我们把\(n^y\)提前。看到有个\(gcd\)，于是就开始喜闻乐见地枚举\(gcd\)：
由\(\sum_{d|n}\mu(d)=[n=1]\)，得\(\sum_{p|i,p|\frac nd}\mu(p)=1\)，然后枚举\(p\)得到：
在这个式子中，你看到了什么？没错，是它，是它，就是它，\(\sum_{i=1}^{\frac n{dp}}i^y\)！
它来了，自然数幂和，它来了！
在上面给出的那道题目中，同样需要自然数幂和。但由于它的次数\(d\)比较小，可以\(O(d^3)\)高斯消元。
但这道题\(y\le3000\)，我觉得它就是在难为我高斯消元，需要采用一些更高级的东西。
自然数幂和
考虑我们设\(S(n)=\sum_{i=1}^ni^y\)。
这里引入一个叫伯努利数的东西，它的定义是\(B_0=1\)，且满足：
由它的性质实际上我们只要把\(B_n\)移到不等式一边，就可以得到递推式：
然后有一个叫做伯努利公式的东西，长这副样子：（注意，这里\(B_i\)指正伯努利数，而上面的\(B_1<0\)，此处要取绝对值）
你会发现这个公式和前面两个式子长得特别像（感觉就是二者的融合版本，然后乘个带\(n\)的项），不过我也没去深究它的证明。。。
把这个公式套到我们的\(S(n)\)上，得到：（其实就是把\(k\)换成\(y\)，然后把\(\frac1{y+1}\)移到了里面。。。）
我们令\(a_i=\frac1{y+1}C_{y+1}^iB_i\)，简化原式得到：
积性函数
我们把\(S(n)\)代回原式，得到：
然后我们提前枚举\(i\)，并令\(D=dp\)，得到：
显然\(d\cdot(\frac Dd)=D\)，我们把它们拼起来，然后又发现\((\frac nD)\cdot D=n\)，于是接着拼起来，最终得到：
然后我们发现后面的一堆项是积性函数，于是设：
假设\(n=\prod_{q=1}^{cnt}P_q^{k_q}\)（注意，\(n\le10^{18}\)，所以这里需要\(Pollard-Rho\)来对\(n\)质因数分解），那么根据积性函数的性质显然有：
所以我们只需要考虑单个质数幂的函数得到：
我们再先单独考虑后面的项，假设：（其中\(D\)是一个质数的幂）
质数的幂，以及，\(\mu\)？根据经验，当两者碰到一起，显然会发生神奇的事情。
根据\(\mu\)的定义，当\(\frac Dd\ge P_q^{2}\)时，\(\mu(\frac Dd)\)就等于\(0\)了，所以我们只需考虑\(\frac Dd=1\)和\(\frac Dd=P_q\)两种情况：
（注意，当\(D=1\)的时候可能会出现点小问题——\(\frac D{P_q}\)挂了，这里先暂时不管，后面会解决的）
然后我们把这个东西放回原式里得到：
考虑我们干脆令\(D=P_q^t\)，然后枚举\(t\)。由于我们之前提到过的\(D=1\)时的小问题，我们单独处理\(D=1\)的情况（不要代入这个式子，代回原式就能发现显然此时值为\(1\)），并从\(1\)开始枚举\(t\)，得到：
看起来已经很完美了，但实际上即便只是快速幂的\(log\)还是会被卡。。。（当然如果你是常数之神就当我没说）
所以我们继续往下，先提取出括号中的\(P_q^{(t-1)(x-y)}=P_q^{t(x-y)-(x-y)}\)得到：
再进一步化简得到：
然后我们就发现，只要先计算出\(P_q^{i-1+x-y},P_q^{x-y}\)，并在枚举\(t\)同时累乘出\(P_q^{t(i-1+x-y)}\)，就不需要再快速幂了。
这样一来这题就真正做完了。
最后再放一下最终的式子：
代码
#include<bits/stdc++.h>
#define Tp template<typename Ty>
#define Ts template<typename Ty,typename... Ar>
#define Reg register
#define RI Reg int
#define Con const
#define CI Con int&
#define I inline
#define W while
#define X 1000000007
#define LL long long
#define RL Reg LL
#define CL Con LL&
using namespace std;
LL n;int x,y;
namespace PollardRho//质因数分解模板
{
I LL gcd(CL x,CL y) {return y?gcd(y,x%y):x;}
I LL QM(CL x,CL y,CL P) {RL k=(1.0L*x*y)/P,t=x*y-k*P;t-=P;W(t<0) t+=P;return t;}
I LL QP(RL x,RL y,CL P) {RL t=1;W(y) y&1&&(t=QM(t,x,P)),x=QM(x,x,P),y>>=1;return t;}
namespace MR//MillerRabin质数测试
{
#define Pt 12
const int P[Pt]={2,3,5,7,11,13,17,19,61,2333,4567,24251};
I bool Check(CL x,CI p)
{
if(!(x%p)||QP(p%x,x-1,x)^1) return 0;
RL k=x-1,t;W(!(k&1)) {if((t=QP(p%x,k>>=1,x))^1&&t^(x-1)) return 0;if(t^1) return 1;}
return 1;
}
I bool IsP(CL x)
{
if(x==1) return 0;for(RI i=0;i^Pt;++i)
{if(x==P[i]) return 1;if(!Check(x,P[i])) return 0;}
return 1;
}
}
int cnt;struct data
{
LL p;int k;I data(CL a=0,CI b=0):p(a),k(b){}
I bool operator < (Con data& o) Con {return p<o.p;}
I bool operator == (Con data& o) Con {return p==o.p;}
}P[60];
I LL Work(CL x,CI y)
{
#define R(x) (1LL*rand()*rand()*rand()*rand()%(x)+1)
RI t=0,k=1;RL v0=R(x-1),v=v0,d,s=1;W(true)
{
if(v=(QM(v,v,x)+y)%x,s=QM(s,abs(v-v0),x),!s) return x;
if(++t==k) {if((d=gcd(s,x))^1) return d;v0=v,k<<=1;}
}
}
I void Resolve(RL x,CI w,RI t)
{
if(x==1) return;if(MR::IsP(x)) return (void)(P[++cnt]=data(x,w));
RL y;W((y=Work(x,t--))==x);RI u=0;W(!(x%y)) x/=y,++u;Resolve(x,w,t-1),Resolve(y,u*w,t-1);
}
I void Resolve(CL x)
{
cnt=0,srand(20050521),Resolve(x,1,302627441),sort(P+1,P+cnt+1);
RI i,tmp=cnt;for(i=2,cnt=1;i<=tmp;++i) P[i]==P[cnt]?(P[cnt].k+=P[i].k):(P[++cnt]=P[i]);//排序合并同一质数
}
}using namespace PollardRho;
I int QP(RI x,RI y) {RI t=1;y<0&&(y+=X-1);W(y) y&1&&(t=1LL*t*x%X),x=1LL*x*x%X,y>>=1;return t;}
class Bernoulli//伯努利数
{
private:
#define SZ 3000
int B[SZ+5],C[SZ+5][SZ+5];
public:
I int operator [] (CI i) Con {return 1LL*QP(y+1,X-2)*C[y+1][i]%X*B[i]%X;}//求出第i项系数
I Bernoulli()//预处理
{
RI i,j;for(C[0][0]=i=1;i<=SZ+1;++i)//组合数
for(C[i][0]=j=1;j<=i;++j) C[i][j]=(C[i-1][j-1]+C[i-1][j])%X;
for(B[0]=i=1;i<=SZ;++i)//伯努利数
{
for(j=0;j^i;++j) B[i]=(1LL*C[i+1][j]*B[j]+B[i])%X;
B[i]=1LL*(X-1)*QP(i+1,X-2)%X*B[i]%X;
}B[1]=X-B[1];//注意是正伯努利数
}
}A;
I int f(CI i,CI p,CI k)//求质数幂的f
{
RI t,res=0,Pnow=1,Pbase=QP(p,i-1+x-y),P0=QP(p,x-y);//先计算好需要的幂
for(t=1;t<=k;++t) res=(1LL*(Pnow=1LL*Pnow*Pbase%X)*(P0-1)+res)%X;//累乘避免快速幂
return 1LL*res*QP(P0,X-2)%X+1;//返回答案
}
I int F(CI i)//求n的F
{
RI q,res=1;for(q=1;q<=cnt;++q) res=1LL*f(i,P[q].p%X,P[q].k)*res%X;//枚举质因数把答案乘起来
return res;
}
int main()
{
RI Tt,i,ans;scanf("%d",&Tt);W(Tt--)
{
scanf("%lld%d%d",&n,&x,&y),Resolve(n);
for(ans=i=0;i<=y;++i) ans=(1LL*A[i]*QP(n%X,y-i+1)%X*F(i)+ans)%X;//根据最终式子计算
printf("%d\n",1LL*QP(n%X,y)*ans%X);
}return 0;
}
