[BZOJ3601] 一个人的数论

[BZOJ3601] 一个人的数论

试题分析

即求：$$f_d(n)=\sum_{gcd(i,n)=1}^n i^d$$。
肯定要将\([gcd(i,n)=1]\)移到里面去，得到：$$f_d(n)=\sum_{i=1}^n [gcd(i,n)=1] i^d $$。
然后由于\(\sum_{d|n} \mu(d) =[n=1]\)，所以可以变形为$$f_d(n)=\sum_{i=1}^n \sum_{x|i,x|j} \mu(x) i^d$$
转而枚举\(x\)：$$f_d(n)=\sum_{x|n} \mu(x) x^d \sum_{i=1}^{\lfloor \frac{n}{x} \rfloor} i^d$$
现在我们需要解决的只有自然数幂求和的问题，也就是解决\(g_d(n)=\sum_{i=1}^n i^d\)
这里有一个神奇的结论：$$g_d(n)=\sum_{i=1}^n i^d=\sum_{i=1} a_i n^i$$
然后我们枚举\(1\to d+1\)列出\(d+1\)个方程以后高斯消元解得\(a_i\)，然后直接数论分块就可以了。

#include<iostream>
#include<cstring>
#include<cstdio>
#include<vector>
#include<algorithm>
 
using namespace std;
#define LL long long
 
inline LL read(){
    LL x=0,f=1; char c=getchar();
    for(;!isdigit(c);c=getchar()) if(c=='-') f=-1;
    for(;isdigit(c);c=getchar()) x=x*10+c-'0';
    return x*f;
}
const LL MAXN = 100010;
const LL INF = 2147483600;
const LL Mod = 1000000007;
 
LL D,W; LL P[MAXN+1],A[MAXN+1];
LL b[MAXN+1],c[MAXN+1],a[113][113];
 
inline LL Pow(LL A,LL B){
    A=(A%Mod+Mod)%Mod; //B=(B%Mod+Mod)%Mod;
    LL res=1; for(; B; B>>=1,A=A*A%Mod) if(B&1) res=res*A%Mod; return res;
}
inline void Gauss(LL n){
    for(LL i=1;i<=n;i++){
        LL r=0; for(LL j=i;j<=n;j++) if(a[j][i]) {r=j; break;}
        if(!r) continue; if(r!=i){
            for(LL j=i;j<=n;j++) swap(a[r][j],a[i][j]);
            swap(b[r],b[i]);
        }
        for(LL j=i+1;j<=n;j++){
            LL tmp=Pow(a[i][i],Mod-2)%Mod*a[j][i]%Mod;
            for(LL k=1;k<=n;k++) a[j][k]=(a[j][k]-tmp*a[i][k]%Mod+Mod)%Mod;
            b[j]=(b[j]-tmp*b[i]%Mod+Mod)%Mod;
        }
    }
    for(LL i=n;i>=1;i--){
        c[i]=b[i]*Pow(a[i][i],Mod-2)%Mod;
        for(LL j=i-1;j>=1;j--)
            b[j]=(b[j]-c[i]*a[j][i]%Mod+Mod)%Mod;
    } 
    return ;
}
 
int main(){
    //freopen(".in","r",stdin);
    //freopen(".out","w",stdout);
    D=read(),W=read(); 
    for(LL i=1;i<=W;i++) P[i]=read(),A[i]=read();
    for(LL i=1;i<=D+1;i++){
        for(LL j=1;j<=D+1;j++) a[i][j]=Pow(i,j)%Mod;
        for(LL j=1;j<=i;j++) b[i]+=Pow(j,D),b[i]%=Mod;
    } Gauss(D+1); LL ans=0;
    for(LL i=1;i<=D+1;i++){
        LL T=1LL;
        for(LL j=1;j<=W;j++)
            T=T*(Pow(P[j],A[j]*i)-Pow(P[j],D+(A[j]-1)*i)%Mod+Mod)%Mod;
        ans=(ans+c[i]*T%Mod)%Mod; ans%=Mod;
    } printf("%lld\n",ans);
    return 0;
}