[JZOJ5261]【GDOI2018模拟8.12】求和

Description

Solution

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

∑d=1k∑i=1nfd(i)∑x=1⌊ni⌋∑y=1⌊ni⌋[gcd(x,y)=1]

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

∑d=1k∑i=1nfd(i)((2∑x=1⌊ni⌋φ(x))−1)

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

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

Fd(n)=∑i=1nμ(i)∑j=1⌊nid+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);
}