BZOJ 3601

BZOJ 3601

文章来自我的新博客

Description：

$~~~~~~~$给定一个非负整数 $d , (d<=100)$ 和正整数 $N$ ,令 $f_d(N)$ 为所有小于 $N$ 且与 $N$ 互质的数的 $d$ 次方和。 对于给定的 $d,N$ 求 $f_d(N)~Mod~10^9+7$的值。
$~~~~~~~$由于 $N$ 很大，所以给出 $N$ 的质因数分解式。

$$N=\prod_{i=1}^{w} p_i^{a_i}(2<=p_i<=10^9~,~1<=a_i<=10^9)$$

$~~~~~~~$
$~~~~~~~$
$~~~~~~~$

Solution：

$~~~~~~~$莫比乌斯反演之后有： $f_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\}$ 使得 $\sum_{i=1}^{T} i^d=\sum_{i=0}^{d+1} a_i*T^d$

$~~~~~~~$那么假设我们构造出了数组 $\{a\}$ ，那么 $f_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}$

$~~~~~~~$$\mu(k)*k^{d-i}$ 显然这是个积性函数， 所以 $\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})$

$~~~~~~~$所以 $f_d(N)=\sum_{i=0}^{d+1} a_i* N^i \prod_{i=1}^{w} (1-*p_i^{d-i})$ ，这样可以轻松做到 $O(d*w)$ 的时间复杂度。但是 $\{a\}$ 怎么求呢，这很简单，高斯消元即可。所以总时间复杂度 $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;
}