洛谷P3768 简单的数学题

本文详细解析了一道涉及杜教筛算法的题目,通过数学推导简化了原问题,并最终给出了具体的实现代码。

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

简述

  这是道比较有意思的题。
  首先用莫比乌斯反演推出答案:
  令g(x)=(x(x+1)2)2g(x)=(x(x+1)2)2
  

ans=x=1nx2g(nx)d|xdμ(xd)ans=∑x=1nx2g(⌊nx⌋)∑d|xdμ(xd)

  后面那个卷积直接线性筛,然后直接做是O(n)O(n)的,能过6060分。良心出题人部分分给的好足啊。
  注意到一个事情,前面的那个g(nx)g(⌊nx⌋)的取值是根号的,因此继续化下式子。
  
ans=x=1ng(nx)x2d|xdμ(xd)ans=∑x=1ng(⌊nx⌋)x2∑d|xdμ(xd)

  令
f(x)=x2d|xdμ(xd)f(x)=x2∑d|xdμ(xd)

  显然我们要预处理这个东西的前缀和才能做这道题,而且肯定是需要低于线性的算法(那就很可能是杜教筛了)。
  这个东西很难看出什么门道,我们显然需要用狄利克雷卷积的运算律给他化简下。
  令h(i)=ih(i)=i
  
f(x)=x2(hμ)(x)f(x)=x2(h∗μ)(x)

  为了好看先把x2x2拿掉,我们看到里面有个μμ,于是非常好奇地卷个11
  
hμ1=h(μ1)=hϵ=h

  咦?
  也就是说(hμ)1=h(h∗μ)∗1=h
  这不就是φφ吗,φ1=hφ∗1=h
  太神奇了(是我知道的太少了)
  所以
f(x)=x2φ(x)f(x)=x2φ(x)

  这个东西怎么杜教筛啊,卷啥呀
  冥思苦想了半天,还是搞出来了,令h2(x)=x2h2(x)=x2
  那么
(fh2)(x)=d|xf(x)h2(xd)=d|xd2φ(d)x2d2=x2d|xφ(x)=x3(f∗h2)(x)=∑d|xf(x)h2(xd)=∑d|xd2φ(d)x2d2=x2∑d|xφ(x)=x3
     这样一个优美的结论着实令我惊讶。
                  ————引用自pty

  令s(n)=ni=1f(x)s(n)=∑i=1nf(x)
  则

h2(1)s(n)=i=1n(fh2)(i)i=2nh2(i)s(ni)s(n)=i=1ni3i=2nh2(i)s(ni)h2(1)s(n)=∑i=1n(f∗h2)(i)−∑i=2nh2(i)s(⌊ni⌋)s(n)=∑i=1ni3−∑i=2nh2(i)s(⌊ni⌋)

  然后你就不能挡我用杜教筛了。
  你问我复杂度?
  大概
O(N12×N12×23)=O(N56)O(N12×N12×23)=O(N56)

  确实比较玄学。

代码

#include <cstdio>
#include <algorithm>
#define maxn 4700000
#define ll long long
using namespace std;
ll N, p, prime[maxn], phi[maxn], f[maxn], s2[maxn], inv;
bool mark[maxn];
inline ll pow(ll a, ll b, ll p)
{
    ll ans=1, t=a;
    for(;b;b>>=1,t=t*t%p)if(b&1)ans=ans*t%p;
    return ans;
}
inline ll sqr(ll x){x%=p;return x*x%p;}
inline ll g(ll x){x%=p;return sqr((x*(1+x)>>1)%p);}
inline ll gets2(ll x)
{
    if(x<maxn)return s2[x];
    x%=p;
    return x*(x+1)%p*(2*x+1)%p*inv%p;
}
void init()
{
    ll i, j, d, x;
    phi[1]=1;
    for(i=2;i<maxn;i++)
    {
        if(!mark[i])prime[++*prime]=i, phi[i]=i-1;
        for(j=1;i*prime[j]<maxn;j++)
        {
            mark[i*prime[j]]=1;
            if(i%prime[j]==0){phi[i*prime[j]]=phi[i]*prime[j]%p;break;}
            phi[i*prime[j]]=phi[i]*phi[prime[j]]%p;
        }
    }
    for(i=1;i<maxn;i++)phi[i]=(phi[i]*sqr(i)%p+phi[i-1])%p;
    for(i=1;i<maxn;i++)s2[i]=(s2[i-1]+i*i)%p;
}
ll getf(ll n){return n<maxn?phi[n]:f[N/n];}
void calc(ll n)
{
    if(!n)return;
    ll i, last, t=N/n;
    if(n<maxn or f[t])return;
    f[t]=g(n);
    for(i=2;i<=n;i=last+1)
    {
        last=n/(n/i);
        calc(n/i);
        f[t]=(f[t]-(gets2(last)-gets2(i-1))*getf(n/i)%p)%p;
    }
}
void solve()
{
    ll ans=0, i, last, t;
    for(i=1;i<=N;i=last+1)
    {
        last=N/(N/i);
        calc(last);
        calc(i-1);
        ans=(ans+g(N/i)%p*((getf(last)-getf(i-1))%p)%p)%p;
    }
    printf("%lld",(ans+p)%p);
}
int main()
{
    ll ans=0, x;
    scanf("%lld%lld",&p,&N);
    inv=pow(6,p-2,p);
    init();
    solve();
    return 0;
}
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值