Description：

求： $$ \sum_{i=1}^n\sum_{j=1}^m\varphi(ij) $$ $1\leqslant n\leqslant 10^5,1\leqslant m\leqslant 10^9$

Solution：

看数据范围应该是枚举 $i$ 然后用一个方法计算 $m$ ，一个新的套路： $$ \begin{align} S(n,m)=&\sum_{i=1}^m\varphi(ni)\\ \end{align} $$ 由于 $\varphi$ 的性质，所以可以把 $n$ 中所有质因子提出来 $p_i^{q_i-1}$ ，即设 $x=\prod p_i^{q_i-1}w=\prod p_i$ ， $\varphi(n)=x\varphi(w)$ ，所以有： $$ \begin{align} S(n,m)=&x\sum_{i=1}^m\varphi(wi)\\ \end{align} $$ 这个时候由于 $w$ 是无平方因子数，因此一定满足若 $a\times b=w$ ， $\varphi(a)\times\varphi(b)=\varphi(w)$ 。因此有： $$ \begin{align} S(n,m)&=x\sum_{i=1}^m\varphi(\frac w {\gcd(w,i)}\times\gcd(w,i)\times i)\\ \because \gcd(&\frac w {\gcd(w,i)},i)=1,\gcd(\frac w {\gcd(w,i)},\gcd(w,i))=1\\ \therefore\gcd(&\frac w {\gcd(w,i)},\gcd(w,i)\times i)=1\\ S(n,m)&=x\sum_{i=1}^m\varphi(\frac w {\gcd(w,i)})\times\varphi(\gcd(w,i)\times i)\\ &=x\sum_{i=1}^m\varphi(\frac w {\gcd(w,i)})\times\varphi(i)\times \gcd(w,i)\\ &=x\sum_{i=1}^m\varphi(\frac w {\gcd(w,i)})\times\varphi(i)\times \sum_{d|w,d|i}\varphi(\frac{\gcd(w,i)}d)\\ &=x\sum_{i=1}^m\varphi(i)\times \sum_{d|w,d|i}\varphi(\frac w d)\\ &=x\sum_{d|w}\varphi(\frac w d)\sum_{i=1}^{\lfloor \frac m d\rfloor}\varphi(id)\\ &=x\sum_{d|w}\varphi(\frac w d)S(d,\lfloor\frac m d\rfloor) \end{align} $$ 然后就可以递归的去做了， $S(n,1)=\varphi(n),S(1,n)$ 可以用杜教筛，题解里用状压枚举约数很妙啊。

分析一下复杂度： $S(n,m)$ 中第一维一定是 $n$ 的约数，上限是 $O(\sqrt n)$ ，实际上远远不到，第二维由于 $\lfloor\frac{\lfloor\frac n a\rfloor}{b}\rfloor=\lfloor\frac n {ab}\rfloor$ ，所以只有 $\sqrt m$ 个，实际也跑不满，最后杜教筛的也是所有的 $\lfloor\frac m d\rfloor$ ，所以总复杂度 $O(n^{\frac 2 3})$ ，由于最外层枚举 $+S(n,m)$ 记忆化第一维可以看成 $1\sim n$ ，所以时间复杂度 $O(n\sqrt m+m^{\frac 2 3})$ ，而且跑不满。

Code：

#include<algorithm>

#include<iostream>

#include<cstdlib>

#include<cstdio>

#include<cmath>

#include<map>

#include<cctype>

#include<cstring>

using namespace std;

int n,m;

#define MAX 1000010

bool isprime[MAX];

int prime[MAX],tot = 0,mindiv[MAX];

typedef long long ll;

ll phi[MAX],sum[MAX];

map<int,ll> phi_;

#define MOD 1000000007

ll calc(int k)

{

	if(k == 0)return 0;

	if(k == 1)return 1;

	if(k < MAX)return sum[k];

	if(phi_.find(k) != phi_.end())return phi_[k];

	ll res = (1ll * k * (k + 1) / 2) % MOD;

	for(int l = 2,r;l <= k;l = r + 1)

	{

		r = k / (k / l);

		res = (res - (r - l + 1) * calc(k / l) % MOD + MOD) % MOD;

	}

	phi_[k] = res;

	return res;

}

#define MAXN 100010

map<int,ll> s[MAXN];

ll S(int n,int m)

{

	if(m == 0)return 0;

	if(m == 1)return phi[n];

	if(n == 1)return calc(m);

	if(s[n].find(m) != s[n].end())return s[n][m];

	ll x = 1,w = 1;

	int fac[20];

	fac[0] = 0;

	ll v = n;

	while(v != 1)

	{

		ll t = mindiv[v];

		w *= t;v /= t;

		fac[++fac[0]] = t;

		while(v % t == 0)

		{

			x *= t;v /= t;

		}

	}

	ll sum = 0;

	for(int i = 0;i < (1 << fac[0]);++i)

	{

		int d = 1;

		for(int t = 1;t <= fac[0];++t)if(i & (1 << (t - 1)))d *= fac[t];

		sum = (sum + phi[w / d] * S(d,m / d) % MOD) % MOD;

	}

	s[n][m] = sum * x % MOD;

	return s[n][m];

}

int main()

{

	for(int i = 2;i < MAX;++i)isprime[i] = true;

	phi[1] = sum[1] = 1;

	for(int i = 2;i < MAX;++i)

	{

		if(isprime[i])

		{

			phi[i] = i - 1;

			prime[++tot] = i;

			mindiv[i] = i;

		}

		for(int j = 1;j <= tot && i * prime[j] < MAX;++j)

		{

			isprime[i * prime[j]] = false;

			mindiv[i * prime[j]] = prime[j];

			if(i % prime[j] == 0)

			{

				phi[i * prime[j]] = phi[i] * prime[j];

				break;

			}

			else phi[i * prime[j]] = phi[i] * (prime[j] - 1);

		}

		sum[i] = (sum[i - 1] + phi[i]) % MOD;

	}

	scanf("%d%d",&n,&m);

	ll ans = 0;

	for(int i = 1;i <= n;++i)

	{

		ans = (ans + S(i,m)) % MOD;

	}

	cout << ans << endl;

	return 0;

}