[hdu6833]A Very Easy Math Problem(莫比乌斯反演)

本文深入解析了一种复杂的数学问题求解算法,通过数论分块和反演技巧,针对特定的数学表达式进行高效计算。算法首先通过引入gcd的概念简化原始问题,随后利用反演原理进一步优化计算过程,最终通过数论分块实现O(nlogn)的时间复杂度。文章详细介绍了算法的推导过程,并提供了完整的代码实现。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >

题意:

给定正整数T,k,xT,k,xT,k,x,接下来TTT组样例,每组样例给出一个正整数nnn,对于每个nnn,要求出
∑a1=1n∑a2=1n…∑ax=1n(∏j=1xajk)f(gcd(a1,a2,…,ax))gcd(a1,a2,…,ax) \sum_{a_{1}=1}^{n}\sum_{a_{2}=1}^{n}…\sum_{a_{x}=1}^{n}(\prod_{j=1}^{x}a_{j}^{k})f(gcd(a_{1},a_{2},…,a_{x}))gcd(a_{1},a_{2},…,a_{x}) a1=1na2=1nax=1n(j=1xajk)f(gcd(a1,a2,,ax))gcd(a1,a2,,ax)
其中f(x)f(x)f(x)的定义是,当xxx存在平方因子时,f(x)=0f(x)=0f(x)=0,否则f(x)=1f(x)=1f(x)=1
(1≤T≤104,1≤k,x≤109,1≤n≤2×105)(1\leq T \leq 10^{4},1\leq k,x \leq 10^{9},1\leq n \leq 2\times 10^{5})(1T104,1k,x109,1n2×105)

推导:

套路提出gcdgcdgcd,变成
∑g=1n∑a1=1n∑a2=1n…∑ax=1n(∏j=1xajk)f(g)g[g=gcd(a1,a2,…,ax)] \sum_{g=1}^{n}\sum_{a_{1}=1}^{n}\sum_{a_{2}=1}^{n}…\sum_{a_{x}=1}^{n}(\prod_{j=1}^{x}a_{j}^{k})f(g)g[g=gcd(a_{1},a_{2},…,a_{x})] g=1na1=1na2=1nax=1n(j=1xajk)f(g)g[g=gcd(a1,a2,,ax)]
然后调整一下式子,
∑g=1ng[g不存在平方因子]∑a1=1n∑a2=1n…∑ax=1n(∏j=1xajk)[g=gcd(a1,a2,…,ax)] \sum_{g=1}^{n}g[g不存在平方因子]\sum_{a_{1}=1}^{n}\sum_{a_{2}=1}^{n}…\sum_{a_{x}=1}^{n}(\prod_{j=1}^{x}a_{j}^{k})[g=gcd(a_{1},a_{2},…,a_{x})] g=1ng[g]a1=1na2=1nax=1n(j=1xajk)[g=gcd(a1,a2,,ax)]
我们令f(g)f(g)f(g)
f(g)=∑a1=1n∑a2=1n…∑ax=1n(∏j=1xajk)[g=gcd(a1,a2,…,ax)] f(g)=\sum_{a_{1}=1}^{n}\sum_{a_{2}=1}^{n}…\sum_{a_{x}=1}^{n}(\prod_{j=1}^{x}a_{j}^{k})[g=gcd(a_{1},a_{2},…,a_{x})] f(g)=a1=1na2=1nax=1n(j=1xajk)[g=gcd(a1,a2,,ax)]
然后F(g)=∑g∣df(d)F(g)=\sum_{g|d}f(d)F(g)=gdf(d)
F(g)=∑g∣df(d)=∑a1=1n∑a2=1n…∑ax=1n(∏j=1xajk)[g∣gcd(a1,a2,…,ax)]=∑g∣a1∑g∣a2…∑g∣ax(∏j=1xajk)=∑a1=1⌊ng⌋∑a2=1⌊ng⌋…∑ax=1⌊ng⌋(∏j=1xajkgk)=∑a1=1⌊ng⌋a1kgk∑a2=1⌊ng⌋a2kgk…∑ax=1⌊ng⌋axkgk=gxk(∑a=1⌊ng⌋ak)x F(g)=\sum_{g|d}f(d)=\sum_{a_{1}=1}^{n}\sum_{a_{2}=1}^{n}…\sum_{a_{x}=1}^{n}(\prod_{j=1}^{x}a_{j}^{k})[g|gcd(a_{1},a_{2},…,a_{x})]\\ =\sum_{g|a_{1}}\sum_{g|a_{2}}…\sum_{g|a_{x}}(\prod_{j=1}^{x}a_{j}^{k})\\ =\sum_{a_{1}=1}^{\lfloor \frac{n}{g} \rfloor}\sum_{a_{2}=1}^{\lfloor \frac{n}{g} \rfloor}…\sum_{a_{x}=1}^{\lfloor \frac{n}{g} \rfloor}(\prod_{j=1}^{x}a_{j}^{k}g^{k})\\ =\sum_{a_{1}=1}^{\lfloor \frac{n}{g} \rfloor}a_{1}^{k}g^{k}\sum_{a_{2}=1}^{\lfloor \frac{n}{g} \rfloor}a_{2}^{k}g^{k}…\sum_{a_{x}=1}^{\lfloor \frac{n}{g} \rfloor}a_{x}^{k}g^{k}\\ =g^{xk}(\sum_{a=1}^{\lfloor \frac{n}{g} \rfloor}a^{k})^{x} F(g)=gdf(d)=a1=1na2=1nax=1n(j=1xajk)[ggcd(a1,a2,,ax)]=ga1ga2gax(j=1xajk)=a1=1gna2=1gnax=1gn(j=1xajkgk)=a1=1gna1kgka2=1gna2kgkax=1gnaxkgk=gxk(a=1gnak)x
F(g)F(g)F(g)式子化简出来后,还是比较简单的,先留着备用。

然后对着原来的式子套反演,
∑g=1ng[g没有平方因子]f(g)=∑g=1n[g没有平方因子]∑g∣dμ(dg)F(d)=∑g=1ng[g没有平方因子]∑g∣dμ(dg)dxk(∑a=1⌊nd⌋ak)x \sum_{g=1}^{n}g[g没有平方因子]f(g)=\sum_{g=1}^{n}[g没有平方因子]\sum_{g|d}\mu(\frac{d}{g}) F(d) \\ =\sum_{g=1}^{n}g[g没有平方因子]\sum_{g|d}\mu(\frac{d}{g})d^{xk}(\sum_{a=1}^{\lfloor \frac{n}{d} \rfloor}a^{k})^{x} g=1ng[g]f(g)=g=1n[g]gdμ(gd)F(d)=g=1ng[g]gdμ(gd)dxk(a=1dnak)x
现在已经可以预处理(∑a=1⌊nd⌋ak)x(\sum_{a=1}^{\lfloor \frac{n}{d} \rfloor}a^{k})^{x}(a=1dnak)x然后o(nlogn)o(nlogn)o(nlogn)做了,但是这样复杂度不够,我们考虑数论分块。

我们用一些函数把里面的一些东西代掉,让它看起来更优雅一些。

我们用is(t)is(t)is(t)表示ttt有没有平方因子,有的时候为000,否则为111

再用A(t)=(∑a=1tak)xA(t)=(\sum_{a=1}^{t}a^{k})^{x}A(t)=(a=1tak)x

原式变成,
∑g=1nis(g)g∑g∣dμ(dg)dxkA(⌊nd⌋) \sum_{g=1}^{n}is(g)g\sum_{g|d}\mu(\frac{d}{g})d^{xk}A(\lfloor \frac{n}{d} \rfloor) g=1nis(g)ggdμ(gd)dxkA(dn)
现在的式子还是不能数论分块。注意,我们现在是枚举ggg,再用ddd枚举ggg的倍数。我们转化成枚举ddd,再枚举ddd的因子ggg的形式,就变成,
∑d=1ndxkA(⌊nd⌋)∑g∣dμ(dg)is(g)g \sum_{d=1}^{n}d^{xk}A(\lfloor \frac{n}{d} \rfloor)\sum_{g|d}\mu(\frac{d}{g})is(g)g d=1ndxkA(dn)gdμ(gd)is(g)g
然后令B(d)=dxk∑g∣dμ(dg)is(g)gB(d)=d^{xk}\sum_{g|d}\mu(\frac{d}{g})is(g)gB(d)=dxkgdμ(gd)is(g)g,也是可以o(nlogn)o(nlogn)o(nlogn)预处理的,式子最后变成,
∑d=1nA(⌊nd⌋)B(d) \sum_{d=1}^{n}A(\lfloor \frac{n}{d} \rfloor)B(d) d=1nA(dn)B(d)
啊!这不就是某种显然的数论分块吗。
所以我们最后只要预处理出B(d)B(d)B(d)的前缀和SB(d)SB(d)SB(d),再预处理出A(d)A(d)A(d),然后数论分块一下就好了。
代码常数巨大。

代码:

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll mod=1000000007;
const ll MAXN=400005;

ll fpow(ll a,ll n)
{
    ll sum=1,base=a%mod;
    while(n!=0)
    {
        if(n%2)sum=sum*base%mod;
        base=base*base%mod;
        n/=2;
    }
    return sum;
}
ll inv(ll a)
{
    return fpow(a,mod-2);
}

ll prime[MAXN+10],notPrime[MAXN+10],mu[MAXN+10],tot;
void initMu(ll n)
{
    notPrime[1]=mu[1]=1;
    for(ll i=2;i<=n;i++)
    {
        if(!notPrime[i])prime[tot++]=i,mu[i]=-1;
        for(ll j=0;j<tot&&i*prime[j]<=n;j++)
        {
            notPrime[i*prime[j]]=1;
            if(i%prime[j])mu[i*prime[j]]=-mu[i];
            else {mu[i*prime[j]]=0;break;}
        }
    }
}
ll T,k,x,n;

ll A[400005];
ll is[400005];
ll B[400005];
ll SB[400005];
int main()
{
    scanf("%lld%lld%lld",&T,&k,&x);
    initMu(200000);
    for(ll i=2;i*i<=200000;i++)
    {
        for(ll j=1;i*i*j<=200000;j++)is[i*i*j]=1;
    }
    for(ll g=1;g<=200000;g++)
    {
        for(ll d=g;d<=200000;d+=g)
        {
            B[d]=(B[d]+g*mu[d/g]*(1-is[g])%mod)%mod;
        }
    }
    for(ll d=1;d<=200000;d++)B[d]=B[d]*fpow(d,x*k)%mod;
    for(ll d=1;d<=200000;d++)SB[d]=(SB[d-1]+B[d])%mod;
    for(ll d=1;d<=200000;d++)
    {
        A[d]=(A[d-1]+fpow(d,k))%mod;
    }
    for(ll d=1;d<=200000;d++)A[d]=fpow(A[d],x);
    while(T--)
    {
        scanf("%lld",&n);
        ll ans=0;
        for(ll l=1,r;l<=n;l=r+1)
        {
            if(n/l)r=n/(n/l);
            else r=n;
            r=min(r,n);
            ans=(ans+A[n/l]*(SB[r]-SB[l-1])%mod)%mod;
        }
        ans=ans%mod;
        printf("%lld\n",(ans%mod+mod)%mod);
    }
    return 0;
}

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值