問題
整数 $N$ が与えられる。座標平面上の格子点 $(x,y)$（$x,y$ は整数、原点 $(0,0)$ は除く）のうち、$x^2+y^2\le N$ かつ $\gcd(|x|,|y|)=1$ を満たすものの個数を求めよ。ただし $\gcd(0,y)=\lvert y\rvert$ とする（よって $(\pm1,0),(0,\pm1)$ は条件を満たし、$(\pm2,0)$ は満たさない）。
入力形式
N制約
| # | 条件 |
|---|---|
| 1 | $1\le N\le10^{12}$ |
| 2 | 入力は整数 |
入出力例
入力例1
1出力例1
4（$(\pm1,0),(0,\pm1)$ の 4 点）
入力例2
10出力例2
24入力例3
100出力例3
192背景と出典
円内の格子点を数える問題は Gauss の円問題として知られ、半径 $r$ の円内の格子点数は $\pi r^2$ に近い。「互いに素な $(x,y)$」は原点から遮られずに見える格子点にあたり、密度が $6/\pi^2$ になる。この $6/\pi^2$ は $1/\zeta(2)=\sum_{d\ge1}\mu(d)/d^2$ から来る。私の記憶では、この個数の漸近が $\frac{6}{\pi^2}\cdot\pi N=\frac{6N}{\pi}$ となる（下の実測で $N=10^{12}$ のとき比が $0.9999999984$）。本セッションでは文献を開いておらず、正しさの根拠は総当たりとの照合と、この漸近の比。
概念図
ヒント（クリックで段階的に開示）
ヒント1（方向性）
$\gcd(x,y)=1$ を直接数えるのは難しい。「$d$ が $x$ と $y$ の両方を割る点」なら数えやすい: $(x,y)=d(x',y')$ と書けて、$x'^2+y'^2\le N/d^2$ の点の数に等しい。
ヒント2（アプローチ）
$\sum_{d\mid g}\mu(d)=[g=1]$（$g=\gcd$）を使う。点の集合を「$\gcd$ がちょうど 1」と「$d\mid\gcd$ の $d$ を $\mu(d)$ 倍して足す」に分けると、答えは $\sum_{d\ge1}\mu(d)\cdot(\text{原点を除く、}x^2+y^2\le N/d^2\text{ の点の数})$。$d\le\sqrt N$ だけ見ればよい。
ヒント3（誘導）
円内の格子点数 $C(m)=\#\{(x,y): x^2+y^2\le m\}$ は、$x$ を固定して $y$ の個数 $2\lfloor\sqrt{m-x^2}\rfloor+1$ を足す $O(\sqrt m)$。各 $d$ で $m=\lfloor N/d^2\rfloor$ に対して $C(m)-1$（原点を除く）を $\mu(d)$ 倍して足す。Python では $\lfloor\sqrt{\cdot}\rfloor$ に math.isqrt を使う。
模範解答 (Python)
模範解答を表示する（まず自力で解いてから）
import sys
from math import isqrt
def main():
n = int(sys.stdin.readline())
r = isqrt(n) # d^2 <= n となる d の上限
# メビウス関数 mu(1..r) を線形篩で求める
mu = [1] * (r + 1)
is_comp = bytearray(r + 1)
primes = []
for i in range(2, r + 1):
if not is_comp[i]:
primes.append(i)
mu[i] = -1
for p in primes:
if i * p > r:
break
is_comp[i * p] = 1
if i % p == 0:
mu[i * p] = 0
break
mu[i * p] = -mu[i]
def disk(m):
"""x^2 + y^2 <= m を満たす整数点 (x, y) の個数（原点を含む）"""
s = isqrt(m)
tot = 0
for x in range(1, s + 1):
tot += isqrt(m - x * x)
return 1 + 4 * s + 4 * tot # x=0 の列 (2s+1 点) と ±、軸上の点 の分をまとめる
ans = 0
for d in range(1, r + 1):
if mu[d]:
ans += mu[d] * (disk(n // (d * d)) - 1) # 原点を除く
print(ans)
main()計算量: $\sum_{d\le\sqrt N}\sqrt{N/d^2}=\sqrt N\sum 1/d=O(\sqrt N\log N)$。実測（Python 3.14.7、M1 Mac）で $N=10^{10}$ が 0.15 s、$N=10^{12}$ が 1.42 s。
Step-by-Step 解説
$x$ を $-s..s$（$s=\lfloor\sqrt m\rfloor$）で動かし、各 $x$ で $|y|\le\lfloor\sqrt{m-x^2}\rfloor$ の $y$ が $2\lfloor\sqrt{m-x^2}\rfloor+1$ 個。$x$ と $-x$ は同じなので $x\ge1$ だけ数えて 2 倍し、$x=0$ の列を別に足すと $C(m)=(2s+1)+2s+4\sum_{x=1}^{s}\lfloor\sqrt{m-x^2}\rfloor=1+4s+4\sum_{x=1}^{s}\lfloor\sqrt{m-x^2}\rfloor$。
原点以外の点では $g=\gcd(|x|,|y|)\ge1$。$\sum_{d\mid g}\mu(d)=[g=1]$ より、 $$V(N)=\sum_{(x,y)\ne0,\ x^2+y^2\le N}\sum_{d\mid\gcd}\mu(d)=\sum_{d\ge1}\mu(d)\cdot\#\{(x,y)\ne0:\ d\mid x,\ d\mid y,\ x^2+y^2\le N\}.$$ $d\mid x,\ d\mid y$ の点は $d(x',y')$ と書け、条件は $x'^2+y'^2\le N/d^2$、すなわち（整数なので）$\lfloor N/d^2\rfloor$。点の数は $C(\lfloor N/d^2\rfloor)-1$。
$d>\sqrt N$ では $N/d^2<1$ で、原点以外の点が無く寄与が 0。$d\le\sqrt N$ だけ見れば足りる。各 $d$ の $C$ の計算は $O(\sqrt{N}/d)$。$\mu(d)=0$ の $d$（約 4 割）は飛ばせる。
$\lfloor\sqrt m\rfloor$ を int(math.sqrt(m)) で書いても、$N\le10^{12}$ では isqrt と一致した（$N=10^9,\ 10^9+7,\ 3456789012,\ 10^{10}$ で全体の出力が一致、実測）。ただし $m=(10^8+1)^2-1$（約 $10^{16}$）では int(math.sqrt(m))=100000001、isqrt(m)=100000000 と 食い違う（実測）。この制約の外（$m\gtrsim2^{52}\approx4.5\times10^{15}$ のあたりから。倍精度の仮数 52 ビットからの導出で、境目そのものは未測定）では float の平方根は使えない。**本問は制約内なので float でも通るが、isqrt を使えば制約を広げても壊れない**。
$C(m)$ は原点を含む。$\mu(d)\cdot C$ のまま足すと $\sum_d\mu(d)=M(r)$ の分だけずれる（Mertens 関数）。126 件の検算中 118 件で誤る（実測）。
迷った他案と採らなかった理由
| 案 | 計算量 | 採否 | 理由 |
|---|---|---|---|
| 全点を調べて $\gcd$ を取る | $O(N)$ | 検算用のみ | $N\le6000$ まで。$N=10^{12}$ には不可 |
| 円内の格子点数 $C(N)$ だけで近似（$6/\pi^2$ を掛ける） | $O(\sqrt N)$ | 不採用 | 整数の答えが出ない。漸近の比の確認には使った |
| $\gcd$ ごとに点を数える（$d$ ごとに互いに素な点の数を再帰） | $O(\sqrt N\log N)$ | 不採用 | 同じ形の再帰になり、Möbius 反転に帰着するだけ |
| Möbius 反転＋$C(\lfloor N/d^2\rfloor)$（本解） | $O(\sqrt N\log N)$ | 採用 | 総当たりと 126 件全件一致、$N=10^{12}$ で 1.42 s |
| $C(m)$ を $O(m^{1/3})$ 程度で求める高速な方法 | 未確認 | 扱わない | 存在すると記憶している（円問題の幾何的な高速計算）が、実装・検証していない。$N\le10^{18}$ に広げるなら必要 |
前提が変わったら:
- **$N\le10^{18}$** になると、$\sqrt N=10^9$ で本解は届かない。$C(m)$ を $m^{1/3}$ 程度で求める方法（凸な曲線上の格子点の数え上げ）が要る
- 正方形・楕円など別の領域でも、「$d$ の倍数の格子点は領域を $1/d$ に縮めた格子点」という構造が成り立てば同じ Möbius 反転が使える。円でなくなると $C(m)$ の計算だけ替える
- **$k$ 次元（互いに素な $k$ 組）**なら $\mu(d)$ の重みは同じで、$C$ が $k$ 次元の球の格子点数になる
よくあるミス
| ミス | 原因 | 正しい書き方 |
|---|---|---|
| 原点を引かない | $C(m)$ に原点が含まれていることを忘れる | disk(...) - 1（118/126 件で誤答） |
| $d$ の範囲を $N$ まで回す | $N/d^2<1$ なら寄与が 0 になることを使わない | d を isqrt(n) まで |
int(math.sqrt(m)) を大きい $m$ で使う | 浮動小数点の丸め（$m\gtrsim4.5\times10^{15}$ で破れる） | math.isqrt(m)（本問の制約内では float でも一致したが、広げると破れる） |
| $\gcd(0,y)$ を 0 や 1 と扱う | 軸上の点の条件を間違える | $\gcd(0,y)=\lvert y\rvert$ より $(0,\pm1),(\pm1,0)$ のみ条件を満たす |
次のステップ
- 発展1: $\gcd(x,y)=1$ かつ $x,y>0$（第 1 象限の互いに素な点）の個数を求める。軸上の点の扱いが変わる
- 発展2: $x^2+y^2=n$ ちょうどの点の個数 $r_2(n)=4(d_1(n)-d_3(n))$（$4k+1$ 型と $4k+3$ 型の約数の個数の差。記憶による式で未確認）を使って、円内の点数を別の方法で数えて突き合わせる
- 次回予告: 焦点帯 M1（DP・ゲーム）に戻る
検証ログ（本セッションで実行した数値）
- 総当たりとの一致（実測）: $x^2+y^2\le N$ の全点で $\gcd$ を取る総当たりと、$N=1..59$、$60..6000$ から無作為 60 個、および $1,2,5,10,100,1000,4096$ の計 126 件でプログラム全体の出力が全件一致。入力例（実測）: $N=1,10,100$ で 4、24、192
- 誤った方法（実測）: 原点を引き忘れると 118/126 件で誤答
- float の平方根（実測）: $N=10^9,\ 10^9+7,\ 3456789012,\ 10^{10}$ で float 版と
isqrt版の出力が一致。$m=(10^8+1)^2-1$ ではint(math.sqrt(m))=100000001、isqrt(m)=100000000 - 漸近との比（実測。導出した期待値 $6N/\pi$ との比）: $N=10^8$ で 1.00000027、$N=10^{10}$ で 1.00000004、$N=10^{12}$ で 0.9999999984（答えは $1909859313984$）
- 実行時間（実測）: $N=10^8$ 0.04 s、$N=10^{10}$ 0.15 s、$N=10^{12}$ 1.42 s
- 扱っていない範囲: $N=10^{12}$ の厳密値を外部の値と突き合わせていない（総当たりは $N\le6000$ まで、それ以上は漸近の比のみ）、$N\le10^{18}$ 向けの高速な円内格子点数、$r_2(n)$ の式を使った別解