BZOJ 3601

BZOJ 3601

文章来自我的新博客

Description：

~~~~~~~给定一个非负整数 d,(d<=100)d , (d<=100) 和正整数 NN ,令 fd(N)f_d(N) 为所有小于 NN 且与 NN 互质的数的 dd 次方和。 对于给定的 d,Nd,N 求 fd(N) Mod 109+7f_d(N)~Mod~10^9+7的值。

~~~~~~~由于 NN 很大，所以给出 NN 的质因数分解式。

N=∏i=1wpaii(2<=pi<=109 , 1<=ai<=109)N=\prod_{i=1}^{w} p_i^{a_i}(2<=p_i<=10^9~,~1<=a_i<=10^9)

~~~~~~~

~~~~~~~

~~~~~~~

Solution：

~~~~~~~莫比乌斯反演之后有： fd(N)=∑k|Nμ(k)∑Nki=1(i∗k)d=∑k|Nμ(k)∗kd∑Nki=1idf_d(N)=\sum_{k|N} \mu(k) \sum_{i=1}^{\frac{N}{k}} (i*k)^d=\sum_{k|N} \mu(k)*k^d \sum_{i=1}^{\frac{N}{k}} i^d

~~~~~~~显然，可以构造数组 {a}\{a\} 使得 ∑Ti=1id=∑d+1i=0ai∗Td\sum_{i=1}^{T} i^d=\sum_{i=0}^{d+1} a_i*T^d

~~~~~~~那么假设我们构造出了数组 {a}\{a\} ，那么 fd(N)=∑k|Nμ(k)∗kd∑d+1i=0ai∗(Nk)i=∑k|Nμ(k)∑d+1i=0ai∗Ni∗kd−i=∑d+1i=0ai∗Ni∑k|Nμ(k)∗kd−if_d(N)=\sum_{k|N} \mu(k)*k^d \sum_{i=0}^{d+1} a_i* (\frac{N}{k})^i=\sum_{k|N} \mu(k) \sum_{i=0}^{d+1} a_i* N^i*k^{d-i}=\sum_{i=0}^{d+1} a_i* N^i \sum_{k|N} \mu(k)*k^{d-i}

~~~~~~~μ(k)∗kd−i\mu(k)*k^{d-i} 显然这是个积性函数， 所以 ∑k|Nμ(k)∗kd−i=∏wi=1(1+μ(pi)∗pd−ii)=∏wi=1(1−pd−ii)\sum_{k|N} \mu(k)*k^{d-i}=\prod_{i=1}^{w} (1+\mu(p_i)*p_i^{d-i})=\prod_{i=1}^{w} (1-p_i^{d-i})

~~~~~~~所以 fd(N)=∑d+1i=0ai∗Ni∏wi=1(1−∗pd−ii)f_d(N)=\sum_{i=0}^{d+1} a_i* N^i \prod_{i=1}^{w} (1-*p_i^{d-i}) ，这样可以轻松做到 O(d∗w)O(d*w) 的时间复杂度。但是 {a}\{a\} 怎么求呢，这很简单，高斯消元即可。所以总时间复杂度 O(d∗w+d3)O(d*w+d^3)。

~~~~~~~

~~~~~~~

~~~~~~~

Code：

#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <algorithm>
#include <iostream>

using namespace std;

const long long Mod=1e9+7;

long long fc[110][110]={{0}};
long long a[110]={0};
int d,W;
long long p_N[110]={0};
int D[1010]={0};
long long temp[1010]={0};

long long power(long long a,long long k)
{
    long long o=1;
    if(k<0) k+=Mod-1;
    for(;k>0;k>>=1)
    {
        if(k&1) o=o*a%Mod;
        a=a*a%Mod;
    }
    return o;
}

void Gauss()
{
    for(int i=1;i<=d+2;i++)
    {
        int g=i;
        for(int j=i;j<=d+2;j++)
            if(fc[j][i]!=0)
            {
                g=j;
                break;
            }
        if(g!=i) swap(fc[i],fc[g]);
        for(int j=1;j<=d+2;j++)
            if(j!=i && fc[j][i]!=0)
            {
                long long tmp=fc[j][i]*power(fc[i][i],Mod-2)%Mod;
                for(int p=0;p<=d+2;p++)
                    fc[j][p]=(fc[j][p]-tmp*fc[i][p]%Mod+Mod)%Mod;
            }
    }
    for(int i=1;i<=d+2;i++)
        a[i-1]=fc[i][0]*power(fc[i][i],Mod-2)%Mod;
    return;
}

long long N=1;

int main()
{
    scanf("%d%d",&d,&W);
    long long sum=0;
    for(int i=1;i<=d+2;i++)
    {
        sum=(sum+power(i,d))%Mod;
        fc[i][0]=sum;
        for(int j=1;j<=d+2;j++)
            fc[i][j]=power(i,j-1);
    }
    Gauss();
    for(int i=1;i<=W;i++)
    {
        int p,q;
        scanf("%d%d",&p,&q);
        N=N*power(p,q)%Mod;
        D[i]=p;
        temp[i]=power(p,Mod-2);
    }
    p_N[0]=1;
    for(int i=1;i<=d+1;i++)
        p_N[i]=p_N[i-1]*N%Mod;
    long long ans=0;
    for(int i=d+1;i>=0;i--)
    {
        long long cnt=1;
        for(int j=1;j<=W;j++)
        {
            cnt=cnt*(1-temp[j]+Mod)%Mod;
            temp[j]=temp[j]*D[j]%Mod;
        }
        ans=(ans+a[i]*p_N[i]%Mod*cnt)%Mod;
    }
    cout<<ans<<endl;
    return 0;
}