謙虚な数

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

定義

7より大きな素因数を持たない正整数。

$i,j,k,m \geq 0$ によって $2^i \cdot 3^j \cdot 5^k \cdot 7^m$ と表せる数。

タスク

• 最初の50個の humble 数を求めよ。
• 1桁～n桁（できるところまで）の humble 数の個数をそれぞれ数えよ。

考える

らしく

フィボナッチ数列をつくるときのようにして、humble数の昇順のリストを遅延計算で構築することを考える。

• 先頭が 1 なのは確実。
• ある humble 数に対して、それに2,3,5,7を掛けたものもまたhumble数。
• 昇順に生成されると仮定すると、それらを昇順にマージすれば、結果も（仮定においたように）昇順になる。

humbleNumbers = 1 : merge (merge hns2 hns3) (merge hns5 hns7)
  where
    hns2 = map (2 *) humbleNumbers
    hns3 = map (3 *) humbleNumbers
    hns5 = map (5 *) humbleNumbers
    hns7 = map (7 *) humbleNumbers
-- インライン化
humbleNumbers =
  merge (merge (map (2 *) humbleNumbers) (map (3 *) humbleNumbers))
        (merge (map (5 *) humbleNumbers) (map (7 *) humbleNumbers))

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

4つの無限ジェネレータ hns2,hns3,hns5,hns7 からの出力を統合して流すイメージ。merge も3つのインスタンスが持続的に存在する。

ジェネレータによる生成を有限で停止させる代わりに、無限にジェネレータが発生するスタイルも考えられる。Haskellにおけるナイーブなエラトステネスのふるいの実装のようなイメージ。

humbleNumbers1 = 1 : gen humbleNumbers1
  where
    gen (x:xs) = merge (map (x *) [2,3,5,7]) $ gen xs

-- 第1引数が空になるので定義を追加
merge [] ys = ys

ジェネレータが有限で停止するので一見こちらの方が良く感じられるが、大きな数 $\times 5, \times 7$ がなかなかハケずにその間 merge が多段に残り続けるので効率よくないのもエラトステネスのふるいと同じ構造。

動かしてみる。

ghci> humbleNumbers1
[1*** Exception: stack overflow

Ctrl+Cすら効かずに自爆するまで何もできなかった。
gen の呼び出しが何もcarを生成しないまま再帰呼び出しするので、

humbleNumbers1 =>
1 : gen (1 : xs) =>
1 : merge [2,3,5,7] (gen xs) =>

の次に xs の先頭が何かを調べようとすると、xs とは merge [2,3,5,7] (gen xs) 自体なので、同じ計算の繰り返しになってしまった。gen が少なくとも何か一つ結果を生成する必要がある。

humbleNumbers1 が昇順であると仮定すると、gen が生成するべき数列の中で最小の要素は x * 2 でよさそうなので（厳密な証明ができるか知らない）、最初の結果としてそれを先に出す。

humbleNumbers1 = 1 : gen humbleNumbers1
  where
    gen (x:xs) = x * 2 : merge (map (x *) [3,5,7]) $ gen xs

ghci> take 20 humbleNumbers
[1,2,3,4,5,6,7,8,9,10,12,14,15,16,18,20,21,24,25,27]
ghci> take 20 humbleNumbers1
[1,2,3,4,5,6,7,8,9,10,12,14,15,16,18,20,21,24,25,27]

どちらも問題なさそうだ。
タスクを実行してみる。

import Control.Monad
import Text.Printf

test1 :: Int -> [Int] -> IO ()
test1 lim xs = do
  print $ take 50 xs
  forM_ (take lim $ countf xs) (\(l,c) -> printf "%5d have %2d digits\n" c l)

countf :: [Int] -> [(Int,Int)]
countf xs = loop 1 xs
  where
    loop lim hns = (lim, length hns1) : loop (succ lim) hns2
      where
        (hns1, hns2) = span (10^lim >) hns

ghci> test1 12 humbleNumbers
[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,40,42,45,48,49,50,54,56,60,63,64,70,72,75,80,81,84,90,96,98,100,105,108,112,120]
    9 have  1 digits
   36 have  2 digits
   95 have  3 digits
  197 have  4 digits
  356 have  5 digits
  579 have  6 digits
  882 have  7 digits
 1272 have  8 digits
 1767 have  9 digits
 2381 have 10 digits
 3113 have 11 digits
 3984 have 12 digits
(0.05 secs, 21,666,336 bytes)
ghci> test1 12 humbleNumbers1
[出力省略]
(5.42 secs, 7,182,667,072 bytes)

結果は圧倒的だった。

判定器による解法

与えられた整数がhumble数であるかどうか、2,3,5,7で割れるだけ割って、他の素因数がなければそう、というロジックで判定するやり方。

Kotlin, C++, Lua などのコード例がこのやり方でしている。

isHumble :: Int -> Bool
isHumble 1 = True
isHumble x
  | mod x 2 == 0 = isHumble $ div x 2
  | mod x 3 == 0 = isHumble $ div x 3
  | mod x 5 == 0 = isHumble $ div x 5
  | mod x 7 == 0 = isHumble $ div x 7
isHumble _ = False

7で一度割れたとき、もう一度2,3,5で割ってみる無駄が気になるが、シンプルなのは確か。

test2 :: Int -> (Int -> Bool) -> IO ()
test2 lim prop = do
  print $ take 50 $ filter prop [1 ..]
  forM_ [1 .. lim] (\l -> do
    let c = length $ filter prop [10^pred l .. pred $ 10^l]
    printf "%5d have %2d digits\n" c (l :: Int)
    )

ghci> test2 7 isHumble
[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,40,42,45,48,49,50,54,56,60,63,64,70,72,75,80,81,84,90,96,98,100,105,108,112,120]
    9 have  1 digits
   36 have  2 digits
   95 have  3 digits
  197 have  4 digits
  356 have  5 digits
  579 have  6 digits
  882 have  7 digits
(41.37 secs, 14,026,415,744 bytes)

桁数を増やすと判定器で調べる対象が10倍10倍と指数的に増えるので、当然どんどん重くなる。

humble数はその定義から指数的にまばらになっていくので、生成器の出力を数えるやり方だと、桁数が大きくなってもあまり数が増えないためこれだけ差が付いた。

もっと真面目に、生成も判定もすることなく、個数を数えるのは面倒な問題になっていそうだ。
11以上の素数どれかの倍数、という集合の要素数を包除原理で数える感じかな？

掲載されているHaskell版

少し直して：

import qualified Data.IntSet as IS

humbleNumbers2 :: [Int]
humbleNumbers2 = go $ IS.singleton 1
  where
    go sofar = x : go (IS.union pruned $ IS.fromList $ map (x *) [2, 3, 5, 7])
      where
        (x, pruned) = IS.deleteFindMin sofar

「humble数だと判明していて、まだ公開していない数の袋」sofar を用意する。
その中の最小値 x をとりだし、humble数として出力する。
xの2,3,5,7倍の値を、今後出力されるべきhumble数として袋にしまう。重複は袋の機能で消え去る。

上の humberNumbers1 が失敗した原因は merge のthunkが残ることだった。
これを IS.union でその場で完了させてしまうことで、問題を解決している。つまり大方針が同じ。

ghci> test1 12 humbleNumbers2
[出力省略]
(0.10 secs, 55,978,848 bytes)

性能も遜色ない。

「らしく」再訪

他のとまた違うスタイルのプログラムが、ALGOL Wあたりからの系譜にあるようだ。

色々とロジックが混ざって読みづらいが、核の部分は多分こんな感じ：

• humble数を入れる配列 H[] を用意する。
• H[0] = 1 とする。
• 既出のhumble数のk倍を次のhumble数として推薦するエージェントを、k=2,3,5,7の4名つくる。
  エージェントは、次に推薦する数の元になるhumble数の添え字を状態として持つ。初期値は全員0
• 以下を繰り返す
  • 4名の推薦する数の最小値を次のhumble数としてH[i++]に入れる。
  • 値が採用されたエージェントは添え字を1つ次に進める。不採用の人は変化しない

過去の値をとっておく配列を設置し、その添え字で値にアクセスする、という泥臭いことに目をつぶれば、このプログラムの振る舞いは、この記事の最初のHaskell解と完全に一致している。
つまり、「らしく」解は単なる車輪の再発明だった？

それより、このロジックを、制御の構造化がせいぜいの ALGOL で実現すること、 JavaScript などのgenerator 関数や stream 指向プログラミング（やその 先祖のobserver パターン）あるいはイテレータでもっとスマートに実現すること、そして Haskell のストリームとしての無限リストで10行足らずで簡潔に書くこと、の苦労の差を想像して、プログラミングの進化を感じる機会と受け止めておく。