N-smooth数とハミング数

この記事は ひとりアドベントカレンダーRosettaCodeで楽しむプログラミング Advent Calendar 2025の9日めの記事です。

定義

n-smooth数とは、nより大きい素因数を持たない数。（nは素数）

7-smooth数はhumble数ともいう。
5-smooth数はHamming数ともいう。

タスク

• $n = 2$ から $n = 29$ までについて、最初の25個の n-smooth数を求めよ。
• $n = 3$ から $n = 29$ までについて、3,000個めからの n-smooth数を3つ求めよ。
• (stretch goal) $n = 503$ から $n = 521$ までについて、30,000 個めからの n-smooth数を20個求めよ。

ヒント：ストレッチゴールの n は結局 503,509,521 の3つである。

やる

Humble数のコードが流用できるから簡単ッッッ!!簡単ッッ!!!
使う素数が少し増えたので、Set を使う版でやる。

import Data.Numbers.Primes
import qualified Data.Set as S

nSmooth n = gen (S.singleton 1)
  where
    ps = takeWhile (n >=) primes
    gen s = x : gen (S.union s1 $ S.fromList $ map (x *) ps)
      where
        (x, s1) = S.deleteFindMin s

• タスク1

task1 f = [(n, take 25 $ f n) | n <- filter isPrime [2 .. 29]]

ghci> task1
ghci> task1 nSmooth
[(2,[1,2,4,8,16,32,64,128,256,512,1024,2048,4096,8192,16384,32768,65536,131072
,262144,524288,1048576,2097152,4194304,8388608,16777216])
,(3,[1,2,3,4,6,8,9,12,16,18,24,27,32,36,48,54,64,72,81,96,108,128,144,162,192])
,(5,[1,2,3,4,5,6,8,9,10,12,15,16,18,20,24,25,27,30,32,36,40,45,48,50,54])
,(7,[1,2,3,4,5,6,7,8,9,10,12,14,15,16,18,20,21,24,25,27,28,30,32,35,36])
,(11,[1,2,3,4,5,6,7,8,9,10,11,12,14,15,16,18,20,21,22,24,25,27,28,30,32])
,(13,[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,18,20,21,22,24,25,26,27,28])
,(17,[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,20,21,22,24,25,26,27])
,(19,[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,24,25,26])
,(23,[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25])
,(29,[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25])]
(0.04 secs, 2,604,400 bytes)

• タスク2

task2 f = [(n, take 3 $ drop 2999 $ f n) | n <- filter isPrime [3 .. 29]]

ghci> task2 nSmooth
[(3,[91580367978306252441724649472,92829823186414819915547541504,94096325042746502515294076928])
,(5,[278942752080,279936000000,281250000000])
,(7,[50176000,50331648,50388480])
,(11,[2112880,2116800,2117016])
,(13,[390000,390390,390625])
,(17,[145800,145860,146016])
,(19,[74256,74358,74360])
,(23,[46552,46575,46585])
,(29,[33516,33524,33534])]
(0.20 secs, 112,774,472 bytes)

• タスク3

task3 f = [(n, take 20 $ drop 29999 $ f n) | n <- filter isPrime [503 .. 521]]

ghci> task3 nSmooth
[(503,[62913,62914,62916,62918,62920,62923,62926,62928,62930,62933,62935,62937,62944,62946,62951,62952,62953,62957,62959,62964])
,(509,[62601,62602,62604,62607,62608,62609,62611,62618,62620,62622,62624,62625,62626,62628,62629,62634,62640,62643,62645,62646])
,(521,[62287,62288,62291,62292,62300,62304,62307,62308,62310,62315,62320,62321,62322,62325,62328,62329,62330,62331,62335,62336])]
(28.57 secs, 6,299,432,504 bytes)
ghci>

RosettaCodeの解

nで括られたものの一つ、5-smooth数、Hamming数の解に基づいている、と。ふむふむ。

ハミング数のページに行くと、The classic version と称して、自分のhumble数の解のやり方は古くさいと一蹴されている。

確かにこれだと、x2 → x3 と x3 → x2 の両方が生成され、先に進むに従って重複が増える。
Avoiding generation of duplicates では、構造を少し変えている。
これで漏れなく生成できると思いつけるのすごいな。

nSmoothA n = foldr step [1] $ takeWhile (n >=) primes
  where
    step p (x:xs) = ys
      where
        ys = x : merge xs (map (p *) ys)

    merge xxs@(x:xs) yys@(y:ys) =
      case compare x y of
        LT -> x : merge xs yys
--      EQ -> x : merge xs  ys -- 不要と示すためあえて定義しない
        GT -> y : merge xxs ys
    merge [] ys = ys
    merge xs [] = xs

なんか

    step p xs = ys
      where
        ys = merge xs (map (p *) ys)

こうでも動いちゃって、しかもこっちの方がすこし速いんだけど、これ ys で loop して返って来ちゃ駄目なやつでは？どうして動くん？

ghci> task1 nSmoothA
[略]
(0.04 secs, 2,069,392 bytes)
ghci> task2 nSmoothA
[略]
(0.06 secs, 15,277,136 bytes)
ghci> task3 nSmoothA
[略]
(0.31 secs, 122,634,984 bytes)

重複を排除しない classic 版もやってみよう。
と思ったら、むしろ classic 版を任意の n でするの難しくないか。

nSmoothB n = xs
  where
    xs = 1 : foldr step [] (takeWhile (n >=) primes)
    step p ys = merge ys $ map (p *) xs

    merge xxs@(x:xs) yys@(y:ys) =
      case compare x y of
        LT -> x : merge xs yys
        EQ -> x : merge xs  ys
        GT -> y : merge xxs ys
    merge [] ys = ys
    merge xs [] = xs

ghci> task1 nSmoothB
[略]
(0.03 secs, 2,119,304 bytes)
ghci> task2 nSmoothB
[略]
(0.11 secs, 46,274,152 bytes)
ghci> task3 nSmoothB
[略]
(1.79 secs, 1,107,841,288 bytes)

確かに重複を削除する必要があり、かつその必要がない版より遅い。

次の Explicit multiples reinserting は、最初に示した Set を使う版の nSmooth と同じ発想のようだ。
ただし Data.List.Ordered.union でも使っているのか、リストで整列を保ってマージしている。
特にnが大きいn-smooth数では Set や IntSet を使うべきだろう。

次の Enumeration by a chain of folded merges は、
Avoiding generation of duplicates の foldr を手で展開した変種で、
さらに merge でフン詰まりにならないように、左の先頭は最小として通す変種を使っている。

Direct calculation through triples enumertion は、Hamming 数の方の最終タスクが「100万個めのハミング数を求めよ」であるために、列挙するのではなく直接取り出す方法を考えている。
logとか出てきて雲行きが怪しい上に、コードもかなりゴチャついている。
さらに Using loops for a faster code ... では大変な事になっている。

DDJでの議論が元になっているらしいが、これは一体何をしているんだ？

ひゃくまん個めのハミング数

今ある、列挙するタイプので計算するとこんな感じ。

ghci> nSmooth 5 !! 999999
519312780448388736089589843750000000000000000000000000000000000000000000000000000000
(15.53 secs, 3,176,170,024 bytes)
ghci> nSmoothA 5 !! 999999
519312780448388736089589843750000000000000000000000000000000000000000000000000000000
(1.24 secs, 489,350,176 bytes)
ghci> nSmoothB 5 !! 999999
519312780448388736089589843750000000000000000000000000000000000000000000000000000000
(2.58 secs, 1,089,587,144 bytes)

掲載されているコードを直しながら引き写す。

• foldl' に foldl_ という別名付けて何が楽しいの。
• nthHamのガード、n <= 0 は利用側のミス、ここまでフォローしなくともよいのでは。
  n < 2 はつまり $n = 1$ のとき。等式を分けておく。
  w >= 1, m < 0, m >= nb は、内部での推定値が外れたので計算が失敗しました、とエラーを吐いて終わる。つまりかなり身勝手な理由での部分関数になっているようだ。
  とりあえず見にくいだけなので消す。
  これらからだけ参照されている変数 nb も消す。
  otherwiseだけになった。
• インデント記法と ; 混ぜんなよ…
• v がこのコードの肝の一つ、n から n 番目のハミング数の log 2 を推定した値。多分魔法。
• estval というヘンな関数は一度だけ使われ、hi = v + 1/v, hi - w = v - 1/v と、おそらく何となくで幅を持たせている。
  w の代わりに lo = v - 1/v を追加してこの関数は消す。
• properFraction もわざわざ別名にしなくていい。
• (c,b) を作るリスト内包表記がこのプログラムの中核。その前に foldl'のステップ関数を見てみると

(\(c,b) (i,t)-> let c2=c+i in c2 `seq`
                     case t of []-> (c2,b);[v]->(c2,v:b) )

名前を付けて整理すると

step (c,b) (i,t) = c2 `seq`
  case t of
    [ ] -> (c2,   b)
    [v] -> (c2, v:b)
  where
    c2 = c + i

タプルの左は整数で足し、右はたかだか要素1のリストで、逆順に連結する。
c2 `seq を無視すれば、step (c,b) (i,t) = (c + i, t ++ b) でよい。
さらにこの後 foldl' を抜けた b は sortBy されてしまうので順不同。
unzip を使えばもっと簡単に書けるとわかる。

ここまででこうなった。

import Data.List
import Data.Function

main = let (r,t) = nthHam 1000000 in print t >> print (trival t)

trival (i,j,k) = 2^i * 3^j * 5^k

-- n個めのハミング数X=2^i*3^j*5^kのlog2 Xと(i,j,k)を求める
nthHam :: Int -> (Double, (Int, Int, Int))
nthHam 1 = (0, (0,0,0))
nthHam n = b !! m
  where
    lb3 = logBase 2 3
    lb5 = logBase 2 5
    lb30_2 = logBase 2 30 / 2
    v = (6 * lb3 * lb5 * fromIntegral n) ** (1/3) - lb30_2
    hi = v + 1/v
    lo = v - 1/v
    m = fromIntegral (sum cs - n)
    b = sortBy (flip compare `on` fst) $ concat bs
    (cs,bs) = unzip
--     (リスト内包表記)

リスト内包表記の中では、
$2^i \cdot 3^j \cdot 5^k \leq 2^v$ を満たす $(i,j,k)$ の組を
$k,j$ を全ての組み合わせに関して二重ループで探している。
このとき、不等式を満たす最大の $i$ が、それより小さいものも含めて $i+1$ 個のハミング数が推定値以下にあることを c でカウントする。
また見つけた $(i,j,k)$ を、そのハミング数の2底の対数値 r とともに b に取り出す。
という計算をしているようだ。

$j,k,p,q$ の計算は計算誤差を考えると

(k,p) <- zip [0 ..] $ takeWhile (hi >=) $ map ((lb5 *) . fromIntegral) [0 ..]

のようにした方が良さそうだが、わかりやすさ優先で直す。

結局 $p = \log_2 (5^k), q = \log_2 (3^j \cdot 5^k)$ なので、
上限 hi までの残り hi - q の整数部が、hi を超えないハミング数の最大の $i$ になる。
このときそのハミング数の対数値は $i + q$ である。
これを r = hi - frac と作らなくても i + q でよい。

リスト内包表記の結果値の左側は $i + 1$ で確かに個数、右側はまたリスト内包表記になっているが、
ガード frac < w を満たすとき1要素、満たさないとき空リスト、という条件付けをしているだけである。
このガードは無くしても sortBy が遅くなるだけで結果に悪影響はないと思うというのはともかく、
i + frac = hi - q より両辺に i を足して i + frac = hi - q < i + w となる。
さらに移項すると hi - w < i + q = r となる。hi - w = v + 1/v - 2/v = v - 1/v = lo なので結局 lo < r のことで、
properFraction も frac もいらなかった。

        [ (succ i, [(r,(i,j,k)) | lo < r] )
        | (k,p) <- zip [0 ..] [0, lb5 .. hi]
        , (j,q) <- zip [0 ..] [q, q + lb3 .. hi]
        , let i = floor $ hi - q
        , let r = fromIntegral i + q
        ]

リスト内包表記もここまですっきりした。

ghci> trival $ snd $ nthHam 1000000
519312780448388736089589843750000000000000000000000000000000000000000000000000000000
(0.13 secs, 17,481,440 bytes)

とんでもなく速い。というか、こんなに対数で計算して誤差でコケないものなんだな、という感想。
（本当はちょくちょくズレているのかもしれない。未検証。）

整数化

Double で計算するのはやっぱり不安なので、その精神を保ったまま Integer で現物を扱う形に直したい。

まず、何らかの上限が与えられたとき、それ以下のハミング数の個数を同様に数えつつ、上限近辺の値を $(j,k)$ の異なる全てのものについて挙げる関数を書いてみる。

hamsLE v = (sum cs, sortBy (flip compare) rs, length rs)
  where
    (cs,rs) = unzip
      [ (length rc, last rc)                        -- rcの長さがiの可能な個数、末尾がrの最大値
      | p <- takeWhile (v >=) $ iterate (5 *) 1     -- p <- 5 ^ k が v 以下のものについて
      , q <- takeWhile (v >=) $ iterate (3 *) p     -- q <- p * 3 ^ j が v 以下のものについて
      , let rc = takeWhile (v >=) $ iterate (2 *) q -- r <- q * 2 ^ i で v 以下のものを列挙し
      ]

ghci> hamsLE 200
(46,[108,120,125,128,135,144,150,160,162,180,192,200],12)
ghci> drop 34 $ take 46 $ nSmooth 5
    [108,120,125,128,135,144,150,160,162,180,192,200]

うまくいっているようだ。では、このロジックを使って N 番目のハミング数を見つける計算を作る。
適当な上限値の初期値で hamsLE を試し、

• reverse (sort rs) !! sum cs - n が取り出せるなら、それが答え
• sum cs が N 未満なら、もっと上限を引き上げてやり直し
• sum cs - n が length rs 以上なら、もっと上限を下げてやり直し
  と反復をかける。

nthHam1 :: Int -> Integer
nthHam1 1 = 1
nthHam1 n = loop v0
  where
    l0 = (logBase 2 3 * logBase 2 5 * 6 * fromIntegral n)**(1/3) - logBase 2 30 / 2
    v0 = floor (2 ** (l0 + 1/l0))
    loop v
      | cnt < n = traceShow (n,'u',v) $ loop $ v * 30
      | length rs <= cnt - n = traceShow (n,'d',v) $ loop $ minimum rs
      | otherwise = sortBy (flip compare) rs !! (cnt - n)
      where
        (cs,rs) = unzip
          [ (length rc, last rc)
          | p <- takeWhile (v >=) $ iterate (5 *) 1
          , q <- takeWhile (v >=) $ iterate (3 *) p
          , let rc = takeWhile (v >=) $ iterate (2 *) q
          ]
        cnt = sum cs

初期値は Rosetta Code 版の謎の式を利用した。
やり直しがかかったときはデバッグメッセージでわかるようにしてみた。

ghci> nthHam1 1000000
519312780448388736089589843750000000000000000000000000000000000000000000000000000000
(0.77 secs, 227,279,704 bytes)

ほう。

ghci> filter id $ zipWith (/=) (nSmooth 5) (map nthHam1 [1 ..])
(2,'d',4)
(849,'u',18653879)
(849,'d',559616370)
(849,'d',279936000)
(849,'d',140625000)
(849,'d',70778880)
(1092,'u',93234427)
(1092,'d',2797032810)
(1092,'d',1399680000)
(1092,'d',703125000)
(1092,'d',353894400)
(2951,'u',241538221885)
(2951,'d',7246146656550)
(2951,'d',3623878656000)
(2951,'d',1813985280000)
(2951,'d',911250000000)
(5061,'u',58588070198025)
(5061,'d',1757642105940750)
(5061,'d',878906250000000)
(5061,'d',439804651110400)
(5061,'d',220150628352000)
(7335,'u',4745135246135350)
(7335,'d',142354057384060500)
(7335,'d',71191406250000000)
(7335,'d',35624176739942400)
(7335,'d',17832200896512000)
Interrupted.

これ以降ログが出る気配はないので、かなり精度は高いようだが、完璧でもないようだ。
実際、これらの値でオリジナルを実行すると、内部事情での異常停止を書いたガードが発動する。

この後さらに
Using loops for a faster code, and a narrower band to save space
Using "roll-your-own" extended precision logarithm values in the error band to extend range
でコードは下手なCソースのようにインデントを深めていくが、その前にこの不完全さを黙っているのは不誠実ではないのかな。
ガードを書いたからそれでいい、ということなのかな？

おしまい。

今日のネタは簡単だったはずなのに。