# 4815: [Cqoi2017]小Q的表格

$ans = \sum_{d=1}^n f(d) \sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{n}{d}}ij[(i,j)=1]$

$ans = \sum_{d=1}^n f(d) S(\frac{n}{d}) \\ S(n) = \sum_{i=1}^n i \sum_{j=1}^n j [(i,j)=1] = \sum_{i=1}^n \varphi(i)*i^2$

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
using namespace std;
typedef long long ll;
const int N = 4e6+5, mo = 1e9+7, M = 3e3+5, inv2 = (mo+1)/2;
char c=getchar(); int x=0,f=1;
while(c<'0' || c>'9') {if(c=='-')f=-1; c=getchar();}
while(c>='0' && c<='9') {x=x*10+c-'0'; c=getchar();}
return x*f;
}

bool notp[N]; int p[N/10], phi[N]; ll s[N];
void sieve(int n) {
phi[1] = 1;
for(int i=2; i<=n; i++) {
if(!notp[i]) p[++p[0]] = i, phi[i] = i-1;
for(int j=1; j <= p[0] && i*p[j] <= n; j++) {
notp[i*p[j]] = 1;
if(i%p[j] == 0) {
phi[i*p[j]] = phi[i] * p[j];
break;
}
phi[i*p[j]] = phi[i] * (p[j] - 1);
}
}
for(int i=1; i<=n; i++) s[i] = (s[i-1] + (ll) i * i %mo * phi[i] %mo) %mo;
}

int gcd(int a, int b) {return !b ? a : gcd(b, a%b);}
ll Pow(ll a, int b) {
ll ans=1;
for(; b; b>>=1, a=a*a%mo)
if(b&1) ans=ans*a%mo;
return ans;
}

int Q, n, a, b, k; ll x;
namespace B {
int pos[N], block, m;
struct _blo{int l, r;} b[N];
void init() {
block = sqrt(n); m = (n-1)/block+1;
for(int i=1; i<=n; i++) pos[i] = (i-1)/block+1, a[i] = (ll)i*i%mo, f[i] = (a[i] + f[i-1]) %mo;
for(int i=1; i<=m; i++) b[i].l = (i-1)*block+1, b[i].r = i*block;
}
int que(int x) { return (f[x] + add[pos[x]]) %mo; }
void cha(int x, int v) {
int d = (v - a[x] + mo) %mo, r = b[pos[x]].r; a[x] = v;
if(d==0) return;
for(int i=x; i<=r; i++) f[i] = (f[i] + d) %mo;
}
} using B::cha; using B::que;

void solve(int n) {
int ans = 0, r, last = 0, now;
for(int i=1; i<=n; i=r+1, last=now) {
r = n/(n/i); now = que(r);
ans = (ans + (ll) (now - last + mo) * s[n/i] %mo) %mo;
}
printf("%d\n", (ans + mo) %mo);
}

int main() {
freopen("in", "r", stdin);