求

(∑i=1n∑j=1nijgcd(i,j))mod p (\sum_{i=1}^{n}\sum_{j=1}^{n}ij\gcd(i,j)) mod~p

n≤1010,p≤1.1×109n\le 10^{10}, p\le 1.1\times10^9

Luogu P3768

Solution

首先化下式子

∑i=1n∑j=1nijgcd(i,j)=∑i=1n∑j=1nij∑d∣gcd(i,j)φ(d)=∑d=1nφ(d)(∑i=1⌊nd⌋id)2=∑d=1nφ(d)d2S(⌊nd⌋)2 \begin{aligned} &\sum_{i=1}^{n}\sum_{j=1}^{n}ij\gcd(i,j)\\ =&\sum_{i=1}^{n}\sum_{j=1}^{n}ij\sum_{d|\gcd(i,j)}\varphi(d)\\ =&\sum_{d=1}^{n}\varphi(d)(\sum_{i=1}^{\lfloor \frac{n}{d} \rfloor }id)^2\\ =&\sum_{d=1}^{n}\varphi(d)d^2S(\lfloor \frac{n}{d} \rfloor )^2 \end{aligned}

其中S(n)=∑i=1niS(n)=\sum_{i=1}^ni

可以发现，后半部分S(⌊nd⌋)2S(\lfloor \frac{n}{d} \rfloor)^2显然是可以整除分块的

那么，我们就是要快速计算某一段φ(d)d2\varphi(d)d^2的和

即需要快速计算f(i)=φ(d)d2f(i)=\varphi(d)d^2的前缀和

考虑如何构造杜教筛

∑d∣nφ(d)d2(nd)2=n2∑d∣nφ(d)=n3⇒f∗Id2=Id3 \begin{aligned} \sum_{d|n}\varphi(d)d^2(\frac{n}{d})^2&=n^2\sum_{d|n}\varphi(d)=n^3\\ \Rightarrow f*\mathrm{Id}^2&=\mathrm{Id}^3 \end{aligned}

由于∑i=1ni3=(∑i=1ni)2\sum_{i=1}^ni^3=(\sum_{i=1}^ni)^2可以快速计算

所以杜教筛的式子长这样：

F(n)=(n(n+1)2)2−∑i=2ni2F(⌊ni⌋) F(n)=(\frac{n(n+1)}{2})^2-\sum_{i=2}^ni^2F(\lfloor\frac{n}{i}\rfloor)

Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
#include <bits/stdc++.h>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/hash_policy.hpp>

#define x first
#define y second
#define y1 Y1
#define y2 Y2
#define mp make_pair
#define pb push_back

using namespace std;

typedef long long LL;
typedef pair<int, int> pii;

template <typename T> inline int Chkmax (T &a, T b) { return a < b ? a = b, 1 : 0; }
template <typename T> inline int Chkmin (T &a, T b) { return a > b ? a = b, 1 : 0; }

inline void proc_status()
{
	ifstream t ("/proc/self/status");
	cerr << string (istreambuf_iterator <char> (t), istreambuf_iterator <char> ()) <<endl;
}

inline LL read ()
{
	LL sum = 0, fl = 1; char ch = getchar();
	for (; !isdigit(ch); ch = getchar()) if (ch == '-') fl = -1;
	for (; isdigit(ch); ch = getchar()) sum = (sum << 3) + (sum << 1) + ch - '0';
	return sum * fl;
}

const int Maxn = 1e7 + 100;

LL N, Mod, inv2, inv6;
int Phi[Maxn], Not_Prime[Maxn], Prime[Maxn], prime_cnt;
__gnu_pbds :: gp_hash_table<int, LL> F;
LL FF[Maxn];

inline void Add (LL &a, LL b) { a += b; if (a >= Mod) a -= Mod; }

inline LL Pow (LL a, int b)
{
	LL ans = 1;
	while (b) { if (b & 1) ans = 1ll * ans * a % Mod; a = 1ll * a * a % Mod; b >>= 1; }
	return ans;
}

inline LL Calc_S (LL x)
{
	x %= Mod;
	LL sum = 1ll * x * (x + 1) % Mod * inv2 % Mod;
	return 1ll * sum * sum % Mod;
}

inline LL Calc_G (LL x)
{
	x %= Mod;
	return x * (x + 1) % Mod * (2 * x + 1) % Mod * inv6 % Mod;
}

inline LL Calc_F (LL n)
{
	if (n <= Maxn - 5) return FF[n];
	if (F[n]) return F[n];
	LL ans = Calc_S(n), sum1, sum2, sum;
	for (LL l = 2, r; l <= n; l = r + 1)
	{
		r = n / (n / l);
		sum1 = Calc_G(l - 1), sum2 = Calc_G(r);
		sum = Calc_F(n / l);
/**/	Add (ans, Mod - sum * (sum2 - sum1 + Mod) % Mod);
	}
	return F[n] = ans;
}

inline void Solve ()
{
	LL sum1, sum2, sum, ans = 0;
	for (LL l = 1, r; l <= N; l = r + 1)
	{
		r = N / (N / l);
		sum1 = Calc_F(l - 1), sum2 = Calc_F(r);
		sum = Calc_S(N / l);
		Add (ans, sum * (sum2 - sum1 + Mod) % Mod);
	}
	cout<<ans<<endl;
}

inline void Input ()
{
	Mod = read(), N = read();
}

inline void Init ()
{
	inv2 = Pow(2, Mod - 2), inv6 = Pow(6, Mod - 2);

	Phi[1] = 1;
	for (int i = 2; i <= min(N, 1ll * Maxn - 5); ++i)
	{
		if (!Not_Prime[i])
		{
			Prime[++prime_cnt] = i;
			Phi[i] = i - 1;
		}
		for (int j = 1; j <= prime_cnt && Prime[j] * i <= min(N, 1ll * Maxn - 5); ++j)
		{
			Not_Prime[i * Prime[j]] = 1;
/**/		if (!(i % Prime[j])) { Phi[i * Prime[j]] = Phi[i] * Prime[j]; break; }
			else Phi[i * Prime[j]] = Phi[i] * (Prime[j] - 1);
		}
	}
	
/**/for (LL i = 1; i <= min(N, 1ll * Maxn - 5); ++i) FF[i] = Phi[i] * i % Mod * i % Mod;
/**/for (LL i = 1; i <= min(N, 1ll * Maxn - 5); ++i) Add (FF[i], FF[i - 1]);
}

int main()
{
#ifndef ONLINE_JUDGE
	freopen("A.in", "r", stdin);
	freopen("A.out", "w", stdout);
#endif
	Input();
	Init();
	Solve();
	return 0;
}

Debug

  • 73L：把减号写成加号
  • 111L：把Prime[j]Prime[j]写成jj
  • 116、117L：求这个函数前缀和的时候要分两步算，一开始直接一步算的。。。