Jordan-Pólya numbers

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

定義

一つ以上の階乗（同じものを複数回使ってもよい）を掛け合わせて作られる数。

タスク

• 最初の50個のJP数
• 100 million未満の最大のJP数
• ボーナス
  • 800,1800,2800,3800番目のJP数と、階乗の構成、あと1050番目
  • JP数を階乗に分割するとき、複数のやり方があるときは、最も大きな階乗を使うものを選ぶこと。
  • ヒント：これらは2^53未満である

調べる

派生関係

Rosetta Code のコード例は、他の言語のものを移植した派生版が多い。
その関係を調べてみる。（明記してあるもののみ）

他の項目でもよくWrenが参考にされているが、聞いたことない言語だ。

Think Smalltalk in a Lua-sized package with a dash of Erlang and wrapped up in a familiar, modern syntax.

面白そう。

OEIS Haskell

OEIS A001013 : Jordan-Polya numbers

Haskell版がある。

import Data.Set (empty, fromList, deleteFindMin, union)
import qualified Data.Set as Set (null)

a001013 n = a001013_list !! (n-1)
a001013_list = 1 : h 0 empty [1] (drop 2 a000142_list) where
   h z s mcs xs'@(x:xs)
    | Set.null s || x < m = h z (union s (fromList $ map (* x) mcs)) mcs xs
    | m == z = h m s' mcs xs'
    | otherwise = m : h m (union s' (fromList (map (* m) $ init (m:mcs)))) (m:mcs) xs'
    where (m, s') = deleteFindMin s

-- A000142 : Factorial numbers
a000142 :: (Enum a, Num a, Integral t) => t -> a
a000142 n = product [1 .. fromIntegral n]
a000142_list = 1 : zipWith (*) [1..] a000142_list

ていうかまたオヌシか、Reinhard Zumkeller サン。
自分風に書き直して解釈する。

import qualified Data.Set as S

jpns0 :: [Integer]
jpns0 = 1 : h 0 S.empty [1] (scanl (*) 2 [3 ..])
  where
    h :: Integer       -- z : 直前に出力した値、重複の削除用
      -> S.Set Integer -- s : これまでに生成した数から作った、以降に生成するべきjpn
      -> [Integer]     -- mcs : 全ての生成した値を貯め込むリスト、増える一方
      -> [Integer]     -- f:fs : 2!, 3!, ... の順に消費されるリスト
      -> [Integer]     -- 結果
    h z s mcs ffs@(f:fs)
      | S.null s ||      -- 最初の起動のとき
        f < m            -- 次に生成する予定の数より、次の階乗の方が小さいとき
                  = h z s2 mcs fs -- 階乗fを一つ消費してs2にmcsとの積を投入してやり直し
      | m == z    = h z s1 mcs ffs -- ダブったmはスルーしてやり直し
      | otherwise = m : h m s3 (m:mcs) ffs -- mがひとつ出力される
      where
        (m, s1) = S.deleteFindMin s -- sの最小値mを取り出す 残りはs1
        s2 = S.union s  (S.fromList $ map (f *) mcs) -- mcsに次の階乗fを掛けたものを追加
        s3 = S.union s1 (S.fromList $ map (m *) $ init $ m:mcs) -- 末尾にある1を除いた
                        -- これまでの全ての出力にmを掛けたもの全てを、mを除いたs1に追加

値を一つ生成するたびに s3 で追加するものが増えていき、大変なことになっていそうだ。

test1 jps = do
  print $ take 50 jps
  print $ last $ takeWhile (100_000_000 >) jps
  mapM_ (\i -> print (i, jps !! pred i)) [800,1050,1800,2800,3800]

ghci> test1 jpns0
[1,2,4,6,8,12,16,24,32,36
,48,64,72,96,120,128,144,192,216,240
,256,288,384,432,480,512,576,720,768,864
,960,1024,1152,1296,1440,1536,1728,1920,2048,2304
,2592,2880,3072,3456,3840,4096,4320,4608,5040,5184]
99532800
(800,18345885696)
(1050,139345920000)
(1800,9784472371200)
(2800,439378587648000)
(3800,7213895789838336)
(5.42 secs, 6,833,487,904 bytes)

後半がだんだん遅くなっていくが、最後まで動いた。

OEIS Python

def aupto(lim, mx=None):
    if lim < 2: return [1]
    v, t = [1], 1
    if mx == None: mx = lim
    for k in range(2, mx+1):
        t *= k
        if t > lim: break
        v += [t*rest for rest in aupto(lim//t, t)]
    return sorted(set(v))
print(aupto(5760))

名前は A001013 upto lim mx (lim: 結果の上限 mx:使う階乗の上限) という感じか。
mxを渡さないのが外部API、使うのは内部で再帰呼び出しするとき、という使い分けのようだ。げろげろ。
v が今までに生成した結果を順不同で貯め込むリスト。
for k で回して、t *= k で階乗の値を順に作って、t が上限を超えるまで、再帰呼び出し aupto(lim//t, t) で作ったもう少し小さいものに t を掛けて v に追加。
全て終わったら set(v) で重複を消して sorted() で昇順にして返す。

と、計算の流れがわかったのでHaskell化する。

import qualified Data.IntSet as IS
import Data.List

jps1 :: Int -> [Int]
jps1 lim0 = IS.elems $ go lim0 lim0
  where
    one = IS.singleton 1
    go lim mx
      | lim < 2   = one
      | otherwise = IS.unions $ one :
          [ IS.map (t *) $ go (div lim t) t
          | t <- takeWhile (lim >=) $ scanl' (*) 2 [3 .. mx] ]

v は初めから IntSet で扱う。再帰呼び出しする内部関数は go として隠す。
それは IntSet を返し、最後にリストにする。

test2 :: (Int -> [Int]) -> IO ()
test2 jf = do
  print $ take 50 $ jf 5760
  print $ last $ jf 100_000_000
  let jps = jf $ 2^53
  mapM_ (\i -> print (i, jps !! pred i)) [800,1050,1800,2800,3800]

test20 :: (Int -> [Int]) -> IO ()
test20 jf = do
  print $ take 50 $ jf 5760
  print $ last $ jf 100_000_000

上限を与えないとリストにならず、ボーナス問題は $2^{53}$ を指定することになるので、かなり時間がかかるためオミットする。

ghci> test20 jpns1
[略]
(0.04 secs, 13,766,904 bytes)

Phix

次は Wren v2 の元になった Phix のを見てみる。
なんかゴチャゴチャやってるなぁ。

function factorials_le(atom limit) // ようは、値がlimit以下の階乗のリストを作って返す
    sequence res = {}
    while true do
        atom nf = factorial(length(res)+1)
        if nf>limit then exit end if
        res &= nf
    end while
    return res
end function

function jp(atom limit)                  // limit以下のJP数の昇順のリストを作る
    sequence res = factorials_le(limit)  // 最初はlimitまでの階乗を入れておく
    integer k=2
    while k<=length(res) do              // 結果とするresの要素を小さい順に見ていく。
// ループ中で追加されるが、それは res[k] より大きく、後ろに追加されるので順序が狂うことはない
        atom rk = res[k]
        for l=2 to length(res) do        // resにある既出の要素全て res[l] について
            atom kl = res[l]*rk          // res[l] * rk^1,2,3... をlimitまで追加
            if kl>limit then exit end if
            do
                integer p = binary_search(kl,res) // resが昇順を保つようにklを挿入
                if p<0 then
                    p = abs(p)
                    res = res[1..p-1] & kl & res[p..$]
                end if                            // つまり Data.List.insert
                kl *= rk
            until kl>limit
        end for
        k += 1
    end while
    return res
end function

だいたい読めたのでHaskellに翻訳する。
res は初めから Data.Set で管理する。

jpns2 :: (Ord a, Num a, Enum a) => a -> [a]
jpns2 limit = 1 : loop 1 s0
  where
-- 階乗をlimit以下で入れる
    s0 = S.fromAscList $ takeWhile (limit >=) $ scanl' (*) 1 [2 .. ]
-- arg1が前回の結果
-- それの次に大きい値が今回の成果、まずそれを探す
    loop jp0 s =
      case S.lookupGT jp0 s of
        Nothing  -> [] -- もうない、つまりlimitまでにはないとき、終了。
        Just jp1 -> jp1 : next jp1 s -- 見つかったなら、それ jp1 (was rk)を出力し、
-- さらに、sの1以外の(自身を含む）全ての要素 res[l] に、rk^1,rk^2,… をlimitまで全て追加して再開
    next jp1 s = loop jp1 $ S.unions $ s :
      [S.fromAscList $ takeWhile (limit >) $ tail $ iterate (jp1 *) rl  | rl <- tail $ S.elems s]

ghci> test2 (map fromIntegral . jpns2 . fromIntegral)
[略]
(21.86 secs, 13,538,045,120 bytes)

ボーナス込みで動いた。

Java も派生が多い。読む気も起こらないからわからないけど、多分これと似たようなものでは？

オリジナル解法

色々見ているうちにひらめいた。
JP数の値でなく、$\prod (k!)^{e_k}$ の $e_k$ の列の方で管理するべき。

基本の考え方

まず、説明のために上限ありで考える。

上限 $l$ を単一で超えないギリギリの階乗 $k ! \leq l < (k+1)!$ を限界として、
それぞれのJP数 $(2!)^{e_2} \cdot \ldots \cdot (k!)^{e_k}$ を作る乗数の組み合わせを長さ $k-1$ のベクトル $\langle e_2, e_3, \dots, e_k \rangle$ で表す。
$e_k$ は0と1しかありえない。その他の変数の上限は他の値との兼ね合いで変化するので、工夫して列挙する。

あるJP数を出力したとき、この値から派生するJP数のベクトル表記は
$\langle e_2+1,e_3,\dots,e_k \rangle$ から $\langle e_2,\dots,e_k+1 \rangle$までの$k-1$とおりある。
そこで、今後出力するべきJP数を優先度付きキューに格納する。優先度にJP数の値、本体にベクトル表記を持たせておく。
キュー先頭の値を出力したら、ベクトル表記から$k-1$個の次の値を作り、キューに投入する。
ここで、ベクトルが同一なもの（$\langle 1,1 \rangle \to \langle 2,1 2,2 \rangle$ と $\langle 1,1 \rangle \to \langle 1,2 \rangle \to \langle 2,2 \rangle$ など）、ベクトルが異なっていても値が同一なものが存在しうるので、出力が重複する。
これを回避するため、値が同一のエントリは、ベクトルが最も長いものに統合する。このため、優先度付きキューとしては本職の Data.Heap ではなく Data.Map を流用する。

拡張

次にこれを、上限なしの構造に拡張する。

長さ固定のベクトルの代わりに、可変長のリスト $[e_2, e_3, \dots]$ で表現する。
ここから派生するJP数は、上限がないため、ベクトル表記のものに加えて
$[e_2,\dots,e_k,1], [e_2,\dots,e_k,0,1], [e_2,\dots,e_k,0,0,1], \cdots$
と無限に出てきてしまって困る。
そこでそうする代わりに、単一の階乗ひとつからなるJP数を特別扱いする。
通常のJP数からは、これまで通り $k-1$ 通りのJP数のみを派生させる。
特別なJP数 $[0, \cdots, 0, 1]$ からは、通常のものに加えて、長さが1増えた次の特別なJP数も派生させる。
新たな階乗が必要になったときに、またその次の階乗が必要になったときに発動する番兵を再度仕込む感じ。

この方式なら、ボーナスタスクの「decomposition も示せ」にも容易に対応できる。

実装

import qualified Data.Map as M

-- Jordan-Polya numbers and decompositions
jpns3 :: [(Integer,[Int])]
jpns3 = ((1,[]) :) $ go (M.singleton 2 [1])
  where
    facts = scanl1 (*) [2 ..]        -- facts !! (n-2) == n !
    ones = iterate (0 :) (1 : zeros) -- [1,0,..],[0,1,0..],[0,0,1,0,..],...
    zeros = 0 : zeros

    go q = res : go (qc $ M.unionWith unifun q1 q2) -- q キューとして使うMap
      where
        (res@(jpn, es), q1) = M.deleteFindMin q     -- キュー先頭を取り出す
        es1 = 0 : es                                -- 番兵用、1 を一つ後ろにずらす
        jpn1 = fst $ last $ zip facts es1           -- その 1 がある位置のfactsの要素が値
        qc = if sum es == 1 then M.insert jpn1 es1 else id  -- 番兵を追加する操作
        q2 = M.fromList [(jpn * f, zipWith (+) one es) | (f, one, _) <- zip3 facts ones es] -- キューに追加する、jpnからの派生要素

    unifun as bs = snd $ max        -- M.unionWith で重複時に残すものを選ぶ
      ((length as, reverse as), as) -- 長い方、長さが一致したら後ろがより重い方
      ((length bs, reverse bs), bs)

test3 jpdcs = do
  print $ take 50 $ map fst jpdcs
  print $ last $ takeWhile (100_000_000 >) $ map fst jpdcs
  mapM_ (\i -> print (i, jpdcs !! pred i)) [800,1050,1800,2800,3800]

ghci> test3 jpns3
[1,2,4,6,8,12,16,24,32,36
,48,64,72,96,120,128,144,192,216,240
,256,288,384,432,480,512,576,720,768,864
,960,1024,1152,1296,1440,1536,1728,1920,2048,2304
,2592,2880,3072,3456,3840,4096,4320,4608,5040,5184]
99532800
(800,(18345885696,[2,0,7]))
(1050,(139345920000,[1,0,0,3,0,0,1]))
(1800,(9784472371200,[15,0,2,0,2]))
(2800,(439378587648000,[0,0,0,0,0,1,0,0,0,0,0,0,1]))
(3800,(7213895789838336,[16,0,8]))
(0.16 secs, 63,289,752 bytes)

圧倒的じゃないか。
お化粧してフィードバックかけましょうかね。

import Text.Printf
import Data.List.Split
import Control.Monad

main = do
  putStrLn "First 50 Jordan-Pólya numbers:"
  putStrLn $ unlines $ chunksOf 50 $ concatMap (printf "%5d") $ take 50 $ map fst jpns3
  printf "The largest Jordan-Pólya number before 100 millon: %d\n" $ last $ takeWhile (100_000_000 >) $ map fst jpns3
  forM_ [800,1050,1800,2800,3800] (\i -> do
    let (n,dc) = jpns3 !! pred i
    printf "\nThe %dth Jordan-Pólya number is : %d\n" i n
    putStrLn $ intercalate " * " $ map (uncurry (printf "(%d!)^%d")) $ filter ((0 <) . snd) $ zip [2 :: Int ..] dc
    )

ghci> main
First 50 Jordan-Polya numbers:
    1    2    4    6    8   12   16   24   32   36
   48   64   72   96  120  128  144  192  216  240
  256  288  384  432  480  512  576  720  768  864
  960 1024 1152 1296 1440 1536 1728 1920 2048 2304
 2592 2880 3072 3456 3840 4096 4320 4608 5040 5184

The largest Jordan-Polya number before 100 millon: 99532800

The 800th Jordan-Polya number is : 18345885696
(2!)^2 * (4!)^7

The 1050th Jordan-Polya number is : 139345920000
(2!)^1 * (5!)^3 * (8!)^1

The 1800th Jordan-Polya number is : 9784472371200
(2!)^15 * (4!)^2 * (6!)^2

The 2800th Jordan-Polya number is : 439378587648000
(7!)^1 * (14!)^1

The 3800th Jordan-Polya number is : 7213895789838336
(2!)^16 * (4!)^8
(0.36 secs, 64,316,568 bytes)