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