题解 P4464 【[国家集训队]JZPKIL】

· · 题解

多次查询，每次给定 n,x,y，求：

\sum_{i=1}^n \gcd(i,n)^x\textrm{lcm}(i,n)^y \pmod {10^9+7} \rm Sol:

答案是：

\sum_{i=1}^n i^y\cdot n^y\gcd(i,n)^{x-y}

枚举 \gcd=d，忽略 n^y 则有：

&\sum_{d|n} d^x\sum_{j=1}^{\frac{n}{d}}[\gcd(j,\frac{n}{d})=1]j^{y} \\&=\sum_{d|n} d^x\sum_{j=1}^{n/d}\sum_{k|j,k|\frac{n}{d}} \mu(k)j^y \\&=\sum_{d|n} d^x\sum_{k|\frac{n}{d}}\mu(k) \sum_{k|j} j^y \\&=\sum_{d|n} d^x\sum_{k|\frac{n}{d}}\mu(k)k^y \sum_{j=1}^{n/(kd)} j^y \\&=\sum_{d|n} d^x \sum_{k|\frac{n}{d}}\mu(k)k^y S_y(\frac{n}{kd}) \end{aligned}

我们将 S_y(\frac{n}{kd}) 视为关于其的 y+1 次多项式，那么此处有：

&\sum_{d|n} d^x \sum_{k|\frac{n}{d}}\mu(k)k^y \sum_{i=0}^{y+1}f_i (\frac{n}{kd})^i \\&=\sum_{i=0}^{y+1}f_i\sum_{d|n} d^x \sum_{k|\frac{n}{d}}\mu(k)k^y(\frac{n}{kd})^i \end{aligned}

注意到后者可以视为 f,g,h 的迪利克雷卷积，令 f(d)=d^x,g(d)=\mu(d)d^y,h^c(d)=d^c，这个卷积满足分配律和交换律以及结合律，同时注意到三者均为积性函数，所以卷积结果也是积性函数。

设结果为 \Z(x)，于是我们所求为 \Z(n)，那么显然有：

  1. \Z(1)=1
  2. \Z(p)=p^x+p^c-p^y
  3. \Z(p^k)=f\cdot h^c(p^k)-f\cdot h^c(p^{k-1})\times p^y

于是只需要考虑 f\cdot h^c(p^k)，其为：

\sum_{i+j=k} f(p^i)\cdot h^c(p^j)

所以通过 MR 进行暴力质因数分解，然后我们 \mathcal O(\log n) 的枚举每个质因子（include 次幂）然后暴力卷积并计算答案即可。

还有一个问题，如何计算自然数幂的和的多项式？

自然数幂和多项式

构造多项式 G_n(x)=\frac{e^{nx}-1}{1-e^{-x}}，一堆神仙分析可以得到 G_n(x)[x^k] 就是 \sum_{i=1}^n i^k（其实应该是根据 G_n(x)[x^k] 是 \sum_{i=1}^n i^k 推出 G_n(x) 的表达式的...）

将 G_n(x) 裂项成 \frac{xe^x}{e^x-1}\times \frac{e^{nx}-1}{x}，前者的 EGF 序列记为 \{B_0,B_1,B_2...\}，后者显然是 EGF 序列 \{\frac{n}{1},\frac{n^2}{2},\frac{n^3}{3}...\}

所以我们会发现其实有 G_n(x)[x^k] 是这两个数列的二项卷积的结果，所以有：

&G_n[x^k]=\sum_{i+j=k} B_i \frac{n^{j+1}}{j+1}\binom{k}{i} \\&= \frac{1}{k+1}\sum_{i=0}^{k} B_i \binom{k+1}{i}n^{k-i+1} \end{aligned}

换而言之如果能够得到伯努利数那么就可以知道自然数幂和多项式

然后只需要递推伯努利数即可。

注意到 \frac{xe^x}{e^x-1}\times \frac{e^x-1}{x}=e^x 所以 \mathbf{EGF}\{B_0,B_1...\}\times \mathbf{EGF}\{\frac{1}{1},\frac{1}{2}...\}=\mathbf{EGF}\{1,1,1...\}

所以有：

&\sum_{i+j=k} B_i\times \frac{1}{j+1}\times \binom{k}{i}=1 \\& \sum_{i}^k B_i\times\binom{k+1}{i}=k+1 \end{aligned}

于是得到 B_k=1-\sum_{i<k} \frac{B_i}{k+1}\binom{k+1}{i}

综上，本题复杂度为 \mathcal O(T\cdot (n^{\frac{1}{4}}\log n+y\log^2 n)+y^2)

然后这道题卡常，你需要尽可能的省略快速幂，如果你是 \log^3 然后还是我这种大常数玩家就会死得非常惨。

然后千万不要写龟速乘，然后挑一个夜深人静的时候交一发，多洗一下脸，没准就过了。

不过我在 darkbzoj 这道题过得非常稳...但是咕咕就死得非常惨

Code: 
#include<bits/stdc++.h>
using namespace std ;
#define Next( i, x ) for( register int i = head[x]; i; i = e[i].next )
#define rep( i, s, t ) for( register int i = (s); i <= (t); ++ i )
#define drep( i, s, t ) for( register int i = (t); i >= (s); -- i )
#define re register
#define int __int128
int gi() {
    char cc = getchar() ; int cn = 0, flus = 1 ;
    while( cc < '0' || cc > '9' ) {  if( cc == '-' ) flus = - flus ; cc = getchar() ; }
    while( cc >= '0' && cc <= '9' )  cn = cn * 10 + cc - '0', cc = getchar() ;
    return cn * flus ;
}
const int M = 3000 + 5 ; 
const int P = 1e9 + 7 ;
int n, fx, fy, fc, maxn, top, cnt, num, tp[M], st[M], w[M], fac[M], inv[M], B[M], f[M] ; 
int mul( int a, int b, int p ) {
    return a * b % p ; 
}
int fpow( int x, int k, int p ) {
    int ans = 1, base = x ; 
    while(k) {
        if(k & 1) ans = mul(ans, base, p) ;
        base = mul(base, base, p), k >>= 1 ; 
    } return ans ; 
}
int sw( int x ) { return ( ( 1ll * rand() << 16ll ) | rand() ) % x + 1 ; }
bool M_R( int p ) {
    if( p == 2 || p == 3 ) return 1 ; 
    if( p == 1 || ( p % 2 == 0 ) ) return 0 ; 
    re int d = p - 1, s = 0 ; while( !( d & 1 ) ) ++ s, d >>= 1 ; 
    rep( i, 0, 5 ) {
        re int a = sw( p - 3 ) + 2 ; 
        re int x = fpow( a, d, p ), y = 0 ;
        for( register int j = 0; j < s; ++ j ) {
            y = mul( x, x, p ) ; if( y == 1 && ( x != 1 ) && x != ( p - 1 ) ) return 0 ;
            x = y ;
        } if( y != 1 ) return 0 ;
    } return 1 ; 
}
int gcd( int x, int y ) {
    return ( x == 0 ) ? y : gcd( y % x, x ) ;
}
int work( int p ) {
    re int k = 2, x = sw(p - 1), y = x, d = 1, c = sw(99997) ;
    for( re int i = 1; d == 1; ++ i ) {
        x = ( mul( x, x, p ) + c ) % p ;
        d = gcd( (x > y) ? x - y : y - x, p ) ;
        if( i == k ) y = x, k <<= 1 ; 
    } return d ; 
}
void Pollard_Rho(int p) {
    if( p == 1 ) return ; 
    if( M_R(p) ) { st[++ top] = p; return ; }
    int x = p ; while( x == p ) x = work(x) ;
    Pollard_Rho(x), Pollard_Rho(p / x) ;
}
int g[M], gg[M], fg[M], fgg[M], fff[M] ; 
int c( int x, int y ) {
    if( x < y ) return 0 ;
    return fac[x] * inv[y] % P * inv[x - y] % P ;  
}
signed main()
{
    srand(time(NULL)) ;
    int T = gi() ; maxn = 3002 ; B[0] = 1, fac[0] = inv[0] = 1 ; 
    rep( i, 1, maxn ) fac[i] = fac[i - 1] * i % P, inv[i] = fpow( fac[i], P - 2, P ) ;
    rep( i, 1, maxn ) {
        for( re int j = 0; j < i; ++ j ) B[i] = ( B[i] - c(i + 1, j) * B[j] % P + P ) % P ; 
        B[i] = B[i] * fpow( i + 1, P - 2, P ) % P, ++ B[i] ;  
    }
    while( T-- ) {
        n = gi(), fx = gi(), fy = gi() ; int k = fy ; top = 0 ; 
        for( re int i = 0; i <= k; ++ i ) f[k - i + 1] = B[i] * c(k + 1, i) % P * fpow( k + 1, P - 2, P ) % P ; 
        Pollard_Rho(n), sort(st + 1, st + top + 1) ; cnt = 0 ; 
        rep( i, 1, top ) {
            if( st[i] != st[i - 1] ) tp[++ cnt] = st[i] ; 
            ++ w[cnt] ; 
        } int Ans = 0 ; 
        rep( j, 1, cnt ) fg[j] = fpow( tp[j], fx, P ), fff[j] = fpow(tp[j], fy, P) ; 
        for( re int ic = 0; ic <= k + 1; ++ ic ) {
            int ans = 1 ;
            for( re int j = 1; j <= cnt; ++ j ) {
                int tk = w[j], qwq = 0, de = 0 ; 
                g[0] = 1, gg[0] = 1 ; fgg[j] = ( !ic ) ? 1 : fgg[j] * tp[j] % P ;  
                int dg = fg[j], dgg = fgg[j] ; 
                for( re int o = 1; o <= w[j]; ++ o ) 
                    g[o] = g[o - 1] * dg % P, gg[o] = gg[o - 1] * dgg % P ; 
                for( re int o = 0; o <= w[j]; ++ o )
                    qwq = ( qwq + g[o] * gg[(w[j] - o)] % P ) % P ; 
                for( re int o = 0; o < w[j]; ++ o ) 
                    de = ( de + g[o] * gg[(w[j] - o - 1)] % P ) % P ; 
                ans = ans * ( qwq - de * fff[j] % P + P ) % P ; 
            }
            Ans = ( Ans + ans * f[ic] % P + P ) % P ; 
        } rep( i, 1, cnt ) w[i] = tp[i] = 0 ;  
        cout << (long long)(Ans * fpow( n, fy, P ) % P) << endl ; 
    }
    return 0 ;
}