【XSY2731】Div 数论 杜教筛 莫比乌斯反演

题目大意

  定义复数a+bi为整数k的约数,当且仅当ab为整数且存在整数cd满足(a+bi)(c+di)=k

  定义复数a+bi的实部为a,虚部为b

  定义f(n)为整数n的所有实部大于0的约数的实部之和。

  给定正整数n,求出ni=1f(i)1004535809取模后得到的值。

  n1010

题解

  以前看到一个数论题就是反演预处理。

  现在看到一个数论题就是反演杜教筛。

  记s(n)=i|nin的因数和,S(n)=ni=1s(i)

  当b=0时答案就是S(n)。以下仅考虑b>0的情况(b<0也是一样的)

  设n=(a+bi)(c+di),那么

{acbdad+bc=n=0ab=cd

  因为这是一道数论题,设
abcdgcd(p,q)=px=qx=py=qy=1

  这样一组x,y,p,q就唯一确定了一组a,b,c,d

  记

g(n)G(n)f(n)F(n)=p2+q2=n[gcd(p,q)=1]p=i=1ng(i)=p2+q2=np=i=1nf(i)

  问题转化为求
==x,y,p,q>0,[(p,q)=1][xy(p2+q2)n]pxi=1n(p2+q2=i[gcd(p,q)=1]p)(xynix)i=1ng(i)S(ni)

  那么怎么求F,G,S呢?
S(n)G(n)F(n)=i=1nj|ij=i=1nini=p2+q2np=i=1nini2=p2+q2n[gcd(i,j)=1]p=i=1niμ(i)j=1ni2p2+q2ni2p=i=1niμ(i)G(ni2)

  这些东西求一次是O(n)的,预处理一下,总的复杂度是O(n23),因为每个n都是题目给的n除以某个东西。

  预处理大家都会,我就不讲了。

  zjt:在O(n23)内求出所有F(n除以某个东西)的一类算法都叫杜教筛。

  时间复杂度:O(n23)

代码

#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long ll;
const ll p=1004535809;
int _gcd[3500][3500];
int gcd(int a,int b)
{
    int &s=_gcd[a][b];
    if(~s)
        return s;
    if(!b)
        return s=a;
    return s=gcd(b,a%b);
}
const int maxn=10000000;
ll n,m;
int vis[10000010];
ll f[10000010];
int miu[10000010];
int b[10000010];
int pri[1000010];
int cnt;
ll g1[10000010];
ll g2[10000010];
ll s[10000010];
int c[10000010];
void init()
{
    int i,j;
    for(i=1;i*i<=maxn;i++)
        for(j=1;i*i+j*j<=maxn;j++)
        {
            if(gcd(i,j)==1)
                (f[i*i+j*j]+=i)%=p;
            (g1[i*i+j*j]+=i)%=p;
        }
    for(i=1;i<=maxn;i++)
    {
        f[i]=(f[i]+f[i-1])%p;
        g1[i]=(g1[i]+g1[i-1])%p;
    }
    miu[1]=1;
    c[1]=1;
    s[1]=1;
    for(i=2;i<=maxn;i++)
    {
        if(!b[i])
        {
            pri[++cnt]=i;
            miu[i]=-1;
            c[i]=i;
            s[i]=i+1;
        }
        for(j=1;j<=cnt&&i*pri[j]<=maxn;j++)
        {
            b[i*pri[j]]=1;
            if(i%pri[j]==0)
            {
                miu[i*pri[j]]=0;
                c[i*pri[j]]=c[i]*pri[j];
                if(c[i]==i)
                    s[i*pri[j]]=(s[i]*pri[j]+1)%p;
                else
                    s[i*pri[j]]=s[c[i*pri[j]]]*s[i/c[i]]%p;
                break;
            }
            miu[i*pri[j]]=-miu[i];
            c[i*pri[j]]=pri[j];
            s[i*pri[j]]=s[i]*(pri[j]+1)%p;
        }
    }
    for(i=1;i<=maxn;i++)
        s[i]=(s[i]+s[i-1])%p;
}
const ll inv2=502267905;
ll S(ll n)
{
    if(n<=maxn)
        return s[n];
    ll s=0,s2;
    ll i,j;
    for(i=1;i<=n;i=j+1)
    {
        j=n/(n/i);
        s=(s+(i+j)%p*(j-i+1)%p*inv2%p*((n/i)%p))%p;
    }
    return s;
}
ll G(ll n)
{
    if(n<=maxn)
        return g1[n];
    if(vis[m/n]&2)
        return g2[m/n];
    vis[m/n]|=2;
    ll i,s=0,j;
    for(i=1;i*i<=n;i++);
    j=i-1;
    for(i=1;i*i<=n;i++)
    {
        while(j*j>n-i*i)
            j--;
        s=(s+i*j)%p;
    }
    g2[m/n]=s;
    return s;
}
ll F(ll n)
{
    if(n<=maxn)
        return f[n];
    ll i,s=0;
    ll now;
    for(i=1;i*i<=n;i++)
        s=(s+miu[i]*i*G(n/i/i))%p;
    return s;
}
int main()
{
#ifndef ONLINE_JUDGE
    freopen("c.in","r",stdin);
    freopen("c.out","w",stdout);
#endif
    memset(_gcd,-1,sizeof _gcd);
    scanf("%lld",&n);
    m=n;
    init();
    ll ans=0;
    ll i,j;
    ll last=0,now;
    for(i=1;i<=n;i=j+1)
    {
        j=n/(n/i);
        now=F(j);
        ans=(ans+(now-last)%p*S(n/i))%p;
        last=now;
    }
    ans=(ans*2+S(n))%p;
    ans=(ans+p)%p;
    printf("%lld\n",ans);
    return 0;
}
发布了172 篇原创文章 · 获赞 26 · 访问量 6万+
展开阅读全文

没有更多推荐了,返回首页

©️2019 CSDN 皮肤主题: 编程工作室 设计师: CSDN官方博客

分享到微信朋友圈

×

扫一扫,手机浏览