[JZOJ5261]【GDOI2018模拟8.12】求和

Description

这里写图片描述
这里写图片描述

Solution

里面随便变换一下
原式化为

d=1ki=1nfd(i)x=1niy=1ni[gcd(x,y)=1]

后面的式子对于每个数单独考虑,原式可以化为

d=1ki=1nfd(i)((2x=1niφ(x))1)

后面的式子可以杜教筛弄出来

考虑前面
根据容斥原理有
Fd(n)=i=1nfd(i),λ(n)=f(n)

Fd(n)=i=1nμ(i)j=1nid+1λ(id+1j)

这是本题的关键,解释一下这个式子
一开始i=1,那就是所有的总和, μ(1)=1
对于i=2,应该将所有2的指数至少是d+1的数的贡献去掉,枚举它的倍数,因此 μ(2)=1
对于i=3与i=2同理

但是会减重复,因为 i=6 的情况被减了两次,又要加回来
因此 μ(6)=1

此处 μ 刚好是容斥系数

继续
可以发现 λ 是完全积性函数
Λ(n)=i=1nλ(i)

原式 λ(id+1) 提出

Fd(n)=i=1nμ(i)λ(id+1)Λ(nid+1)

i明显到不了n,修改上界

Fd(n)=i=1n(1d+1)μ(i)λ(id+1)Λ(nid+1)

那前面很小,可以暴力枚举
考虑如何算 Λ

根据定义可以发现 i|nλ(i)=[n 为完全平方数 ] <script type="math/tex" id="MathJax-Element-4669">]</script>,n的每个质因子指数均为偶数

单独考虑某一个质因子,如果是奇数它的和刚好为0,偶数才有答案,全为偶数则为1

于是也可以杜教筛

Code

本人被卡常,TLE80分

#include <cstdio>
#include <algorithm>
#include <iostream>
#include <cmath>
#include <cstdlib>
#include <cstring>
#define fo(i,a,b) for(int i=a;i<=b;i++)
#define fod(i,a,b) for(int i=a;i>=b;i--)
#define N 10000005
#define mo 1073741824
#define mo1 1073741823
#define H 30000007
#define LL long long
using namespace std;
LL n,m;
int pr[N],l,mu[N],s[N],phi[N],sp[N],sum[N],f[N];
LL h[H][2],hs[H][2];
int hash(LL k)
{
    int k1=k%H;
    while(h[k1][0]!=0&&h[k1][0]!=k) (k1+=1)%=H;
    return k1;
}
int has(LL k)
{
    int k1=k%H;
    while(hs[k1][0]!=0&&hs[k1][0]!=k) (k1+=1)%=H;
    return k1;
}
bool bz[N];
void prp()
{
    s[1]=mu[1]=1;
    f[1]=sum[1]=sp[1]=1;
    phi[1]=1;
    fo(i,2,10000000)
    {
        if(!bz[i])
        {
            mu[i]=-1;
            pr[++l]=i;
            phi[i]=i-1;
            f[i]=-1;
        }
        for(int j=1;j<=l&&pr[j]*i<=10000000;j++)
        {
            bz[i*pr[j]]=1;
            f[i*pr[j]]=-f[i];
            if(i%pr[j]==0)
            {
                mu[i*pr[j]]=0;
                phi[i*pr[j]]=phi[i]*pr[j];
                break;
            }
            mu[i*pr[j]]=-mu[i];
            phi[i*pr[j]]=phi[i]*phi[pr[j]];
        }
        s[i]=s[i-1]+mu[i];
        sp[i]=(sp[i-1]+phi[i])%mo;
        sum[i]=sum[i-1]+f[i];
    }   
}
LL get(LL k)
{
    if(k<=10000000) return sp[k];
    int p=hash(k);
    if(h[p][0]!=0) return h[p][1];
    h[p][0]=k;
    h[p][1]=((k*(k+1))/2)&mo1;
    LL i=2;
    while(i<=k)
    {
        LL i1=k/(k/i);
        h[p][1]=(h[p][1]-(get(k/i)*(i1-i+1)&mo1)+mo)&mo1;
        i=i1+1;
    } 
    return h[p][1];
}
LL ksm(LL k,LL n)
{
    if(!n) return 1;
    LL s=ksm(k,n/2);
    return (n%2)?s*s*k:s*s;
}
LL fd(LL k)
{
    if(k<=10000000) return sum[k];
    int p=has(k);
    if(hs[p][0]==k) return hs[p][1];
    hs[p][0]=k,hs[p][1]=sqrt(k);
    LL d=2;
    while(d<=k)
    {
        LL d1=k/(k/d);
        hs[p][1]=(hs[p][1]-(fd(k/d)*(d1-d+1)&mo1)+mo)&mo1;
        d=d1+1;
    }
    return hs[p][1];
}
LL gt(LL n,LL d)
{
    LL s=0;
    int lm=(int)(pow(n,1.0/(d+1))+0.001);
    fo(i,1,lm)
    {
        LL v=1;
        if((d+1)%2) v*=f[i];
        s=(s+v*mu[i]*fd(n/ksm(i,d+1)))&mo1;
    }   
    return s;
}
int main()
{
    freopen("sum.in","r",stdin);
    //freopen("sum.out","w",stdout);
    cin>>n>>m;
    prp();
    LL ans=0;
    fo(ld,1,m)
    {
        LL p=1;
        while(p<=n)
        {
            LL m1=n/p,p1=n/m1,s1=0,i=1;
            (ans+=(LL)(gt(p1,ld)-gt(p-1,ld)+mo)%mo*(((LL)2*get(m1)-1+mo)%mo)+mo)%=mo;
            p=p1+1;
        }
    }
    printf("%lld\n",ans);
}
  • 0
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

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

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

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值