LOJ6052 DIV
题面：LOJ
解析
题中所给条件即为：
[ac-bd=k
]
[ad+bc=0
]
那么，不妨设 (b gt 0)，(frac{a}{b}=-frac{c}{d}=frac{p}{q})，其中：((p,q)=1)。
那么有：(ac-bd=p^2frac{ac}{p^2}-q^2frac{bd}{q^2}=p^2frac{ac}{p^2}+q^2frac{ac}{p^2}=frac{ac}{p^2}(p^2+q^2)=xy(p^2+q^2)=k)，其中:(x=frac{a}{p})，(y=frac{c}{p})。
那么(b>0)部分的答案即为：
[=sum_{i=1}^{n}sum_{p^2+q^2|i,(p,q)=1}psum_{x|frac{i}{p^2+q^2}}x
]
[=sum_{i=1}^{n}sum_{p^2+q^2|i,(p,q)=1}p imessigma(frac{i}{p^2+q^2})
]
[=sum_{i=1}^{n}sigma(i)sum_{p^2+q^2leqlfloorfrac{n}{i}
floor}p[(p,q)=1]
]
[=sum_{i=1}^{n}sigma(i)sum_{p^2+q^2leqlfloorfrac{n}{i}
floor}psum_{d|p,d|q}mu(d)
]
[=sum_{d=1}^{lfloorsqrt{n}
floor}imu(i)F(frac{n}{i^2})
]
其中：(F(n)=sum_{i=1}^{n}sigma(i)sum_{p^2+q^2leqlfloorfrac{n}{i} floor}p)。
考虑计算(F(n))，对(lfloorfrac{n}{i}
floor)数论分块，有(sum_{i=1}^{n} sigma(i)=sum_{i=1}^{n}ilfloorfrac{n}{i}
floor)，
又有(G(n)=sum_{p^2+q^2leq n}p=sum_{i=1}^{lfloorsqrt{n}
floor}ilfloorsqrt{k-p^2}
floor)。
类似于杜教筛，预处理加记忆化以后复杂度为(O(n^{frac{2}{3}}))。
而(b<0)的部分与(b>0)的部分对称，答案乘(2)即可。
对于(b=0)，加上：(sum_{acleq n}a=sum_{i=1}^{n}sigma(i)=sum_{i=1}^{n}ilfloorfrac{n}{i} floor)即可。
代码
#include <cmath>
#include <cstdio>
typedef long long ll;
const int mod=1004535809;
const int inv2=502267905;
const int N=4e6+5;
ll n,m,no_pri[N],pri[N],pri_cnt,mu[N],minp[N],sigma[N],sum[N],ans;
void sieve(int n){
for(int i=1;i*i<=n;++i){
sum[i*i]=(sum[i*i]+i)%mod;
for(int j=1;i*i+j*j<=n;++j) sum[i*i+j*j]=(sum[i*i+j*j]+2*i)%mod;
}
no_pri[1]=mu[1]=minp[1]=sigma[1]=1;
for(int i=2;i<=n;++i){
if(!no_pri[i]){
pri[++pri_cnt]=i; mu[i]=-1; minp[i]=i; sigma[i]=i+1;
for(ll j=i;j*i<=n;j*=i) sigma[j*i]=sigma[j]*i+1;
}
for(int j=1;pri[j]*i<=n;++j){
no_pri[i*pri[j]]=1;
if(i%pri[j]==0){
minp[i*pri[j]]=minp[i]*pri[j];
sigma[i*pri[j]]=sigma[i/minp[i]]*(sigma[minp[i]]*pri[j]+1);
break;
}
mu[i*pri[j]]=-mu[i];
minp[i*pri[j]]=pri[j];
sigma[i*pri[j]]=sigma[i]*(pri[j]+1);
}
}
for(int i=1;i<=n;++i) sum[i]=(sum[i-1]+sum[i])%mod;
for(int i=1;i<=n;++i) sigma[i]=(sigma[i-1]+sigma[i])%mod;
}
ll Sqrt(ll n){
ll x=sqrt(n); if(x*x<=n) ++x; if(x*x>n) --x; return x;
}
struct Array{
ll a[100005],b[100005];
ll& operator [] (ll x){ return x<=100000?a[x]:b[n/x]; }
}_S,_G,_F;
ll S(ll n){
if(n<=m) return sigma[n];
if(_S[n]) return _S[n];
ll res=0;
for(ll l=1,r;l<=n;l=r+1){
r=n/(n/l);
res=(res+(l+r)%mod*((r-l+1)%mod)%mod*inv2%mod*(n/l%mod))%mod;
}
return _S[n]=res;
}
ll G(ll n){
if(n<=m) return sum[n];
if(_G[n]) return _G[n];
ll res=0;
for(ll i=1;i*i<=n;++i) res=(res+i*(Sqrt(n-i*i)<<1|1))%mod;
return _G[n]=res;
}
ll F(ll n){
if(_F[n]) return _F[n];
ll res=0;
for(ll l=1,r;l<=n;l=r+1){
r=n/(n/l);
res=(res+(S(r)-S(l-1))*G(n/l))%mod;
}
return _F[n]=res;
}
int main(){
scanf("%lld",&n); m=pow(n,0.66); sieve(m);
for(ll i=1;i*i<=n;++i) ans=(ans+i*mu[i]*F(n/i/i))%mod;
printf("%d
",(ans+mod)%mod);
return 0;
}