bzoj3601-一个人的数论

题目

定义:
\[ \begin{aligned} f_d(n)=\sum _{gcd(x,n)=1} x^d \end{aligned} \]
给定\(d\)\(n\) (\(n\)以质因数分解的形式给出),求\(f_d(n)\ (mod \ 10^9+7)\)

分析

本来是学莫比乌斯反演时候的题,拖到了高斯消元来做。
看到一个互质,即\(gcd(x,n)=1\),马上就能想到莫比乌斯函数的那个性质:
\[ 令 \ f(n)=\sum _{d|n} \mu(d) \\ 则有 \ f(n)= \left\{\begin{aligned} 1,n=1\\0,n>1\end{aligned} \right. \]
下面我们证明这个结论:
\(n=1\),结论显然。
\(n>1\),那么我们设
\[ \begin{aligned} s=\prod _{i=1} ^k p_i^{a_i},k|n \end{aligned} \]
当任意一个\(a\_i>1\)时,\(\mu(s)=0\)
所以上面的\(f(n)\)等价于\(f(s)\),其中\(s=\prod\_{i=1}^k p\_i\)
根据莫比乌斯函数的定义,当\(k\)为奇数时,\(f(s)=-1\)\(k\)为偶数时,\(f(s)=1\).
所以\(f(n)\)可以写成:
\[ \begin{aligned} f(n)=C_k^0-C_k^1+C_k^2....=\sum _{i=0}^k (-1)^i*C_k^i \end{aligned} \]
所以我们现在要证明的就是:
\[ \begin{aligned} \sum _{i=0}^k (-1)^i*C_k^i=0 \end{aligned} \]
注意到二项式定理:
\[ \begin{aligned} (x+y)^n=\sum _{i=0}^k C_k^i*x^i*y^{k-i} \end{aligned} \]
与我们要证明的结论的相似性,我们令\(x=-1\)\(y=1\),即可得证。

这样,我们可以对原式进行变形:
\[ \begin{aligned} f_d(n)&=\sum _{gcd(x,n)=1} x^d \\ &=\sum _{x=1}^n x^d[gcd(x,n)=1] \\ &=\sum _{x=1}^n x^d\sum _{c|gcd(x,n)}\mu(c) \\ &=\sum _{c|n} \sum _{x=1}^{\frac{n}{c}} (cx)^d\mu(c) \\ &=\sum _{c|n} \mu(c)c^d \sum _{x=1} ^{\frac{n}{c}} x^d \\ \end{aligned} \]
令:
\[ \begin{aligned} h(n)=\sum _{i=1}^n i^d \end{aligned} \]
我们称\(h\)为自然数幂和。
则有:
\[ \begin{aligned} f_d(n)=\sum _{c|n} \mu(c)c^d h(\frac{n}{c}) \\ \end{aligned} \]
这里有一个很奇妙的,并且正确的猜想:
\[ \begin{aligned} h(n)=\sum _{i=1} ^{d+1}a_in^i \end{aligned} \]
也就是说,自然数幂和可以用\(n\)的多项式表示出来。
于是奇妙地用到了高斯消元,直接用\(1\)\(d+1\)为例子,列出\(d+1\)个方程,强行把多项式的系数解出来。其实自然数幂和有通项公式,近期会有相关的文章。于是我们有:
\[ \begin{aligned} f_d(n)&=\sum _{c|n} \mu(c)c^d \sum _{i=1} ^{d+1}a_i(\frac{n}{c})^i \\ &=\sum _{i=1}^{d+1} a_i \sum _{c|n} \mu(c)c^d(\frac{n}{c})^i \end{aligned} \]
令:
\[ \begin{aligned} g_i(n)=\sum _{c|n} \mu(c)c^d(\frac{n}{c})^i \end{aligned} \]
那么有:
\[ \begin{aligned} f_d(n)=\sum _{i=1}^{d+1} a_i *g_i(n) \end{aligned} \]
现在我们可以在\(O(d^3)\)的时间内求出所有的\(a_i\),只要对每个\(i\)求出\(g_i(n)\),我们就能\(O(d)\)算出答案了。
若令\(r(c)=\mu(c)c^d\),那么有:
\[ \begin{aligned} g_i(n)=\sum _{c|n} r(c)(\frac{n}{c})^i \end{aligned} \]
这是一个狄利克雷卷积(即形如 \(t(n)=\sum _{d|n}f(d)g(\frac{n}{d})\)的函数)!
狄利克雷卷积函数有一个性质,如果原来两个函数都是积性函数,那么新函数也是积性函数。证明如下:
如果两个函数不完全积性,那么需要\(gcd(x,y)=1\),否则不需要。
\[ \begin{aligned} t(x)t(y)&=\sum _{d|x}f(d)g(\frac{x}{d}) \sum _{e|y}f(e)g(\frac{y}{e}) \\ &=\sum _{de|xy}f(de)g(\frac{xy}{de}) \\ &=t(xy) \end{aligned} \]
所以我们可以通过\(n\)的质因数分解求出每个\(g_i(n)\)
\[ \begin{aligned} g_i(n)=\prod _{j=1}^w g_i(p_j^{q_j}) \end{aligned} \]
而对于只包含一个质数的\(g\)函数:
\[ \begin{aligned} g_i(p^q) &=\sum _{j=0}^q \mu(p^j)p^{jd}(\frac{p^q}{p^j})i \\ &= \left\{\begin{aligned} p^{qi},j=0\\-p^{qi+d-i},j=1\\ 0, otherwise\end{aligned} \right. \end{aligned} \]
所以有:
\[ \begin{aligned} g_i(p^q) =p^{qi}-p^{qi+d-i} \end{aligned} \]
直接计算即可。还是很简单的嘛!

代码

公式推错了四次~

#include<cstdio>
#include<algorithm>
#include<cctype>
using namespace std;
typedef long long giant;
giant read() {
    giant 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 giant maxn=1e3+5;
const giant maxd=105;
const giant q=1e9+7;
giant p[maxn],s[maxn];
giant a[maxd][maxd],as[maxn];
giant mi(giant x,giant y) {
    giant ret=1;
    while (y) {
        if (y&1) (ret*=x)%=q;
        y>>=1,(x*=x)%=q;
    }
    return ret;
}
void elm(giant n) {
    for (giant i=1;i<n;++i) {
        if (!a[i][i]) {
            for (giant j=i+1;j<=n;++j) if (a[j][i]) {
                for (giant k=1;k<=n+1;++k) swap(a[j][k],a[i][k]);
                break;
            }   
        }   
        giant inv=mi(a[i][i],q-2);
        for (giant j=i+1;j<=n;++j) if (a[j][i]) {
            giant tmp=(a[j][i]*inv)%q;
            for (giant k=1;k<=n+1;++k) a[j][k]=((a[j][k]-a[i][k]*tmp+q)%q+q)%q;
        }
    }
    for (giant i=n;i>1;--i) {
        if (!a[i][i]) {
            for (giant j=i-1;j>0;--j) if (a[j][i]) {
                for (giant k=1;k<=n+1;++k) swap(a[j][k],a[i][k]);
                break;
            }   
        }   
        giant inv=mi(a[i][i],q-2);
        for (giant j=i-1;j>0;--j) if (a[j][i]) {
            giant tmp=(a[j][i]*inv)%q;
            for (giant k=1;k<=n+1;++k) a[j][k]=((a[j][k]-a[i][k]*tmp+q)%q+q)%q; 
        }
    }
}
int main() {
    #ifndef ONLINE_JUDGE
        freopen("test.in","r",stdin);
        freopen("my.out","w",stdout);
    #endif
    giant d=read(),w=read();
    for (giant i=1;i<=w;++i) p[i]=read(),s[i]=read();
    giant sum=0;
    for (giant i=1;i<=d+1;++i) {
        (sum+=mi(i,d))%=q;
        giant tmp=1;
        for (giant j=1;j<=d+1;++j) (tmp*=i)%=q,a[i][j]=tmp;
        a[i][d+2]=sum;
    }
    elm(d+1); 
    for (giant i=1;i<=d+1;++i) as[i]=(a[i][d+2]*mi(a[i][i],q-2))%q;
    giant ans=0;
    for (giant i=1;i<=d+1;++i) {
        giant g=1;
        for (giant j=1;j<=w;++j) {
            giant tmp=(mi(p[j],s[j]*i)-mi(p[j],s[j]*i+d-i)+q)%q;
            (g*=tmp)%=q;
        }
        (ans+=(as[i]*g))%=q;
    }
    printf("%lld\n",ans);
}

转载于:https://www.cnblogs.com/owenyu/p/6724695.html

  • 0
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值