$Description:$
给定$n,x,y$，求$sumlimits_{i=1}^{n} gcd(i,n)^x lcm(i,n)^y$
$x,y le 3000$,$n le 10^{18}$,$mod=10^9+7$
挺久没有为单独一道题写一篇博客了，但是这是真的大神题。（最近总是被各路大神题吊起来爆捶）
这个式子乍一眼看起来挺朴实，化两步就越发绝望。
首先$gcd$和$lcm$乘起来是原数乘积。
$=sumlimits_{i=1}^{n} gcd(i,n)^{x-y} (in)^y$
$=n^ysumlimits_{i=1}^{n} gcd(i,n)^{x-y} i^y$
$=n^y sumlimits_{p|n} sumlimits_{i=1}^{n} [gcd(i,n)=p]p^{x-y} i^y$
$=n^y sumlimits_{p|n} sumlimits_{i=1}^{frac{n}{p}} [gcd(i,frac{n}{p})] p^x i^y $
$=n^y sumlimits_{p|n} p^x sumlimits_{i=1}^{frac{n}{p}} i^y sumlimits_{d|i and d|frac{n}{p} } mu(d)$
$=n^y sumlimits_{p|n} p^x sumlimits_{d|frac{n}{p} } d^y mu(d)sumlimits_{i=1}^{frac{n}{pd}} i^y$
化不动了。。。？但是在最后面我们看到了熟悉的自然数幂和。反正肯定算不了，于是先用伯努利数代换了再说。
这里用到的是正伯努利数$B_i^+$。相较于普通的伯努利数$B_1=-frac{1}{2}$，正伯努利数$B_1^+=frac{1}{2}$
其余项均相同，于是它们使用的公式也不太一样。
$sumlimits_{i=0}^{n} i^k =frac{1}{k+1} sumlimits_{i=0}^{k} C_{k+1}^{i} B_i (n+1)^{k+1-i} =sumlimits_{i=0}^{k} C_{k+1}^{i} B_i^+ n^{k+1-i}$
多数情况下貌似还是后者好用一些。
然而因为这是第一次真正用到伯努利数，所以再记载一下求法。
不少时候用这个递推式$O(n^2)$也是可以接受的：
$B_0=1,sumlimits_{i=0}^{n} B_i C_{n+1}^{i} =0$，也即$B_n=frac{1}{n} sumlimits_{i=0}^{n-1} B_i C_{n+1}^{i}$
然而有时候会被卡，可以用多项式生成函数的知识把上面的式子化一化得到下面这个极其好看的结论：$B(x)=frac{x}{e^x-1}$
写个多项式求逆可以做到$O(n log n)$。这个不会留坑待填。
然后不管用哪种方法求完之后把$B_1$取正就得到了更好用的正伯努利数。
所以回归本题，我们把正伯努利数代入之后：
$=n^y sumlimits_{p|n}p^x sumlimits_{d|frac{n}{p}} d^y mu(d) frac{1}{y+1} sumlimits_{i=0}^{y} C_{y+1}^{i} B_i^+ (frac{n}{pd})^{y+1-i}$
$=frac{n^y}{y+1} sumlimits_{i=0}^{y} C_{y+1}^{i} B_i^+ sumlimits_{p|n} p^x sumlimits_{d|frac{n}{p}} d^y mu(d) (frac{n}{pd})^{y+1-i}$
然后意犹未尽的我又对着这个式子化了半天，结果就越来越不可做。
然而停下来仔细看一看这个式子，最麻烦的还是整除关系，所以看看后两个和式，这种形式是卷积？
再仔细看一看就发现它们几个都是积性函数。这好啊，于是我们就可以对于每种质因子求解然后把贡献相乘就是整个$n$的函数值。
然而为什么要分解？其实主要还是利用了出现在函数内部的$mu$，因为如果只含有一种因子$p$的话，那么整个式子就只有在$p^0$和$p^1$时有值。
$=frac{n^y}{y+1}sumlimits_{i=0}^{y} C_{y+1}^{i} B_i^+ prodlimits_{p^k|n} sumlimits_{t=0}^{k} (p^t)^x sumlimits_{u=0}^{min(k-t),1} (p^u)^y mu(p^u) (p^{k-t-u})^{y+1-i}$
只有在$k=t$时特殊，$u$不能取$1$。剩下的情况都一样，分$u=0$和$u=1$考虑，列出式子：
$=frac{n^y}{y+1}sumlimits_{i=0}^{y} C_{y+1}^{i} B_i^+ prodlimits_{p^k|n} sumlimits_{t=0}^{k-1} p^{tx} (p^{(k-t)(y+1-i)}-p^{y+(k-t-1)(y+1-i)} + p^{kx}) $
这个式子就可以算了。最外层枚举$i$，上界$y$，内层枚举质因子及其次数是$O(log n)$级别。快速幂的$O(log)$算上。复杂度合法，稍微有点卡常。
这先不说，关键是回避了一个问题：怎么质因数分解啊？
这就是$Pollard Rho$的锅了。
这个理解不深刻只能大概口胡。
我们要求$n$的质因子，这不好弄，于是退而求其次，我们求$n$的一个因子$p$然后递归地对$p$和$frac{n}{p}$求因子。
递归边界就是这个数是个质数，那么它就是原数的质因子了。判质数套一个$Miller Rabin$复杂度也是可以接受的。
然而现在的问题是如何快速的找到一个数的一个因子？
运用生日悖论，我们既然直接找到一个因子很困难，那么我们就随机一些数然后两两做差。
这里有一个原理不明的随机数生成：$x_i=x_{i-1}^2 +c(mod n)$。c是随机常数。很多题解都在声明这玩意的确随机，那么为什么不直接用$rand()$呢？
不解，随大溜。然而这种只有一个参数的取模函数一定会进入循环，进入循环节之后要及时跳出。
然后又引出了$Floyd$的龟兔赛跑判环，然而并用不上不再赘述。
最后被广泛采纳的方法是倍增，其实有点迭代加深的意思？
因为$gcd$的复杂度是$O(log)$的，复杂度会很卡，所以我们怎么能少做几次$gcd$？
我们可以把随机数相邻项做差，累乘起来，如果其中一项和$n$的$gcd$非1，那么最终累乘后的$gcd$也不会是1
然后每做多少次进行一次$gcd$比较好？前人的经验告诉我们：$127$。好，我记住了就这样吧。
为什么说有迭代加深/倍增呢？因为对于任意质数我们都跑$127$次太浪费了，我们先试试1次行不行，再2次，再4次。。。如果找到就直接跳出。
邻项做差到底是如何提高概率的，这个我们无从知道。算是实践出真知吧。
然而如果你脸黑随机出的c无论跑多少次，进了环还是找不到因子，那么其实应该要及时跳出的。但是貌似并不存在这种情况（脸不够黑）
所以就这么稀里糊涂的会用了。于是这题就可以做了。
1 #include<bits/stdc++.h> 2 using namespace std; 3 #define ll long long 4 map<ll,int>M; 5 #define mod 1000000007 6 int x,y,pc,k[19];ll C[3333][3333],B[3333],ans,n,p[19]; 7 ll gcd(ll a,ll b){return b?gcd(b,a%b):a;} 8 ll mul(ll a,ll b,ll m){ 9 a%=m;b%=m; 10 if(a<0x7fffffff&&b<0x7fffffff)return a*b%m; 11 ll c=(long double)a/m*b,r=(unsigned ll)a*b-(unsigned ll)c*m;return (r%m+m)%m; 12 } 13 ll pow(ll b,ll t,ll m,ll a=1){for(;t;t>>=1,b=mul(b,b,m))if(t&1)a=mul(a,b,m);return a;} 14 ll pow0(ll b,int t,ll a=1){b%=mod;for(;t;t>>=1,b=b*b%mod)if(t&1)a=a*b%mod;return a;} 15 bool Miller_Rabin(ll x,ll b){ll k=x-1; 16 while(!(k&1)){ 17 ll y=pow(b,k,x); 18 if(y!=x-1&&y!=1)return 0; 19 if(y==x-1)return 1;k>>=1; 20 }return 1; 21 } 22 bool isprime(ll x){ 23 if(x==2||x==3||x==5||x==7||x==11||x==13||x==61)return 1; 24 if(x%2==0||x%3==0||x%5==0||x%7==0||x%11==0||x%13==0||x==1)return 0; 25 return Miller_Rabin(x,2)&&Miller_Rabin(x,11)&&Miller_Rabin(x,61); 26 } 27 ll Pollard_Rho(ll x){ 28 ll s=0,t=0,c=rand()%(x-1)+1,v=1,d; 29 for(int g=1;;g<<=1,s=t,v=1){ 30 for(int stp=1;stp<=g;++stp){ 31 t=(mul(t,t,x)+c)%x; v=mul(v,abs(t-s),x); 32 if(stp%127==0){d=gcd(v,x);if(d>1)return d;} 33 }d=gcd(v,x);if(d>1)return d; 34 } 35 } 36 void Div(ll x){ 37 if(x<2)return; 38 if(isprime(x)){M[x]++;return;} 39 ll p=Pollard_Rho(x);Div(p);Div(x/p); 40 } 41 int main(){ 42 for(int i=0;i<3003;++i)C[i][0]=1;B[0]=1; 43 for(int i=1;i<3003;++i)for(int j=1;j<=i;++j)C[i][j]=(C[i-1][j-1]+C[i-1][j])%mod; 44 for(int i=1;i<3003;++i){ 45 for(int j=0;j<i;++j)B[i]=(B[i]+mod-B[j]*C[i+1][j]%mod)%mod; 46 B[i]=B[i]*pow0(i+1,mod-2)%mod; 47 }B[1]=mod-B[1]; 48 int t;cin>>t;while(t--){ 49 cin>>n>>x>>y;Div(n); 50 for(auto o:M)p[++pc]=o.first,k[pc]=o.second; 51 for(int i=0;i<=y;++i){ 52 ll rate=C[y+1][i]*B[i]%mod,tot=1; 53 for(int pn=1;pn<=pc;++pn){ 54 ll tsum=0,P=p[pn],K=k[pn],bs=pow0(P,(K-1)*(y+1-i)+y),stp=pow0(P,x-y-1+i+mod-1),sub=pow0(P,mod-i)-1; 55 for(int t=0;t<K;++t)tsum=(tsum+bs*sub)%mod,bs=bs*stp%mod; 56 tot=tot*(tsum+pow0(P,K*x))%mod; 57 }ans=(ans+tot*rate)%mod; 58 }cout<<ans*pow0(n,y)%mod*pow0(y+1,mod-2)%mod<<endl; 59 pc=ans=0;M.clear(); 60 } 61 }