HDU 6340 Problem I. Delightful Formulas（莫比乌斯反演+数论+FFT）

Description

给出$n,K$，令$s_i=\sum\limits_{j\le i}j^K$，求$\sum\limits_{1\le i\le n}s_i[gcd(i,n)=1]$

Input

第一行一整数$T$表示用例组数，每组用例首先输入一整数$K$，之后输入一整数$m$表示$n$的素因子数量，最后输入$m$对整数$p_i,a_i$表示$n=\prod\limits_{i=1}^np_i^{a_i}$

$(K\le 10^5,\sum K\le 10^6,m\le 20,p_i,a_i\le 10^9,p_i<p_{i+1})$

Output

输出答案，结果模$10^9+7$

Sample Input

2
1
2
2 1
3 1
233
1
23333 1

Sample Output

16
32727388

Solution

对所求表达式反演有

$$\begin{array}{rcl} ans&=&\sum\limits_{i=1}^n[(i,n)=1]\sum\limits_{j=1}^ij^K\\ &=&\sum\limits_{i=1}^n\sum\limits_{d|(i,n)}\mu(d)S(i,K)\\ &=&\sum\limits_{d|n}\mu(d)\sum\limits_{i=1}^{\frac{n}{d}}S(id,K)\\ \end{array}$$
其中 $S(n,k)=\sum\limits_{i=1}^ni^k=\frac{1}{k+1}\sum\limits_{i=1}^{k+1}C_{k+1}^iB_{k+1-i}(n+1)^i$， $B_0=1,\sum\limits_{i=0}^nC_{n+1}^iB_i=0$为伯努利序列

预处理伯努利序列，由于$\sum\limits_{i=0}^nC_{n+1}^iB_i=(n+1)!\sum\limits_{i=0}^n\frac{B_i}{i!}\frac{1}{(n+1-i)!}=0$，而$B_0=1$，故多项式$f(x)=\sum\limits \frac{B_i}{i!}x^i$即为多项式$g(x)=\sum\limits \frac{(n+1-i)!}x^i$的逆，用多项式求逆即可$O(Klog^2K)$得到伯努利序列

考虑到$S(n,k)$为一个关于$n$的$k+1$阶多项式，设$S(n,k)=\sum\limits_{j=0}^{k+1}a_{k,j}n^j$，先求$a_{k,j}$，由于

$$\begin{array}{rcl} S(n,k)&=&\frac{1}{k+1}\sum\limits_{i=1}^{k+1}C_{k+1}^iB_{k+1-i}(n+1)^i\\ &=&\frac{1}{k+1}\sum\limits_{i=0}^{k+1}C_{k+1}^iB_{k+1-i}(n+1)^i-\frac{1}{k+1}B_{k+1}\\ &=&\frac{1}{k+1}\sum\limits_{i=0}^{k+1}C_{k+1}^iB_{k+1-i}\sum\limits_{j=0}^iC_i^jn^j-\frac{1}{k+1}B_{k+1}\\ &=&\sum\limits_{j=0}^{k+1}n^j\sum\limits_{i=j}^{k+1}C_{k+1}^iC_i^jB_{k+1-i}-\frac{1}{k+1}B_{k+1} \end{array}$$
显然有 $a_{k,0}=0$，当 $j\ge 1$时有
$$\begin{array}{rcl} a_{k,j}&=&\frac{1}{k+1}\sum\limits_{i=j}^{k+1}C_{k+1}^iC_i^jB_{k+1-i}\\ &=&\frac{1}{k+1}\sum\limits_{i=0}^{k+1-j}C_{k+1}^iC_{k+1-i}^jB_i\\ &=&\frac{1}{k+1}\sum\limits_{i=0}^{k+1-j}\frac{(k+1)!(k+1-i)!}{i!(k+1-i)!j!(k+1-i-j)!}B_i\\ &=&\frac{1}{k+1}\sum\limits_{i=0}^{k+1-j}\frac{(k+1)!(k+1-j)!}{i!(k+1-j)!j!(k+1-i-j)!}B_i\\ &=&\frac{1}{k+1}\sum\limits_{i=0}^{k+1-j}C_{k+1-j}^iC_{k+1}^jB_i\\ &=&\frac{1}{k+1}C_{k+1}^j(\sum\limits_{i=0}^{k-j}C_{k+1-j}^iB_i+B_{k+1-j})\\ &=& \left\{\begin{array}{lcl} \frac{1}{k+1}C_{k+1}^jB_{k+1-j}&&1\le j<k\\ B_0+B_1&&j=k\\ \frac{1}{k+1}&&j=k+1 \end{array}\right. \end{array}$$
注意到只要我们令 $B_1=B_0+B_1$，那么有 $a_{k,j}=\frac{1}{k+1}C_{k+1}^jB_{k+1-j},j\ge 0,a_{k,0}=0$

得到$S(n,k)$后我们有

$$\begin{array}{rcl} ans&=&\sum\limits_{d|n}\mu(d)\sum\limits_{i=1}^{\frac{n}{d}}S(id,K)\\ &=&\sum\limits_{d|n}\mu(d)\sum\limits_{i=1}^{\frac{n}{d}}\sum\limits_{j=0}^{K+1}a_{K,j}i^jd^j\\ &=&\sum\limits_{d|n}\mu(d)\sum\limits_{j=0}^{K+1}a_{K,j}d^j\sum\limits_{i=1}^{\frac{n}{d}}i^j\\ &=&\sum\limits_{d|n}\mu(d)\sum\limits_{j=0}^{K+1}a_{K,j}d^j\sum\limits_{k=0}^{j+1}a_{j,k}n^kd^{-k}\\ &=&\sum\limits_{d|n}\mu(d)\sum\limits_{i=-1}^{K+1}d^i\sum\limits_{j=0}^{K+1}\sum\limits_{k=0}^{j+1}[j-k=i]a_{K,j}a_{j,k}n^k\\ &=&\sum\limits_{d|n}\mu(d)\sum\limits_{i=-1}^{K+1}d^iF_{i+1}\\ &=&\sum\limits_{i=-1}^{K+1}F_{i+1}\sum\limits_{d|n}\mu(d)d^i\\ &=&\sum\limits_{i=-1}^{K+1}F_{i+1}\prod\limits_{j=1}^m(1-p_j^i)\\ \end{array}$$
其中 $F_i=\sum\limits_{j=0}^{K+1}\sum\limits_{k=0}^{j+1}[j-k=i-1]a_{K,j}a_{j,k}n^k=\sum\limits_{j=0}^{K+1}a_{K,j}a_{j,j-i+1}n^{j-i+1},0\le i\le K+2$

只要可以快速求出$F_{0},...,F_{K+2}$，则可以$O(mK)$的得到答案，最后考虑求$F_i$，由于$a_{K+1,0}=0$，故$F_{K+2}=0$

$$\begin{array}{rcl} F_i&=&\sum\limits_{j=0}^{K+1}a_{K,j}a_{j,j-i+1}n^{j-i+1}\\ &=&\sum\limits_{j=1}^{K+1}a_{K,j}a_{j,j-i+1}n^{j-i+1}\\ &=&\sum\limits_{j=0}^{K}a_{K,j+1}a_{j+1,j-i+2}n^{j-i+2}\\ &=&\sum\limits_{j=0}^K\frac{1}{K+1}C_{K+1}^{j+1}B_{K-j}\frac{1}{j+2}C_{j+2}^{i}B_in^{j-i+2}\\ &=&\frac{K!B_i}{i!}\sum\limits_{j=0}^K\frac{B_{K-j}}{(K-j)!}\frac{n^{j-i+2}}{(j-i+2)!} \end{array}$$
以上为一个卷积形式，用 $FFT$即可 $O(KlogK)$求出 $F_0,...,F_{K+1}$

Code

#include<cstdio>
#include<algorithm>
#include<cmath>
using namespace std;
typedef long long ll;
#define maxn 100005
#define maxfft 262144+5
#define mod 998244353
int add(int x,int y)
{
    x+=y;
    if(x>=mod)x-=mod;
    return x;
}
int mul(int x,int y)
{
    ll z=1ll*x*y;
    return z-z/mod*mod;
}
int Pow(int a,int b)
{
    int ans=1;
    while(b)
    {
        if(b&1)ans=mul(ans,a);
        a=mul(a,a);
        b>>=1;
    }
    return ans;
}
const double pi=acos(-1.0);
struct cp
{
    double a,b;
    cp operator +(const cp &o)const {return (cp){a+o.a,b+o.b};}
    cp operator -(const cp &o)const {return (cp){a-o.a,b-o.b};}
    cp operator *(const cp &o)const {return (cp){a*o.a-b*o.b,b*o.a+a*o.b};}
    cp operator *(const double &o)const {return (cp){a*o,b*o};}
    cp operator !() const{return (cp){a,-b};}
}w[maxfft];
int pos[maxfft];
void fft_init(int len)
{
    int j=0;
    while((1<<j)<len)j++;
    j--;
    for(int i=0;i<len;i++)
        pos[i]=pos[i>>1]>>1|((i&1)<<j);
}
void fft(cp *x,int len,int sta)
{
    for(int i=0;i<len;i++)
        if(i<pos[i])swap(x[i],x[pos[i]]);
    w[0]=(cp){1,0};
    for(int i=2;i<=len;i<<=1)
    {
        cp g=(cp){cos(2*pi/i),sin(2*pi/i)*sta};
        for(int j=i>>1;j>=0;j-=2)w[j]=w[j>>1];
        for(int j=1;j<i>>1;j+=2)w[j]=w[j-1]*g;
        for(int j=0;j<len;j+=i)
        {
            cp *a=x+j,*b=a+(i>>1);
            for(int l=0;l<i>>1;l++)
            {
                cp o=b[l]*w[l];
                b[l]=a[l]-o;
                a[l]=a[l]+o;
            }
        }
    }
    if(sta==-1)for(int i=0;i<len;i++)x[i].a/=len,x[i].b/=len;
}
cp x[maxfft],y[maxfft],z[maxfft];
void FFT(int *a,int *b,int n,int m,int *c)
{
    int len=1;
    while(len<n+m)len<<=1;
    fft_init(len);
    for(int i=0;i<len;i++)
    {
        int aa=i<n?a[i]:0,bb=i<m?b[i]:0;
        x[i]=(cp){(aa>>15),(aa&32767)},y[i]=(cp){(bb>>15),(bb&32767)};
    }
    fft(x,len,1),fft(y,len,1);
    for(int i=0;i<len;i++)
    {
        int j=len-1&len-i;
        z[i]=((x[i]+!x[j])*(y[i]-!y[j])+(x[i]-!x[j])*(y[i]+!y[j]))*(cp){0,-0.25};
    }
    fft(z,len,-1);
    for(int i=0;i<n+m-1;i++)
    {
        ll ta=(ll)(z[i].a+0.5)%mod;
        ta=(ta<<15)%mod;
        c[i]=ta;
    }
    for(int i=0;i<len;i++)
    {
        int j=len-1&len-i;
        z[i]=(x[i]-!x[j])*(y[i]-!y[j])*(cp){-0.25,0}+(x[i]+!x[j])*(y[i]+!y[j])*(cp){0,0.25};
    }
    fft(z,len,-1);
    for(int i=0;i<n+m-1;i++)
    {
        ll ta=(ll)(z[i].a+0.5)%mod,tb=(ll)(z[i].b+0.5)%mod;
        ta=(ta+(tb<<30))%mod;
        c[i]=(c[i]+ta)%mod;
    }
}
int temp1[maxfft];
void Poly_Inv(int *poly,int n,int *ans)
{
    ans[0]=Pow(poly[0],mod-2);
    for(int i=2;i<=n;i<<=1)
    {
        FFT(poly,ans,i,i/2,temp1);
        FFT(ans,temp1+i/2,i/2,i/2,temp1);
        for(int j=0;j<i/2;j++)ans[j+i/2]=temp1[j]==0?0:mod-temp1[j];
    }
}
int inv[maxn],fact[maxn],B[2*maxn];
void init(int n=100002)
{
    fact[0]=1;
    for(int i=1;i<=n;i++)fact[i]=mul(i,fact[i-1]);
    inv[1]=1;
    for(int i=2;i<=n;i++)inv[i]=mul(mod-mod/i,inv[mod%i]);
    inv[0]=1;
    for(int i=1;i<=n;i++)inv[i]=mul(inv[i-1],inv[i]);
    int len=1;
    while(len<n-1)len<<=1;
    Poly_Inv(inv+1,len,B);
    for(int i=0;i<n;i++)B[i]=mul(B[i],fact[i]);
    B[1]++;
}
int T,K,m,a[22],b[22],f[2*maxn],g[maxn],h[maxn]; 
int main()
{
    init();
    scanf("%d",&T);
    while(T--)
    {
        scanf("%d%d",&K,&m);
        int n=1;
        for(int i=1;i<=m;i++)
        {
            scanf("%d%d",&a[i],&b[i]);
            n=mul(n,Pow(a[i],b[i]));
        }
        for(int i=0;i<=K+2;i++)h[i]=1;
        for(int i=1;i<=m;i++)
        {
            int p=Pow(a[i],mod-2);
            for(int j=0;j<=K+2;j++)
            {
                h[j]=mul(h[j],1+mod-p);
                p=mul(p,a[i]);
            }
        }
        for(int i=0;i<=K;i++)f[i]=mul(B[K-i],inv[K-i]);
        int res=n;
        for(int i=0;i<=K+1;i++)
        {
            g[i]=mul(inv[i+1],res);
            res=mul(res,n);
        }
        reverse(g,g+K+2);
        FFT(f,g,K+1,K+2,f);
        for(int i=0;i<=K+1;i++)
            g[i]=mul(f[i+K],mul(fact[K],mul(inv[i],B[i])));
        int ans=0;
        for(int i=0;i<=K+1;i++)ans=add(ans,mul(g[i],h[i]));
        printf("%d\n",ans);
    }
    return 0;
}