今回は「数」を題材に簡単なプログラムを Haskell で作ってみましょう。
点を多角形の形に並べたとき、その総数を多角数 (polygonal number) といいます。三角形に配置したものを三角数 (triangular number)、四角形に配置したものを四角数 (square number)、五角形に配置したものを五角数 (pentagonal number) といいます。
三角数、四角数、五角数を下図に示します。
１ ３ ６ 10 15
● ● ● ● ●
●● ●● ●● ●●
●●● ●●● ●●●
●●●● ●●●●
●●●●●
図 : 三角数
|
１ ４ ９ 16 25
● ●● ●●● ●●●● ●●●●●
●● ●●● ●●●● ●●●●●
●●● ●●●● ●●●●●
●●●● ●●●●●
●●●●●
図 : 四角数
|
１ ５ 12 22
● ● ● ●
● ● ● ● ● ●
● ● ● ●
● ● ● ● ● ●
● ● ● ● ●
● ● ●
● ● ● ●
● ● ● ●
● ● ● ●
図 : 五角数
多角数の点の増分を表に示すと、次のようになります。
n 三角数 四角数 五角数
---+-----------------------------------------------------------
1 | 1 1 1
2 | 3 = 1+2 4 = 1+3 5 = 1+4
3 | 6 = 1+2+3 9 = 1+3+5 12 = 1+4+7
4 | 10 = 1+2+3+4 16 = 1+3+5+7 22 = 1+4+7+10
5 | 15 = 1+2+3+4+5 25 = 1+3+5+7+9 35 = 1+4+7+10+13
6 | 21 = 1+2+3+4+5+6 36 = 1+3+5+7+9+11 51 = 1+4+7+10+13+16
・・・・・・ ・・・・・・・ ・・・・・
n | n(n + 1) / 2 n^2 n(3n - 1) / 2
表を見ればお分かりのように、三角数は公差 1、四角数は公差 2、五角数は公差 3、p 角数は公差 p - 2 の等差数列の和になります。初項を a, 公差を d とすると、等差数列の和 \(S_n\) は次式で求めることができます。
a = 1, d = p - 2 を代入して計算すると、多角数 \(P_{p,n}\) は次式で求めることができます。
この式を Haskell でプログラムすると、次のようになります。
リスト : 多角数 polygonalNum :: Integer -> Integer -> Integer polygonalNum p n = ((p - 2) * n * n - (p - 4) * n) `div` 2
それでは実行してみましょう。
ghci> map (polygonalNum 3) [1..20] [1,3,6,10,15,21,28,36,45,55,66,78,91,105,120,136,153,171,190,210] ghci> map (polygonalNum 4) [1..20] [1,4,9,16,25,36,49,64,81,100,121,144,169,196,225,256,289,324,361,400] ghci> map (polygonalNum 5) [1..20] [1,5,12,22,35,51,70,92,117,145,176,210,247,287,330,376,425,477,532,590] ghci> map (polygonalNum 6) [1..20] [1,6,15,28,45,66,91,120,153,190,231,276,325,378,435,496,561,630,703,780] ghci> map (polygonalNum 7) [1..20] [1,7,18,34,55,81,112,148,189,235,286,342,403,469,540,616,697,783,874,970] ghci> map (polygonalNum 8) [1..20] [1,8,21,40,65,96,133,176,225,280,341,408,481,560,645,736,833,936,1045,1160]
p 角数の初項は 1 で、第 2 項は p になります。多角数 - Wikipedia によると、多角数にはいろいろな面白い性質があるようです。
整数 m が p 角数か判定するプログラムも簡単です。次式で n が整数値になれば、m は p 角数の第 n 項であることがわかります。
リスト : 多角数の判定
isPolygonalNum :: Integer -> Integer -> Bool
isPolygonalNum p m =
let a = p - 2
b = p - 4
c = b * b + 8 * a * m
d = floor (sqrt (fromIntegral c))
in d * d == c && (b + d) `mod` (2 * a) == 0
p - 2 を変数 a に、p - 4 を変数 b に、ルートの中を計算して変数 c にセットします。sqrt で c の平方根を求め、floor で整数値に変換して d にセットします。あとは、d * d が c に等しく、b + d が 2 * a で割り切れることを確認するだけです。
それでは実行してみましょう。
ghci> isPolygonalNum 3 55 True ghci> isPolygonalNum 3 54 False ghci> isPolygonalNum 4 64 True ghci> isPolygonalNum 4 65 False ghci> isPolygonalNum 5 117 True ghci> isPolygonalNum 5 118 False
もちろん、n の値を求めることもできます。興味のある方はプログラムを修正してみてください。
ピタゴラスの定理は、平面幾何学において直角三角形の辺を \(a, b, c \ (a + b \gt c)\) とすると、次式が成り立つという皆さんお馴染みの有名な定理です。
「ピタゴラス数」または「ピタゴラスの三つ組数 (pythagoras triple)」は、上式を満たす自然数の組 (a, b, c) のことで、とくに a, b, c が互いに素であるとき、(a, b, c) を「原始ピタゴラス数」といいます。たとえば、(3, 4, 5), (5, 12, 13), (8, 15, 17) などがあります。
ピタゴラスの定理 - Wikipedia によると、原始ピタゴラス数は次の方法で簡単に見つけることができるそうです。
本ページでは証明を割愛しますが、ピタゴラス数 - 桃の天然水 の説明がわかりやくて参考になると思います。inamori さんに感謝いたします。
それでは、実際にプログラムを作ってみましょう。三辺の合計値 (a + b + c) が k 以下の原始ピタゴラス数をすべて求めることにします。プログラムは次のようになります。
リスト : 三辺の合計値が k 以下の原始ピタゴラス数を求める
primitivePythagoras :: Integer -> [(Integer, Integer, Integer)]
primitivePythagoras k =
foldr (\m xs -> iter 1 m ++ xs) [] [2 .. mh]
where mh = ceiling (sqrt (fromIntegral k * 0.5))
iter n m
| n >= m || 2 * m * (m + n) > k = []
| odd (m - n) && gcd m n == 1 =
let a = m * m - n * n
b = 2 * m * n
c = m * m + n * n
in (if a < b then (a, b, c) else (b, a, c)) : iter (n + 1) m
| otherwise = iter (n + 1) m
三辺の合計値は次の式で求めることができます。
上式から三辺の合計値は必ず偶数になることがわかります。
変数 m の値を \(\sqrt{\frac{k}{2}}\) とすると、n = 1 のときに三辺の合計値は \(2m^2 + 2m = k + 2 \times \sqrt{\frac{k}{2}}\) となり、k よりも大きくなるので、m の上限値を \(\sqrt{\frac{k}{2}}\) とします。この値を変数 mh にセットします。ceiling で小数点を切り上げていることに注意してください。
foldr で m の値を 2 から mh まで増やしていき、ラムダ式の中の局所関数 iter で原始ピタゴラス数を生成します。n が m 以上になるか、三辺の合計値が k を超えたならば空リストを返します。そうでなければ、m - n が奇数で、m と n が互いに素であることを確認します。そして、辺 a, b, c を計算して iter の返り値のリストに (a, b, c) を追加するだけです。
それでは実行してみましょう。
ghci> primitivePythagoras 100 [(3,4,5),(5,12,13),(8,15,17),(7,24,25),(20,21,29),(9,40,41),(12,35,37)] ghci> primitivePythagoras 500 [(3,4,5),(5,12,13),(8,15,17),(7,24,25),(20,21,29),(9,40,41),(12,35,37),(11,60,61), (28,45,53),(33,56,65),(13,84,85),(16,63,65),(48,55,73),(39,80,89),(15,112,113), (36,77,85),(65,72,97),(17,144,145),(20,99,101),(60,91,109),(51,140,149),(19,180,181), (44,117,125),(88,105,137),(85,132,157),(57,176,185),(21,220,221),(24,143,145), (119,120,169),(95,168,193),(52,165,173),(104,153,185),(133,156,205),(28,195,197), (84,187,205)]
三辺の合計値が k となる全てのピタゴラス数を求めることも簡単です。次のリストを見てください。
リスト : 合計値が k となるピタゴラス数をすべて求める
pythagoras :: Integer -> [(Integer, Integer, Integer)]
pythagoras k =
foldr check [] $ primitivePythagoras k
where check (a, b, c) xs =
let p = k `mod` (a + b + c)
q = k `div` (a + b + c)
in if p /= 0
then xs
else (q * a, q * b, q * c) : xs
primitivePythagoras k で原始ピタゴラス数を生成し、k が三辺の合計値 (a + b + c) で割り切れることを確認するだけです。とても簡単ですね。
それでは実行してみましょう。
ghci> pythagoras 12 [(3,4,5)] ghci> pythagoras 120 [(30,40,50),(20,48,52),(24,45,51)] ghci> pythagoras 240 [(60,80,100),(40,96,104),(48,90,102),(15,112,113)]
もうひとつ、クールな方法を紹介します。この方法は nineties さんのブログ「ブートストラッピングでコンパイラを作る日記」の Project Euler (Problem 3～10) を参考にさせていただきました。nineties さんに感謝いたします。
式 \(a^2 + b^2 = c^2\) と \(a + b + c = k\) を使って変数 c を削除して因数分解すると、次のようになります。
上式より k - a と k - b は \(\frac{k^2}{2}\) の約数であることがわかります。つまり、\(\frac{k^2}{2}\) の約数を求め、条件 \(a \lt b \lt k\) を満たすものを探せばいいわけです。プログラムは次のようになります。
リスト : 合計値が k となるピタゴラス数をすべて求める (2)
pythagoras' :: Integer -> [(Integer, Integer, Integer)]
pythagoras' x = iter xs (reverse xs)
where xs = divisor (x * x `div` 2)
iter (y:ys) (z:zs) =
if z < y
then []
else if z >= x
then iter ys zs
else (a, b, c) : iter ys zs
where a = x - z
b = x - y
c = x - (a + b)
関数 divisor は約数をリストに昇順に格納して返します。拙作のページ Yet Another Haskell Problems 問題 54 で作成したプログラムと同じです。簡単な実行例を示します。
ghci> divisor 123456789 [1,3,9,3607,3803,10821,11409,32463,34227,13717421,41152263,123456789] ghci> divisor 100000 [1,2,4,5,8,10,16,20,25,32,40,50,80,100,125,160,200,250,400,500,625,800,1000,1250, 2000,2500,3125,4000,5000,6250,10000,12500,20000,25000,50000,100000] ghci> divisor 1000000 [1,2,4,5,8,10,16,20,25,32,40,50,64,80,100,125,160,200,250,320,400,500,625,800,1000, 1250,1600,2000,2500,3125,4000,5000,6250,8000,10000,12500,15625,20000,25000,31250, 40000,50000,62500,100000,125000,200000,250000,500000,1000000]
約数の一つを y とすると、もう一つの約数は z になります。z が y よりも大きい場合はすべての組み合わせを調べたので空リストを返します。z が x よりも大きい場合は条件を満たさないので次の組み合わせを調べます。そうでなければ条件を満たしているので、a, b, c を計算してリストに追加します。
それでは実行してみましょう。
ghci> :set +s ghci> pythagoras 1000000 [(200000,375000,425000),(218750,360000,421250)] (0.55 secs, 261,519,288 bytes) ghci> pythagoras' 1000000 [(200000,375000,425000),(218750,360000,421250)] (0.00 secs, 250,296 bytes) 実行環境 : GHCi, version 9.6.6, Ubunts 22.04 (WSL2), Intel Core i5-6200U 2.30GHz
k が大きな値の場合、約数を高速に求めることができれば pythagoras' のほうが速くなるようです。興味のある方はいろいろ試してみてください。
実数 \(a\) の平方根 \(\sqrt a\) の値を求める場合、方程式 \(x^2 - a = 0\) を Newton (ニュートン) 法で解くことが多いと思います。方程式を f(x), その導関数を f'(x) とすると、ニュートン法は次の漸化式の値が収束するまで繰り返す方法です。
平方根を求める場合、導関数は \(f'(x) = 2x\) になるので、漸化式は次のようになります。
参考文献『Ｃ言語による最新アルゴリズム事典』によると、\(\sqrt a\) より大きめの初期値から出発し、置き換え x <- (x + a / x) / 2 を減少が止まるまで繰り返すことで \(\sqrt a\) の正確な値を求めることができるそうです。
Haskell でプログラムすると、次のようになります。
リスト : 平方根を求める
sqrt' :: Double -> Double
sqrt' x
| x < 0 = error "sqrt': domain error"
| x == 0 = 0
| otherwise = iter (if x > 1 then init 1 x else 1)
where init s x
| s >= x = s
| otherwise = init (s * 2) (x / 2)
iter p = let q = (p + x / p) / 2
in if q >= p then p else iter q
局所関数 init は \(\sqrt x\) よりも大きめの初期値を求めます。たとえば、\(\sqrt {123456}\) を求める場合、初期値の計算は次のようになります。
s x ------------------- 1.0 123456.0 2.0 61728.0 4.0 30864.0 8.0 15432.0 16.0 7716.0 32.0 3858.0 64.0 1929.0 128.0 964.5 256.0 482.25 512.0 241.125 √123456 = 351.363060095964
s を 2 倍、x を 1 / 2 していき、s >= x となったときの s が初期値 (512) となります。\(4, 16, 64, 256, \cdots, 2^{2n} \) の平方根はこれだけで求めることができます。
あとは漸化式を計算して変数 q にセットし、q がひとつ前の値 p 以上になったら p を返すだけです。\(\sqrt {123456}\) を求めたときの p と q の値を示します。
p q -------------------------------------- 512.0 376.5625 376.5625 352.20622925311204 352.20622925311204 351.3640693544162 351.3640693544162 351.3630600974135 351.3630600974135 351.363060095964 351.363060095964 351.363060095964 √123456 = 351.363060095964
6 回の繰り返しで \(\sqrt {123456}\) を求めることができます。
整数 n の平方根の整数部分を求めることも簡単です。次のリストを見てください。
リスト : 整数 n の平方根を求める
iSqrt :: Integer -> Integer
iSqrt n
| n < 0 = error "iSqrt: domain error"
| n == 0 = 0
| otherwise = iter (init 1 n)
where init s x
| s >= x = s
| otherwise = init (s * 2) (x `div` 2)
iter p = let q = (n `div` p + p) `div` 2
in if q >= p then p else iter q
除算を / から `div` に変えただけです。\(\sqrt {123456}\) と \(\sqrt {123456789}\) のときの p, q の値の変化を示します。
p q (√123456 = 351) ------- 512 376 376 352 352 351 351 351 |
p q (√123456789 = 11111) ------------ 16384 11959 11959 11141 11141 11111 11111 11111 |
もうひとつ、面白い方法を紹介しましょう。次の公式を使って平方根の整数部分を求めます。
式 (1) は、奇数 1 から 2n - 1 の総和は \(n^2\) になることを表しています。式 (2) のように、整数 m の値が \(n^2\) より大きくて \((n + 1)^2\) より小さいのであれば、m の平方根の整数部分は n であることがわかります。これは m から奇数 \(1, 3, 5, \cdots, (2n - 1), (2n + 1)\) を順番に引き算していき、引き算できなくなった時点の (2n + 1) / 2 = n が m の平方根になります。参考文献『平方根計算法 (PDF)』によると、この方法を「めのこ平方」と呼ぶそうです。
プログラムは次のようになります。
リスト : めのこ平方 iSqrt'' :: Integer -> Integer iSqrt'' n m | n < m = m `div` 2 | otherwise = iSqrt'' (n - m) (m + 2)
アルゴリズムをそのままプログラムしただけなので、とくに難しいところは無いと思います。
それでは実行してみましょう。
ghci> iSqrt'' 4 1 2 ghci> iSqrt'' 16 1 4 ghci> iSqrt'' 64 1 8 ghci> iSqrt'' 80 1 8 ghci> iSqrt'' 81 1 9 ghci> iSqrt'' 82 1 9 ghci> iSqrt'' 100 1 10
この方法はとても簡単ですが、数が大きくなると時間がかかるようになります。そこで、整数を 2 桁ずつ分けて計算していくことにします。次の図を見てください。
整数 6789 を 67 と 89 に分ける 1 + 3 + ... + 15 = 82 < 67 両辺を 100 倍すると 802 < 6700 < 6789 802 = 1 + 3 + ... + 159 (= 2 * 80 - 1) 161 + 163 < (6789 - 6400 = 389) < 161 + 163 + 165
整数 6789 を 67 と 89 に分けます。最初に 67 の平方根を求めます。この場合は 8 になり、82 < 67 を満たします。次に、この式を 100 倍します。すると、802 < 6700 になり、6700 に 89 を足した 6789 も 802 より大きくなります。802 は 1 から 159 までの奇数の総和であることはすぐにわかるので、6789 - 6400 = 389 から奇数 161, 163, ... を順番に引き算していけば 6789 の平方根を求めることができます。
プログラムは次のようになります。
リスト : めのこ平方 (改良版)
iSqrt' :: Integer -> Integer
iSqrt' n
| n < 100 = iSqrt'' n 1
| otherwise = let m = 10 * iSqrt' (n `div` 100)
in iSqrt'' (n - m * m) (2 * m + 1)
iSqrt' は n の平方根の整数部分を求めます。n が 100 未満の場合は iSqrt'' で平方根を求めます。これが再帰呼び出しの停止条件になります。n が 100 以上の場合は、n の下位 2 桁を取り除いた値 (n `div` 100) の平方根を iSqrt' で求め、その値を 10 倍して変数 m にセットします。そして、iSqrt'' で n - m * m から奇数 2 * m + 1, 2 * m + 3 ... を順番に引き算していって n の平方根を求めます。
それでは実行してみましょう。
ghci> iSqrt' 6789 82 ghci> iSqrt' 123456789 11111 ghci> iSqrt' 1234567890 35136 ghci> :set +s ghci> iSqrt $ 2 ^ 1000 3273390607896141870013189696827599152216642046043064789483291368096133796404674 554883270092325904157150886684127560071009217256545885393053328527589376 (0.05 secs, 462,816 bytes) ghci> iSqrt' $ 2 ^ 1000 3273390607896141870013189696827599152216642046043064789483291368096133796404674 554883270092325904157150886684127560071009217256545885393053328527589376 (0.01 secs, 587,880 bytes) 実行環境 : GHCi, version 9.6.6, Ubunts 22.04 (WSL2), Intel Core i5-6200U 2.30GHz
iSqrt' は数が大きくなると遅くなりますが、思っていたよりも実行速度が速くて驚きました。実装がとても簡単なので、数が大きくなければ実用的に使えるかもしれません。興味のある方はいろいろ試してみてください。
N 個の要素の順列は N! 通りあります。この順列に 0 から N! - 1 までの番号をつけることを考えます。拙作のページ「お気楽Ｃ言語プログラミング超入門」 Yet Another Clang Problems (6) では、配列を使った方法を紹介しましたが、リストを使ってプログラムすることもできます。
たとえば、0 から 8 までの 9 個の整数の順列で、番号の振り方を考えてみましょう。最初が 0 で始まるパターンは 8! = 40320 通りありますね。このパターンには 0 - 40319 までの番号を割り当てます。そして、1 で始まるパターンには 40320 - 80639 までの番号を割り当てます。残りのパターンも同じです。
次に 2 番目の数字を考えましょう。01 で始まるパターンは 7! = 5040 通りあります。 したがって、このパターンには 0 - 5039 までの番号を割り当てます。10 で始まるパターンには 40320 - 45359 までの番号を、12 で始まるパターンには 45360 - 50399 までの番号を割り当てます。あとはこれを 9 番目までの数字まで続ければ、すべてパターンに番号を割り当てることができます。
では実際に 867254301 というパターンで試してみましょう。次の図を見てください。
８ 8 * 8! ６ [0 1 2 3 4 5 6 7] : 8*8! + 6*7! ７ [0 1 2 3 4 5 7] : 8*8! + 6*7! + 6*6! ２ [0 1 2 3 4 5] : 8*8! + 6*7! + 6*6! + 2*5! ５ [0 1 3 4 5] : 8*8! + 6*7! + 6*6! + 2*5! + 4*4! ４ [0 1 3 4] : 8*8! + 6*7! + 6*6! + 2*5! + 4*4! + 3*3! ３ [0 1 3] : 8*8! + 6*7! + 6*6! + 2*5! + 4*4! + 3*3! + 2*2! ０ [0 1] : 8*8! + 6*7! + 6*6! + 2*5! + 4*4! + 3*3! + 2*2! + 0*1! １ [1] : 番号：357478
注意すべき点は、数字をそのまま掛け算してはいけないところです。たとえば、7 に注目してください。このとき、残されている数字は 0, 1, 2, 3, 4, 5, 7 がありますね。番号は順番に振っていくので、867 は 86 で始まるパターンの 6*6! 番目から始まるのです。つまり、残っている数字の中で何番目に位置しているのかを求める必要があります。
それでは、Haskell でプログラムを作りましょう。次のリストを見てください。
リスト : 順列を番号に変換する
import Data.List
import Data.Maybe
fromPerm :: Eq a => [a] -> [a] -> Integer
fromPerm xs ys = iter xs ys 0
where
iter [_] _ n = n
iter zs@(x:xs) ys n = iter xs (delete x ys) (n + m * (fromIntegral i))
where m = fact $ fromIntegral (length zs - 1)
i = fromJust $ findIndex (== x) ys
関数 fromPerm の第 1 引数が番号に変換する順列、第 2 引数が要素を昇順に並べたリストです。実際の処理は局所関数 iter で行います。iter の第 3 引数は累積変数で、ここに順列の番号が求まります。
第 1 引数の要素が一つになったら順列を番号に変換できたので、第 3 引数の n を返します。そうでなければ、第 1 引数 zs の先頭要素 x に番号を割り当てます。最初に、zs の長さから 1 を引いた値の階乗を求めて、変数 m にセットします。次に、findIndex で ys にある x の位置を求めて、変数 i にセットします。findIndex の返り値は Maybe 型です。この場合、x は必ず見つかるので formJust で Just から値を取り出しています。あとは n に i * m を加算するだけです。iter を再帰呼び出しするときは、delete で ys から x を削除することをお忘れなく。
それでは実行してみましょう。
ghci> fromPerm [1,2,3,4] [1,2,3,4] 0 ghci> fromPerm [2,1,3,4] [1,2,3,4] 6 ghci> fromPerm [3,2,1,4] [1,2,3,4] 14 ghci> fromPerm [4,3,2,1] [1,2,3,4] 23 ghci> map (\xs -> fromPerm xs [1..4]) $ permutation 4 [1,2,3,4] [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23] ghci> fromPerm [0,1,2,3,4,5,6,7,8] [0..8] 0 ghci> fromPerm [4,5,6,7,8,0,1,2,3] [0..8] 184896 ghci> fromPerm [8,7,6,5,4,3,2,1,0] [0..8] 362879
permutation は拙作のページ「順列と組み合わせ」で作成した順列を生成する関数です。fromPerm は正常に動作していますね。
番号を順列に変換することも簡単です。次のリストを見てください。
リスト : 番号を順列に変換する (番号は 0 から開始)
toPerm :: Eq a => Integer -> [a] -> [a]
toPerm n xs = iter n xs []
where
iter _ [] a = reverse a
iter n xs a =
iter (n - m * p) (delete x xs) (x : a)
where m = fact (fromIntegral(length xs) - 1)
p = n `div` m
x = xs !! (fromIntegral p)
toPerm の第 1 引数が順列を表す番号、第 2 引数が要素を昇順に並べたリストです。実際の処理は局所関数 iter で行います。iter の第 3 引数は累積変数で、ここに順列が求まります。
iter の第 2 引数が空リストになったら番号を順列に変換できたので、第 3 引数 a を reverse で反転して返します。そうでなければ、xs の長さから 1 を引いた値の階乗を求めて変数 m にセットします。n `div` m で xs にある要素の位置がわかるので、それを演算子 !! で取り出して、変数 x にセットします。あとは、iter を再帰呼び出しするとき、n から m * p を引き算し、delete で x を xs から取り除き、x を累積変数 a に追加します。これで番号に対応した順列を求めることができます。
それでは実行してみましょう。
ghci> toPerm 0 [1,2,3,4] [1,2,3,4] ghci> toPerm 12 [1,2,3,4] [3,1,2,4] ghci> toPerm 23 [1,2,3,4] [4,3,2,1] ghci> map (\x -> toPerm x [1,2,3,4]) [0..23] [[1,2,3,4],[1,2,4,3],[1,3,2,4],[1,3,4,2],[1,4,2,3],[1,4,3,2],[2,1,3,4],[2,1,4,3], [2,3,1,4],[2,3,4,1],[2,4,1,3],[2,4,3,1],[3,1,2,4],[3,1,4,2],[3,2,1,4],[3,2,4,1], [3,4,1,2],[3,4,2,1],[4,1,2,3],[4,1,3,2],[4,2,1,3],[4,2,3,1],[4,3,1,2],[4,3,2,1]] ghci> toPerm 0 [0..8] [0,1,2,3,4,5,6,7,8] ghci> toPerm 362879 [0..8] [8,7,6,5,4,3,2,1,0]
正常に動作していますね。
--
-- number.hs : 数で遊ぼう
--
-- Copyright (C) 2022 Makoto Hiroi
--
import Data.List
import Data.Maybe
-- 多角数
polygonalNum :: Integer -> Integer -> Integer
polygonalNum p n = ((p - 2) * n * n - (p - 4) * n) `div` 2
isPolygonalNum :: Integer -> Integer -> Bool
isPolygonalNum p m =
let a = p - 2
b = p - 4
c = b * b + 8 * a * m
d = floor (sqrt (fromIntegral c))
in d * d == c && (b + d) `mod` (2 * a) == 0
-- 素因数分解
factorization :: Integer -> [(Integer, Integer)]
factorization n =
if c > 0 then iter 3 m [(2, c)] else iter 3 n []
where
(c, m) = factor n 2 0
factor n m c =
if n `mod` m /= 0 then (c, n)
else factor (n `div` m) m (c + 1)
iter i n a
| n == 1 = reverse a
| n < i * i = reverse ((n, 1) : a)
| otherwise =
let (c, m) = factor n i 0
in if c == 0 then iter (i + 2) n a
else iter (i + 2) m ((i, c):a)
sigma :: Integer -> Integer -> [Integer]
sigma p n = map (p^) [0 .. n]
-- 約数
divisor :: Integer -> [Integer]
divisor n =
let ((p1, x1):xs) = factorization n
product xs ys = [x * y | x <- xs, y <- ys]
in sort $ foldl (\a (p, x) -> product (sigma p x) a) (sigma p1 x1) xs
-- ピタゴラス数
primitivePythagoras :: Integer -> [(Integer, Integer, Integer)]
primitivePythagoras k =
foldr (\m xs -> iter 1 m ++ xs) [] [2 .. mh]
where mh = ceiling (sqrt (fromIntegral k * 0.5))
iter n m
| n >= m || 2 * m * (m + n) > k = []
| odd (m - n) && gcd m n == 1 =
let a = m * m - n * n
b = 2 * m * n
c = m * m + n * n
in (if a < b then (a, b, c) else (b, a, c)) : iter (n + 1) m
| otherwise = iter (n + 1) m
pythagoras :: Integer -> [(Integer, Integer, Integer)]
pythagoras k =
foldr check [] $ primitivePythagoras k
where check (a, b, c) xs =
let p = k `mod` (a + b + c)
q = k `div` (a + b + c)
in if p /= 0
then xs
else (q * a, q * b, q * c) : xs
pythagoras' :: Integer -> [(Integer, Integer, Integer)]
pythagoras' x = iter xs (reverse xs)
where xs = divisor (x * x `div` 2)
iter (y:ys) (z:zs) =
if z < y
then []
else if z >= x
then iter ys zs
else (a, b, c) : iter ys zs
where a = x - z
b = x - y
c = x - (a + b)
-- 平方根
sqrt' :: Double -> Double
sqrt' x
| x < 0 = error "sqrt': domain error"
| x == 0 = 0
| otherwise = iter (if x > 1 then init 1 x else 1)
where init s x
| s >= x = s
| otherwise = init (s * 2) (x / 2)
iter p = let q = (p + x / p) / 2
in if q >= p then p else iter q
iSqrt :: Integer -> Integer
iSqrt n
| n < 0 = error "iSqrt: domain error"
| n == 0 = 0
| otherwise = iter (init 1 n)
where init s x
| s >= x = s
| otherwise = init (s * 2) (x `div` 2)
iter p = let q = (n `div` p + p) `div` 2
in if q >= p then p else iter q
-- めのこ平方
iSqrt'' :: Integer -> Integer -> Integer
iSqrt'' n m
| n < m = m `div` 2
| otherwise = iSqrt'' (n - m) (m + 2)
iSqrt' :: Integer -> Integer
iSqrt' n
| n < 100 = iSqrt'' n 1
| otherwise = let m = 10 * iSqrt' (n `div` 100)
in iSqrt'' (n - m * m) (2 * m + 1)
-- 要素の選択
select :: [a] -> [(a, [a])]
select [x] = [(x, [])]
select (x:xs) = (x, xs) : map (\(y, ys) -> (y, x:ys)) (select xs)
-- 順列の生成
permutation :: Int -> [a] -> [[a]]
permutation 0 _ = [[]]
permutation n xs =
concatMap (\(y, ys) -> map (y:) (permutation (n - 1) ys)) $ select xs
-- 階乗
fact :: Integer -> Integer
fact n = if n == 0 then 1 else n * fact (n - 1)
-- 順列に番号をつける
fromPerm :: Eq a => [a] -> [a] -> Integer
fromPerm xs ys = iter xs ys 0
where
iter [_] _ n = n
iter zs@(x:xs) ys n = iter xs (delete x ys) (n + m * (fromIntegral i))
where m = fact $ fromIntegral (length zs - 1)
i = fromJust $ findIndex (== x) ys
-- 番号を順列に変換する (番号は 0 から開始)
toPerm :: Eq a => Integer -> [a] -> [a]
toPerm n xs = iter n xs []
where
iter _ [] a = reverse a
iter n xs a =
iter (n - m * p) (delete x xs) (x : a)
where m = fact (fromIntegral(length xs) - 1)
p = n `div` m
x = xs !! (fromIntegral p)