Description
输入一个整数n和一个整数p，你需要求出$(sum_{i=1}^nsum_{j=1}^n ijgcd(i,j))~mod~p$，其中gcd(a,b)表示a与b的最大公约数。
Input
一行两个整数p、n。
Output
一行一个整数$(sum_{i=1}^nsum_{j=1}^n ijgcd(i,j))~mod~p$。
Sample Input
998244353 2000
Sample Output
883968974
HINT
对于20%的数据，$n leq 1000$。
对于30%的数据，$n leq 5000$。
对于60%的数据，$n leq 10^6$，时限1s。
对于另外20%的数据，$n leq 10^9$，时限3s。
对于最后20%的数据，$n leq 10^{10}$，时限6s。
对于100%的数据，$5 imes 10^8 leq p leq 1.1 imes 10^9$且p为质数。
题解
题目要求 $$sum_{i=1}^nsum_{j=1}^n ijcdot gcd(i,j)$$
提出 $gcd$ egin{aligned}&sum_{d=1}^ndsum_{i=1}^{n}sum_{j=1}^nij[gcd(i,j)=d]\=&sum_{d=1}^nd^3sum_{i=1}^{leftlfloorfrac{n}{d} ight floor}sum_{j=1}^{leftlfloorfrac{n}{d} ight floor}ij[gcd(i,j)=1]\=&sum_{d=1}^nd^3sum_{i=1}^{leftlfloorfrac{n}{d} ight floor}sum_{j=1}^{leftlfloorfrac{n}{d} ight floor}ijsum_{kmid gcd(i,j)}mu(k)\=&sum_{d=1}^nd^3sum_{k=1}^{leftlfloorfrac{n}{d} ight floor}mu(k)k^2sum_{i=1}^{leftlfloorfrac{n}{kd} ight floor}sum_{j=1}^{leftlfloorfrac{n}{kd} ight floor}ijend{aligned}
令 $F(x)=sumlimits_{i=1}^xi=frac{x(x+1)}{2}$ ， $T=kd$ egin{aligned}Rightarrow&sum_{d=1}^nd^3sum_{k=1}^{leftlfloorfrac{n}{d} ight floor}mu(k)k^2F^2left(frac{n}{kd} ight)\=&sum_{T=1}^nF^2left(frac{n}{T} ight)sum_{dmid T}d^3left(frac{T}{d} ight)^2muleft(frac{T}{d} ight)\=&sum_{T=1}^nF^2left(frac{n}{T} ight)T^2sum_{dmid T}dcdotmuleft(frac{T}{d} ight)end{aligned}
前面的很好处理，但那个狄利克雷卷积是什么鬼...我们把它拎出来调教一下: $$sum_{dmid T}dcdotmuleft(frac{T}{d} ight)$$
这个形式似乎不好看，我们让它女装变个样子： $$sum_{dmid T}frac{T}{d}cdotmu(d)=Tsum_{dmid T}frac{mu(d)}{d}$$
这玩意不就是 $varphi(T)$ 么。带入原柿 $$sum_{T=1}^nF^2left(frac{n}{T} ight)T^2varphi(T)$$
现在就好搞♂了，美滋滋。记 $f(T)=T^2varphi(T)$ 显然他是个积性函数，可以杜教筛了。
考虑求 $S(n)=sumlimits_{i=1}^nf(i)$
上述式子 $$g(1)S(n)=sum_{i=1}^n(g*f)(i)-sum_{i=2}^ng(i)Sleft(leftlfloorfrac{n}{i} ight floor ight)$$
考虑到 $sumlimits_{dmid n}varphi(d)=n$ ，又由于 $(g*f)(n)=sumlimits_{dmid n}varphi(d)d^2cdot gleft(frac{n}{d} ight)$ 。我们考虑让 $g(n)=id^2(n)$ ，那么 $(id*f)(n)=sumlimits_{dmid n}n^2cdotvarphi(d)=n^3$ 。由于 $sumlimits_{i=1}^ni^3=frac{n^2(n+1)^2}{4}=left(frac{n(n+1)}{2} ight)^2=F^2(n)$ 。显然这个卷积的前缀为 $sumlimits_{i=1}^n(g*f)(i)=F^2(n)$ 。
故对于 $f$ $$S(n)=F^2(n)-sum_{i=2}^ni^2cdot Sleft(leftlfloorfrac{n}{i} ight floor ight)$$
由公式 $sumlimits_{i=1}^ni^2=frac{n(n+1)(2n+1)}{6}$ 现在整个柿子都好算了。
1 //It is made by Awson on 2018.1.25 2 #include <set> 3 #include <map> 4 #include <cmath> 5 #include <ctime> 6 #include <queue> 7 #include <stack> 8 #include <cstdio> 9 #include <string> 10 #include <vector> 11 #include <cstdlib> 12 #include <cstring> 13 #include <iostream> 14 #include <algorithm> 15 #define LL long long 16 #define Abs(a) ((a) < 0 ? (-(a)) : (a)) 17 #define Max(a, b) ((a) > (b) ? (a) : (b)) 18 #define Min(a, b) ((a) < (b) ? (a) : (b)) 19 #define Swap(a, b) ((a) ^= (b), (b) ^= (a), (a) ^= (b)) 20 #define writeln(x) (write(x), putchar(' ')) 21 #define lowbit(x) ((x)&(-(x))) 22 using namespace std; 23 const int N = 2333333; 24 void read(LL &x) { 25 char ch; bool flag = 0; 26 for (ch = getchar(); !isdigit(ch) && ((flag |= (ch == '-')) || 1); ch = getchar()); 27 for (x = 0; isdigit(ch); x = (x<<1)+(x<<3)+ch-48, ch = getchar()); 28 x *= 1-2*flag; 29 } 30 void write(LL x) { 31 if (x > 9) write(x/10); 32 putchar(x%10+48); 33 } 34 35 LL n, p, phi[N+5], inv2, inv6; 36 LL prime[N+5], isprime[N+5], tot; 37 map<LL,LL>mp; 38 39 LL quick_pow(LL a, LL b) { 40 LL ans = 1; a %= p; 41 while (b) { 42 if (b&1) ans = ans*a%p; 43 a = a*a%p; b >>= 1; 44 } 45 return ans; 46 } 47 void get_pre() { 48 memset(isprime, 1, sizeof(isprime)); isprime[1] = 0, phi[1] = 1; 49 for (int i = 2; i <= N; i++) { 50 if (isprime[i]) prime[++tot] = i, phi[i] = 1ll*(i-1)*i%p*i%p; 51 for (int j = 1; j <= tot && i*prime[j] <= N; j++) { 52 isprime[i*prime[j]] = 0; 53 if (i%prime[j]) phi[i*prime[j]] = phi[i]*(prime[j]-1)%p*prime[j]%p*prime[j]%p; 54 else {phi[i*prime[j]] = phi[i]*prime[j]%p*prime[j]%p*prime[j]%p; break; } 55 } 56 (phi[i] += phi[i-1]) %= p; 57 } 58 } 59 LL sum(LL n) {n %= p; return n*(n+1)%p*inv2%p; } 60 LL sum2(LL n) {n %= p; return n*(n+1)%p*(n*2+1)%p*inv6%p; } 61 LL get_phi(LL n) { 62 if (n <= N) return phi[n]; 63 if (mp.count(n)) return mp[n]; 64 LL ans = sum(n)*sum(n)%p; 65 for (LL i = 2, last; i <= n; i = last+1) { 66 last = n/(n/i); (ans -= get_phi(n/i)*((sum2(last)-sum2(i-1))%p)%p) %= p; 67 } 68 return mp[n] = ans; 69 } 70 void work() { 71 read(p), read(n); get_pre(); inv2 = quick_pow(2, p-2), inv6 = quick_pow(6, p-2); LL ans = 0; 72 for (LL i = 1, last; i <= n; i = last+1) { 73 last = n/(n/i); LL s = sum(n/i); 74 (ans += s*s%p*((get_phi(last)-get_phi(i-1))%p)%p) %= p; 75 } 76 writeln((ans+p)%p); 77 } 78 int main() { 79 work(); 80 return 0; 81 }