[数论记录]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)的预处理 。

不会的话请见 多项式计数杂谈

得到\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。

贡献=\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很小,暴力求就好,下面同理。

-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)

#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;
}