给出 $n$ 和 $p$ ，求 $(sumlimits_{i=1}^nsumlimits_{j=1}^nijgcd(i,j))mod p$ 。
$nle 10^{10}$ 。
题解
欧拉函数(欧拉反演)+杜教筛
推式子：
$$egin{align}&sumlimits_{i=1}^nsumlimits_{j=1}^nijgcd(i,j)\=&sumlimits_{i=1}^nsumlimits_{j=1}^nijsumlimits_{d|gcd(i,j)}varphi(d)\=&sumlimits_{i=1}^nsumlimits_{j=1}^nijsumlimits_{d|i,d|j}varphi(d)\=&sumlimits_{d=1}^nvarphi(d)sumlimits_{i=1}^{lfloorfrac nd floor}idsumlimits_{j=1}^{lfloorfrac nd floor}jd\=&sumlimits_{d=1}^nd^2varphi(d)(sumlimits_{i=1}^{lfloorfrac nd floor}i)^2\=&sumlimits_{d=1}^nd^2varphi(d)(frac{lfloorfrac nd floor(lfloorfrac nd floor+1)}2)^2end{align}$$
对 $lfloorfrac nd floor$ 分块处理，只需要求出 $f(n)=n^2varphi(n)$ 的前缀和即可。
显然这是一个积性函数，然而 $n$ 有 $10^{10}$ 之大，不能线性筛。
考虑杜教筛。设 $g(n)=n^2$ ，那么有：
$$egin{align}&(f·g)(n)\=&sumlimits_{d|n}f(d)g(frac nd)\=&sumlimits_{d|n}d^2varphi(d)·(frac nd)^2\=&n^2sumlimits_{d|n}varphi(d)\=&n^3end{align}$$
所以：
$$egin{align}&sumlimits_{i=1}^ni^3\=&sumlimits_{i=1}^n(f·g)(i)\=&sumlimits_{i=1}^nsumlimits_{d|i}f(d)·g(frac id)\=&sumlimits_{i=1}^nsumlimits_{d|i}f(frac id)g(d)\=&sumlimits_{d=1}^ng(d)sumlimits_{i=1}^{lfloorfrac nd floor}f(i)\=&sumlimits_{d=1}^nd^2S(lfloorfrac nd floor)end{align}$$
故有：
$$S(n)=sumlimits_{i=1}^ni^3-sumlimits_{i=2}^ni^2S(lfloorfrac nd floor)$$
线性筛预处理出 $n^{frac 23}$ 以内的 $S(i)$ ，对超过 $n^{frac 23}$ 的部分进行杜教筛即可。
可能需要用到的公式：
$$sumlimits_{i=1}^ni^2=frac{n(n+1)(2n+1)}6\sumlimits_{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
" , ans);
return 0;
}