# [中等题] Project Euler 608 Divisor Sums

$\begin{array}{}\text{(29)}& D\left(m,n\right)& =& \sum _{d|m}\sum _{k=1}^{n}{\sigma }_{0}\left(kd\right)\text{(30)}& & =& \sum _{d|m}\sum _{k=1}^{n}\sum _{a|k}\sum _{b|d}\left[\left(a,b\right)=1\right]\text{(31)}& & =& \sum _{a=1}^{n}⌊\frac{n}{a}⌋\sum _{d|m}\sum _{b|d}\left[\left(a,b\right)=1\right]\text{(32)}& & =& \sum _{a=1}^{n}⌊\frac{n}{a}⌋\sum _{d|m}\sum _{b|d}\sum _{i|a,i|b}\mu \left(i\right)\text{(33)}& & =& \sum _{i|m}\mu \left(i\right)×\left({\sigma }_{0}\ast 1\right)\left(\frac{m}{i}\right)\sum _{id\le n}⌊\frac{n}{id}⌋\text{(34)}& & =& \sum _{i|m}\mu \left(i\right)×\left({\sigma }_{0}\ast 1\right)\left(\frac{m}{i}\right)\sum _{d=1}^{⌊\frac{n}{i}⌋}⌊\frac{⌊\frac{n}{i}⌋}{d}⌋\text{(35)}& & =& \sum _{i|m}\mu \left(i\right)×\left({\sigma }_{0}\ast 1\right)\left(\frac{m}{i}\right)×f\left(⌊\frac{n}{i}⌋\right)\end{array}$

$f$$f$ 这个东西全部预处理好，总计算量也就杜教筛复杂度 ${10}^{9}$$10^9$ 级别，还是能接受的。

#include<cstdio>
#include<cstdlib>
#include<map>
#include<cmath>
#include<algorithm>
using namespace std;
typedef long long ll;

const ll n=1e12,m=200;
const int P=1e9+7;

int p[205],num;
int vst[205];

int q[205];

inline void Pre(int n){
for (int i=2;i<=n;i++){
if (!vst[i]) p[++num]=i;
for (int j=1;j<=num && i*p[j]<=n;j++){
vst[i*p[j]]=1;
if (i%p[j]==0)
break;
}
}
}

inline ll C(ll x){
return (x+2)*(x+1)/2;
}

ll Map1[1000005],Map2[1000005];

inline ll Calc(ll x){
ll ret=0;
for (ll i=1,j,t;i<=x;i=j+1){
j=x/(t=x/i);
ret+=t*(j-i+1);
}
return (x<=1e6?Map1[x]:Map2[n/x])=ret%P;
}
inline ll Map(ll x){
return x<=1e6?Map1[x]:Map2[n/x];
}

ll ans=0;

inline void dfs(ll t,ll cur,ll mu,ll phi){
static int Cnt=0;
if (t==num+1){
if ((++Cnt)%1000000==0)
printf("%d\n",Cnt);
if (mu==1)
ans+=phi*Map(n/cur)%P;
else
ans+=P-phi*Map(n/cur)%P;
return;
}
if (cur*p[t]<=n)
dfs(t+1,cur*p[t],-mu,phi*C(q[t]-1)%P);
dfs(t+1,cur,mu,phi*C(q[t])%P);
}

inline ll Brute(){
ll mm=1;
for (int i=1;i<=m;i++) mm*=i;
int ret=0;
for (int i=1;i<=mm;i++)
if (mm%i==0)
for (int j=1;j<=n;j++)
for (int k=1;k<=i*j;k++)
if ((i*j)%k==0)
ret++;
printf("%d\n",ret);
return ret;
}
inline void GG(){
//ll ret=0;
ll Cnt=0;
for (ll i=1,j;i<=n;i=j+1){
j=n/(n/i);
Calc(n/i);
if ((++Cnt)%10000==0)
printf("%lld %lld\n",i,Cnt);
}
//printf("%lld\n",ret);
//return ret;
}

int main(){
Pre(m);
for (int i=1;i<=num;i++){
q[i]=0;
for (ll t=p[i];t<=m;t*=p[i])
q[i]+=m/t;
}
GG();
dfs(1,1,1,1);
printf("%lld\n",ans%P);
//Brute();
return 0;
}