题目大意
定义复数\(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\)为\(n\)的因数和，\(S(n)=\sum_{i=1}^ns(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#include #include 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;}