この記事は ひとりアドベントカレンダーRosettaCodeで楽しむプログラミング Advent Calendar 2025の1日めの記事です。
定義
正整数 n で、その正の約数の平均が整数であるものを算術数（訳語？）という。
タスク
- 算術数の最初の100個を示せ。
- 1,000番目、10,000番目の算術数はいくつか。
- 最初の1,000個、10,000個の算術数のうち、合成数はいくつあるか。
- (stretch goal) タスク2,3を100,000、1,000,000で行え。
1は算術数であり、素数でなく、合成数でもないことに注意せよ。
考える
$n$ の全ての約数を取り出すには、$1 \leq p \leq \sqrt n$ の範囲で $n$ を割り切れる $p$ とその商を組み合わせればよい。
他にないかと検索したら、高校数学のテクニックが見つかった。
nが素因数分解されて $n = \prod_i {p_i}^{k_i}$ のとき、
- $n$ の約数の総和は $\displaystyle \prod_i \sum_{j=0}^{k_i} {p_i}^j = \prod_i \frac{p_i^{k_i+1} - 1}{p_i - 1}$
- $n$ の約数の個数は $\prod_i (k_i + 1)$
素因数分解の方が被除数がどんどん小さくなるのでもっと早く終わりそうな気がする。のでこれを使ってみよう。
import Data.Numbers.Primes
import Data.List
import Data.Bool
-- 算術数か、合成数か、両方を報告する
isArithCompo :: Int -> (Bool, Bool)
isArithCompo n = (mod sumFs numFs == 0, numFs > 2)
where
pks = map (\xs -> (head xs, length xs)) $ group $ primeFactors n -- 素因数分解
sumFs = product [div (pred $ p ^ succ k) (pred p) | (p,k) <- pks] -- 約数の総和
numFs = product $ map (succ . snd) pks -- 約数の個数
task1 :: [Int]
task1 = take 100 [i | i <- [1 ..], fst $ isArithCompo i]
task23 :: [(Int, Int)]
task23 = map sub [10^3,10^4,10^5,10^6]
where
sub cnt = (cnt, fst $ ics !! pred cnt, length $ filter snd $ take cnt ics)
ics = [(i, c) | i <- [1 ..], let (a, c) = isArithCompo i, a]
ghci> task1
[1,3,5,6,7,11,13,14,15,17
,19,20,21,22,23,27,29,30,31,33
,35,37,38,39,41,42,43,44,45,46
,47,49,51,53,54,55,56,57,59,60
,61,62,65,66,67,68,69,70,71,73
,77,78,79,83,85,86,87,89,91,92
,93,94,95,96,97,99,101,102,103,105
,107,109,110,111,113,114,115,116,118,119
,123,125,126,127,129,131,132,133,134,135
,137,138,139,140,141,142,143,145,147,149]
(0.02 secs, 2,107,056 bytes)
ghci> task23
[(1000,1361,782),(10000,12953,8458),(100000,125587,88219),(1000000,1228663,905043)]
(42.79 secs, 119,859,878,744 bytes)
Rosetta CodeにHaskell版がないので、お化粧して出しておきましょう。
main = do
forM_ (chunksOf 10 $ map (printf "%4d") task1) (\as -> sequence_ as >> putChar '\n')
putChar '\n'
forM_ task23 (\(i,a,c) -> do
printf "%dth arithmetic number is %d\n" i a
printf "Number of composite arithmetic numbers <= %d: %d\n\n" a c
)
OEIS
OEIS A00361にはHaskellのコードが出ているのだが、他の数列を参照していて自己完結していない。
動かすためには全部持ってこないといけない。
-- Numbers j such that the average of the divisors of j is an integer: sigma_0(j) divides sigma_1(j). Alternatively, numbers j such that tau(j) (A000005(j)) divides sigma(j) (A000203(j)).
a003601 n = a003601_list !! (n-1)
a003601_list = filter ((== 1) . a245656) [1..]
-- Characteristic function of arithmetic numbers
a245656 = (0 ^) . a054025 :: (Integral a, Integral t) => a -> t
-- Sum of divisors of n read modulo (number of divisors of n).
import Data.List (genericIndex)
a054025 n = genericIndex a054025_list (n-1)
a054025_list = zipWith mod a000203_list a000005_list
-- the number of divisors of n.
a000005 = product . map (+ 1) . a124010_row
-- Triangle in which first row is 0, n-th row (n>1) lists the exponents of distinct prime factors ("ordered prime signature") in the prime factorization of n.
a124010 n k = a124010_tabf !! (n-1) !! (k-1)
a124010_row 1 = [0]
a124010_row n = f n a000040_list where
f 1 _ = []
f u (p:ps) = h u 0 where
h v e | m == 0 = h v' (e + 1)
| m /= 0 = if e > 0 then e : f v ps else f v ps
where (v', m) = divMod v p
a124010_tabf = map a124010_row [1..]
-- The prime numbers.
import Data.List (genericIndex)
a000040 n = genericIndex a000040_list (n - 1)
a000040_list = base ++ larger where
base = [2, 3, 5, 7, 11, 13, 17]
larger = p : filter prime more
prime n = all ((> 0) . mod n) $ takeWhile (\x -> x*x <= n) larger
_ : p : more = roll $ makeWheels base
roll (Wheel n rs) = [n * k + r | k <- [0..], r <- rs]
makeWheels = foldl nextSize (Wheel 1 [1])
nextSize (Wheel size bs) p = Wheel (size * p)
[r | k <- [0..p-1], b <- bs, let r = size*k+b, mod r p > 0]
data Wheel = Wheel Integer [Integer]
-- a(n) = sigma(n), the sum of the divisors of n. Also called sigma_1(n).
a000203 n = product $ zipWith (\p e -> (p^(e+1)-1) `div` (p-1)) (a027748_row n) (a124010_row n)
-- Irregular triangle in which first row is 1, n-th row (n > 1) lists distinct prime factors of n.
import Data.List (unfoldr)
a027748 n k = a027748_tabl !! (n-1) !! (k-1)
a027748_tabl = map a027748_row [1..]
a027748_row 1 = [1]
a027748_row n = unfoldr fact n where
fact 1 = Nothing
fact x = Just (p, until ((> 0) . (`mod` p)) (`div` p) x)
where p = a020639 x -- smallest prime factor of x
-- Lpf(n): least prime dividing n (when n > 1); a(1) = 1. Or, smallest prime factor of n, or smallest prime divisor of n.
a020639 n = spf a000040_list where
spf (p:ps) | n < p^2 = n
| mod n p == 0 = p
| otherwise = spf ps
A000040「素数」はインデントが潰れていたのを復元した。
これは自己完結しているのだが、これはまたこれで、説明が欲しいコードだ。
どうにか集めきったと思いきや、
a054025_list = zipWith mod a000203_list a000005_list
の右辺二つがない。a000203 a000005 は n 番目を取り出す関数として定義され、_list がない。
a054025 n = genericIndex a054025_list (n-1)
の方の定義と合わせて
a054025 n = mod (a000203 (n-1)) (a000005 (n-1))
としてみたが、では算術数を求めてみようとすると
> take 100 a003601_list
Interrupted.
循環参照でもしているのか、ウンともスンとも言わない。
1つずれてるとかかなと
a054025 n = mod (a000203 n) (a000005 n)
とすると
ghci> take 100 a003601_list
*** Exception: divide by zero
ハイクを詠め、Reinhard Zumkeller＝サン。
A000040「素数」のHaskellコードを読み解く
自己完結しているのはいいが、やたらとややこしい。
読み解いた結果、反対側から説明した方がいいと思うのでそうする。
非常によく知られた、素数リストを無限に生成するコンパクトなコードはこんな感じ。
primes1 :: [Integer]
primes1 = 2 : filter isPrime1 [3, 5 ..]
isPrime1 :: Integer -> Bool
isPrime1 x = all ((0 /=) . mod x) $ takeWhile ((x >=) . (^ 2)) primes1
まず 2 は素数。その後、素数になりうるのは奇数だけなのでそういう候補 [3, 5 ..] で、素数判定に成功した値をリストに連ねる。
とその素数判定 isPrime は、今考えている値 x の平方根までの素数について、どれでも割り切れないこと、という感じに primes と isPrime が相互参照している。
これをもう少し効率化したいと思うと、「5以上の素数は6の倍数±1だけ」という知識が使える。
primes2 :: [Integer]
primes2 = 2 : 3 : filter isPrime2 [x | w <- [5,11..], x <- [w, w+2]]
isPrime2 :: Integer -> Bool
isPrime2 x = all ((0 /=) . mod x) $ takeWhile ((x >=) . (^ 2)) primes2
isPrime2 の定義は1と変わっていない。prime2 の方は、最初の素数のタネに3が追加され、候補の生成器が「奇数」から「6の倍数±1」に置き換えられている。
システマティックに見ると、これは、既出の素数2,3の最小公倍数（つまり積）である6を周期として、そのなかの数 0,1,2,3,4,5 から、2の倍数 を除くと 1,3,5 となり、3の倍数も除くと 1,5 が残る、という、周期の幅でだけエラトステネスのふるいをしている。完全なエラトステネスのふるいならこのパターンが篩全体にスタンプされるうえにさらに大きな素数の倍数についてもマークされるのを、ここまでで止めている、というハイブリッドな計算をしていると理解できる。
さて、OEISのプログラムに戻って、この「区間までの篩（に残った素数候補）」を表すのが Wheel データ型である。
data Wheel = Wheel Integer [Integer]
大仰な名前が付いているけど、type Wheel = (Integer, [Integer]) で十分なのでそうしてもよい。
完成した Wheel を回すことで、素数の候補がどんどん出てくる。
これが roll 関数の仕事。
roll :: Wheel -> [Integer]
roll (Wheel n rs) = [n * k + r | k <- [0..], r <- rs]
-- 掛け算など不要。
roll (Wheel s rs) = [b + r | b <- [0, s ..], r <- rs]
ある素数までについて構築された Wheel s rs があるとき、次の素数 p に関してこれを拡大するには、周期は単に掛け算 s * p で、そしてそこまでの数について roll で生成してみて、p で割りきれるものを捨てる、エラトステネスのふるいそのものを行えばできる。
（元のコードの nextSize のこと、その定義は roll を使わずにコードクローンしているので意味不鮮明）
expand w@(size, _) p = (size * p, [r | r <- takeWhile (size * p >) $ roll w, mod r p > 0])
小さな素数を好きなだけ用意して、なるべく大きな Wheel を作っておくと効率がいいだろう。
base = [2, 3, 5, 7, 11, 13, 17]
makeWheels = foldl expand (Wheel 1 [1])
これで役者が揃った。素数列の先頭は base で、それ以降は、makeWheels base を roll して得られる候補のうち、平方根までの素数で割りきれないものだけ残したリストになる。これは冒頭で示したスタイルの、wheel の大きさに関する一般形となっている。
primes = base ++ larger
where
larger = filter isPrime $ roll $ makeWheels base
isPrime n = all ((> 0) . mod n) $ takeWhile (\x -> x*x <= n) larger
Wheelを回して出てくる数が base に対して素であることは自明なので、isPrime では larger についてだけ調べる。
…動かない。roll $ makeWheels base の先頭 19 が isPrime に素数と認められるためには、（primes からではなく）larger から $\sqrt 19$ 超でない範囲の全ての素数を取ってきて調べる必要がある。が、19はまだその先頭なので、判定が終わるまで larger に出て行けない。デッドロック。
先頭は素数と保証できるとして、このチェックをバイパスするquick hackで切り抜ける。
あと、roll の結果の先頭にある 1 も邪魔なので取り除く。
primes = base ++ larger
where
larger = p : filter isPrime more
1 : p : more = roll $ makeWheels base
-- isPrime n = all ((> 0) . mod n) $ takeWhile (\x -> x*x <= n) larger -- これは変更なし
以上で、a000040_list の内容は完全に解明された。