Haskell 超巨大実数ライブラリ
やけくそにでかい実数を適当な精度で扱えるライブラリを作ってみた。
きっかけは小数小野小町算のツイート。うん、屁理屈なんだ、すまない。
でもどんくらいでっかい数がでてくるものか、興味があったんだ。
バグとかあったら教えてください。
追記:案の定バグってたんで、テストつけてバグとって更新しました。
module Main ( DoublePlus(..), doublePlus, logB, main ) where import Control.Monad import Data.List import System.Random base = 10 data DoublePlus = Imm Double | Exp DoublePlus deriving (Eq, Ord, Show, Read) logB a = log a / log base doublePlus a = if a < 10 then Imm a else Exp $ doublePlus $ logB a eval (Imm a) = a eval (Exp a) = base ** eval a add small large = let ans = eval small + eval large in if ans < 1e10 then doublePlus ans else add' small large add' small large@(Exp b) = let increase = case small/large of Imm a -> log (1+a) in Exp (b + Imm increase) sub large@(Imm a) small@(Imm b) = Imm (a-b) sub large@(Exp a) small = let increase = case small/large of Imm b -> log (1-b) in Exp (a + Imm increase) instance Num DoublePlus where Imm a + Imm b = doublePlus (a+b) a + b = if a < b then add a b else add b a Imm a * Imm b = doublePlus (a*b) Imm a * Exp b = Exp (Imm (logB a) + b) Exp a * Imm b = Imm b * Exp a instance Fractional DoublePlus where {- assert a < b -} Imm a / Imm b = Imm (a/b) Imm a / Exp b = Imm $ a / eval (Exp b) Exp a / Exp b = Imm $ 1 / base**eval (sub b a) fromRational x = doublePlus (fromRational x) instance Floating DoublePlus where Imm a ** Imm b = doublePlus (a**b) Imm a ** Exp b = Exp (Imm (logB a) * Exp b) Exp a ** b = Exp (a * b) main = main2 main2 = do putStrLn "hi" print $ ans where ans = foldr1 (**) $ map Imm [ 1.1 , 1.2 , 1.3 , 1.4 , 1.5 , 1.6 , 1.7 , 1.8 , 1.9 , 2.0] ans2 = foldr1 (**) $ map Imm [ 1.4, 1.5, 1.6 , 1.7, 1.8, 1.9, 2.0, 1.3/(1.2-1.1)] maintest = do testSets <- fmap concat $ replicateM 1000 makeTestSet mapM_ printTest $ sort $ map badness $ testSets makeTestSet = do n' <- (randomIO) let n = 3 + mod (div n' 10000) 10 rands <- replicateM n (randomIO) let xs = map ((1.0::Double)+) rands let infinity = 2**100000 let ans = foldr1 (**) xs return $ if ans < infinity then [(ans,xs)] else [] badness (ans, xs) = let numAns = eval $ foldr1 (**) $ map Imm xs gosa = abs(numAns-ans) / ans msg = " or " ++ show ans ++ " vs " ++ show numAns ++ " for " ++ show xs in (gosa,msg) printTest (gosa,msg) = do putStrLn $ "Error was " ++ show gosa ++ " " ++ msg