给出 $n$ 和 $p$ ，求 $(\sum\limits_{i=1}^n\sum\limits_{j=1}^nij\gcd(i,j))\mod p$ 。
$n\le 10^{10}$ 。
题解
欧拉函数(欧拉反演)+杜教筛
推式子：
$$\begin{align}&\sum\limits_{i=1}^n\sum\limits_{j=1}^nij\gcd(i,j)\\=&\sum\limits_{i=1}^n\sum\limits_{j=1}^nij\sum\limits_{d|\gcd(i,j)}\varphi(d)\\=&\sum\limits_{i=1}^n\sum\limits_{j=1}^nij\sum\limits_{d|i,d|j}\varphi(d)\\=&\sum\limits_{d=1}^n\varphi(d)\sum\limits_{i=1}^{\lfloor\frac nd\rfloor}id\sum\limits_{j=1}^{\lfloor\frac nd\rfloor}jd\\=&\sum\limits_{d=1}^nd^2\varphi(d)(\sum\limits_{i=1}^{\lfloor\frac nd\rfloor}i)^2\\=&\sum\limits_{d=1}^nd^2\varphi(d)(\frac{\lfloor\frac nd\rfloor(\lfloor\frac nd\rfloor+1)}2)^2\end{align}$$
对 $\lfloor\frac nd\rfloor$ 分块处理，只需要求出 $f(n)=n^2\varphi(n)$ 的前缀和即可。
显然这是一个积性函数，然而 $n$ 有 $10^{10}$ 之大，不能线性筛。
考虑杜教筛。设 $g(n)=n^2$ ，那么有：
$$\begin{align}&(f·g)(n)\\=&\sum\limits_{d|n}f(d)g(\frac nd)\\=&\sum\limits_{d|n}d^2\varphi(d)·(\frac nd)^2\\=&n^2\sum\limits_{d|n}\varphi(d)\\=&n^3\end{align}$$
所以：
$$\begin{align}&\sum\limits_{i=1}^ni^3\\=&\sum\limits_{i=1}^n(f·g)(i)\\=&\sum\limits_{i=1}^n\sum\limits_{d|i}f(d)·g(\frac id)\\=&\sum\limits_{i=1}^n\sum\limits_{d|i}f(\frac id)g(d)\\=&\sum\limits_{d=1}^ng(d)\sum\limits_{i=1}^{\lfloor\frac nd\rfloor}f(i)\\=&\sum\limits_{d=1}^nd^2S(\lfloor\frac nd\rfloor)\end{align}$$
故有：
$$S(n)=\sum\limits_{i=1}^ni^3-\sum\limits_{i=2}^ni^2S(\lfloor\frac nd\rfloor)$$
线性筛预处理出 $n^{\frac 23}$ 以内的 $S(i)$ ，对超过 $n^{\frac 23}$ 的部分进行杜教筛即可。
可能需要用到的公式：
$$\sum\limits_{i=1}^ni^2=\frac{n(n+1)(2n+1)}6\\\sum\limits_{i=1}^ni^3=\frac{n^2(n+1)^2}4$$
时间复杂度 $O(n^{\frac 23})$ ，这里偷懒使用map，复杂度多一个 $\log$ 。
#include <map>
#include <cstdio>
#define N 10000010
using namespace std;
typedef long long ll;
const int m = 10000000;
map<ll , ll> mp;
int phi[N] , prime[N] , tot , np[N];
ll sum[N] , p , inv4 , inv6;
void init()
{
int i , j;
phi[1] = 1;
for(i = 2 ; i <= m ; i ++ )
{
if(!np[i]) phi[i] = i - 1 , prime[++tot] = i;
for(j = 1 ; j <= tot && i * prime[j] <= m ; j ++ )
{
np[i * prime[j]] = 1;
if(i % prime[j] == 0)
{
phi[i * prime[j]] = phi[i] * prime[j];
break;
}
else phi[i * prime[j]] = phi[i] * phi[prime[j]];
}
}
for(i = 1 ; i <= m ; i ++ ) sum[i] = (sum[i - 1] + 1ll * phi[i] * i % p * i) % p;
inv4 = (p + 1) / 2;
if(p % 3 == 1) inv6 = (2 * p + 1) / 3;
else inv6 = (p + 1) / 3;
inv6 = inv6 * inv4 % p , inv4 = inv4 * inv4 % p;
}
ll calc2(ll n) {n %= p; return n * (n + 1) % p * (2 * n + 1) % p * inv6 % p;}
ll calc3(ll n) {n %= p; return n * n % p * (n + 1) % p * (n + 1) % p * inv4 % p;}
ll solve(ll n)
{
if(n <= m) return sum[n];
if(mp.find(n) != mp.end()) return mp[n];
ll i , last , ans = calc3(n);
for(i = 2 ; i <= n ; i = last + 1) last = n / (n / i) , ans = (ans - (calc2(last) - calc2(i - 1) + p) * solve(n / i) % p + p) % p;
return mp[n] = ans;
}
int main()
{
ll n , i , last , ans = 0;
scanf("%lld%lld" , &p , &n);
init();
for(i = 1 ; i <= n ; i = last + 1)
last = n / (n / i) , ans = (ans + (solve(last) - solve(i - 1) + p) * calc3(n / i)) % p;
printf("%lld\n" , ans);
return 0;
}