【XSY2731】Div 數論 杜教篩 莫比烏斯反演
題目大意
定義復數\(a+bi\)為整數\(k\)的約數，當且僅當\(a\)和\(b\)為整數且存在整數\(c\)和\(d\)滿足\((a+bi)(c+di)=k\)。
定義復數\(a+bi\)的實部為\(a\)，虛部為\(b\)。
定義\(f(n)\)為整數\(n\)的所有實部大於\(0\)的約數的實部之和。
給定正整數\(n\)，求出\(\sum_{i=1}^nf(i)\)對\(1004535809\)取模後得到的值。
\(n\leq {10}^{10}\)
題解
以前看到一個數論題就是反演預處理。
現在看到一個數論題就是反演杜教篩。
記\(s(n)=\sum_{i|n}i\)
當\(b=0\)時答案就是\(S(n)\)。以下僅考慮\(b>0\)的情況（\(b<0\)也是一樣的）
設\(n=(a+bi)(c+di)\)，那麽
\[
\begin{cases}
ac-bd&=n\ad+bc&=0
\end{cases}\\frac{a}{b}=-\frac{c}{d}
\]
因為這是一道數論題，設
\[
\begin{align}
a&=px\b&=qx\c&=py\d&=-qy\\gcd(p,q)&=1\\end{align}
\]
這樣一組\(x,y,p,q\)就唯一確定了一組\(a,b,c,d\)
記
\[
\begin{align}
g(n)&=\sum_{p^2+q^2=n}[\gcd(p,q)=1]p\G(n)&=\sum_{i=1}^ng(i)\f(n)&=\sum_{p^2+q^2=n}p\F(n)&=\sum_{i=1}^nf(i)
\end{align}
\]
問題轉化為求
\[
\begin{align}
&\sum_{x,y,p,q>0,[(p,q)=1]}[xy(p^2+q^2)\leq n]px\=&\sum_{i=1}^n(\sum_{p^2+q^2=i}[\gcd(p,q)=1]p)(\sum_{xy\leq\lfloor\frac{n}{i}\rfloor}x)\=&\sum_{i=1}^ng(i)S(\lfloor\frac{n}{i}\rfloor)
\end{align}
\]
那麽怎麽求\(F,G,S\)呢？
\[ \begin{align} S(n)&=\sum_{i=1}^n\sum_{j|i}j\&=\sum_{i=1}^ni\lfloor\frac{n}{i}\rfloor\G(n)&=\sum_{p^2+q^2\leq n}p\&=\sum_{i=1}^\sqrt ni\lfloor\sqrt{n-i^2}\rfloor\F(n)&=\sum_{p^2+q^2\leq n}[\gcd(i,j)=1]p\&=\sum_{i=1}^\sqrt{n}i\mu(i)\sum_{j=1}^\frac{n}{i^2}\sum_{p^2+q^2\leq\frac{n}{i^2}}p\&=\sum_{i=1}^\sqrt{n}i\mu(i)G(\lfloor\frac{n}{i^2}\rfloor) \end{align} \]
這些東西求一次是\(O(\sqrt{n})\)的，預處理一下，總的復雜度是\(O(n^\frac{2}{3})\)，因為每個\(n\)都是題目給的\(n\)除以某個東西。
預處理大家都會，我就不講了。
zjt：在\(O(n^\frac{2}{3})\)內求出所有\(F(n\)除以某個東西\()\)的一類算法都叫杜教篩。
時間復雜度：\(O(n^\frac{2}{3})\)
代碼
#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;
}【XSY2731】Div 數論 杜教篩 莫比烏斯反演