[国家集训队] JZPKIL

题目链接

洛谷：https://www.luogu.org/problemnew/show/P4464

Solution

这题是真的毒....数论大杂烩，窝断断续续写了两天。

众所周知：

\[{\rm lcm}(x,y)=\frac{xy}{\gcd(x,y)} \]

带进去，顺便枚举\(\gcd\)值套\(mobius\)反演公式然后换一下求和符号：

\[\begin{align} ans=&n^y\sum_{d|n}d^{x-y}\sum_{i=1}^{n}i^y[\gcd(i,n)=1]\\ =&n^y\sum_{d|n}d^x\sum_{i=1}^{n/d}i^y[\gcd(id,n)=1]\\ =&n^y\sum_{d|n}d^x\sum_{i=1}^{n/d}i^y\sum_{t|i,t|\frac{n}{d}}\mu(t)\\ =&n^y\sum_{d|n}d^x\sum_{t|\frac{n}{d}}\mu(t)t^y\sum_{i=1}^{n/dt}i^y\\ \end{align} \]

众所周知，自然数\(k\)次幂和可以表示为一个\(k+1\)次多项式，设多项式系数为\(s\)，即：

\[\sum_{i=1}^{n}i^y=\sum_{i=0}^{y+1}s_in^i \]

那么带进去顺便把\(s_i\)丢最前面：

\[\begin{align} ans=&n^y\sum_{d|n}d^x\sum_{t|\frac{n}{d}}\mu(t)t^y\sum_{i=1}^{y+1}s_i(\frac{n}{dt})^i\\ =&n^y\sum_{i=1}^{y+1}s_i\sum_{d|n}d^{x}\sum_{t|\frac{n}{d}}\mu(t)t^{y}(\frac{n}{dt})^i\\ \end{align} \]

可以注意到后面是一个狄利克雷卷积的形式，设：

\[f_i(n)=\sum_{d|n}d^{x}\sum_{t|\frac{n}{d}}\mu(t)t^{y}(\frac{n}{dt})^i\\ a(n)=d^n,b(n)=\mu(n)n^y,c_i(n)=n^i \]

显然：

\[f_i=a*b*c_i \]

也就是说\(f_i\)是一个积性函数，我们只需要关系\(f(p^w)\)的值就好了，然后用\(\rm pollard\ rho\)分解\(n\)，暴力算就好了。

显然根据\(\mu\)函数的性质，\(f_i\)的定义式中枚举的\(t\)只有等于\(1\ or \ p\)的时候才有贡献，我们分类讨论：

• \(t=1\)：

  \[f(n)=n^i\sum_{d|n}d^{x-i}\\ f(p^w)=p^{wi}\sum_{j=0}^{w}p^{j(x-i)} \]

• \(t=p\)：

  \[f(n)=-n^i\sum_{d|n}d^{x-i}p^{y-i}\\f(p^w)=-p^{wi+y-i}\sum_{j=0}^{w-1}p^{j(x-i)} \]

这里可以用等比数列求和，当然暴力枚举复杂度也是对的，本题复杂度瓶颈在质因数分解上，但是这题卡pollard rho常就很不友好。

那个自然数幂和的系数可以用伯努利数求，具体上百度，伯努利数直接大力预处理就好了。

然后就做完了，复杂度\(O(T\sqrt[4]{n}+Ty\log n+y^2)\)。

// #pragma GCC optimize(3)
#include<bits/stdc++.h>
using namespace std;

template<typename T> void read(T &x) {
    x=0;T f=1;char ch=getchar();
    for(;!isdigit(ch);ch=getchar()) if(ch=='-') f=-f;
    for(;isdigit(ch);ch=getchar()) x=x*10+ch-'0';x*=f;
}

void print(int x) {
    if(x<0) putchar('-'),x=-x;
    if(!x) return ;print(x/10),putchar(x%10+48);
}
void write(int x) {if(!x) putchar('0');else print(x);putchar('\n');}

#define lf double
#define ll long long 

#define pii pair<int,int >
#define vec vector<int >

#define pb push_back
#define mp make_pair
#define fr first
#define sc second

#define FOR(i,l,r) for(int i=l,i##_r=r;i<=i##_r;i++) 

const int maxn = 4e3+10;
const int N = 3e3+10;
const int inf = 1e9;
const lf eps = 1e-8;

ll r[maxn],cnt;

namespace Pollard_rho {
    const int pri[] = {0,2,3,5,7,11,13,23,31};

    ll gcd(ll a,ll b) {return !b?a:gcd(b,a%b);}
    
    ll mul(ll x,ll y,ll p) {
        ll t=(long double)x*y/p;
        ll res=x*y-t*p;res%=p;if(res<0) res+=p;
        return res;
    }
    
    ll qpow(ll a,ll x,ll p) {
        ll res=1;
        for(;x;x>>=1,a=mul(a,a,p)) if(x&1) res=mul(res,a,p);
        return res;
    }
    
    bool MR(ll n) {
        if(n<2) return 1;
        ll e=(n-1)>>__builtin_ctz(n-1);
        for(int i=1;i<=8;i++) {
            if(pri[i]>=n) break;
            for(ll d=e,lst=qpow(pri[i],d,n),now;d!=n-1;d<<=1,lst=now) {
                now=mul(lst,lst,n);
                if(lst!=1&&lst!=n-1&&now==1) return 0;
            }if(qpow(pri[i],n-1,n)!=1) return 0;
        }return 1;
    }
    
    ll calc(ll n,ll v) {
        ll a=rand()%n,b=a,g=1;
        while(g==1) {
            a=(mul(a,a,n)+v)%n;
            b=(mul(b,b,n)+v)%n;
            b=(mul(b,b,n)+v)%n;
            g=gcd(abs(a-b),n);
        }return g;
    }
    
    void PR(ll n) {
        if(MR(n)) return n<2?0:r[++cnt]=n,void();
        ll d=n;
        while(d==n) d=calc(n,1ull*rand()*rand()*rand()*rand()%(n-1)+1);
        PR(d),PR(n/d);
    }
    
    void divide(ll n) {
        for(int i=1;i<=8;i++) while(n%pri[i]==0) n/=pri[i],r[++cnt]=pri[i];
        PR(n);
    }
}

const int mod = 1e9+7;

int fac[maxn],ifac[maxn],inv[maxn],b[maxn];
    
int add(int x,int y) {return x+y>=mod?x+y-mod:x+y;}
int del(int x,int y) {return x-y<0?x-y+mod:x-y;}
int mul(int x,int y) {return 1ll*x*y-1ll*x*y/mod*mod;}

int c(int x,int y) {return mul(fac[x],mul(ifac[y],ifac[x-y]));}

int qpow(int a,int x) {
    int res=1,f=0;if(x<0) x=-x,f=1;
    for(;x;x>>=1,a=mul(a,a)) if(x&1) res=mul(res,a);
    return f?qpow(res,mod-2):res;
}

namespace Bernoulli {
    void sieve() {
        fac[0]=ifac[0]=inv[0]=inv[1]=1;
        for(int i=1;i<maxn;i++) fac[i]=mul(fac[i-1],i);
        for(int i=2;i<maxn;i++) inv[i]=mul(mod-mod/i,inv[mod%i]);
        for(int i=1;i<maxn;i++) ifac[i]=mul(ifac[i-1],inv[i]);
    }

    void get_B() {
        b[0]=1;
        for(int i=1;i<N;i++) {
            for(int j=0;j<i;j++) b[i]=add(b[i],mul(b[j],c(i+1,j)));
            b[i]=del(0,mul(b[i],inv[i+1]));
        }b[1]++;
    }
}

int x,y,s[maxn],tot;ll n;
struct data {ll p;int x;}a[maxn];

int calc(int p,int w,int i) {
    int res=0;
    for(int j=0;j<=w;j++) res=add(res,qpow(p,j*(x-i)+w*i));
    for(int j=0;j<=w-1;j++) res=del(res,qpow(p,j*(x-i)+w*i-i+y));
    return res;
}

void solve() {
    read(n),read(x),read(y);cnt=0;s[0]=0;tot=0;
    if(!n) return puts("0"),void();
    for(int i=0;i<=y;i++) s[y+1-i]=mul(inv[y+1],mul(b[i],c(y+1,i)));
    Pollard_rho :: divide(n);
    sort(r+1,r+cnt+1);
    for(int i=1;i<=cnt;i++)
        if(r[i]!=r[i-1]) a[++tot].p=r[i],a[tot].x=1;
        else a[tot].x++;
    int ans=0;
    for(int i=0;i<=y+1;i++) {
        int res=s[i];
        for(int j=1;j<=tot;j++) {
            int p=a[j].p%mod,w=a[j].x;
            res=mul(res,calc(p,w,i));
        }ans=add(ans,res);
    }ans=mul(ans,qpow(n%mod,y));
    write(ans);
    for(int i=0;i<=cnt+5;i++) a[i].p=a[i].x=r[i]=0;
    for(int i=0;i<=y+1+5;i++) s[i]=0;
}

int main() {
    srand(19280817);
    Bernoulli :: sieve(),Bernoulli :: get_B();
    int t;read(t);while(t--) solve();
    return 0;
}