ICode9

精准搜索请尝试: 精确搜索
首页 > 其他分享> 文章详细

一个人的数论[莫比乌斯反演+拉格朗日插值]

2022-01-01 19:36:20  阅读:152  来源: 互联网

标签:插值 res sum mu 反演 int 莫比 void mod


看了半天巨神 \(\color{red}{C}\color{black}{Yjian}\) 的题解才终于想明白。

规定变量:

  • 题中\(d\to k\)
  • 题中\(n\to m\) (下文的 \(n\) 制作未知数

令 \(F(n)=\sum_{i=1}^ni^k\) , \(g(n)=\sum_{i=1}^ni^k[\gcd(n,i)=1]\) ,那么莫比乌斯反演,有

\[F(n)=\sum_{d|n}d^kg(\frac{n}{d})\Leftrightarrow g(n)=\sum_{d|n}\mu(d)d^kF(\frac{n}{d}) \]

且最终答案即为 \(g(m)\) 。

发现 \(F(n)\) 为一个 \(k+1\) 次多项式,即有 \(F(n)=\sum_{i=0}^{k+1}f_in^i\) 。把他代入上式。

\[g(n)=\sum_{d|n}\mu(d)d^k\sum_{i=0}^{k+1}f_i(\frac{n}{d})^i \]

\[g(n)=\sum_{i=0}^{k+1}f_in^i\sum_{d|n}\mu(d)d^{k-i} \]

\(\mu(d)d^{k-i}\) 在 \(d\) 做参数时为积性函数,所以将 \(n\) 分解为 \(\prod_{p_i}p_i^{a_i}\) 后,式子可以写成

\[g(n)=\sum_{i=0}^{k+1}f_in^i\prod_{p_j}\sum_{l=0}^{a_i}\mu(p_j^{l})p_j^{l(k-i)} \]

而 \(\mu(p_i^{a_i})\) 在 \(a_i>1\) 时取 \(0\) ,因此只考虑 \(l=0\) 与 \(l=1\) 即可。因此

\[g(n)=\sum_{i=0}^{k+1}f_in^i\prod_{p_j}(1-p_j^{k-i}) \]

题目给出的 \(m\) 就是以质因子的形式输入的。直接拉格朗日暴力插出多项式系数计算就好。

Talk is Code(?
#include<bits/stdc++.h>
using namespace std;

namespace IO{
    typedef long long LL;
    typedef double DB;
    #define int LL
    int read(){
        int x=0,f=0; char ch=getchar();
        while(ch>'9'||ch<'0'){ f|=(ch=='-'); ch=getchar(); }
        while(ch>='0'&&ch<='9'){ x=(x<<1)+(x<<3)+(ch^48); ch=getchar(); }
        return f?-x:x;
    } char output[50];
    void write(int x,char sp){
        int len=0;
        if(x<0) putchar('-'), x=-x;
        do{ output[len++]=x%10+'0'; x/=10; }while(x);
        for(int i=len-1;~i;i--) putchar(output[i]); putchar(sp);
    }
    void ckmax(int& x,int y){ x=x>y?x:y; }
    void ckmin(int& x,int y){ x=x<y?x:y; }
} using namespace IO;

const int NN=1010,mod=1e9+7;
int n,k,m,ans,p[NN];

namespace Mathematics{
    int len,x[NN],y[NN],a[NN],f[NN];
    int qpow(int a,int b,int res=1){
        for(;b;b>>=1,a=a*a%mod)
            if(b&1) res=res*a%mod;
        return res;
    }
    void poly_mul(int x,int y){
        ++len;
        for(int i=len;i;i--)
            a[i]=(a[i]*y+a[i-1]*x)%mod;
        a[0]=a[0]*y%mod;
    }
    void poly_pls(){
        for(int i=0;i<=len;i++)
            f[i]=(f[i]+a[i])%mod, a[i]=0;
        len=0; a[0]=1;
    }
    void numb_mul(int x){
        for(int i=0;i<=len;i++)
            a[i]=a[i]*x%mod;
    }
    void poly_calc(int k){
        a[0]=1;
        for(int mul,i=1;i<=k+2;i++){
            mul=qpow(y[i],mod-2);
            for(int j=1;j<=k+2;j++) if(i^j){
                poly_mul(1,mod-x[j]);
                mul=mul*(mod+x[i]-x[j])%mod;
            }
            numb_mul(qpow(mul,mod-2)); poly_pls();
        }
    }
} using namespace Mathematics;

signed main(){
    k=read(); n=read(); m=1;
    for(int i=1;i<=n;i++)
        p[i]=read(), m=m*qpow(p[i],read())%mod;
    for(int i=1;i<=k+2;i++)
        x[i]=i, y[i]=(y[i-1]+qpow(x[i],k))%mod;
    poly_calc(k);
    for(int tmp=m,i=1;i<=k+1;i++){
        int res=tmp*f[i]%mod;
        tmp=tmp*m%mod;
        for(int j=1;j<=n;j++)
            res=res*(1+mod-qpow(p[j],mod-1+k-i))%mod;
        ans=(ans+res)%mod;
    }
    write(ans,'\n');
    return 0;
}

标签:插值,res,sum,mu,反演,int,莫比,void,mod
来源: https://www.cnblogs.com/keeeeeeen/p/15755921.html

本站声明: 1. iCode9 技术分享网(下文简称本站)提供的所有内容,仅供技术学习、探讨和分享;
2. 关于本站的所有留言、评论、转载及引用,纯属内容发起人的个人观点,与本站观点和立场无关;
3. 关于本站的所有言论和文字,纯属内容发起人的个人观点,与本站观点和立场无关;
4. 本站文章均是网友提供,不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属;如您发现该文章侵犯了您的权益,可联系我们第一时间进行删除;
5. 本站为非盈利性的个人网站,所有内容不会用来进行牟利,也不会利用任何形式的广告来间接获益,纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

专注分享技术,共同学习,共同进步。侵权联系[81616952@qq.com]

Copyright (C)ICode9.com, All Rights Reserved.

ICode9版权所有