【XSY2721】求和 杜教筛

题目描述

  设 n=apii ,那么定义 fd(n)=(1)pi[pid] 。特别的, f1(n)=μ(n)

  给你 n,k ,求

i=1nj=1nd=1kfd(gcd(i,j))

   n1010,k40

题解

  先做一些简单的处理

ans=i=1nj=1nd=1kfd(gcd(i,j))=i=1nd=1kfd(i)(2(j=1niφ(j))1)

  后面 φ 用杜教筛可以在 O(n23) 内搞出来。

  设 λ(n)=f(n)=(1)pi

  考虑容斥,有:

fd(n)=λ(n)id|nμ(i)

Fd(n)=i=1nfd(i)=i=1nλ(i)jd|iμ(j)=i=1nμ(i)j=1nid+1λ(id+1j)=i=1nd+1λd+1(i)μ(i)Λ(nid+1)

   n107 的部分预处理,其他的每次枚举。这部分每次枚举是 n 的。总的是 O(n23) 的。(和杜教筛的分析方法一样。)
j|iλ(j)i=1nj|iλ(j)n=i=1nj[j|i]λ(i)Λ(n)=[i]=n=ij=1nj=1nijλ(j)=i=1nΛ(ni)=ni=2nΛ(ni)

  后面 Λ(n) 用杜教筛可以在 O(n23) 内搞出来

  反正总的是 O(n23) 的就对了。。。

  时间复杂度: O(n23)

代码

#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long ll;
int getsqrt(ll n)
{
    int l=1,r=1000000;
    while(l<r)
    {
        int mid=(l+r+1)>>1;
        if((ll)mid*mid>n)
            r=mid-1;
        else
            l=mid;
    }
    return l;
}
ll n;
ll _sqrt;
namespace hashset
{
    int getnum(ll x)
    {
        return n/x;
    }
}
using hashset::getnum;
int miu[10000010];
int phi[10000010];
int c[10000010];
int cs[10000010];
const int maxn=10000000;
int b[10000010];
int pri[1000010];
int cnt;
int d[10000010];
int e[10000010];
int f[10000010];
int k;
int vis[10000010];
int qp[10000010];
int qc[10000010];
void init()
{
    c[1]=phi[1]=miu[1]=f[1]=e[1]=1;
    d[1]=f[1]=0;
    int i,j;
    for(i=2;i<=maxn;i++)
    {
        if(!b[i])
        {
            miu[i]=-1;
            phi[i]=i-1;
            c[i]=-1;
            pri[++cnt]=i;
            d[i]=e[i]=1;
            f[i]=1;
        }
        for(j=1;j<=cnt&&i*pri[j]<=maxn;j++)
        {
            b[i*pri[j]]=1;
            c[i*pri[j]]=-c[i];
            f[i*pri[j]]=f[i]+1;
            if(i%pri[j]==0)
            {
                miu[i*pri[j]]=0;
                phi[i*pri[j]]=phi[i]*pri[j];
                d[i*pri[j]]=d[i]+1;
                e[i*pri[j]]=max(d[i*pri[j]],e[i]);
                break;
            }
            d[i*pri[j]]=1;
            e[i*pri[j]]=e[i];
            miu[i*pri[j]]=-miu[i];
            phi[i*pri[j]]=phi[i]*phi[pri[j]];
        }
    }
    for(i=1;i<=maxn;i++)
    {
        if(e[i]>k)
            f[i]=0;
        else
            f[i]=(f[i]&1?-1:1)*(k-e[i]+1);
        f[i]+=f[i-1];
//      miu[i]+=miu[i-1];
        phi[i]+=phi[i-1];
        cs[i]=cs[i-1]+c[i];
    }
}
int getphi(ll n)
{
    if(n<=maxn)
        return phi[n];
    int x=getnum(n);
    if(vis[x]&1)
        return qp[x];
    vis[x]|=1;
    ll i,j;
    int s=n*(n+1)>>1;
    for(i=2;i<=n;i=j+1)
    {
        j=n/(n/i);
        s-=(j-i+1)*getphi(n/i);
    }
    qp[x]=s;
    return s;
}
int getc(ll n)
{
    if(n<=maxn)
        return cs[n];
    int x=getnum(n);
    if(vis[x]&2)
        return qc[x];   
    vis[x]|=2;
    int s=getsqrt(n);
    ll i,j;
    for(i=2;i<=n;i=j+1)
    {
        j=n/(n/i);
        s-=(j-i+1)*getc(n/i);
    }
    qc[x]=s;
    return s;
}
ll pw[1000010];
int pw2[1000010];
int pw3[1000010];
int getfd(ll n)
{
    if(n<=maxn)
        return f[n];
    int i,j;
    for(i=1;(ll)i*i<=n;i++)
    {
        pw[i]=i;
        pw2[i]=pw3[i]=c[i];
    }
    int m=i-1;
    int s=0;
    for(j=1;j<=k;j++)
    {
        for(i=1;i<=m;i++)
        {
            pw[i]*=i;
            if(pw[i]>n)
                break;
            pw2[i]*=pw3[i];
            s+=miu[i]*pw2[i]*getc(n/pw[i]);
        }
    }
    return s;
}
int main()
{
#ifndef ONLINE_JUDGE
    freopen("b.in","r",stdin);
    freopen("b.out","w",stdout);
#endif
    scanf("%lld%d",&n,&k);
    _sqrt=getsqrt(n);
    init();
    int s=0;
    ll i,j;
    int now,last=0;
    int ans=0;
    for(i=1;i<=n;i=j+1)
    {
        j=n/(n/i);
        now=getfd(j);
        ans+=(now-last)*(2*getphi(n/i)-1);
        last=now;
    }
    ans&=(1<<30)-1;
    printf("%d\n",ans);
    return 0;
}
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值