1. 程式人生 > >HDU 6340 2018HDU多校賽 第四場 Delightful Formulas(莫比烏斯反演+伯努利數+NTT+積性)

HDU 6340 2018HDU多校賽 第四場 Delightful Formulas(莫比烏斯反演+伯努利數+NTT+積性)

大致題意:給你k和m,還有n分解質因子之後的質因子及其對應的指數,讓你求 \sum_{i=1}^{N}[gcd(i,N)==1]\sum_{j=1}^{i}j^k

首先,這種含有gcd的式子,第一步肯定是進行莫比烏斯反演,這裡由於前面好幾篇都由類似的反演形式,所以我就不展開了,直接就得出反演之後的結果:

                                            \large f(1)=ans=\sum_{d|N}^{N}\mu(d)*\sum_{i=1}^{\frac{N}{d}}\sum_{j=1}^{i*d}j^k

對於最右邊的式子 \sum_{j=1}^{i*d}j^k,我們把i*d看作定值,這就是關於i*d的一個k次冪和。對於這個k次冪和,我們可以用伯努利數進行展開。有公式:\sum_{i=1}^{N}i^k=\frac{1}{k+1}\sum_{i=1}^{k+1}C_{k+1}^{i}*B_i*N^{k+1-i},即 \sum_{i=1}^{N}i^k 一定是一個k階多項式,那麼可以改寫成這樣的形式,而這個多項式的係數可以證明與伯努利數有關。於是我們可以確定那麼我們令a_j=\frac{1}{k+1}C_{k+1}^{k+1-j}*B_{k+1-j},那麼上面的式子可以寫成:

                                           \large f(1)=ans=\sum_{d|N}^{N}\mu(d)*\sum_{i=1}^{\frac{N}{d}}\sum_{j=1}^{k+1}a_j*(i*d)^j

交換一下求和次序:          \large f(1)=ans=\sum_{d|N}^{N}\mu(d)*\sum_{j=1}^{k+1}d^j*a_j*\sum_{i=1}^{\frac{N}{d}}i^j

我們發現,這個最右邊的東西如果把j看作是不變的,那麼它還是一個冪和的形式,於是我們考慮再次用伯努利數對它進行展開。我們令 b_i=\frac{1}{j+1}C_{j+1}^{j+1-i}*B_{j+1-i},那麼上面的式子可以寫成:

                                          \large f(1)=ans=\sum_{d|N}^{N}\mu(d)*\sum_{j=1}^{k+1}d^j*a_j*\sum_{i=1}^{j+1}b_i*(\frac{N}{d})^i

接下來,用一個特別騷的替換操作,把式子簡化。不妨設p=j-i,那麼p的取值範圍是[-1,k],那麼原本的i就可以替換成 j-p,然後再交換求和次序整理一下,那麼原式可以寫成:

                                         \large f(1)=ans=\sum_{d|N}^{N}\mu(d)*\sum_{p=-1}^{k}d^p*\sum_{j=1}^{k+1}a_j*b_{j-p}*N^{j-p}

我們不妨令\large g(p)=\sum_{j=1}^{k+1}a_j*b_{j-p}*N^{j-p},對於這個g(p)我們很高興的發現,它是一個卷積的形式,因此我們可以用NTT在 \large O(k\log k) 的時間複雜度內預處理出 g(p) ,那麼現在,原式就是這樣的:

                                       \large f(1)=ans=\sum_{d|N}^{N}\mu(d)*\sum_{p=-1}^{k}d^p*g(p)

再次交換求和次序:      \large f(1)=ans=\sum_{p=-1}^{k}g(p)*\sum_{d|N}^{N}\mu(d)*d^p

注意到 \large d^p 是冪和的某一項,根據 dls 的論文,我們知道 \large d^p 是具有積性的。然後莫比烏斯函式 \large \mu(d) 也是一個積性函式,因此這兩個東西對應項相乘也是一個積性函式。於是可以用積性的性質去優化這個過程。不妨令h(N)=\sum_{d|N}^{N}\mu(d)*d^p,由於具有積性,所以h(N)=h(prime_1^{a_1})*h(prime_2^{a_2})*...*h(prime_m^{a_m}),因此考慮每一個質因子即可。              

                                      \large h(prime_i^{a_i})=\sum_{j=0}^{ai}\mu(prime^j)*prime^p

注意到,根據莫比烏斯函式的性質,當自變數是某一個質數的2次方及以上的時候,其函式值為0,那麼只有當 j==1或2的時候式子才有意義。那麼我們就可以很自然的寫出 h(prime_i^{a_i}) 的通項:

                                       \large h(prime_i^{a_i})=1-prime_i^p

那麼最後的答案就是: \large f(1)=ans=\sum_{p=-1}^{k}g(p)*\prod_{i}^{m}(1-prime_i^p)

終於,我們推導完畢。我們發現,這個式子是O(KM)的,很愉快的可以解出這道題目。

最後呢,關於這個具體做法,還有一些細節對於第一次接觸這種題的人來說需要交代一下。

a_j=\frac{1}{k+1}C_{k+1}^{k+1-j}*B_{k+1-j}                               b_{j-p}=\frac{1}{j+1}C_{j+1}^{p+1}*B_{p+1}=\frac{j!*B_{p+1}}{(j-p)!*(p+1)!}

但我們的式子是 \large g(p)=\sum_{j=1}^{k+1}a_j*b_{j-p}*N^{j-p} ,實際的卷積形式應該是\large s(n)=\sum_{i=1}^{n}A(n-j)*B(j),我們正好反過來了,因此實際用的時候要把下標統一,再反回去。首先的話把右邊兩項整理一下,變成c_{j-p}=\frac{j!*N^{j-p}}{(j-p)!}\frac{B_{p+1}}{(p+1)!} 。後面哪一個只和p有關所以不用放到NTT中去求。然後是把下標換成 a_{k+1-j} 以及 c_{k+2-j}。那麼有:

                                         \large a_{k+1-j}=\frac{1}{k+1}C_{k+1}^j*B_j=\frac{k!*B_j}{j!*(k+1-j)!}

                                         \large c_{k+2-j}=\frac{j!*N^{k+2-j}}{(k+2-j)!}

然後你發現,當兩個相乘的時候,一個是要除以 j! 一個是要乘以 j!,因此這兩個可以抵消,所以我們又可以少幾項。這樣,我們列舉k+1-j和k+2-j,分別構造出a和c,然後這兩個做一個NTT的卷積運算,最後每一項再乘以一個\frac{B_{p+1}}{(p+1)!},就可以得到 \large g(p) 。最後再積性搞搞即可。 

由於這個k可以到1e5這個級別,因此暴力的O(N^2)預處理伯努利數是不行的,因此還得加上一個多項式求逆。這個也是板子咯,反正還要用到NTT的卷積運算。具體見程式碼:

#include<bits/stdc++.h>
#define file(x) freopen(#x".in","r",stdin),freopen(#x".out","w",stdout)
#define IO ios::sync_with_stdio(0);cin.tie(0);cout.tie(0);
#define mod 998244353
#define LL long long
#define N 550010
using namespace std;

int inv[N],fac[N],ifac[N],b[N],tmp[N];
int w,c[N],d[N],p[N],a[N],pw[N];

int qpow(int a,int b)
{
    int ans=1;
    while(b)
    {
        if(b&1)ans=(LL)ans*a%mod;
        a=(LL)a*a%mod; b>>=1;
    }
    return ans;
}

namespace NTT
{
    #define g 3
    int x[N<<2],y[N<<2],wn[N<<2];

    void init()
    {
        for(int i=0;i<21;i++)
        {
            int t=1<<i; wn[i]=qpow(g,(mod-1)/t);
        }
    }

    void brc(int *F,int len)
    {
        int j=len/2;
        for(int i=1;i<len-1;i++)
        {
            if(i<j)swap(F[i],F[j]);
            int k=len/2;
            while(j>=k) j-=k,k>>=1;
            if(j<k)j+=k;
        }
    }

    void NTT(int *F,int len,int t)
    {
        int id=0; brc(F,len);
        for(int h=2;h<=len;h<<=1)                   
        {
            id++;
            for(int j=0;j<len;j+=h)
            {
                int E=1;
                for(int k=j;k<j+h/2;k++)
                {
                    int u=F[k],v=(LL)E*F[k+h/2]%mod;
                    F[k]=(u+v)%mod,F[k+h/2]=((u-v)%mod+mod)%mod;
                    E=(LL)E*wn[id]%mod;
                }
            }
        }
        if(t==-1)
        {
            for(int i=1;i<len/2;i++)swap(F[i],F[len-i]);
            LL inv=qpow(len,mod-2);
            for(int i=0;i<len;i++)F[i]=(LL)F[i]%mod*inv%mod;
        }
    }

    void multiply(int *a,int len1,int *b,int len2)
    {
        int len=1;
        while(len<len1+len2)len<<=1;
        for (int i = len1; i < len; i++) a[i] = 0;
        for (int i = len2; i < len; i++) b[i] = 0;
        NTT(a,len,1); NTT(b,len,1);
        for(int i=0;i<len;i++) a[i]=(LL)a[i]*b[i]%mod;
        NTT(a,len,-1);
    }
}

void get_inv(int A[],int A0[],int t)
{
    if(t==1) 
    {
        A0[0]=qpow(A[0],mod-2);
        return;
    }
    get_inv(A,A0,t/2);
    for(int i=0;i<t;i++) tmp[i]=A[i];
    for(int i=t;i<2*t;i++) tmp[i]=0;
    for(int i=t/2;i<2*t;i++) A0[i]=0;
    NTT::NTT(tmp,t<<1,1); NTT::NTT(A0,t<<1,1);
    for(int i=0;i<2*t;i++)
    {
        tmp[i]=(2-1LL*tmp[i]*A0[i]%mod)%mod;
        if(tmp[i]<0) tmp[i]+=mod;
        A0[i]=1LL*A0[i]*tmp[i]%mod;
    }
    NTT::NTT(A0,t<<1,-1);
}

void init()
{
    fac[0]=ifac[0]=inv[0]=1;
    fac[1]=ifac[1]=inv[1]=1;
    for(int i=2;i<N;i++)
    {
        fac[i]=fac[i-1]*(LL)i%mod;
        inv[i]=(mod-mod/i)*(LL)inv[mod%i]%mod;
        ifac[i]=ifac[i-1]*(LL)inv[i]%mod;
    }
    for(int i=0;i<N-1;i++) a[i]=ifac[i+1];
    int len=1<<17; NTT::init(); get_inv(a,b,len);
    for(int i=0;i<len;i++) 
        b[i]=b[i]*(LL)fac[i]%mod;
    b[1]=mod-b[1];
}

int main()
{
    // file(g);
    IO; init();
    int T; cin>>T;
    while(T--)
    {
        int k,w,n=1;
        cin>>k>>w;
        int len=1;
        while(len<k*2+10) len<<=1;
        for(int i=0;i<=len;i++) c[i]=d[i]=0;
        for(int i=1;i<=w;i++)
        {
            cin>>p[i]>>a[i];
            n=n*(LL)qpow(p[i],a[i])%mod;
        }
        //cout<<n<<endl;
        for(int i=0;i<=k;i++)
            c[k+1-i]=fac[k]*(LL)ifac[i]%mod*b[i]%mod;
        int pow=1;
        for(int i=1;i<=k+2;i++)
        {
            pow=pow*(LL)n%mod;
            d[k+2-i]=pow*(LL)ifac[i]%mod;
        }
        NTT::multiply(c,len,d,len);
        for(int i=1;i<=w;i++)
            pw[i]=qpow(p[i],mod-2);
        int ans=0;
        for(int d=-1;d<=k;d++)
        {
            int gg=c[d+k+2]*(LL)b[d+1]%mod*ifac[d+1]%mod;
            for(int i=1;i<=w;i++)
            {
                gg=gg*(LL)(1-pw[i]+mod)%mod;
                pw[i]=pw[i]*(LL)p[i]%mod;
            }
            ans=(ans+gg)%mod;
        }
        cout<<ans<<endl;
    }
    return 0;
}