bzoj 3601 一个人的数论 - 莫比乌斯反演 - 拉格朗日插值 - 高斯消元

题目大意：定义$f_k(n)=\sum_{1\le i\le n,(i,n)=1}i^k，求f_k(\prod_{i=1}^wp_i^{a_i})$，其中$k\le1e2,\ w\le1e3$。

题解（前半部分在胡扯可以跳过）：

首先互质不好限制，考虑从总数中减去。

$$f_k(n)=\sum_{i=1}^ni^k-\sum_{d|n,d>1}\sum_{1\le i\le n,(n,i)=d}i^k$$

$$f_k(n)=\sum_{i=1}^ni^k-\sum_{d|n,d>1}\sum_{1\le i\le \frac{n}{d}}(id)^k[(i,\frac{n}{d})==1]$$

$$f_k(n)=\sum_{i=1}^ni^k-\sum_{d|n,d>1}d^k\sum_{1\le i\le\frac{n}{d},(\frac{n}{d},i)}i^k$$

$$f_k(n)=\sum_{i=1}^ni^k-\sum_{d|n,d>1}d^k f_k(\frac{n}{d})$$

注意到涉及的状态只会是n的因数，那么对于n<=1e9是随便跑的，但这个题却不行。

首先要解决的第一个问题是如何求$\sum_{i=1}^n i^k$，其中给出n的因式分解；

把这个东西拉格朗日差值出来，高斯消元一下求出系数表示（反正只有100这么大）。

这样因为要对大质数取模所以直接把n取模即可O(wlgw)的时间对任意n求出答案。

但是没有办法处理后半部分。

向莫比乌斯反演方向思索：

$$f_k(n)=\sum_{i=1}^n i^k [(n,i)==1]=\sum_{i=1}^n i^k\sum_{d|(n,i)}\mu(d)$$

$$f_k(n)=\sum_{d|n}\mu(d)\sum_{d|i}i^k=\sum_{d|n}\mu(d)d^k\sum_{i=1}^{\frac{n}{d}}i^k$$

注意到后面的求和是一个k+1次多项式，显然可以高斯消元求出它的系数$\{a_{k+1}\}$，然后可以考虑对g的每一项求和；

$$f_k(n)=\sum_{i=0}^{k+1}a_i\sum_{d|n}\mu(d)d^k(\frac{n}{d})^i=\sum_{i=0}^{k+1}a_ig_i(n)$$

后面那个$g_i(n)$显然是一个积性函数啦！

因此枚举次数之后就可以对每一个$p_j^{a_j}$单独计算答案：

$g_i(p^c)=(p^c)^i-p^k(p^{c-1})^i$，然后乘起来即可。

显然g函数可以log时间求出，这样问题再$O(k^3+kwlgc)$级别内做出这个题。

代码：

#include<iostream>
#include<cstring>
#include<cstdio>
#include<algorithm>
#define N 110
#define K N
#define W 1010
#define gc getchar()
#define lint long long
#define inv(x) fast_pow(x,mod-2)
#define mod 1000000007
#define debug(x) cerr<<#x<<"="<<x
#define sp <<" "
#define ln <<endl
using namespace std;
inline int inn()
{
    int x,ch;while((ch=gc)<'0'||ch>'9');
    x=ch^'0';while((ch=gc)>='0'&&ch<='9')
        x=(x<<1)+(x<<3)+(ch^'0');return x;
}
int k,w,a[N][N],b[N],p[W],c[W];
inline int fast_pow(int x,int k,int ans=1)
{   for(;k;k>>=1,x=(lint)x*x%mod) (k&1)?ans=(lint)ans*x%mod:0;return ans;   }
inline int f(int n)
{
    int ans=0;
    for(int i=1;i<=n;i++)
        (ans+=fast_pow(i,k))%=mod;
    return ans;
}
inline int g(int i,int p,int c)
{
    return (fast_pow(fast_pow(p,c),i)-(lint)fast_pow(p,k)*fast_pow(fast_pow(p,c-1),i)%mod+mod)%mod;
}
inline int Gauss(int n)
{
    for(int i=1;i<=n;i++)
    {
        int x=inv(a[i][i]);b[i]=(lint)b[i]*x%mod;
        for(int j=i;j<=n;j++) a[i][j]=(lint)a[i][j]*x%mod;
        for(int j=1,y;j<=n;j++)
            if(i^j)
            {
                y=a[j][i],b[j]=(b[j]-(lint)y*b[i]%mod+mod)%mod;
                for(int k=i;k<=n;k++)
                    a[j][k]=(a[j][k]-(lint)y*a[i][k]%mod+mod)%mod;
            }
    }
    return 0;
}
int main()
{
    k=inn(),w=inn();
    for(int i=1;i<=w;i++) p[i]=inn(),c[i]=inn();
    for(int i=1;i<=k+2;b[i]=f(i),i++)
        for(int j=1;j<=k+2;j++)
            a[i][j]=fast_pow(i,j-1);
    Gauss(k+2);int ans=0;
    for(int i=0;i<=k+1;i++)
    {
        int res=b[i+1];
        for(int j=1;j<=w;j++)
            res=(lint)res*g(i,p[j],c[j])%mod;
        (ans+=res)%=mod;
    }
    return !printf("%d\n",ans);
}