問題
解法
$$ f ( x , y ) : = \# \{ ( i , j ) : 1 \leq i \leq j \leq N , \; \gcd ( i , j ) = x , \; \gcd ( P _ i , P _ j ) = y \} $$
とおくと、$x\gt N$ または $y\gt N$ のとき $f(x,y)=0$ なので、答えは $(1)$ で表されます。
$$ \sum _ {x = 2} ^ N \sum _ { y = 2 } ^ N f ( x , y ) \tag{1} $$
$\gcd ( a , b ) = x$ などの等号条件よりも $x \mid \gcd ( a , b ) \iff x\mid a \land x\mid b$ のような倍数条件の方が扱いやすいです。そこで、次のように $g$ を定義します。
$$ g ( x , y ) : = \# \{ ( i , j ) : 1 \leq i \leq j \leq N , \; x \mid \gcd ( i , j ), \; y \mid \gcd ( P _ i , P _ j ) \} $$
このとき、$g(x,y)$ は次のように表されます。
$$ g ( x , y ) = \# \left \{ ( i , j ) \; :\; 1 \leq i \leq j \leq \left\lfloor\frac { N }{ x } \right\rfloor, y \mid P _ { i x }, y \mid P _ { j x} \right \} $$
$c(x,y)$ を $c ( x, y ) := \# \left \{ i \; :\; 1 \leq i \leq \lfloor N / x \rfloor,\; y \mid P _ { i x } \right \}$ で定義すれば、$g(x,y)$ は $c$ を用いて以下のように計算できます。
$$ g(x,y)= \dfrac{c(x,y) \cdot (c(x,y) + 1)}{2} $$
$g$ の計算方法は分かったので、次は $f(x,y)$ を $g$ を使って表します。$\displaystyle g( x , y ) = \sum _ { x \mid x ' } \sum _ { y \mid y ' } f ( x ' , y ' ) $ より、$\displaystyle f ( x , y ) = \sum _ { x \mid a ' } \sum _ { y \mid b ' } \mu ( a ' / x ) \cdot \mu (b ' / y ) \cdot g ( a' , b' )$ です。
これは以下のようにして確かめられますが、倍数系メビウス変換の式なので覚えてしまった方が速いと思います。
証明
$$
\begin{aligned}
f ( x , y )
&= \sum _ { x \mid x ' } f ( x' , y ) \cdot \delta _ { x , x ' } \\
&= \sum _ { x \mid x ' } f ( x' , y ) \sum _ { a \mid (x ' / x )} \mu ( a ) \\
&= \sum _ { x \mid x ' } \sum _ { ax \mid x ' } f ( x' , y ) \cdot \mu ( a ) \\
&= \sum _ { x \mid a ' \land a' \mid x ' } f ( x' , y ) \cdot \mu ( a ' / x ) \\
&= \sum _ { x \mid a ' } \mu ( a ' / x ) \sum _ { a' \mid x ' } f ( x' , y ) \\
&= \sum _ { x \mid a ' } \sum _ { y \mid b ' } \mu ( a ' / x ) \cdot \mu (b ' / y ) \sum _ { a' \mid x ' } \sum _ { b' \mid y ' } f ( x' , y' ) \\
&= \sum _ { x \mid a ' } \sum _ { y \mid b ' } \mu ( a ' / x ) \cdot \mu (b ' / y ) \cdot g ( a' , b' )
\end{aligned}
$$
(書いてから気づいたんですが、和を展開するのが素直ですね)
以上より、答えは次のように表されます。
$$ \begin{aligned} \sum _ { x = 2 } ^ N \sum _ { y = 2 } ^ N f ( x , y ) &= \sum _ { x = 2 } ^ N \sum _ { y = 2 } ^ N \sum _ { x \mid a ' } \sum _ { y \mid b ' } \mu ( a ' / x ) \cdot \mu (b ' / y ) \cdot g ( a' , b' ) \\ &= \sum _ { a' = 2 } ^ N \sum _ { b' = 2 } ^ N g ( a' , b' ) \sum _ { x \mid a ' , x \neq 1 } \mu ( a ' / x ) \sum _ { y \mid b ', y \neq 1 } \mu (b ' / y ) \\ &= \sum _ { a' = 2 } ^ N \sum _ { b' = 2 } ^ N g ( a' , b' ) \cdot ( \delta _ { a ', 1 } - \mu ( a ' ) ) \cdot ( \delta _ { b ', 1 } - \mu ( b ' ) ) \\ &= \sum _ { a' = 2 } ^ N \sum _ { b' = 2 } ^ N g ( a' , b' ) \cdot \mu ( a ' ) \cdot \mu ( b ' ) \\ &= \sum _ { a' = 2 } ^ N \mu ( a ' ) \sum _ { b' = 2 } ^ N g ( a' , b' ) \cdot \mu ( b ' ) \\ &= \sum _ { a' = 2 } ^ N \mu ( a ' ) \sum _ { b' = 2 } ^ N \dfrac{c(a', b')\cdot (c(a', b') + 1)}{2} \cdot \mu ( b ' ) \end{aligned} $$
まず、外側の和で $\mu (a') = 0$ となる部分は無視してよいです。そして、内側の和でも $\mu (b') = 0$ となる部分は無視してよいです。従って $c(x,y)$ は $\mu(x)\neq 0\land \mu(y)\neq 0$ となる部分しか必要なく、$y \mid P _ { i x } \land \mu(y) \neq 0$ を満たす $y$ の個数は $2 ^ { \omega(P _ { i x }) } - 1$ 個 ($\omega(n)$ は $n$ の異なる素因数の個数) です。$n \leq 2\cdot 10 ^ 5$ では $\omega(n)\leq 6$ なので、これは高々 $63$ 個です。
以下のようなコードを書いて最大ケースでどれくらいの計算回数になるかを見積もると $90809523 \lt 10 ^ 8$ 程度であり、常に $\omega(n)=6$ となる訳ではないことを加味すれば十分高速であることが確かめられます。
実験コード
#include <array> #include <iostream> #include <vector> template <uint32_t N> struct MobiusFunction { constexpr MobiusFunction() : mpf{}, dat{} { dat.fill(1); for (uint64_t p = 2; p <= N; ++p) { if (mpf[p]) continue; mpf[p] = p; dat[p] = -1; for (uint64_t q = p * 2; q <= N; q += p) { if (not mpf[q]) mpf[q] = p; dat[q] = q % (p * p) ? -dat[q] : 0; } } } constexpr int32_t operator()(const uint32_t n) const { return dat[n]; } private: std::array<uint32_t, N + 1> mpf; std::array<int32_t, N + 1> dat; }; constexpr uint32_t N = 200000; MobiusFunction<N> mobius; int main() { int64_t sum = 0; for (uint32_t a = 2; a <= N; ++a) { if (not mobius(a)) continue; sum += (N / a) * 63; } std::cout << sum << std::endl; return 0; }
実装
コンパイル時に素因数列挙やメビウス関数の計算をすることで 123 ms で通すことが出来ました (2021/12/15 時点で Fastest)。
感想
メビウス関数で和の形にするのは典型手法なんでしょうか
