この記事は ひとりアドベントカレンダーRosettaCodeで楽しむプログラミング Advent Calendar 2025の21日めの記事です。
定義
$u[1] = 1, u[2] = 2$ とする。
$u[n] \ (n > 2)$ は、$u[n-1] < u[n]$, $u[n] = u[i] + u[j] \ (i < j)$ となる表し方が一通りだけある最小のもの
タスク
$u[n]$ を求める関数を作れ。
考える
自作
既出の数ふたつで作れる数を全て作って持っておく。
またそれは、作り方が唯一か、二つ以上あるかも記録しておく。
これは、IntMap Bool でできるだろう。
次の要素を生成するには、この中から、作り方が唯一な最小値を取り出せばよい。
次の要素が決まると、既出の全ての値とこの値とを足した値を全て、作れる数として追加する。
このとき、二つめならそのように印を設定する。
既出の全ての値は、自身の出力結果リストを前から何個まで使ってよいかのカウンタだけ持って回ればできる。
import qualified Data.IntMap as IM
ulams1 :: [Int]
ulams1 = 1 : 2 : gen 2 im0
where
im0 = IM.singleton 3 True
gen :: Int -> IM.IntMap Bool -> [Int]
gen k im = u : gen (succ k) im1
where
-- imのTrueな最小値が次の要素
u = fst $ head $ filter snd $ IM.assocs im
-- uと既出の要素との和が新規の要素候補。これをimに反映させ、またu以下は消去する
im1 = foldl' step (snd $ IM.split u im) $ map (u +) $ take k ulams1
step im v =
case IM.lookup v im of
Nothing -> IM.insert v True im -- 新規ならTrueを割り当て
Just True -> IM.insert v False im -- 2度めならFalseを割り当て
_ -> im -- 3度め以降は何もしない
ghci> take 50 ulams1
[1,2,3,4,6,8,11,13,16,18,26,28,36,38,47,48,53,57,62,69,72,77,82,87,97,99,102,106,114
,126,131,138,145,148,155,175,177,180,182,189,197,206,209,219,221,236,238,241,243,253]
ghci> ulams1 !! 9
18
(0.01 secs, 87,376 bytes)
ghci> ulams1 !! 99
690
(0.01 secs, 2,084,904 bytes)
ghci> ulams1 !! 999
12294
(0.22 secs, 147,519,392 bytes)
ghci> ulams1 !! 9999
132788
(24.79 secs, 13,278,612,288 bytes)
Rosetta Code 版
遅延リストなのはいいけど、色々慣れてない感じ。
import Data.List
ulam
:: Integral i =>
Int -> i
ulam 1 = 1
ulam 2 = 2
ulam n
| n > 2 = ulams !! (n-1)
ulams
:: Integral n =>
[n]
ulams = 1:2:(nexts [2,1]) -- 1:2:あるなら上いらなくない？
nexts us = u: (nexts (u:us))
where
n = length us
[u] = head . filter isSingleton . group . sort $
[v | i <- [0 .. n-2], j <- [i+1 .. n-1] -- 配列じゃないんだindex使うな
, let s = us !! i
, let t = us !! j
, let v = s+t
, v > head us
]
isSingleton :: [a] -> Bool
isSingleton as
| length as == 1 = True -- length数えんな
| otherwise = False
勝手に赤ペン。
import Data.List
ulam2 :: Int -> Int
ulam2 n = ulams !! pred n
where
ulams = 1 : 2 : nexts [2,1]
nexts us = u : nexts (u:us)
where
w = head us
u = head $ head $ filter isSingleton $ group $ sort
[v | s:ts <- tails us, t <- ts, let v = s + t, w < v]
isSingleton [_] = True
isSingleton _ = False
ghci> ulam2 10
18
(0.00 secs, 113,216 bytes)
ghci> ulam2 100
690
(0.07 secs, 29,100,056 bytes)
ghci> ulam2 1000
12294
(69.70 secs, 28,325,519,472 bytes)
音を上げるのが一桁早いぞ。
OEIS 版
型は性能比較のため Int に統一した。
ulam 関数は局所定義に隠した。
a002858 n = a002858_list !! (n-1)
where
a002858_list = 1 : 2 : ulam 2 2 a002858_list
ulam :: Int -> Int -> [Int] -> [Int]
ulam n u us = u' : ulam (n + 1) u' us where
u' = f 0 (u+1) us'
f 2 z _ = f 0 (z + 1) us'
f e z (v:vs) | z - v <= v = if e == 1 then z else f 0 (z + 1) us'
| z - v `elem` us' = f (e + 1) z vs
| otherwise = f e z vs
us' = take n us
ghci> a002858 10
18
(0.00 secs, 96,904 bytes)
ghci> a002858 100
690
(0.02 secs, 5,028,328 bytes)
ghci> a002858 1000
12294
(5.73 secs, 438,580,624 bytes)
ulam が何をしているのか見ていく。
引数 us は変更されることなく引き渡され続けるだけで、その実体は結果のリスト全体である。これは名前参照で外せる。
引数 n で、これまで出力済みの値だけを us'に取り出して使っている。
引数 u は、直前に出力した値。次に作る値はこれより大きいものなので、探索の始点として使っている。
次に出力する値 u' は、内部関数 f で作っている。
引数 z が現在確認中の出力候補の値で、成功するとこの関数の返り値になる。
引数 v:vs は既出の値のリストで、これを順に消費して計算を進めている。
z - v とは、z を作る $u[i]$ が v と仮定したときの相方 $u[j]$ の値。
ここでは捜索範囲を $u[i] < u[j]$ と限定しているので、z - v <= v とは、この範囲を逸脱したことを意味する。
このとき、引数 e が、z を作れるペアの見つかった個数のカウンタで、
それが 1 なら答えを見つけることができたので返している。
そうでないとき、z を1増やして最初からやり直しになる。
z - v が us に含まれるとは、つまり相方 $u[j]$ が見つかったということなので、見つけたカウンタを1増やして vs を次に進める。
含まれないとき、単にvsを次に進める。
vs を使い切った場合が書かれていない。次の答えまでに、us を組み合わせて作れない数はなく、除外される理由は、作り方が複数回あることだけということになる。
見つけたカウンタが2になったとき、それ以上探すことを放棄して次の候補で仕切り直す。
これを別の行に分ける必要はなく、 e + 1 を行うときに同時に判断すればよい。
すると、引数 e は Int でなく Bool で済む。1つ発見済みのとき True とする。
以上の検討から、このように変わった。
ulam3 n = a002858_list !! pred n
where
a002858_list = 1 : 2 : ulam 2 2
ulam :: Int -> Int -> [Int]
ulam n u = w : ulam (succ n) w
where
us = take n a002858_list
w = f False (succ u) us
-- f _ z [] = f False (succ z) us
f e z (v:vs) | z - v <= v = if e then z else next
| elem (z - v) us = if e then next else f True z vs
| otherwise = f e z vs
where
next = f False (succ z) us
ghci> ulam3 10
18
(0.00 secs, 92,152 bytes)
ghci> ulam3 100
690
(0.02 secs, 4,814,688 bytes)
ghci> ulam3 1000
12294
(5.65 secs, 434,769,160 bytes)
さらに、us を n から作る昇順リストにしておく代わりに、
既出の要素を全て貯め込んだ IntSet にすることで elem を高速化する。
import qualified Data.IntSet as IS
ulam4 :: Int -> Int
ulam4 n = a002858_list !! pred n
where
a002858_list = 1 : 2 : ulam (IS.fromList [1,2]) 2
ulam :: IS.IntSet -> Int -> [Int]
ulam us u = w : ulam (IS.insert w us) w
where
w = f False (succ u) (IS.elems us)
f e z (v:vs) | z - v <= v = if e then z else next
| IS.member (z - v) us = if e then next else f True z vs
| otherwise = f e z vs
where
next = f False (succ z) (IS.elems us)
ghci> ulam4 10
18
(0.01 secs, 94,344 bytes)
ghci> ulam4 100
690
(0.02 secs, 5,324,936 bytes)
ghci> ulam4 1000
12294
(0.65 secs, 463,304,760 bytes)
ghci> ulam4 10000
Interrupted.
動く範囲では速いが、やはり爆発的に遅くなる。
Rosetta Code の状況
他言語の解では 100,000 まで結果を出しているようなものもある。
それらの内容と由来を辿ってみると、XPL0版が集大成になっているような雰囲気。
sieveだとか独自の用語が使われているが、基本方針は自分の ulams1 と同じで、IntMap 的なデータ構造で管理する代わりに巨大な配列に記録することで性能をたたき出しているだけのようだ。
トンデモ解
同じ事を Data.Array.ST で再現しても芸がないので、もう少しひねってみる。
整数一つに対して、
- 既出のウラム数かどうか
- 作り方がひとつ見つかったウラム数候補であるか
- 作り方がもう一つ見つかったウラム数失格候補であるか
の情報を持ち回す必要がある。
Integerのビット演算を用いて、無限サイズのビット要素配列を実装すれば、たった3つのIntegerでこれらを追跡できる。
ulams5 :: [Int]
ulams5 = 1 : 2 : go 3 im0 jm0 km0
where
im0 = bit 1 .|. bit 2
jm0 = bit 0
km0 = 0
go :: Int -- 確認するべき候補の値 ＝ jm, km のオフセット
-> Integer -- im 生成済みUlam数のビット配列 0-origin
-> Integer -- jm Ulam数候補のビット配列 u-origin
-> Integer -- km 失格候補のビット配列 u-origin
-> [Int] -- 答え
go u im jm km
| isUlam = u : go (succ u) im1 (jm1 .>>. 1) (km1 .>>. 1)
| otherwise = go (succ u) im (jm .>>. 1) (km .>>. 1)
where
isUlam = testBit jm 0 && not (testBit km 0)
im1 = setBit im u
jm1 = jm .|. im
km1 = km .|. (im .&. jm)
jm km は随時右シフトすることで、常に今回の候補がビット0に位置するようにしている上、発見したときにこれらを更新する計算（ulams1 の im1 を作る計算）のときに、im .<<. u とすることもなしにそのまま重ねるだけでよくしている。
ghci> ulams5 !! 9
18
(0.00 secs, 74,936 bytes)
ghci> ulams5 !! 99
690
(0.00 secs, 563,776 bytes)
ghci> ulams5 !! 999
12294
(0.02 secs, 29,414,904 bytes)
ghci> ulams5 !! 9999
132788
(0.40 secs, 2,662,570,360 bytes)
ghci> ulams5 !! 99999
1351223
(28.78 secs, 268,520,146,296 bytes)
$10^5$ についても答えが得られたので満足です。
遅延リストらしい解法その2
Haskellにおける、ダメなエラトステネスのふるいの実装：
primes = 2 : sieve [3, 5 ..]
where
sieve (p:xs) = p : sieve (filter (\x -> mod x p /= 0) xs)
に寄せるスタイルでも作れる。
ulams6 :: [Int]
ulams6 = 1 : 2 : pull [3]
where
pull (x:xs)
| x < 0 = pull xs
| otherwise = x : push x xs
push u stream = pull $ merge stream $ map (u +) $ takeWhile (u >) ulams6
merge xxs@(x:xs) yys@(y:ys) =
case compare (abs x) (abs y) of
EQ -> negate (abs x) : merge xs ys
LT -> x : merge xs yys
GT -> y : merge xxs ys
merge xs [] = xs
merge [] ys = ys
ストリームには次の値の候補が並ぶが、merge で同じ値が重複したとき負の値に変えられる。それは pull で除外される。
pull は正の値である有効な候補を見つけたとき、それを種にストリームを更新して次の値を探す push を起動する。
ghci> ulams6b !! 9
18
(0.01 secs, 83,296 bytes)
ghci> ulams6b !! 99
690
(0.01 secs, 5,839,960 bytes)
ghci> ulams6b !! 999
12294
(1.11 secs, 998,813,432 bytes)
ghci> ulams6b !! 9999
132788
(144.38 secs, 129,941,917,576 bytes)
性能的にはふるわないが仕方ない。