この記事は ひとりアドベントカレンダーRosettaCodeで楽しむプログラミング Advent Calendar 2025の11日めの記事です。
定義
正整数 $n$ に対して $S(n)$ とは、$k!$ が $n$ で割りきれる最小の $k$
便宜上 $S(1) = 1$ とし、$S(0)$ は未定義とする。
タスク
- 1から50のケンプナー数を求めよ
- (stretch goal) 77135679311 から 77135679321 のケンプナー数を求めよ
考える
素朴に
順に $k!$ を求めては $n$ で割ってみる、というやり方はさすがに naive すぎるので、素因数分解によって少しだけ効率的に実行しよう。
$S(n) = k$ とは、$n, k!$ の素因数分解 $n = \prod p_j^{e_j}$, $k! = \prod p_j^{d_{k,j}}$ について、$e_j \leq d_j$ となる（そして $k-1$ ではならない）ということ。
これを調べるためには、$d_{k,j}$ の表があればよいだろう。
表を構築するには、$d_{1,j} = 0$ から始めて、$k$ を素因数分解した結果に $d_{k-1,j}$ を足し合わせていけばよい。
素数 $p_j$ をキー、個数 $d_{k,j}$ を値とする IntMap を、$k$ の順に無限リストとして生成する。
この表があれば、前から順に探せば $k$ がわかる。
import Data.List
import Data.Numbers.Primes
import qualified Data.IntMap as IM
kempner :: Int -> Int
kempner 1 = 1
kempner n = length $ takeWhile (cmp $ pf n) fpfs
where
pf = IM.fromDistinctAscList . map (\xs -> (head xs, length xs)) . group . primeFactors
fpfs = scanl (IM.unionWith (+)) IM.empty $ map pf [1 ..]
cmp im1 im2 = or [e > d | (p,e) <- IM.assocs im1, let d = IM.findWithDefault 0 p im2]
theRange :: [Int]
theRange = [77135679311 .. 77135679321]
実行。
ghci> map kempner [1 .. 50]
[1,2,3,4,5,3,7,4,6,5,11,4,13,7,5,6,17,6,19,5,7,11,23,4,10,13,9,7,29,5,31,8,11,17,7,6,37,19,13,5,41,7,43,11,6,23,47,6,14,10]
(0.02 secs, 3,726,512 bytes)
ghci> map kempner [77135679311 .. 77135679321]
[18191,Interrupted.
タスク1は順調に動いたけど、タスク2が二つめで沈黙してしまった。
どうなっているのか見てみる。
ghci> map primeFactors theRange
[[599,7079,18191],[2,2,2,2,3,41,39194959],[7,7,11,11,61,271,787],[2,89,449,809,1193]
,[3,3,5,13723,124909],[2,2,79,2753,88667],[17,4537392901],[2,3,12855946553]
,[29,2659851011],[2,2,2,5,7,31,8886599],[3,199,1289,100237]]
なるほどやらしい区間ですね。
うまくやる
階乗の中の素因数として $p$ は $p!$ で初めて出現し、以降、$2p!, 3p!, \cdots$ に出現する。
ならば $n = \prod p_j^{e_j}$ に対して $S(n) = \max(e_j \cdot p_j)$ でいいか、というとそうでもない。
例えば $p = 3$ のとき、素因数分解の中の3の個数は、$(1 \cdot p)! = 3!$ には1つ、$(2 \cdot p)! = 6!$ には2つ、$(3 \cdot p)! = 9!$ には3つでなく4つになる。
階乗の素因数分解の様子を $3 < p_j$ に注目して考えてみると、
| $k$ | 1 | 2 | … | $p$ | … | $2p$ | … | $3p$ | … | $mp$ |
|---|---|---|---|---|---|---|---|---|---|---|
| $k$の素因数分解 | 1 | 2 | … | $p$ | … | $2, p$ | … | $3, p$ | … | $m$ の素因数分解と $p$ |
| $d_{k,j} - d_{k-1,j}$ | 0 | 0 | … | 1 | … | 1 | … | 1 | … | $1 + {}$（$m$の素因数分解に出現する個数） |
$d_{k,j}$ が増加するのは $k$ が $p$ の倍数のときだけ ($k = mp$) で、具体的には $m$ の素因数分解に $p$ が出現する個数＋１だけ増える。
つまり、$e_j$ に対して、$m = 1,2,\dots$ について順に $m \cdot p_j$ を $p_j$ で割れる回数を数えて累積和をとる。
この値が $e_j$ 以上になる最小の $m$ を見つけたら、素因数 $p_j$ に関する $S(n)$ の候補が $m \cdot p_j$ となる。
それらの最大値が実際の答えになる。
kempner2 :: Int -> Int
kempner2 1 = 1
kempner2 n = maximum
[ m * p
| ps <- group $ primeFactors n, let p = head ps, let e = length ps
, let m = length $ takeWhile (e >) $ scanl (+) 0 $ map (succ . divtimes p) [1 ..]
]
where
divtimes p x = if mod x p /= 0 then 0 else succ (divtimes p (div x p))
ghci> map kempner2 [1 .. 50]
[1,2,3,4,5,3,7,4,6,5,11,4,13,7,5,6,17,6,19,5,7,11,23,4,10,13,9
,7,29,5,31,8,11,17,7,6,37,19,13,5,41,7,43,11,6,23,47,6,14,10]
(0.01 secs, 751,192 bytes)
ghci> map kempner2 theRange
[18191,39194959,787,1193,124909,88667,4537392901,12855946553,2659851011,8886599,100237]
(0.03 secs, 83,095,568 bytes)
うまくできました。
OEIS Haskell
OEIS A002034 にHaskellコードが掲載されている。
-- OEIS Reinhard Zumkeller
a002034 1 = 1
a002034 n = fromJust (a092495 n `elemIndex` a000142_list)
-- それが階乗のリストのいくつめだったか、つまり何の階乗か
-- A092495 Least factorial multiple of n.
a092495 n = fromJust $ find ((== 0) . (`mod` n)) $ a000142_list
-- 階乗のリストを前から見て、最初のnの倍数
-- A000142 Factorial numbers
a000142_list = 1 : zipWith (*) [1..] a000142_list
-- 階乗のリスト
またオヌシか。
そしてこれは定義通りの素朴すぎるコードで、タスク2にはどうにもならない。
OEIS Python (2nd)
Wren 版が、OEISのPythonコードのふたつめを参考にしたと書いてある。
from sympy import factorint
def valp(n, p):
s = 0
while n: n //= p; s += n
return s
def K(p, e):
if e <= p: return e*p
t = e*(p-1)//p*p
while valp(t, p) < e: t += p
return t
def A002034(n):
return 1 if n == 1 else max(K(p, e) for p, e in factorint(n).items())
K(p,e) は、自分の kempner2 のリスト内包表記で m を求める計算に相当していると思うのだけど、
t の初期値とかロジックがよくわからない。
下請けの valp(n,p) と合わせて、(p,e) ごとに二重ループを回している。
…よくわからない。
延長戦
と思っていたら、解説があった。
上のPythonプログラムのしていることは
- 素因数分解 $n = \prod {p_j}^{e_j}$ を求める
- $j$ ごとに ${p_j}^{e_j}$ のケンプナー数を求める
- その最大値が答え
と、自分の kempner2 と大方針は同じようだが、途中からなんでそんなことをしているのかよくわからない。
pedia にある
1918年、オーブリー・J・ケンプナーが最初に正しい計算アルゴリズムを与えた[1]。
のアルゴリズムを notebookLM さんに要約してもらった結果とも何だか合わない。
くだんのコードそのものを notebookLM さんに突っ込んだら、
valp(n, p)関数は、ルジャンドルの公式に基づいて$n!$（$n$の階乗）を割り切る素数$p$の最大の指数を計算します。
と、これまた違う方針なようだ。
乗りかかった船なので、ケンプナーさん本人によるアルゴリズムを忠実にコードにしてみる。
kempner4 :: Int -> Int
kempner4 1 = 1
kempner4 n = maximum [key (head ps) (length ps) | ps <- group $ primeFactors n]
where
key p a = a * pred p + sum ks
where
as = takeWhile (a >=) $ iterate (succ . (p *)) 1
ks = map fst $ scanr (\ai (_, r) -> divMod r ai) (0, a) as
実行速度は kemper2 と変わらないのだけど、「試行を排除するというルールに則ったアルゴリズム」なのだそうだ。
動作原理が理解できないのは、波動エンジンを設計図から組み立てているみたいで気持ち悪いのだけど、notebookLM さんに
Haskell版の
kempner4は、Kempnerが論文で提示した、指数 $a$ の特殊な $p$ 進展開を利用して $\mu(p^a)$ を直接的、かつ試行を排除して計算するRULEを、正確に実装しています。
と太鼓判を頂いたので、これで完了にする。