Description：

求： $$ \sum_{i=1}^n\sum_{j=1}^nij\gcd(i,j) \mod p $$ $1\leqslant n\leqslant 10^{10},5\times10^8\leqslant p\leqslant1.1\times10^9$ 且 $p$ 为质数。

Solution：

用 $\sum_{d|i,d|j}\varphi(d)$ 来化 $\gcd(i,j)$ ，然后交换求和号就是： $$ \sum_{d=1}^n\varphi(d)d^2(\sum_1^{\lfloor\frac n d\rfloor}i)^2 $$ 前面 $\varphi(d)d^2$ 可以用杜教筛求，后面除法分块。

Code：

#include<algorithm>

#include<iostream>

#include<cstdlib>

#include<cstdio>

#include<cmath>

#include<map>

#include<cctype>

#include<cstring>

using namespace std;

typedef long long ll;

ll p,n;

#define MAX 3000010

bool isprime[MAX];

ll prime[MAX],tot = 0;

ll phi[MAX],sum[MAX];

ll inv2,inv6;

inline ll power(ll a,ll b)

{

	register ll res = 1;

	while(b > 0)

	{

		if(b & 1)res = res * a % p;

		a = a * a % p;

		b = b >> 1;

	}

	return res;

}

inline ll calc1(ll k){k %= p;return k * (k + 1) % p * inv2 % p;}

inline ll calc2(ll k){k %= p;return k * (k + 1) % p * (2 * k + 1) % p * inv6 % p;}

inline ll calc3(ll k){return calc1(k) * calc1(k) % p;}

map<ll,ll> sum_;

ll calc_sum(ll k)

{

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

	if(k == 1)return 1;

	if(k == 0)return 0;

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

	register ll res = calc3(k);

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

	{

		r = k / (k / l);

		res = (res - (calc2(r) - calc2(l - 1) + p) % p * calc_sum(k / l) % p + p) % p;

	}

	sum_[k] = res;

	return res;

} 

int main()

{

	scanf("%lld%lld",&p,&n);

	inv2 = power(2,p - 2);inv6 = power(6,p - 2);

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

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

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

	{

		if(isprime[i])phi[i] = i - 1,prime[++tot] = i;

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

		{

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

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

			else

			{

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

				break;

			}

		}

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

	}

	register ll ans = 0;

	for(register ll l = 1,r;l <= n;l = r + 1)

	{

		r = n / (n / l); 

		ans = (ans + (calc_sum(r) - calc_sum(l - 1) + p) % p * calc1(n / l) % p * calc1(n / l) % p) % p;

	}

	cout << ans << endl;

	return 0;

}