Day 175-Q4 — 2026-10-11 | min_25 篩(f(p^e) = p xor e の乗法的関数の和を N=10^9 まで)

2026-10-11 赤色 Master (Phase 8+) ★★★★★★★★★ 数論(焦点帯 M3)。過去回の Lucy の篩(Lucy_Hedgehog 法)は「$N$ 以下の素数の個数・和」までだった。今回はその表を第 1 段に使い、第 2 段で最小素因数ごとに分けて乗法的関数の和を組み立てる(min_25 篩)。$f(p)$ が素数 $p$ の多項式で書け、$f(p^e)$ がすぐ計算できれば、$f$ の具体形に頼らずに $O(N^{3/4}/\log N)$ 程度で和が出る

問題

乗法的関数 $f$ を $f(1)=1$、素数 $p$ と $e\ge1$ について

$$f(p^e)=p\oplus e$$

で定める($\oplus$ はビットごとの排他的論理和)。$\displaystyle\sum_{i=1}^{N}f(i)$ を $10^9+7$ で割った余りを求めよ。

入力形式

N

制約

  • $1\le N\le10^9$

入出力例

入力例1

6

出力例1

16

($f(1..6)=1,3,2,0,4,6$。$f(4)=2\oplus2=0$、$f(6)=f(2)f(3)=3\cdot2$)

入力例2

100

出力例2

3106

入力例3

1000000000

出力例3

338258534

(例 1・2 は篩による全数計算と一致。例 3 は本解の出力で、$N\le10^7$ までは全数計算と一致を確認した)

---

背景と出典

  • 関数 $f(p^e)=p\oplus e$ は LibreOJ 6053「简单的函数」($N\le10^{10}$、C++ 前提)の関数と同じだと記憶している(問題文は本セッションで未確認。本問の入出力と制約は Python 向けに作り直した)
  • min_25 篩: 中国の競技プログラミングで min_25 氏が広めた方法と記憶している(一次資料は未確認)。2 段からなる
    • 第 1 段(Lucy の篩と同じ): $v\in\{\lfloor N/i\rfloor\}$ のすべてについて、$v$ 以下の素数の個数 $\pi(v)$ と和 $\Sigma_p(v)$ を求める
    • 第 2 段: $S(v,k)=\sum_{2\le i\le v,\ \mathrm{lpf}(i)\ge p_k}f(i)$(lpf は最小素因数)を再帰で求める
  • 素数での値: 奇素数 $p$ は $p\oplus1=p-1$、$p=2$ は $2\oplus1=3$。よって $\sum_{p\le v}f(p)=\Sigma_p(v)-\pi(v)+2$($v\ge2$)。**$f(p)$ が $p$ の多項式で書ける**ことが第 1 段を使う条件で、xor 全体は多項式でないが $e=1$ のときだけ多項式になる

概念図

N = 100 の値の集合 {⌊N/i⌋}: 小さい側 lo[v](v ≤ 10)と大きい側 hi[i] = ⌊100/i⌋(i ≤ 10)12345678910lo[1..10]100503325201614121110hi[1..10] = 100, 50, 33, …, 10第 1 段(Lucy の篩): 素数 p ごとに v ≥ p² の g(v) から「最小素因数が p の合成数」を引く → π(v) と Σp(v) F_prime(v) = Σp(v) − π(v) + 2(奇素数は p−1、p=2 だけ 3)第 2 段: 2..v を最小素因数 p と指数 e で分けるS(v, k)素数 p ≥ p_k の部分p^e · j(j ≥ 2)p^(e+1)(j = 1)表から O(1)f(p^e)·S(⌊v/p^e⌋, t+1)。⌊v/p^e⌋ < p_{t+1}² なら表から(N=10⁹ で再帰 118 万回 → 5.2 万回)

ヒント(クリックで段階的に開示)

ヒント1(方向性)

$N=10^9$ まで全部篩うのは無理($10^7$ で 20 s)。値の種類 $\lfloor N/i\rfloor$ は $2\sqrt N$ 個しかない。素数の部分と合成数の部分に分ける。

ヒント2(アプローチ)

2 以上の整数 $i\le v$ を最小素因数 $p$ と、その指数 $e$ で分ける: $i=p^e\cdot j$($j$ の最小素因数 $>p$)。$f(i)=f(p^e)f(j)$。$j=1$ の場合($i=p^e$)と $j\ge2$ の場合を分けて足すと、$S(v,k)$ の漸化式になる。

ヒント3(誘導)

$$S(v,k)=\Bigl(F_{\text{prime}}(v)-F_{\text{prime}}(p_{k-1})\Bigr)+\sum_{t\ge k,\ p_t^2\le v}\ \sum_{e\ge1,\ p_t^{e+1}\le v}\Bigl(f(p_t^e)\,S(\lfloor v/p_t^e\rfloor,t+1)+f(p_t^{e+1})\Bigr)$$

答えは $S(N,1)+1$。$\lfloor v/p^e\rfloor<p_{t+1}^2$ なら $S$ の中身は素数だけなので、再帰せずに表から引ける($N=10^9$ で呼び出し回数が約 23 分の 1 になる)。

模範解答 (Python)

模範解答を表示する(まず自力で解いてから)
import sys
from math import isqrt

MOD = 10**9 + 7


def main():
    n = int(sys.stdin.readline())
    if n == 1:
        print(1)
        return
    r = isqrt(n)
    # 値の集合 { n // i } を、小さい側 lo[v](v <= r)と大きい側 hi[i] = 値 n // i(i <= r)に分けて持つ
    lo_c = [v - 1 for v in range(r + 1)]                  # v 以下の 2 以上の整数の個数(篩う前)
    lo_s = [v * (v + 1) // 2 - 1 for v in range(r + 1)]   # その和
    hi_c = [0] + [n // i - 1 for i in range(1, r + 1)]
    hi_s = [0] + [(n // i) * (n // i + 1) // 2 - 1 for i in range(1, r + 1)]
    primes = []
    for p in range(2, r + 1):
        if lo_c[p] == lo_c[p - 1]:
            continue                                      # 合成数
        primes.append(p)
        pc, ps = lo_c[p - 1], lo_s[p - 1]                 # p 未満の素数の個数と和
        p2 = p * p
        lim = min(r, n // p2)
        # 大きい側: 値 n//i >= p^2 のもの
        for i in range(1, lim + 1):
            d = i * p
            if d <= r:
                c, s = hi_c[d], hi_s[d]
            else:
                v = n // d
                c, s = lo_c[v], lo_s[v]
            hi_c[i] -= c - pc
            hi_s[i] -= (s - ps) * p
        # 小さい側: v >= p^2 のもの(大きい v から)
        for v in range(r, p2 - 1, -1):
            w = v // p
            lo_c[v] -= lo_c[w] - pc
            lo_s[v] -= (lo_s[w] - ps) * p
    # f(p) = p xor 1: 奇素数は p - 1、p = 2 は 3。素数についての和 F を値の集合の上で表にする
    lo_f = [lo_s[v] - lo_c[v] + (2 if v >= 2 else 0) for v in range(r + 1)]
    hi_f = [0] + [hi_s[i] - hi_c[i] + 2 for i in range(1, r + 1)]
    np_ = len(primes)
    pref = [0] * (np_ + 1)                                 # 最初の k 個の素数での f(p) の和
    for k, p in enumerate(primes):
        pref[k + 1] = pref[k] + (p ^ 1)
    primes.append(r + 1)                                   # 番兵: (r+1)^2 > n

    def S(v, k):
        # 2..v のうち最小素因数が primes[k] 以上の整数 i の f(i) の和
        res = (lo_f[v] if v <= r else hi_f[n // v]) - pref[k]
        t = k
        p = primes[t]
        while p * p <= v:
            q = primes[t + 1]
            qq = q * q
            pe = p
            e = 1
            while pe * p <= v:
                w = v // pe
                if qq > w:                                 # この先は素数だけ: 再帰せずに表から
                    sub = (lo_f[w] if w <= r else hi_f[n // w]) - pref[t + 1]
                else:
                    sub = S(w, t + 1)
                res += (p ^ e) * sub + (p ^ (e + 1))
                pe *= p
                e += 1
            t += 1
            p = q
        return res

    print((S(n, 0) + 1) % MOD)


main()

計算量: 第 1 段 $O(N^{3/4}/\log N)$、第 2 段は $O(N^{1-\varepsilon})$ と言われる(厳密な評価は本セッションで未確認)。実測で $N=10^9$ が 1.24〜1.32 s(第 2 段 0.48 s)、$3\times10^9$ が 2.75 s、$10^{10}$ が 6.80 s(Python 3.14.7、M1 MacBook Pro 8GB、負荷平均 10〜12、起動込み)。制約を $10^9$ にしたのはこの差による。

Step-by-Step 解説

1値の集合を 2 つの配列で持つ

$r=\lfloor\sqrt N\rfloor$。$v\le r$ は lo[v]、$v=\lfloor N/i\rfloor>r$ は hi[i]($i\le r$)に置く。これで $2r$ 個の値をハッシュなしで引ける。

2第 1 段(Lucy の篩)

最初は「2 以上の整数をすべて素数とみなした」個数 $v-1$ と和 $\frac{v(v+1)}2-1$。素数 $p$ を小さい順に処理し、$v\ge p^2$ について

$$g(v)\leftarrow g(v)-p^{k}\bigl(g(\lfloor v/p\rfloor)-g(p-1)\bigr)\quad(k=0\text{ は個数},\ k=1\text{ は和})$$

最小素因数が $p$ の合成数を取り除く。lo は大きい $v$ から更新すれば、右辺の古い値を壊さない。

3素数での $f$ の和

$F_{\text{prime}}(v)=\Sigma_p(v)-\pi(v)+2$($v\ge2$)。$p=2$ だけ $f(2)=3=(2-1)+2$ で、補正の $+2$ が要る。

4第 2 段(再帰)

ヒント3 の式。$f(p^e)\cdot S(\lfloor v/p^e\rfloor,t+1)$ は「$p^e\cdot j$($j\ge2$、$j$ の最小素因数 $>p$)」、$f(p^{e+1})$ は「$j=1$ で指数 $e+1$」の分。$e$ の上限を $p^{e+1}\le v$ にしているのは、$j\ge2$ の項が意味を持つのは $p^e\cdot p_{t+1}\le v$ のときだけで、$j=1$ の $p^{e+1}$ もこのループで拾うため($p^1$ 自体は素数の部分 $F_{\text{prime}}$ に入っている)。

5再帰を減らす

$w=\lfloor v/p^e\rfloor$ が $p_{t+1}^2$ 未満なら、$S(w,t+1)$ は素数の和だけ: $F_{\text{prime}}(w)-F_{\text{prime}}(p_t)$。表から引けば関数呼び出しが要らない。$N=10^9$ で呼び出しは 118 万 5818 回 → 5 万 1925 回になり、1.50 s → 1.24 s(実測。負荷の高い時間帯に 1.88 s → 1.19 s)。

迷った他案と採らなかった理由

案計算量採否理由
最小素因数の篩で全数計算$O(N)$検算用のみ$N=10^7$ で 20.1 s(実測)
min_25、第 2 段の再帰をそのまま$O(N^{3/4}/\log N)$ 程度不採用$N=10^9$ で 1.50〜1.88 s、呼び出し 118 万回
第 1 段をリストの内包表記でまとめて更新同じ不採用$N=10^9$ で 1.86〜2.00 s と遅くなった(実測)。素数ごとに作る一時リストの分が上回る
Powerful Number 法($f=g*h$、$g(p)=f(p)$ の $g$ の和を杜教篩で)$O(N^{2/3})$ 程度不採用$g$ を $p=2$ だけ補正した $\varphi$ にする必要があり、和の式が増える(本セッションでは未実装)
min_25、素数だけの部分は表から(本解)同上採用1.24〜1.32 s

前提が変わったら

  • **$N=10^{10}$** なら本解で 6.80 s(実測)。C++ なら 1 s 未満が普通と記憶している。Python では第 1 段を numpy の配列演算にするのが次の手(本セッションの環境には numpy が無く未計測)
  • **$f(p)$ が $p$ の多項式でない**(例: $f(p)=p\oplus3$)なら、第 1 段の「素数の $p^k$ の和」に分けられない。$p\bmod 4$ ごとの素数の個数・和など、別の表が要る
  • **和ではなく $Q$ 個の $N$ について答える**なら、第 1 段は $N$ ごとにやり直し。$N$ が共通なら値の集合も共通なので、第 2 段だけ繰り返す

よくあるミス

ミス原因正しい書き方
$f(2)$ を $2-1=1$ にする$p\oplus1=p-1$ を偶素数にも使う$2\oplus1=3$。素数の和に $+2$
答えに $f(1)=1$ を足し忘れる$S$ は 2 以上の和$S(N,1)+1$
$f(p^e)$ を $f(p)^e$ にする完全乗法的と混同$f(4)=2\oplus2=0$、$f(2)^2=9$。$N=4$ で 6 にならない
lo を小さい $v$ から更新古い値 $g(\lfloor v/p\rfloor)$ を先に書き換える大きい $v$ から
$v=\lfloor N/i\rfloor$ を浮動小数の N / i で求める$10^9$ 近くで切り捨てがずれる整数の //

次のステップ

  • 発展1: $f(p^e)=p\oplus e$ を $N=10^{10}$ で 2 s 以内に(第 1 段を配列演算に、第 2 段を非再帰に)
  • 発展2: $f(p)=p\oplus3$ のように $f(p)$ が多項式でない関数で、$p\bmod4$ ごとの表を作る
  • 発展3: Powerful Number 法で同じ和を求め、min_25 と速さを比べる

検証ログ(本セッションで実行した数値)

  • 全数計算との照合(実測): 最小素因数の篩で $f(1..2\times10^5)$ の累積和を作り、$N=1..199$、$[200,2\times10^5]$ の乱数 150 個、素数の平方 21 個、素数の平方 −1 の 9 個、計 379 個で一致。さらに $N=10^6$(518881640)、$N=10^7$(202772081)でも全数計算と一致
  • 改良前後の一致(実測): 再帰そのままの版と本解で $N=1,2,6,99991,10^6$ の出力が一致、$10^9$ も同じ 338258534
  • 実行時間(実測、起動込み): $N=10^7$ 0.11 s、$10^8$ 0.28 s、$10^9$ 1.24〜1.32 s、$3\times10^9$ 2.75 s、$10^{10}$ 6.80 s。全数計算は $10^7$ で 20.1 s
  • 計測時の負荷: 負荷平均 10〜18。±30% 程度ぶれる
  • 扱っていない範囲: $N=10^{10}$ を 2 s で、Powerful Number 法、LibreOJ 6053 の問題文・公式の答えとの照合

自己評価

自分の回答

気づき・メモ