Word32 (無符号 32 bit 整数) とリストを使ってシンプルな無符号多倍長整数 (Bignum) を実装します。基数を 65536 とし、リストの先頭が下位の桁、末尾を上位の桁とします。簡単な例を示します。
0 => [0] 1 => [1] 65535 => [0xFFFF] 65536 => [0,1] 4294967295 => [0xFFFF,0xFFFF] 4294967296 => [0,0,1]
整数 (Integer) を多倍長整数に変換する関数 toBignum n を定義してください。
type Bignum = [Word32] toBignum :: Integer -> Bignum
ghci> toBignum 0 [0] ghci> toBignum 1 [1] ghci> toBignum 65535 [65535] ghci> toBignum 65536 [0,1] ghci> toBignum 0xffffffff [65535,65535] ghci> toBignum 0x100000000 [0,0,1]
多倍長整数を整数 (Integer) に変換する関数 bignumToInteger xs を定義してください。
bignumToInteger :: Bignum -> Integer
ghci> bignumToInteger [0] 0 ghci> bignumToInteger [1] 1 ghci> bignumToInteger [65535] 65535 ghci> bignumToInteger [0,1] 65536 ghci> bignumToInteger [0,0,1] 4294967296
2 つの多倍長整数 xs, ys の論理積を求める関数 bignumAND xs ys を定義してください。
bignumAND :: Bignum -> Bignum -> Bignum
ghci> bignumAND [1,1,1] [0] [0] ghci> bignumAND [1,1,1] [1] [1] ghci> bignumAND [1,2,1] [3,3,3] [1,2,1] ghci> bignumAND [1,2,1] [3,3,2] [1,2] ghci> bignumAND [1,2,1] [3,0,2] [1]
2 つの多倍長整数 xs, ys の論理和を求める関数 bignumOR xs ys を定義してください。
bignumOR :: Bignum -> Bignum -> Bignum
ghci> bignumOR [1,1,1] [0] [1,1,1] ghci> bignumOR [1,1,1] [2,2,2] [3,3,3] ghci> bignumOR [1,1,1] [2,2,2,2,2,2] [3,3,3,2,2,2]
2 つの多倍長整数 xs, ys の排他的論理和を求める関数 bignumXOR xs ys を定義してください。
bignumXOR :: Bignum -> Bignum -> Bignum
ghci> bignumXOR [0,1,2,3] [0] [0,1,2,3] ghci> bignumXOR [0,1,2,3] [1] [1,1,2,3] ghci> bignumXOR [0,1,2,3] [1,2,3,4] [1,3,1,7] ghci> bignumXOR [0,1,2,3] [0,1,2,3] [0]
多倍長整数 xs を左へ n bit シフトする関数 bignumShiftLeft xs n を定義してください。
bignumShiftLeft :: Bignum -> Int -> Bignum
ghci> bignumShiftLeft [1] 0 [1] ghci> bignumShiftLeft [1] 1 [2] ghci> bignumShiftLeft [1,1,1,1] 15 [32768,32768,32768,32768] ghci> bignumShiftLeft [1,1,1,1] 16 [0,1,1,1,1] ghci> bignumShiftLeft [1,1,1,1] 32 [0,0,1,1,1,1]
多倍長整数 xs を右へ n bit シフトする関数 bignumShiftRight xs n を定義してください。
bignumShiftRight :: Bignum -> Int -> Bignum
ghci> bignumShiftRight [0,1] 0 [0,1] ghci> bignumShiftRight [0,1] 1 [32768] ghci> bignumShiftRight [0,1] 15 [2] ghci> bignumShiftRight [0,1] 16 [1] ghci> bignumShiftRight [0,1] 17 [0] ghci> bignumShiftRight [0,0,1] 32 [1]
多倍長整数を比較する関数 bignumEQ, bignumNE, bignumGT, bignumGE, bignumLT, bignumLE を定義してください。
bignumEQ :: Bignum -> Bignum -> Bool -- == bignumNE :: Bignum -> Bignum -> Bool -- /= bignumGT :: Bignum -> Bignum -> Bool -- > bignumGE :: Bignum -> Bignum -> Bool -- >= bignumLT :: Bignum -> Bignum -> Bool -- < bignumLE :: Bignum -> Bignum -> Bool -- <=
ghci> bignumEQ [1,1,1] [1,1,1] True ghci> bignumEQ [1,1,1] [1,0,1] False ghci> bignumNE [1,1,1] [1,1,1] False ghci> bignumNE [1,1,1] [0,1,1] True ghci> bignumGT [0,0,1] [0,1] True ghci> bignumGT [0,0,1] [0,0,2] False ghci> bignumGE [0,0,1] [0,0,1] True ghci> bignumGE [0,0,2] [0,0,1] True ghci> bignumGE [0,0,1] [0,0,2] False ghci> bignumLT [0,0,2] [0,0,1] False ghci> bignumLT [0,2] [0,0,1] True ghci> bignumLE [0,0,1] [0,0,1] True ghci> bignumLE [0,0,1] [0,0,2] True ghci> bignumLE [0,0,1] [0,2] False
多倍長整数 xs と基数未満の整数 n を加算する関数 bignumAddInt xs n を定義してください。
bignumAddInt :: Bignum -> Word32 -> Bignum
ghci> bignumAddInt [65535,65535] 0 [65535,65535] ghci> bignumAddInt [65535,65535] 1 [0,0,1] ghci> bignumAddInt [65535,65535] 65535 [65534,0,1]
多倍長整数 xs. ys を加算する関数 bignumAdd xs ys を定義してください。
bignumAdd :: Bignum -> Bignum -> Bignum
ghci> bignumAdd [0,65535] [65535,1] [65535,0,1] ghci> bignumAdd [65535,65535] [65535,65535] [65534,65535,1] ghci> bignumAdd [65535,65535] [1,0,1] [0,0,2]
多倍長整数 xs と基数未満の整数 n を減算する関数 bignumSubInt xs n を定義してください。
bignumSubInt :: Bignum -> Word32 -> Bignum
ghci> bignumSubInt [1] 0 [1] ghci> bignumSubInt [1] 1 [0] ghci> bignumSubInt [1] 2 [65535*** Exception: bignum: underflow ghci> bignumSubInt [0,1] 1 [65535] ghci> bignumSubInt [0,0,0,0,1] 1 [65535,65535,65535,65535] ghci> bignumSubInt [0,0,0,0,1] 65535 [1,65535,65535,65535]
多倍長整数 xs. ys を減算する関数 bignumSub xs ys を定義してください。
bignumSub :: Bignum -> Bignum -> Bignum
ghci> bignumSub [0,0,1] [65535,65535] [1] ghci> bignumSub [0,0,1] [0,0,1] [0] ghci> bignumSub [0,0,1] [0,1] [0,65535] ghci> bignumSub [0,0,1] [1,0,1] [65535,65535,65535*** Exception: bignum: underflow
多倍長整数 xs と基数未満の整数 n を乗算する関数 bignumMulInt xs n を定義してください。
bignumMulInt :: Bignum -> Word32 -> Bignum
ghci> bignumMulInt [1,2,3,4,5] 0 [0] ghci> bignumMulInt [1,2,3,4,5] 1 [1,2,3,4,5] ghci> bignumMulInt [1,2,3,4,5] 2 [2,4,6,8,10] ghci> bignumMulInt [65535, 65535] 65535 [1,65535,65534]
多倍長整数 xs. ys を乗算する関数 bignumMul xs ys を定義してください。
bignumMul :: Bignum -> Bignum -> Bignum
ghci> bignumMul [1,1,1] [1,1,1] [1,2,3,2,1] ghci> bignumMul [65535,65535,65535] [1,1,1] [65535,65534,65534,0,1,1] ghci> bignumMul [65535,65535,65535] [65535,65535,65535] [1,0,0,65534,65535,65535]
多倍長整数 xs と基数未満の整数 n を除算する関数 bignumDivInt xs n を定義してください。返り値は商と剰余をタプルで返すものとします。
bignumDivInt :: Bignum -> Word32 -> (Bignum, Bignum)
ghci> bignumDivInt [2,4,6,8,10] 1 ([2,4,6,8,10],[0]) ghci> bignumDivInt [2,4,6,8,10] 2 ([1,2,3,4,5],[0]) ghci> bignumDivInt [2,4,6,8,10] 3 ([21846,1,2,21848,3],[0]) ghci> bignumDivInt [2,4,6,8,10] 4 ([0,32769,1,32770,2],[2]) ghci> bignumDivInt [2,4,6,8,10] 65535 ([28,24,18,10],[30]) ghci> bignumDivInt [2,4,6,8,10] 0 *** Exception: bignum: divide by zero
多倍長整数 xs. ys を除算する関数 bignumDiv xs ys を定義してください。返り値は商と剰余をタプルに格納して返すものとします。
bignumDiv :: Bignum -> Bignum -> (Bignum, Bignum)
ghci> bignumDiv [0,0,0,1] [0,0,0,1] ([1],[0]) ghci> bignumDiv [0,0,0,1] [0,0,1] ([0,1],[0]) ghci> bignumDiv [1,0,0,1] [0,0,1] ([0,1],[1]) ghci> bignumDiv [0,0,0,1] [65535,65535] ([0,1],[0,1])
多倍長整数 xs を文字列に変換する関数 bignumToString xs r を定義してください。r は基数を表す整数値 (r <= 16) です。
bignumToString :: Bignum -> Int -> String
ghci> bignumToString [0,0,0,1] 10 "281474976710656" ghci> bignumToString [0,0,0,1] 16 "1000000000000" ghci> bignumToString [65535,65535,65535] 10 "281474976710655" ghci> bignumToString [65535,65535,65535] 16 "FFFFFFFFFFFF"
文字列 str を多倍長整数に変換する関数 stringToBignum str r を定義してください。r は基数を表す整数値 (r <= 16) です。
stringToBignum :: String -> Int -> Bignum
ghci> stringToBignum "1234567890" 10 [722,18838] ghci> stringToBignum "FFFFFFFFFFFFFFFF" 16 [65535,65535,65535,65535] ghci> stringToBignum "1234567890ABCDEF" 16 [52719,37035,22136,4660] ghci> stringToBignum "12345678901234567890" 10 [2770,60191,43404,43860]
多倍長整数の平方根の整数部分を求める関数 bignumSqrt xs を定義してください。
bignumSqrt :: bignum -> bignum
ghci> bignumToInteger $ bignumSqrt $ toBignum $ 2 * 10^20 14142135623 ghci> bignumToInteger $ bignumSqrt $ toBignum $ 2 * 10^40 141421356237309504880 ghci> bignumToInteger $ bignumSqrt $ toBignum $ 2 * 10^80 14142135623730950488016887242096980785696 ghci> bignumToInteger $ bignumSqrt $ toBignum $ 3 * 10^20 17320508075 ghci> bignumToInteger $ bignumSqrt $ toBignum $ 3 * 10^40 173205080756887729352 ghci> bignumToInteger $ bignumSqrt $ toBignum $ 3 * 10^80 17320508075688772935274463415058723669428
リスト : 整数を多倍長整数に変換する base :: Word32 base = 65536 baseL :: Integer baseL = 65536 toBignum :: Integer -> Bignum toBignum x | x < 0 = error "toBignum : domain error" | x < baseL = [fromIntegral x] | otherwise = fromIntegral (x `mod` baseL) : toBignum (x `div` baseL)
toBignum は簡単です。x が baseL よりも小さければ x を Wrod32 に変換してからリストに格納して返します。あとは、整数 n を基数 baseL で割り算していき、剰余を toBignum の返り値 (リスト) に追加していくだけです。
リスト : 多倍長整数を整数に変換する bignumToInteger :: Bignum -> Integer bignumToInteger xs = foldr (\x a -> a * baseL + fromIntegral x) 0 xs
bignumToInteger も簡単です。引数 xs を foldr で末尾から畳み込むだけです。
リスト : 多倍長整数の論理積 bignumAND :: Bignum -> Bignum -> Bignum bignumAND xs ys = if null zs then [0] else zs where zs = dropWhileEnd (== 0) $ zipWith (.&.) xs ys
bignumAND は xs と ys の要素の論理積を zipWith で求め、末尾から連続する 0 を dropWhileEnd で取り除きます。その結果、zs が空リストになったならば [0] を返します。そうでなければ zs をそのまま返します。
リスト : 多倍長整数の論理和 bignumLogical :: (Word32 -> Word32 -> Word32) -> Bignum -> Bignum -> Bignum bignumLogical _ [] [] = [] bignumLogical f xs [] = xs bignumLogical f [] ys = ys bignumLogical f (x:xs) (y:ys) = f x y : bignumLogical f xs ys bignumOR :: Bignum -> Bignum -> Bignum bignumOR = bignumLogical (.|.)
bignumOR と bignumXOR の処理は関数 bignumLogical で行います。bignumLogical は xs と ys の要素を順番に取り出して関数 f を適用し、その結果を bignumLogical の返り値 (リスト) に追加します。ys が空リストになった場合は xs を返し、xs が空リストになった場合は ys を返します。両方とも空リストの場合は空リストを返します。
リスト : 多倍長整数の排他的論理和 bignumXOR :: Bignum -> Bignum -> Bignum bignumXOR xs ys = if null zs then [0] else zs where zs = dropWhileEnd (==0) $ bignumLogical xor xs ys
bignumXOR は bignumLogical を呼び出して、リストの要素を排他的論理和を求め、dropWhileEnd で末尾から連続する 0 を取り除きます。その結果、zs が空リストになったら [0] を返します。そうでなければ zs をそのまま返します。
リスト : 多倍長整数の左シフト
baseBit :: Int
baseBit = 16
mask :: Word32
mask = 0xffff
bignumShiftLeft16 :: Bignum -> Int -> Bignum
bignumShiftLeft16 xs b = iter xs 0
where iter [] c = if c == 0 then [] else [c]
iter (x:xs) c = z : iter xs (shiftR y baseBit')
where y = shiftL x b
z = (y .|. c) .&. mask
bignumShiftLeft :: Bignum -> Int -> Bignum
bignumShiftLeft xs n
| n == 0 || xs == [0] = xs
| n < baseBit = bignumShiftLeft16 xs n
| otherwise = replicate a 0 ++ bignumShiftLeft16 xs b
where a = n `div` baseBit
b = n `mod` baseBit
bignumShiftLeft は引数 n が 0 または xs が [0] の場合は xs をそのまま返します。baseBit 未満の場合は関数 bignumShiftLeft16 を呼び出します。そうでなければ、n を baseBit で除算し、商を a に、剰余を b にセットします。そして、xs を b ビットシフトした結果に a 個の 0 を先頭に追加します。
実際のビットシフトは bignumShiftLeft16 で行います。局所関数 iter でリストの要素を順番に取り出します。引数 c には左ビットシフトしたときに溢れたビットをセットします。x を b ビット左シフトした値を y にセットし、y と c の論理和を求め、それと mask の論理積の結果を z にセットします。この z を iter の返り値 (リスト) に追加します。溢れたビットは y を右へ 16 ビットシフトすれば求めることができます。
リスト : 多倍長整数の右シフト
bignumShiftRight16 :: Bignum -> Int -> Bignum
bignumShiftRight16 xs b = if null zs then [0] else zs
where zs = dropWhileEnd (==0) $ snd $ iter xs
mask = shiftL 1 b - 1
iter [x] = (x .&. mask, [shiftR x b])
iter (x:xs) = let (c, ys) = iter xs
a = (shiftR x b) .|. (shiftL c (baseBit - b))
in (x .&. mask, a:ys)
bignumShiftRight :: Bignum -> Int -> Bignum
bignumShiftRight xs n
| n == 0 || xs == [0] = xs
| n < baseBit = bignumShiftRight16 xs n
| otherwise = let ys = drop (n `div` baseBit) xs
b = n `mod` baseBit
in if null ys
then [0]
else bignumShiftRight16 ys b
bignumShiftRight は引数 n が 0 または xs が [0] の場合は xs をそのまま返します。baseBit 未満の場合は関数 bignumShiftRight16 を呼び出します。
そうでなければ、n を base_bit で除算し、xs の先頭から商の数だけ要素を取り除き、それを ys にセットします。剰余は b にセットします。もし ys が空リストならば [0] を返します。そうでなければ、bignumShiftRigth16 を呼び出して、ys を b ビット右へシフトします。
実際のビットシフトは関数 bignumShiftRight16 で行います。局所関数 iter で xs の要素を順番に取り出します。最初に iter を再帰呼び出しして、溢れたビット c とシフトした結果 ys を求めます。それから、x を右に b ビットシフトし、それと c を左へ baseBit - b ビットシフトした値の論理和を求めます。
この値 a を ys に追加し、溢れるビットを x と mask の論理積で求めます。mask の値は 1 を b ビット左シフトした値から 1 を引き算するだけです。これで下位 b ビットを 1 にすることができます。
リスト : 多倍長整数の比較
bignumCompare :: Bignum -> Bignum -> Ordering
bignumCompare [] [] = EQ
bignumCompare [] _ = LT
bignumCompare _ [] = GT
bignumCompare (x:xs) (y:ys) =
case bignumCompare xs ys of
EQ -> compare x y
r -> r
bignumEQ, bignumNE, bignumGT, bignumLT, bignumGE, bignumLE :: Bignum -> Bignum -> Bool
bignumEQ xs ys = bignumCompare xs ys == EQ
bignumNE xs ys = bignumCompare xs ys /= EQ
bignumGT xs ys = bignumCompare xs ys == GT
bignumLT xs ys = bignumCompare xs ys == LT
bignumGE xs ys = bignumCompare xs ys >= EQ
bignumLE xs ys = bignumCompare xs ys <= EQ
関数 bignumCompare は xs と ys が等しい場合は EQ を、xs のほうが大きい場合は GT を、xs のほうが小さい場合は LT を返します。
bignumCompare は 2 つのリストを順番にたどっていき、xs が先に空リストになったら LT を、ys が先に空リストになったら GT を返します。両方とも空リストになった場合は EQ を返します。bignumCompare を再帰呼び出しして、その結果が EQ の場合は関数 compare で x と y を比較します。そうでなければ、返り値 r をそのまま返します。あとは、各関数から bignumCompare を呼び出して結果を比較するだけです。
リスト : 多倍長整数と整数の加算 integerAdd :: Word32 -> Word32 -> Word32 -> (Word32, Word32) integerAdd x y c = let n = x + y + c in if n < base then (n, 0) else (n - base, 1) bignumAddInt :: Bignum -> Word32 -> Bignum bignumAddInt [] n = if n == 0 then [] else [n] bignumAddInt (x:xs) n = y : bignumAddInt xs c where (y, c) = integerAdd x 0 n
bignumAddInt は最下位の桁と引数 c を加算し、桁上がりがあればそれを上位の桁に加算します。あとは、桁上げの処理を繰り返すだけです。整数同士の加算は関数 integerAdd で行います。引数 x, y, c を加算し、その値 n が base 未満であれば n と 0 をタプルで返します。そうでなければ、n - base と 1 をタプルで返します。
リスト : 多倍長整数の加算
bignumAdd :: Bignum -> Bignum -> Bignum
bignumAdd xs ys = iter xs ys 0
where iter [] ys c = bignumAddInt ys c
iter xs [] c = bignumAddInt xs c
iter (x:xs) (y:ys) c =
let (n, c') = integerAdd x y c
in n : iter xs ys c'
bignumAdd の処理は局所関数 iter で行います。第 3 引数 c は桁上がりを表します。xs または ys が空リストの場合は bignumAddInt を呼び出して、もう一方の値と c を加算します。
そうでなければ、xs と ys の要素と c を integerAdd で加算します。あとは、iter に c' を渡してを再帰呼び出しし、その返り値に n を追加するだけです。
リスト : 多倍長整数と整数の減算
integerSub :: Word32 -> Word32 -> Word32 -> (Word32, Word32)
integerSub x y c =
let n = base + x - y - c
in if n < base then (n, 1) else (n - base, 0)
bignumSubInt :: Bignum -> Word32 -> Bignum
bignumSubInt xs n = if null zs then [0] else zs
where zs = dropWhileEnd (== 0) $ iter xs n
iter [] n =
if n > 0 then error "bignum: underflow" else []
iter (x:xs) n =
let (y, c) = integerSub x 0 n
in y : iter xs c
bignumSubInt は最下位の桁と引数 n を減算し、桁借りがあればそれを上位の桁から減算します。あとは、桁借りの処理を繰り返すだけです。実際の処理は局所関数 iter で行います。引数が空リストで、桁借りの値 n が正であれば、計算結果は負になるのでエラーを送出します。
そうでなければ、空リストを返します。あとは、integerSub で x - n を計算し、桁借り c を iter に渡して再帰呼び出しし、その返り値に y を追加するだけです。最後に、末尾から連続する 0 を dropWhileEnd で取り除き、その結果 zs が空リストであれば [0] を、そうでなければ zs を返します。
integerSub は base + x から y と c を減算して変数 n にセットします。もし、n が base よりも小さい場合、n と桁借りを表す 1 をタプルで返します。そうでなければ、n - base と 0 をタプルで返します。
リスト : 多倍長整数の減算
bignumSub :: Bignum -> Bignum -> Bignum
bignumSub xs ys = if null zs then [0] else zs
where zs = dropWhileEnd (==0) $ iter xs ys 0
iter [] ys c = if null ys && c == 0
then []
else error "bignum: underflow"
iter xs [] c = bignumSubInt xs c
iter (x:xs) (y:ys) c =
let (n, c') = integerSub x y c
in n : iter xs ys c'
bignumSub の処理は局所関数 iter で行います。iter は xs と ys の要素と桁借りを表す引数 c を integerSub で減算し、その結果を iter の返り値 (リスト) に格納していきます。xs が空リストの場合、ys が空リストで c が 0 であれば空リストを返します。そうでなければ、結果は負になるのでエラーを送出します。ys が空リストで xs が空リストでない場合、bignumSubInt で xs から c を減算します。最後に、dropWhileEnd で末尾から連続する 0 を取り除き、その結果 zs が空リストであれば [0] を、そうでなければ zs を返します。
リスト : 多倍長整数と整数の乗算
integerMul :: Word32 -> Word32 -> Word32 -> (Word32, Word32)
integerMul x y c =
let n = x * y + c
in if n < base
then (n, 0)
else (n `mod` base, n `div` base)
bignumMulInt :: Bignum -> Word32 -> Bignum
bignumMulInt xs n
| n == 0 = [0]
| n == 1 = xs
| otherwise = iter xs 0
where iter [] c =
if c == 0 then [] else [c]
iter (x:xs) c = let (y, c') = integerMul x n c
in y : iter xs c'
bignumMulInt は引数 n が 0 ならば [0] を、1 ならば xs をそのまま返します。それ以外の場合、最下位の桁から順番に x と掛け算して、iter の返り値 (リスト) に追加します。このとき、桁上がり c' を iter に渡して、上位の桁に足し算します。整数の乗算は関数 integerMul で行います。引数 x, y が乗算する整数、c が桁上がりで加算する値です。x * y + c を n にセットし、値が base 未満であれば n と 0 をタプルで返します。そうでなければ、n と base の剰余と商をタプルで返します。
リスト : 多倍長整数の乗算
bignumMul :: Bignum -> Bignum -> Bignum
bignumMul [x] ys = bignumMulInt ys x
bignumMul xs [y] = bignumMulInt xs y
bignumMul xs ys = iter xs ys [0]
where iter _ [] a = a
iter xs (y:ys) a =
iter (0:xs) ys (bignumAdd (bignumMulInt xs y) a)
多倍長整数同士の乗算は筆算と同じ方法で行います。簡単な例を示しましょう。
xs : (4 3 2 1)
ys : (7 6 5)
1 2 3 4
* 5 6 7
----------------------
7 14 21 28
6 12 18 24 0
5 10 15 20 0 0
----------------------
5 16 34 52 45 28
図 : 多倍長整数の乗算
上図のように、xs を 16 bit 左シフトしながら ys の要素を掛け算し、その値を加算していけばいいわけです。
bignumMul は引数 xs, ys が base 未満の整数であれば、bignumMulInt を呼び出して計算します。そうでなければ、xs と y の乗算を bignumMulInt で求め、累積変数 a にその値を bignumAdd で加算します。ys の次の要素を乗算するとき、xs の先頭に 0 を挿入して 16 bit 左シフトします。
なお、今回の方法は桁数が多くなると遅くなります。これよりも高速な方法として「Karatsuba 法」や「高速フーリエ変換」を使った方法があります。これらのアルゴリズムについては、Fussy さんの「乗算処理の高速化」, 「高速フーリエ変換」、M.Kamada さんの「離散フーリエ変換を用いた多倍長乗算の話」が参考になると思います。
リスト : 多倍長整数と整数の除算
integerDiv :: Word32 -> Word32 -> Word32 -> (Word32, Word32)
integerDiv x y c =
let n = c * base + x
in (n `div` y, n `mod` y)
bignumDivInt :: Bignum -> Word32 -> (Bignum, Bignum)
bignumDivInt xs 0 = error "bignum: divide by zero"
bignumDivInt xs 1 = (xs, [0])
bignumDivInt xs n = iter (reverse xs) 0 []
where iter [] c [] = ([0], [c])
iter [] c a = (a, [c])
iter (x:xs) c a =
let (p, q) = integerDiv x n c
in iter xs q (if p == 0 && null a then a else p:a)
bignumDivInt は引数 x が 0 の場合はエラーを送出し、1 の場合は xs と [0] をタプルで返します。それ以外の場合は、xs の上位の桁から順番に整数 x で除算していきます。このため、xs を reverse で反転しています。局所関数 iter の変数 c には上位の桁の余りをセットします。あとは、関数 integerDiv で xs の要素と n の除算を行います。このとき、c * base を加えてから x で割ることに注意してください。あとは商と剰余をタプルで返します。
bignumDivInt は上位の桁から処理を行うため、リストの末尾に 0 が付加されないように工夫する必要があります。商 p が 0 で、かつ累積変数 a が空リストの場合、p を a に追加しません。それ以外の場合は p を a に追加します。最後に、a が空リストであれば zero と [c] を、そうでなければ a と [c] をタプルで返します。
多倍長整数の除算は筆算と同じ方法で行いますが、かなり複雑な処理になります。ここではアルゴリズムの概略を説明するだけにとどめます。詳細は参考文献をお読みください。
リスト : 多倍長整数の除算 (擬似コード)
xs = (x1 x2 ... xn), ys = (y1 y2 ... ym) とし、xs / ys の商と剰余を求める
*base* / 2 <= ym * d < *base* を満たす d を求め、(xs * d) / (ys * d) を計算する
xs1 = xs * d とする
xs1 と同じ桁数になるよう (ys * d) の下位に 0 を追加たものを ys1 とする
このとき、追加した 0 の個数を s とする
qs = ()
while( s >= 0 ){
xs1 / ys1 の仮の商 q' を求める。
(1) xs1 が ys1 よりも少ない桁数の場合、q' は 0 である
(2) xs1 と ys1 の桁数 (n) が同じ場合、q' = xn / yn とする
(3) xs1 が n 桁, ys1 が n - 1 桁の場合、q' = min( (xn * base + xn-1) / yn-1, base - 1 ) とする
if( q' > 0 ){
ys2 = ys1 * q'
while( xs1 < ys2 ){
q' = q' - 1
ys2 = ys2 - ys1
}
xs1 = xs1 - ys2
}
q' を qs に追加する
ys1 の最下位から 0 を取り除く
s = s - 1
}
商は qs, 剰余は xs1 / d となる。
ポイントは仮の商 q' を求める処理です。ys1 の最上位の桁 ym が条件 (A) base / 2 <= ym < base を満たしている場合、(2) であれば q' は 0 か 1 になります。(3) であれば xs1 の上位 2 桁と ys1 の上位 1 桁 (ym) から仮の商を求めます。このとき、真の商を q とすると、条件 (A) を満たしているならば次の式が成り立ちます。
q <= q' <= q + 2
したがって、q の値は q', q' - 1, q' - 2 のどれかになります。ys2 = ys1 * q' を計算し、xs1 < ys2 であれば q' から 1 を、ys2 から ys1 を引き算します。これを xs1 >= ys2 になるまで繰り返しますが、最悪でも 2 回の繰り返しで済むわけです。
商 q が q' - 1 と q' - 2 になる例を示します。
xs1 = [65535, 65535, 32767] ys1 = [65535, 32768] q' = (32767 * base + 65535) / 32768 = 65535 ys2 = [65535, 32768] * 65535 = [1, 32766, 32768] > xs1 q' = q' - 1 = 65534 ys2 = ys2 - ys1 = [2, 65533, 32767] < xs1 q' = 65534, xs1 = xs1 - ys2 = [65533, 2] ----------------------------------------------------- xs1 = [65535, 0, 32767] ys1 = [65535, 32768] q' = (32767 * base + 0) / 32768 = 65534 ys2 = [65535, 32768] * 65534 = [2, 65533, 32767] > xs1 q' = q' - 1 ys2 = ys2 - ys1 = [3, 32764, 32767] > xs1 q' = q' - 1 ys2 = ys2 - ys1 = [4, 65531, 32766] < xs1 q' = 65532, xs1 = xs1 - ys2 = [65531, 5]
なお、(3) を満たしているとき、より高い精度で仮の商 q' を求める方法があります。有名なクヌース先生のアルゴリズムＤはこの方法を使っています。除算のアルゴリズムについては、参考文献 [2] がわかりやすくまとまっていると思います。また、乗算の処理が高速な場合、ニュートン法で ys の逆数 1 / ys を求め、xs * (1 / ys) を計算することで除算を高速に実行することができます。
擬似コードをそのままプログラムすると次のようになります。
リスト : 多倍長整数の除算
halfBase = 0x8000
-- シフトするビット数を求める
getShiftBit :: Word32 -> Int
getShiftBit n = iter n 0
where iter n c
| n >= halfBase = c
| otherwise = iter (n * 2) (c + 1)
-- 仮の商を求める
getQuot :: Bignum -> Bignum -> Word32
getQuot [] _ = 0
getQuot [x] [y] = x `div` y
getQuot (x1:x2:xs) [y] = (x2 * base + x1) `div` y
getQuot (_:xs) (_:ys) = getQuot xs ys
bignumDiv :: Bignum -> Bignum -> (Bignum, Bignum)
bignumDiv xs [y] = bignumDivInt xs y
bignumDiv xs ys
| bignumLT xs ys = ([0], xs)
| otherwise = iter s xs' ys' []
where d = getShiftBit $ last ys
xs' = bignumShiftLeft xs d
s = length xs' - length ys
ys' = bignumShiftLeft ys (baseBit * (fromIntegral s) + (fromIntegral d))
iter s xs1 ys1 q =
if s < 0
then (q, bignumShiftRight xs1 d)
else iter (s - 1)
xs1'
(tail ys1)
(if quot == 0 && null q then q else quot:q)
where (quot, xs1') = iter2 xs1 ys1
iter2 xs1 ys1 =
let quot = min (getQuot xs1 ys1) (base - 1)
ys2 = bignumMulInt ys1 quot
in if quot == 0
then (quot, xs1)
else iter3 xs1 ys1 ys2 quot
iter3 xs1 ys1 ys2 quot =
if bignumLT xs1 ys2
then iter3 xs1 ys1 (bignumSub ys2 ys1) (quot - 1)
else (quot, bignumSub xs1 ys2)
関数 getShiftBit は ys の最上位の値が base / 2 以上になるよう、左シフトするビット数を求めます。関数 getQuot は仮の商を求めます。xs が空リストならば、xs の桁は ys よりも少ないので 0 を返します。ys が末尾の要素で、かつ xs も末尾の要素であれば、同じ桁数なので x / y を返します。そうでなければ、xs の上位 2 桁を求め、それを y で割り算します。関数 bignumDiv は説明をそのままプログラムしただけなので、とくに難しいところはないと思います。
リスト : 多倍長整数を文字列に変換する
charTable :: String
charTable = "0123456789ABCDEF"
bignumToString :: Bignum -> Word32 -> String
bignumToString xs r = iter xs []
where iter [0] a = a
iter xs a = iter ys ((charTable !! (fromIntegral m)) : a)
where (ys, [m]) = bignumDivInt xs r
bignumToString は簡単です。bignumDivInt で xs を基数 r で割り算し、charTable から m 番目の要素を求め、それを累積変数 a のリストに追加します。この処理を xs が [0] になるまで繰り返し、最後に a を返します。
リスト : 文字列を多倍長整数に変換する
stringToBignum :: String -> Word32 -> Bignum
stringToBignum s r = iter s [0]
where iter [] a = a
iter (x:xs) a =
case findIndex (== toUpper x) charTable of
Nothing -> error "stringToBignum: illegal char"
Just n -> iter xs (bignumAddInt (bignumMulInt a r) (fromIntegral n))
stringToBignum も簡単です。局所関数 iter で 1 文字ずつ順番に取り出し、関数 findIndex で文字を数値 n に変換します。このとき、英小文字をモジュール Data.Char の関数 toUpper で英大文字に変換しています。文字が見つからない場合はエラーを送出します。あとは、bignumMulInt で累積変数 a と基数 r を掛け算し、それに bignumAddInt で n を加算していくだけです。最後に a を返します。
数 a の平方根 √a の値を求める場合、方程式 x2 - a = 0 を Newton (ニュートン) 法で解くことが多いと思います。方程式を f(x), その導関数を f'(x) とすると、ニュートン法は次の漸化式の値が収束するまで繰り返す方法です。
xn+1 = xn - f(xn) / f'(xn)
平方根を求める場合、導関数は f'(x) = 2x になるので、漸化式は次のようになります。
xn+1 = (xn + a / xn) / 2
参考文献『Ｃ言語による最新アルゴリズム事典』によると、√a より大きめの初期値から出発し、置き換え x <- (x + a / x) / 2 を減少が止まるまで繰り返すことで √a の正確な値を求めることができるそうです。
これをプログラムすると次のようになります。
リスト : 平方根の整数部を求める
bignumSqrt :: Bignum -> Bignum
bignumSqrt [0] = [0]
bignumSqrt xs = iter xs
where iter p = let (q, _) = bignumDiv xs p
q' = bignumShiftRight16 (bignumAdd p q) 1
in if bignumGE q' p then p else iter q'
漸化式を計算して変数 q' にセットし、q' がひとつ前の値 p 以上になったら p を返すだけです。
ところで、Newton 法は初期値の設定によって性能が大きく変わります。引数 xs をそのまま初期値として使う場合、xs が大きくなると実行速度はかなり遅くなります。そこで、xs を半分に分けて、上位の桁を初期値として使うことにします。ただし、そのままでは √xs よりも小さくなる場合があるので、16 ビット左シフトすることにします。
プログラムは次のようになります。
リスト : 平方根の整数部を求める (改良版)
bignumSqrt' :: Bignum -> Bignum
bignumSqrt' [0] = [0]
bignumSqrt' xs = iter (init xs)
where init xs = let k = length xs `div` 2
in if k == 0 then xs else 0 : drop k xs
iter p = let (q, _) = bignumDiv xs p
q' = bignumShiftRight16 (bignumAdd p q) 1
in if bignumGE q' p then p else iter q'
init で初期値を求めます。xs を半分の桁を k とし、k が 0 ならば xs をそのまま返します。そうでなければ、drop で xs から k 個の要素を削除し、その先頭に 0 を追加します。
簡単な実行例を示します。
ghci> bignumToInteger $ bignumSqrt $ toBignum $ 2 * 10^100 141421356237309504880168872420969807856967187537694 (0.12 secs, 46,704,752 bytes) ghci> bignumToInteger $ bignumSqrt' $ toBignum $ 2 * 10^100 141421356237309504880168872420969807856967187537694 (0.04 secs, 12,779,384 bytes) ghci> bignumToInteger $ bignumSqrt $ toBignum $ 2 * 10^200 14142135623730950488016887242096980785696718753769480731766797379907324784621070388503875343276415727 (0.66 secs, 287,033,664 bytes) ghci> bignumToInteger $ bignumSqrt' $ toBignum $ 2 * 10^200 14142135623730950488016887242096980785696718753769480731766797379907324784621070388503875343276415727 (0.07 secs, 27,698,400 bytes) 実行環境 : GHC ver 9.6.6, Ubuntu 22.04 (WSL2), Intel Core i5-6200U 2.30GHz
bignumSqrt' のほうが高速ですね。改良の効果は十分に出ていると思います。