[数论记录]P4464 [国家集训队]JZPKIL

神仙题目,做一题学一堆系列…… 题意:求$\sum_{i=1}^{n}gcd(i,n)^x lcm(i,n)^y \bmod (10^9+7)$ $n<=10^{18};\ x,y<3000$ $ANS=\sum_{i=1}^{n}gcd(i,n)^x lcm(i,n)^y$ $=n^y\sum_{i=1}^{n}gcd(i,n)^{x-y}i^y$ 先忽略掉$n^y$。 $=\sum_{d|n}\sum_{i=1}^{n}[(i,n)==d]d^{x-y}i^y$ $=\sum_{d|n}\sum_{i=1}^{n/d}[(i,\frac{n}{d})==1]d^{x-y}(id)^y$ $=\sum_{d|n}d^x\sum_{i=1}^{n/d}[(i,\frac{n}{d})==1]i^y$ $=\sum_{d|n}d^x\sum_{i=1}^{n/d}i^y\sum\limits_{j|i,j|\frac{n}{d}}\mu(j)$ $=\sum_{d|n}d^x\sum\limits_{j|\frac{n}{d}}\mu(j)\sum_{j|i}^{n/d}i^y$ $=\sum_{d|n}d^x\sum\limits_{j|\frac{n}{d}}\mu(j)\sum_{i=1}^{n/dj}(ij)^y$ $=\sum_{d|n}d^x\sum\limits_{j|\frac{n}{d}}\mu(j)j^y\sum_{i=1}^{n/dj}i^y$ 推导到这里,先冷静观察一下式子: 后面$\sum_{i=1}^{n/dj}i^y$是自然数幂和,这个好办。 伯努利数$O(y)$搞出多项式系数就好,需要$O(y^2)$的预处理 。 不会的话请见 [多项式计数杂谈](https://www.luogu.com.cn/blog/command-block/sheng-cheng-han-shuo-za-tan) 得到$\sum_{i=1}^{n}i^y=\sum_{i=0}^{y+1}c_in^i$($c$是多项式系数) $ANS=\sum_{d|n}d^x\sum\limits_{j|\frac{n}{d}}\mu(j)j^y\sum_{i=0}^{y+1}c_i(n/dj)^i$ $=\sum_{i=0}^{y+1}c_i\sum_{d|n}d^x\sum\limits_{j|\frac{n}{d}}\mu(j)j^y(n/dj)^i$ 仔细看看就会发现后面的$\sum_{d|n}d^x\sum\limits_{j|\frac{n}{d}}\mu(j)j^y(n/dj)^i$是三个函数卷积的形式: $F1(n)=n^x;\ \ F2(n)=\mu(n)n^y;\ \ F3(n)=n^i$ 那么这个卷积后的东西$R(n)=F1*F2*F3$也应该是个积性函数。 我们只要设法求出任意的$R(p^k)$($p$指质数),就可以利用积性了。 $R(p^k)=\sum\limits_{j|p^k}\mu(j)j^y\sum_{d|\frac{p^k}{j}}d^x(p^k/dj)^i$ 根据$\mu$很有可能为0,这里的$j=p$或$1$。 - $j=1$ 贡献$=\sum_{d|p^k}d^x(p^k/d)^i$ $=p^{ik}\sum_{d|p^k}d^{x-i}$ $=p^{ik}\sum_{t=0}^{k}p^{t(x-i)}$ 为了躲掉逆元,我们把负指数变一下: $=\sum_{t=0}^{k}p^{tx+(k-t)i}$ 这里$k$很小,暴力求就好,下面同理。 - $j=p$ $-p^y\sum_{d|p^{k-1}}d^x(p^{k-1}/d)^i$ $=-p^{y}($上面$k-1$的情况$)$ 懒得汇总了…… 我们算$y$次这个$R$就能解决问题了。 上面的全部暴力快速幂实现,求一次也就$log^2$。 那么现在还剩下一个问题就是分解,Prho硬上即可,分解的结果要保存下来。 总复杂度$O(T(\sqrt[4]{n}+ylog^2n)+y^2)$ ```cpp #include<algorithm> #include<cstdio> #define mod 1000000007 #define ll long long using namespace std; ll powM(ll a,ll t=mod-2) { ll ans=1;a%=mod; while(t){ if(t&1)ans=ans*a%mod; a=a*a%mod; t>>=1; }return ans; } namespace BNL { #define Maxn 3050 ll B[Maxn],ifac[Maxn],fac[Maxn]; void Init() { fac[0]=ifac[0]=1; for (int i=1;i<=3005;i++){ ifac[i]=ifac[i-1]*powM(i)%mod; fac[i]=fac[i-1]*i%mod; }B[0]=1; for (int j=1;j<=3005;j++){ for (int i=0;i<j;i++) B[j]=(B[j]+B[i]*ifac[j-i+1])%mod; B[j]=mod-B[j]; } } void getf(ll *f,int k) { for (int i=1;i<=k+1;i++) f[i]=B[k-i+1]*ifac[i]%mod*fac[k]%mod; f[0]=0;if (k)f[k]++; } } inline ll mul(ll a,ll b,ll m){ ll d=((long double)a/m*b+0.5); ll r=a*b-d*m; return r<0?r+m:r; } ll powM(ll a,ll t,ll m) { ll ans=1; while(t){ if (t&1)ans=mul(ans,a,m); a=mul(a,a,m); t>>=1; }return ans; } bool mr(ll n,ll a) { ll t=n-1; while(!(t&1))t>>=1; ll buf=powM(a,t,n); if (buf==1||buf==n-1)return 0; while((t<<=1)<n-1){ buf=mul(buf,buf,n); if (buf==n-1)return 0; }return 1; } const int testp[8]={2,3,5,7,13,19,61,233}; bool ptest(ll n) { if (n<2)return 0; for (int i=2;i*i<=min(n,10000ll);i++) if (n%i==0)return 0; if (n<=10000)return 1; for (int i=0;i<8;i++) if (mr(n,testp[i]))return 0; return 1; } inline ll gcd(ll a,ll b){ if(a==0)return b; if(a<0)a=-a; ll t; while(b){t=a%b;a=b;b=t;} return a; } ll sav[150]; const int lim=128; ll prho(ll n,ll c) { ll x1=(c+1)%n,x2=(mul(x1,x1,n)+c)%n,buf; int tot=0; while(1){ buf=mul(x1-x2,buf,n); sav[++tot]=x1-x2; if (tot==lim){ if (gcd(buf,n)>1)break; tot=0; } x1=(mul(x1,x1,n)+c)%n; x2=(mul(x2,x2,n)+c)%n; x2=(mul(x2,x2,n)+c)%n; } for (int i=1;i<=tot;i++){ buf=gcd(sav[i],n); if (buf>1)return buf; }return n; } int tot,savc[127]; ll savf[127]; void factor(ll n) { if (ptest(n)){ for (int i=1;i<=tot;i++) if (savf[i]==n) {savc[i]++;return ;} savf[++tot]=n;savc[tot]=1; return ; } ll sav=prho(n,rand()%(n-1)+1); while(sav==n)sav=prho(n,rand()%(n-1)+1); factor(sav);factor(n/sav); } ll c[3050],n;int x,y; ll calc(ll p,int k,int i) { ll ans=0; for (int t=0;t<=k;t++) ans=(ans+powM(p,t*x+(k-t)*i))%mod; return ans; } void solve() { scanf("%lld%d%d",&n,&x,&y); ll savn=n%mod; BNL::getf(c,y); for (int i=2;i<=10000&&i*i<=n;i++) if (n%i==0){ savf[++tot]=i; while (n%i==0){ n/=i;savc[tot]++; } } if (n>1)factor(n); ll ans=0; for (int i=0;i<=y+1;i++){ ll buf=1; for (int j=1;j<=tot;j++){ ll sav=calc(savf[j],savc[j],i); sav=(sav-powM(savf[j],y) *calc(savf[j],savc[j]-1,i))%mod; buf=buf*sav%mod; } ans=(ans+c[i]*buf)%mod; }printf("%lld\n",(ans*powM(savn,y)%mod+mod)%mod); for (int i=1;i<=tot;i++)savc[i]=0; tot=0; } int main() { BNL::Init(); int T;scanf("%d",&T); while(T--)solve(); return 0; } ```