-- BigFloat: Infinite-precision floating point type for Haskell. -- Recoded from my original "bignum" library written in Miranda. -- -- Copyright (c) Martin Guy, martinwguy@yahoo.it, October 2004. -- -- Email: Web: http://medialab.freaknet.org/bignum -- Free bignum server in Internet: telnet medialab.freaknet.org 31415 -- -- You can do as you please with this code, as long as you neither try to -- prevent anyone else from doing as *they* please with it, nor tell lies! -- -- This is an original work, except for edigits, modified from DAT's example -- script, and the scientific algorithms ripped from the library of bc. -- -- History: -- Version 0, April 1986, Canterbury, UK: addition and multiplication -- of positive numbers only. -- Version 1, October 2001, Raddusa, Sicily, with full arithmetic, e, pi, -- sqrt, sin, cos, atan. -- Version 2, 22 dec 2003, with logarithms and space and speed improvements. -- Version 3, October 2004, Newcastle, England: addition of bn_abs, bn_trunc, -- bn_frac, bn2int and derived scientific fns. bn_rnd removed. -- More speed improvements. First Haskell version (see WARNING below). -- -- Bugs: -- The exp function hangs after giving a few digits on certain critical -- values such as 0.2 module BigFloat (BigFloat,readBigFloat) where -- Exports type BigFloat, a floating point type of unlimited precision and -- member of classes Eq, Ord, Read, Show, Num, Fractional, Floating. import Numeric(lexDigits) import Char (isDigit) import Ratio -- to build fromRational {- WARNING!!! WARNING!!! WARNING!!! WARNING!!! WARNING!!! WARNING!!! WARNING!!! Under nhc98-1.16 and 1.18, programs using sin, cos, exp, sqrt and functions derived from these, will bomb out with Segmentation fault after a number of result digits. In other cases, they may silently give the wrong results! This is to a bug in the runtime of nhc98's Integer div/mod code. Under mira and hugs-Nov2003 they are fine. I have not tested with ghc or hbc. You can now work around this by setting integer_division_is_broken to True, The replacement code is not fast, but stops nhc98 from bombing out. WARNING!!! WARNING!!! WARNING!!! WARNING!!! WARNING!!! WARNING!!! WARNING!!! -} -- Enable workarounds for Integer div/mod bug in nhc98 ? -- The workaround code is roughly 10,000 times slower. integer_division_is_broken = False {- Note: If you want to be able to process infinite-length text input and still obtain results you must use readBigFloat directly instead of "read" (because the definition of "read" in te prelude insists on reading all the decimal places before returning any of them for processing). readBigFloat is in the ReadS family, so you pass it a String and it returns [] if there is not a valid BigFloat at the start of the string, or [(bf,t)] where bf is your BigFloat and t is the rest of the string that followed the converted text. Valid BigFloat strings are: an optional minus sign followed by one or more digits, optionally followed by a dot and one or more digits. Unlike Haskell, BigFloat also accepts ".5" as a valid form. -} -- For testing: check for impossible cases and complain if so. -- You can set this False when you are using this for real work -- for a slight speed increase. None of these traps ever spring. sanity_checks = False ---- ---- End of introductory gumf. Start of BigFloat library ---- ---- -- This source file is in three parts: -- the first section is native Haskell, implementing the new types -- and simple operations on them, including Eq, Ord, Read & Show. -- the second (short) section creates the Haskell bindings to put -- BigFloat into classes Num, Rational and Floating -- the third, by far the longest section, is the guts of my old Miranda code -- implementing the operations on these last three classes but using -- simple type synonyms Bigfloat & UBigfloat instead of the newtypes -- BigFloat and UBigFloat. This saves a lot of tedious mucking about -- converting to and from different types, but means that comparison -- and arithmetic operations on them are all functions, not operators. ---- Types ---- -- Component parts of our BigFloats. -- -- The Sign only ever has value of 1 or -1. This is better -- than a "data" type or Bool, because we can multiply -- signs together and suchlike, which turns out handy. -- -- The Exponent is 0 for values from 0.1 to 0.999.. and is positive -- for values >= 1.0, negative for 0.0999.. down. -- -- The Mantissa is a list of decimal digits which never starts -- with a 0 (i.e. is always normalised) though to save us endless checking -- it may end in one or more zeroes (these are always removed when a value -- is "shown"). It may also be an infinite list. -- -- Zero is always represented with a positive sign, exponent of zero -- and an empty mantissa. -- -- The value of a BigFloat is thus Sign * 0.[Mantissa] * 10 ^ Exponent. type Sign = Int type Exponent = Int type Mantissa = [Int] newtype BigFloat = MakeBigFloat (Sign,Exponent,Mantissa) newtype UBigFloat = MakeUBigFloat (Exponent,Mantissa) newtype M = MakeM (Mantissa) ---- Constructors ---- -- Construct a BigFloat from its component parts, normalising as we go. -- Ensure that the 1st digit of the mantissa is always significant. toBigFloat (s,e,0:x) = toBigFloat (s,e-1,x) -- Always give zero in its normal form toBigFloat (s,e,[]) = bf_0 -- Other cases are already normalised toBigFloat (s,e,m) = MakeBigFloat (s,e,m) -- Construct an unsigned BigFloat from its component parts, similarly. toUBigFloat (e,m@(a:x)) = if a == 0 then toUBigFloat (e-1, x) -- Normalise the abnormal else MakeUBigFloat (e,m) toUBigFloat (e,[]) = ubf_0 -- Canonical zero -- Mantissa type constructor toM m = MakeM m -- Canonical forms of zero bf_0 = MakeBigFloat (1,0,m_0) ubf_0 = MakeUBigFloat (0,m_0) m_0 = [] -- Simple function to test the Sign s_is_negative s = (s < 0) -- Is a Mantissa numerically equivalent to zero? m_is_zero :: Mantissa -> Bool m_is_zero m = and (map ((==) 0) m) ---- Deconstructors ---- fromBigFloat (MakeBigFloat (s,e,m)) = (s,e,m) fromUBigFloat (MakeUBigFloat(e,m)) = (e,m) fromM (MakeM m) = m ---- Equality functions ---- -- Equality: Since BigFloat and UBigFloat are always normalised, we can blindly compare -- sign and exponent. Mantissae require more careful comparison because -- two numerically equal ones may have different amounts of trailing zeroes. instance Eq BigFloat where MakeBigFloat(as,ae,am) == MakeBigFloat(bs,be,bm) = (as == bs) && (MakeUBigFloat(ae,am) == MakeUBigFloat(be,bm)) instance Eq UBigFloat where MakeUBigFloat(ae,am) == MakeUBigFloat (be,bm) = (ae == be) && (toM am == toM bm) instance Eq M where MakeM a == MakeM b = m_eq a b -- Are two mantissae numerically equivalent? m_eq :: Mantissa -> Mantissa -> Bool m_eq (a:x) (b:y) = (a == b) && m_eq x y m_eq a [] = m_is_zero a m_eq [] b = m_is_zero b ---- Ordering functions ---- instance Ord BigFloat where compare (MakeBigFloat(as,ae,am)) (MakeBigFloat(bs,be,bm)) | (as /= bs) = compare as bs | (as == bs) = if not (s_is_negative as) then compare (MakeUBigFloat (ae,am)) (MakeUBigFloat (be,bm)) else compare (MakeUBigFloat (be,bm)) (MakeUBigFloat (ae,am)) instance Ord UBigFloat where compare (MakeUBigFloat (ae,am)) (MakeUBigFloat (be,bm)) | (ae > be) = GT | (ae < be) = LT | (ae == be) = m_cmp am bm instance Ord M where compare (MakeM a) (MakeM b) = m_cmp a b (MakeM a) < (MakeM b) = m_lt a b (MakeM a) > (MakeM b) = m_lt b a (MakeM a) <= (MakeM b) = not (m_lt b a) (MakeM a) >= (MakeM b) = not (m_lt a b) -- Compare two mantissae. m_cmp :: Mantissa -> Mantissa -> Ordering m_cmp (a:x) (b:y) | a > b = GT | a < b = LT | a == b = m_cmp x y m_cmp a [] = if m_is_zero a then EQ else GT m_cmp [] b = if m_is_zero b then EQ else LT -- Optimised less-than m_lt (a:x) (b:y) = if a == b then m_lt x y else (a < b) m_lt a [] = False m_lt [] b = not (m_is_zero b) m_gt a b = m_lt b a m_le a b = not (m_gt a b) m_ge a b = not (m_lt a b) ---- Show functions ---- -- output syntax (as per Haskell): -- [-] digit {digit} [ "." digit {digit} ] -- [x] = optional -- {x} = 0 or more instance Show BigFloat where show (MakeBigFloat (s,e,[])) = "0" show (MakeBigFloat (-1,e,m)) = '-':show (MakeUBigFloat(e,m)) show (MakeBigFloat ( 1,e,m)) = show (MakeUBigFloat(e,m)) instance Show UBigFloat where show (MakeUBigFloat (e,m)) = showWhole (e,m) ++ showFrac (e,m_trim m) where -- show the bit before the point showWhole (e,m) | e <= 0 = "0" | e > 0 = take e (map int2digit m ++ repeat '0') -- show the point and the bit after it (if any) -- Fortunately, drop and replicate on negative values do the -- same as on zero, which simplifies this code. showFrac (e,m) | frac == [] = "" | otherwise = '.' : replicate (-e) '0' ++ showMantissa frac where frac = drop e m -- convert string of Int digits to string of characters showMantissa m = map int2digit m -- convert Int (0..9) to corresponding character int2digit i = toEnum (fromEnum '0' + i) -- Strip trailing zeroes from the end of a mantissa m_trim :: Mantissa -> Mantissa m_trim l = case l of (0:x) -> if m_is_zero x then [] else 0:m_trim x (a:x) -> a : m_trim x [] -> [] ---- Read functions ---- -- Read bigfloats, a modified version of readFloat from the Prelude. -- -- Input syntax is as per Haskell (including "10." being a total failure rather -- than a partial match), except that we don't accept the "E+10" notation since -- this would exclude the possibility of reading an infinite number of decimal -- places. instance Read BigFloat where readsPrec _ = readBigFloat -- readBigFloat is the read funtion for BigFloats. -- We cannot simply use "readSigned readUBigFloat" here cos a) it only works on -- Real, and Real requires a toRational method, which unlimited bigfloats -- can never supply, and b) it uses "lex", which bottoms out on infinite lists -- of decimal places. -- Like readSigned, this tries to match a positive BigFloat and a negative BigFloat and -- concatenates the two resulting lists (since both cannot succeed). readBigFloat :: ReadS BigFloat readBigFloat r = [ (MakeBigFloat(1,e,m),t) | (MakeUBigFloat(e,m),t) <- readUBigFloat r] ++ -- We apply toBigFloat rather than MakeBigFloat to deal with "-0" [ (toBigFloat(-1,e,m),t) | r /= [] && head r == '-', (MakeUBigFloat(e,m),t) <- readUBigFloat (tail r) ] -- Since the exponent is determined by the whole number part, this accepts -- an infinite list of decimal places and starts returning something useful -- as soon as it finds the ".". -- The exponent of the result is the number of digits before the point, and -- toUBigFloat takes care of normalising out any initial zeroes. -- Unlike Haskell, we also accept ".5" as valid syntax bcos it is likely that -- existing sets of data are in that form. readUBigFloat :: ReadS UBigFloat readUBigFloat r = if r /= "" && head r == '.' then [(toUBigFloat(0,map dig2Int f),t) | (f,t) <- readFrac r] else [(toUBigFloat(length w,map dig2Int (w++f)),t) | (w,s) <- readWhole r, (f,t) <- readFrac s] where -- readWhole reads the digits before the optional "." -- and returns exponent, mantissa and "the rest", -- or an empty list if no match. readWhole r = [ (digits, s) | (digits,s) <- lexDigits r ] -- convert a (Dhar) digit to the corresponding Int dig2Int :: Char -> Int dig2Int dig = fromEnum dig - fromEnum '0' -- readFrac gives the digits after the ".", if any. -- If there's a numeric part after a ".", return it and "rest"; -- If not, return an empty mantissa and "the rest" as passed; -- a trailing dot followed by no digits fails the whole match. -- We cannot use lexDigits here bcos it cannot handle endless -- lists of digits (for the whole part it makes no difference -- since we cannot return anything until we are sure that the -- whole part is properly terminated) readFrac ('.':digits) | digits == [] || not (isDigit (head digits)) = [] | otherwise = lexWhile isDigit digits -- no ".123" part at all is a successful match of an integer readFrac s = [("",s)] -- lexWhile is a generic lexing function that matches zero or -- more characters for which the supplied function is True. -- Its return value is in the style of the read* functions: -- [(s,t)] where s is the characters matched and t is the rest -- of the supplied string. Unlike most read* and lex* functions, -- it never fails giving []. lexWhile f x = [(takeWhile f x, dropWhile f x)] ---- Glue to overload Haskell operators using the old Miranda-based code ---- -- We use MakeBigFloat here instead of toBigFloat since all the bn_* functions -- normalise their results. instance Num BigFloat where x + y = MakeBigFloat (bn_add (fromBigFloat x) (fromBigFloat y)) x - y = MakeBigFloat (bn_sub (fromBigFloat x) (fromBigFloat y)) x * y = MakeBigFloat (bn_mul (fromBigFloat x) (fromBigFloat y)) negate = MakeBigFloat . bn_neg . fromBigFloat abs = MakeBigFloat . bn_abs . fromBigFloat signum = MakeBigFloat . bn_signum . fromBigFloat fromInteger = MakeBigFloat . bn . show instance Fractional BigFloat where x / y = MakeBigFloat ((fromBigFloat x) `bn_div` (fromBigFloat y)) recip = MakeBigFloat . (bn_div bn_1) . fromBigFloat fromRational r = MakeBigFloat (((bn.show.toInteger.numerator) r) `bn_quot` (toInteger (denominator r))) instance Floating BigFloat where pi = MakeBigFloat bn_pi exp = MakeBigFloat . bn_exp . fromBigFloat log = MakeBigFloat . bn_ln . fromBigFloat sqrt = MakeBigFloat . bn_sqrt . fromBigFloat x ** y = MakeBigFloat ( (fromBigFloat x) `bn_pow` (fromBigFloat y) ) sin = MakeBigFloat . bn_sin . fromBigFloat cos = MakeBigFloat . bn_cos . fromBigFloat tan = MakeBigFloat . bn_tan . fromBigFloat atan = MakeBigFloat . bn_atan . fromBigFloat sinh = MakeBigFloat . bn_sinh . fromBigFloat cosh = MakeBigFloat . bn_cosh . fromBigFloat tanh = MakeBigFloat . bn_tanh . fromBigFloat asinh = MakeBigFloat . bn_asinh . fromBigFloat acosh = MakeBigFloat . bn_acosh . fromBigFloat atanh = MakeBigFloat . bn_atanh . fromBigFloat -- Reverse glue for testing purposes bn_show bf = show (MakeBigFloat bf) -------------------------- HERE BEGINNETH... ---------------------------------- -- This is just the guts of the old Miranda code converted into -- Haskell and munged a bit. That's why it looks a bit odd. -- From here on, the BigFloat and UBigFloat types are replaced by -- Bigfloat and UBigfloat, which are simple type synonyms. -- The urge in this code is always to be able to return as many digits as -- possible of the result by inspection of as few as possible of the -- operands. Being able to return the sign and exponent independently of the -- mantissa calculation is also a win. -- bignum.m converted to Haskell by mira2hs on Wed Jul 28 02:20:42 CEST 2004 -- base: The number base we are working in. -- It's not guaranteed that anything other than 10 will work - but it might. -- ubn_text (text to bignum conversion) only works up to base 10. -- NOTE: If you want to use odd number bases (3,5,7 etc) you must reinstate -- the commented lines in bn_twice and bn_half. base :: Int base = 10 -- Constants derived from base. base_1 = base - 1 ---- Simple conversion functions ---- -- bn converts a decimal string representation of a number into a Bigfloat. -- It's the easiest way of entering constants into an equation. bn :: [Char] -> Bigfloat bn str = (s,e,m) where [(MakeBigFloat(s,e,m),_)] = readBigFloat str integer2bn :: Integer -> Bigfloat integer2bn = bn.show ---- Declarations for functions defined below ---- -- signed <-> unsigned Bigfloat converters bn2ubn :: Bigfloat -> UBigfloat ubn2bn :: UBigfloat -> Bigfloat -- bn_normalise strips leading zeroes from mantissa, adjusting -- exponent accordingly. It also corrects the "-0" case. bn_normalise :: Bigfloat -> Bigfloat -- Functions dealing directly with Bigfloats: -- bn_make and bn_unmake construct and destruct a Bigfloat. bn_make :: (Int,Int,[Int]) -> Bigfloat bn_unmake :: Bigfloat -> (Int,Int,[Int]) -- signof returns the exponent of a Bigfloat sign_of :: Bigfloat -> Sign -- exponent_of returns the exponent of a Bigfloat exponent_of :: Bigfloat -> Exponent -- mantissa_of returns the mantissa of a Bigfloat mantissa_of :: Bigfloat -> Mantissa -- bn_is_neg tells you whether a Bigfloat is a negative number -- (0 is not negative!) -- bn_is_zero tells whether it has value 0. -- bn_is_integer tells you whether a Bigfloat represents a whole number. bn_is_neg, bn_is_zero, bn_is_integer :: Bigfloat -> Bool -- bn_cmp compares two Bigfloats bn_cmp :: Bigfloat -> Bigfloat -> Ordering -- The six usual comparison functions: -- == /= > < >= <= bn_eq, bn_ne, bn_gt, bn_lt, bn_ge, bn_le :: Bigfloat -> Bigfloat -> Bool -- bn_abs returns absolute (positive) value of a bigfloat -- bn_neg negates a Bigfloat bn_abs, bn_neg, bn_signum :: Bigfloat -> Bigfloat -- bn_trunc returns the integer BF if integral, or the nearest one to it -- between it and zero. -- bn2Int returns the same value as an Int, bn2Integer as an Integer bn_trunc :: Bigfloat -> Bigfloat bn2Integer :: Bigfloat -> Integer bn2Int :: Bigfloat -> Int -- bn_frac is the complement of bn_trunc, returning the part that bn_trunc -- abandons. Some laws: -- -1 < bn_frac x < 1 -- bn_trunc x + bn_frac x = x bn_frac :: Bigfloat -> Bigfloat -- bn_add adds two Bigfloats together, returning their sum. -- bn_sub subtracts its second Bigfloat parameter from its first, -- returning the difference. Since negative numbers are not -- representable (yet), a negative result is a fatal error. bn_add, bn_sub :: Bigfloat -> Bigfloat -> Bigfloat -- bn_twice multiplies by two. -- bn_half divides by two. bn_twice, bn_half :: Bigfloat -> Bigfloat -- bn_times multiplies an Integer by a Bigfloat by repeated addition. bn_times :: Integer -> Bigfloat -> Bigfloat -- bn_mul multiplies two Bigfloats. bn_mul :: Bigfloat -> Bigfloat -> Bigfloat -- bn_raise raises a Bigfloat to an Integer power bn_raise :: Bigfloat -> Integer -> Bigfloat -- bn_quot divides a Bigfloat by an Integer bn_quot :: Bigfloat -> Integer -> Bigfloat -- divide-by-int version for intermal use bn_quotI :: Bigfloat -> Int -> Bigfloat -- bn_div performs division on two Bigfloats bn_div :: Bigfloat -> Bigfloat -> Bigfloat -- bn_sqr returns the square of its argument. -- bn_sqrt returns the square root of its argument. bn_sqr, bn_sqrt :: Bigfloat -> Bigfloat -- Logarithmic and trigonometric functions bn_ln, bn_exp :: Bigfloat -> Bigfloat bn_pow :: Bigfloat -> Bigfloat -> Bigfloat bn_sin, bn_cos, bn_atan :: Bigfloat -> Bigfloat -- Here beginneth... type Bigfloat = ( Sign, Exponent, Mantissa ) type UBigfloat = ( Exponent, Mantissa ) -- Note: herein, when we use type "mantissa", we mean [num] in which -- every element x is 0 <= x < base. -- In all cases dealing with partial results which may contain values -- outside the permitted range for our number base, we say [num]. -- Some simple constants bn_0 = bn_make (1, 0, []) bn_1 = bn_make (1, 1, [1]) bn_2 = bn_make (1, 1, [2]) -- unsigned Bigfloat constants ubn_0 :: UBigfloat ubn_0 = (0, []) -- -- ---------- simple conversion functions ---------- -- -- make a Bigfloat from its parts and vice versa (type-fooling glue) bn_make (s,e,m) = (s,e,m) bn_unmake (s,e,m) = (s,e,m) -- functions to select parts of a Bigfloat sign_of (s,e,m) = s exponent_of (s,e,m) = e mantissa_of (s,e,m) = m -- functions to select parts of a UBigfloat uexponent_of (e,m) = e umantissa_of (e,m) = m -- Convert Bigfloat to unsigned Bigfloat and vice versa bn2ubn (s,e,m) | sanity_checks && s /= 1 = error "Internal error: bn2ubn of negative number" | otherwise = (e,m) ubn2bn (e,m) = bn_normalise (1,e,m) -- bn_normalise converts a Bigfloat into "canonical form", that is to say: -- - the first digit of the mantissa is significant (not 0). -- - convert the value zero to canonical form (eliminating "-0"). -- To do this it strips off leading zeroes from mantissa, adjusting exponent -- accordingly. bn_normalise (s, e, []) = bn_0 bn_normalise (s, e, 0:m) = bn_normalise (s, e-1, m) bn_normalise other = other -- -- ----- utility functions used here and there in the rest of the library ----- -- -- ubn_same_exp takes a pair of UBigfloats and returns the same pair such that -- they have the same exponent. This means prefixing a load of zeroes and -- upping the exponent of the UBigfloat with the smaller exponent. ubn_same_exp :: (UBigfloat, UBigfloat) -> (UBigfloat, UBigfloat) ubn_same_exp (a@(ae,am), b@(be,bm)) | ae == be = ((ae,am), (be,bm)) | ae < be = (ubn_abnormalise be a, b) | ae > be = (a, ubn_abnormalise ae b) -- ubn_abnormalise sets the exponent of its second parameter to its first -- parameter, giving as a result an un-normalised Bigfloat. This simply involves -- prefixing a quantity of zeroes to the mantissa and adjusting the exponent. -- Of course, the new exponent must be >= the old exponent! ubn_abnormalise :: Exponent -> UBigfloat -> UBigfloat ubn_abnormalise new_e (old_e, m) | sanity_checks && new_e < old_e = error "ubn_abnormalise called with duff parameters" | otherwise = (new_e, prefix0s (new_e - old_e) m) -- Add a number of zeroes to the head of a mantissa. prefix0s n x = if sanity_checks && n < 0 then error "Internal error in BigFloat.prefix0s" else replicate n 0 ++ x -- -- ----------- Simple one-argument functions ------------ -- -- Return absolute value of a Bigfloat bn_abs (s,e,m) = (1,e,m) -- Negate a Bigfloat, avoiding the "-0" syndrome. bn_neg (s,e,m) | m_is_zero m = bn_0 | otherwise = (-s,e,m) bn_signum (s,e,m) | m_is_zero m = bn_0 | otherwise = (s,1,[1]) -- +/- bn_1 -- Truncate a Bigfloat towards 0, giving a Bigfloat bn_trunc (s,e,m) | e <= 0 || m_is_zero m = bn_0 | otherwise = (s,e,take e m) -- Truncate a Bigfloat towards 0, giving an Integer bn2Integer (s,e,m) | e <= 0 || m_is_zero m = 0 -- Not strictly necessary. | otherwise = toInteger s * m2Integer (take e (m ++ repeat 0)) -- Truncate a Bigfloat towards 0, giving an Int bn2Int (s,e,m) | e <= 0 || m_is_zero m = 0 -- Not strictly necessary | otherwise = s * m2Int (take e (m ++ repeat 0)) -- Convert an Int or Integer to a Bigfloat i2bn :: Integral i => i -> Bigfloat i2bn = bn.show -- Convert a mantissa to an Integer m2Integer :: [Int] -> Integer m2Integer = foldl f 0 where f a b = toInteger base * a + toInteger b -- Convert a mantissa to an Int m2Int :: [Int] -> Int m2Int m = -- Integers go up to 2^31 = 2,147,XXX,XXX -- (well, 2^29 is guaranteed by the language definition) if sanity_checks && length (take 10 m) == 10 then error "Integer overflow in m2Int" else foldl f 0 m where f a b = base * a + b -- Just give the part after the decimal point bn_frac b@(s,e,m) | e <= 0 = b | otherwise = bn_normalise (s,0,drop e m) -- -- ----------- Stuff for addition/subtraction of Bigfloats ------------ -- -- bn_add gives the arithmetic sum of two Bigfloats. -- 1) both +ve or both -ve -- 2) sign differs and absolute value of a >= abs value of b. -- 3) sign differs and abs(a) < abs(b) bn_add (as,ae,am) (bs,be,bm) = bn_normalise (s,e,m) where (s,(e,m)) | as == bs = (as, ubn_add (ae,am) (be,bm)) | (be,bm) `ubn_lt` (ae,am) = (as, ubn_sub (ae,am) (be,bm)) | otherwise = (bs, ubn_sub (be,bm) (ae,am)) -- bn_sub performs subtraction of 2 Bigfloats. bn_sub a b = bn_add a (bn_neg b) -- ubn_add gives the arithmetic sum of two unsigned Bigfloats. -- Strategy: denormalise the arguments so that the have the same exponent, -- form the sum of the mantissae, then perform the sum after tacking a 0 onto -- the heads to catch possible carry from the first digit. It would be -- slightly faster to do m_carry ( 0: m_add2 ... ) but this is simpler. ubn_add :: UBigfloat -> UBigfloat -> UBigfloat ubn_add = ubn_add2 -- Older, simpler implementation of ubn_add. {- ubn_add1 :: UBigfloat -> UBigfloat -> UBigfloat ubn_add1 a b = (ae+1, m_addcarry0 am bm) where ((ae,am),(be,bm)) = ubn_same_exp (a, b) -} -- ubn_add2 instead takes advantage of unequal exponents in the operands to -- avoid putting a string of 0s on the head of one parameter and then adding -- them in uselessly. To be able to do this -- 1) the exponents must differ by at least two -- 2) there must be a non-9 digit in the first, non-overlapping digits of the -- bigger operand, to prevent possible carry due to the second from affecting -- the first digit of the result. ubn_add2 :: UBigfloat -> UBigfloat -> UBigfloat ubn_add2 (ae,am) (be,bm) -- First, make sure ae >= be | ae == be = (ae+1, m_addcarry0 am bm) | ae > be = ubn_add2' (ae,am) (be,bm) | otherwise = ubn_add2' (be,bm) (ae,am) -- ubn_add2' knows that ae > be ubn_add2' :: UBigfloat -> UBigfloat -> UBigfloat ubn_add2' (ae,am) (be,bm) -- First case: difference is just one | ae == be + 1 = ubn_add2'a (ae,am) (be,bm) -- Second case (the fun one): diff >= 2 | ae > be + 1 = ubn_add2'b (ae,am) (be,bm) -- ubn_add2'a: ae = be + 1 -- It is (just barely) worth separating these two cases, as far as -- speed is concerned, instead of always using m_addcarry0. ubn_add2'a :: UBigfloat -> UBigfloat -> UBigfloat ubn_add2'a (ae,(a:x)) (be,bm) -- a) exponent differs by 1 and first digit < 9: -- there can be no overflow from digit 1. | a < base_1 = (ae, m_addcarry (a:x) (0:bm)) -- b) exponent differs by 1 and first digit of a is 9: -- there may be an overflow from the first digit. | otherwise = (ae+1, m_addcarry0 (a:x) (0:bm)) -- ubn_add2'b: ae >= be + 2 -- Make sure there can be no overflow from the result -- and pass the work on to the lookahead function. ubn_add2'b :: UBigfloat -> UBigfloat -> UBigfloat -- i) a = 0 ubn_add2'b (ae,[]) (be,bm) = (be,bm) -- ii) second digit of a is (implied) 0 ubn_add2'b (ae,[x]) (be,bm) = (ae, x:0:(prefix0s (ae-be-2) bm)) -- iii) general case ubn_add2'b (ae,am) (be,bm) -- There can be no carry into the first digit if the difference -- between the exponents is >= 2 and any of the intervening digits -- are less than (base-1), 'cos they are sure to absorb any carry -- of 1 from the first pair of overlapping digits. | any_lt_base_1 am (ae-be) = (ae, ubn_add2'' (ae,am) (be,bm)) -- All the protruding digits of a are 9! This is hard to believe :-) | otherwise = (ae+1, m_addcarry0 am bm) -- ubn_add2'' knows that ae-be >= 2 and that there can be no overflow -- from the first digit of the result. Since its result always has -- the same exponent as its first parameter, we just return the mantissa. ubn_add2'' :: UBigfloat -> UBigfloat -> Mantissa ubn_add2'' (ae,[]) (be,bm) = prefix0s (ae-be) bm ubn_add2'' (ae,[a]) (be,bm) = a:0:(prefix0s (ae-be-2) bm) ubn_add2'' (ae,am) (be,bm) | ae == be + 2 = ubn_add2''a (ae,am) (be,bm) -- End of recursion | otherwise = ubn_add2''b (ae,am) (be,bm) -- This one may recurse -- ae = be+2: End of recursion ubn_add2''a (ae,(p:q:x)) (be,bm) | q < base_1 = p : m_addcarry (q:x) (0:bm) -- WIN! | otherwise = m_addcarry (p:q:x) (0:0:bm) -- q == base_1 :-( -- ae > be+2: Can recurse ubn_add2''b (ae,(a:x)) (be,bm) | any_lt_base_1 x (ae-be-1) = a : (ubn_add2'' (ae-1,x) (be,bm)) -- WIN AND PLAY AGAIN! | otherwise = m_addcarry (a:x) (prefix0s (ae-be) bm) -- Boring case -- Are any of the first n digits of a mantissa less than base-1 ? -- (if there are less than n digits in the mantissa, then yes, -- because the omitted digits are implicitly 0) any_lt_base_1 x n = or (map (< base_1) (take n x)) -- ubn_sub performs subtraction of 2 unsigned Bigfloats. -- It's up to the caller to ensure that a >= b. -- This uses the version of subcarry that notices when arg2 is shorter than -- arg1, so that "n `bn_sub` bn_1" costs almost nothing instead of running a -- pointless carry down an infinite list. ubn_sub :: UBigfloat -> UBigfloat -> UBigfloat ubn_sub = ubn_sub1 ubn_sub1 :: UBigfloat -> UBigfloat -> UBigfloat ubn_sub1 a@(ae,am) b@(be,bm) | sanity_checks && (ae < be) = error "Internal error in BigFloat.ubn_sub1: exponents are the wrong way round" | otherwise = (ae, ubn_sub' a b) -- ubn_sub2 is an alternative implementation that takes advantage of ae > be -- and returns the initial part of am unaltered rather than subtracting a -- string of zeroes from it. -- As in ubn_add2, we can return the first digit of am without inspecting -- bm if ae is at least 2 greater than be, and there is at least one -- non-zero digit between the first digit of a and the start of b (for a -- possible borrow to happen from). -- Some calculations turn out .5% faster than with ubn_sub1, others are .5% slower. ubn_sub2 :: UBigfloat -> UBigfloat -> UBigfloat ubn_sub2 a@(ae,am) b@(be,bm) | sanity_checks && (ae < be) = error "Internal error in BigFloat.ubn_sub2: exponents are the wrong way round" | otherwise = (ae, ubn_sub2' a b) ubn_sub2' :: UBigfloat -> UBigfloat -> Mantissa ubn_sub2' a@(ae,am) b@(be,bm) | ae-be >= 2 && (not $ m_is_zero $ drop 1 $ take (ae-be) am) = head am : ubn_sub2' (ae-1, tail am) b | otherwise = ubn_sub' a b -- ubn_sub' just equalises the exponents and subtracts the mantissae -- It is guaranteed that a >= b, therefore ae >= be, so the (unnormalised) -- result's exponent will always be ae, and this function just returns the -- mantissa. ubn_sub' :: UBigfloat -> UBigfloat -> Mantissa ubn_sub' (ae,am) (be,bm) = m_subcarry2 am (prefix0s (ae-be) bm) -- m_addcarry performs addition of 2 mantissas and carry simultaneously. -- As usual, it must be guaranteed that there will be no carry from the -- first digit. Used in add_skew_carry. -- It only applies carry1 to the part resulting from an addition. m_addcarry :: Mantissa -> Mantissa -> Mantissa m_addcarry a b = m_carry1a tocarry ++ nocarry where (tocarry,nocarry) = m_add2 a b -- m_addcarry0 is like m_addcarry except that it prefixes a zero to the -- result before performing the carry (thereby guaranteeing no overflow) m_addcarry0 :: Mantissa -> Mantissa -> Mantissa m_addcarry0 a b = m_carry1a (0:tocarry) ++ nocarry where (tocarry,nocarry) = m_add2 a b -- m_subcarry subtracts b from a and propagates carry (ok, borrow) m_subcarry :: Mantissa -> Mantissa -> Mantissa m_subcarry a b = m_borrow1 (m_sub a b) -- If the second operand is shorter that the first, the final part -- of the first is returned "as is" since there's no need to do any -- borrowing from it. m_subcarry2 a b = (m_borrow1a toborrow) ++ noborrow where (toborrow,noborrow) = m_sub2 a b -- "m_add" adds two mantissas without performing carry across the result -- leaving each digit with value from 0 to base*2 m_add :: Mantissa -> Mantissa -> Mantissa m_add (a:x) (b:y) = (a+b):m_add x y m_add a [] = a m_add [] b = b -- m_add2 returns the mantissa in two parts: the first part is the one -- resulting from addition (to which carry must be applied); the second -- part is the part resulting from unequal lengths of a and b. m_add2 :: Mantissa -> Mantissa -> (Mantissa,Mantissa) m_add2 (a:x) (b:y) = ((a+b):tocarry, nocarry) where (tocarry,nocarry) = m_add2 x y m_add2 a [] = ([],a) m_add2 [] b = ([],b) -- "m_sub" subtracts one mantissa from another without performing borrow, -- leaving each digit with value from -(base-1) to +(base-1). m_sub :: Mantissa -> Mantissa -> Mantissa m_sub (a:x) (b:y) = (a-b):m_sub x y m_sub a [] = a m_sub [] b = map negate b -- m_sub2, like m_add2, returns the part that borrowing needs to be performed -- on and the part it doesn't. -- Unlike add2, we only win if the first param is longer than the second. m_sub2 :: Mantissa -> Mantissa -> (Mantissa,Mantissa) m_sub2 (a:x) (b:y) = ((a-b):tocarry, nocarry) where (tocarry,nocarry) = m_sub2 x y m_sub2 a [] = ([], a) m_sub2 [] b = (map negate b, []) -- m_carry1 performs carry throughout a mantissa. -- Each term must be <= (base-1)*2. -- This means that the maximum carry from any digit to the next will be 1, -- To determine the first digit, it inspects at least the first two elements. -- It assumes that the digits are sufficient to hold the result, ie that -- there will be no carry from the first digit. m_carry1 :: [Int] -> [Int] m_carry1 (a:b:x) | b < base_1 = a : m_carry1 (b : x) | b > base_1 = (a + 1) : m_carry1 ((b - base) : x) | otherwise = (a + carry): m_carry1 ((b-carry*base):x) where carry = m_carry1from (b:x) m_carry1 [a] | a > 0 = [a] | otherwise = [] m_carry1 [] = [] -- m_carry1a is like m_carry1 except that it never removes trailing zeroes m_carry1a :: [Int] -> [Int] m_carry1a (a:b:x) | b < base_1 = a : m_carry1a (b : x) | b > base_1 = (a + 1) : m_carry1a ((b - base) : x) | otherwise = (a + carry): m_carry1a ((b-carry*base):x) where carry = m_carry1from (b:x) m_carry1a x = x -- cases [a] and [] -- m_carry1from returns the carry from a list of digits -- where each digit is <= (base-1)*2 -- To determine the carry from a list, it inspects at least the first element. m_carry1from :: [Int] -> Int m_carry1from (a:x) | a < base_1 = 0 | a > base_1 = 1 | otherwise = m_carry1from x m_carry1from [] = 0 -- m_borrow1 takes a list of digits from -9 to +9 and performs the borrowing -- necessary to bring them all in the range 0-9. The first term of the list -- must not be negative, and must be sufficient to provide any borrowing -- required of it by the following digit(s). m_borrow1 :: [Int] -> [Int] m_borrow1 (a:b:x) | b > 0 = a : m_borrow1 (b:x) | b < 0 = (a-1) : m_borrow1 ((b+base):x) | otherwise = (a-borrow) : m_borrow1 ((b + borrow * base) : x) where borrow = m_borrow1from (b:x) m_borrow1 [a] | sanity_checks && a < 0 = error "Internal error: Negative digit in m_borrow1" | a == 0 = [] | otherwise = [a] m_borrow1 [] = [] -- m_borrow1a is like m_borrow1 except that it doesn't strip trailing zeroes. m_borrow1a :: [Int] -> [Int] m_borrow1a (a:b:x) | b > 0 = a : m_borrow1a (b:x) | b < 0 = (a-1) : m_borrow1a ((b+base):x) | otherwise = (a-borrow) : m_borrow1a ((b + borrow * base) : x) where borrow = m_borrow1from (b:x) m_borrow1a [a] = [a] m_borrow1a [] = [] m_borrow1from (a:x) | a > 0 = 0 | a < 0 = 1 | otherwise = m_borrow1from x m_borrow1from [] = 0 -- Simple doubler. It gains in lookahead because we know that, when you double -- a number in an even number base, the value of input digit n+2 cannot affect -- the result digit n. -- This relies on base being an even number. bn_twice (s,e,m) | m_is_zero m = bn_0 -- reinstate this if u want odd number bases -- | odd base = bn_add (s,e,m) (s,e,m) | head m < half_base = (s, e, m_twice m) | otherwise = (s, e+1, m_twice (0:m)) -- m_twice is always applied to a list whose first element is < base/2 m_twice :: [Int] -> [Int] m_twice (a:b:x) | b < half_base = (a+a) : m_twice (b:x) | otherwise = (a+a+1) : m_twice ((b - half_base) : x) m_twice [0] = [] -- Final digit 5 (ok, base/2) generates a trailing 0. m_twice [a] = [a+a] m_twice [] = [] bn_half (s,e,m) = bn_normalise (s, e, m_half m) -- When base can be odd, you need the following instead -- | even base = bn_normalise (s, e, m_half m) -- | otherwise = bn_quotI (s,e,m) 2 m_half (a:x) | even a = (a `div` 2) : m_half x | otherwise = (a `div` 2) : (m_add [half_base] (m_half x)) m_half [] = [] half_base = base `div` 2 -- -- --------- Simple comparison functions ---------- -- -- Is a Bigfloat negative? bn_is_neg (s,e,m) = (s < 0) -- Is it zero? bn_is_zero (s,e,m) = m_is_zero m -- Is it a whole number? Yes, if there are no significant digits after the . bn_is_integer (s,e,m) = m_is_zero (drop e m) bn_cmp a b = compare (MakeBigFloat a) (MakeBigFloat b) bn_eq a b = MakeBigFloat a == MakeBigFloat b bn_ne a b = MakeBigFloat a /= MakeBigFloat b bn_lt a b = MakeBigFloat a < MakeBigFloat b bn_le a b = MakeBigFloat a <= MakeBigFloat b bn_gt a b = MakeBigFloat a > MakeBigFloat b bn_ge a b = MakeBigFloat a >= MakeBigFloat b -- Unsigned Bigfloat comparison; result as per bn_cmp ubn_cmp :: UBigfloat -> UBigfloat -> Ordering ubn_cmp a b = compare (MakeUBigFloat a) (MakeUBigFloat b) ubn_eq, ubn_ne, ubn_lt, ubn_le, ubn_gt, ubn_ge :: UBigfloat -> UBigfloat -> Bool ubn_eq a b = MakeUBigFloat a == MakeUBigFloat b ubn_ne a b = MakeUBigFloat a /= MakeUBigFloat b ubn_lt a b = MakeUBigFloat a < MakeUBigFloat b ubn_le a b = MakeUBigFloat a <= MakeUBigFloat b ubn_gt a b = MakeUBigFloat a > MakeUBigFloat b ubn_ge a b = MakeUBigFloat a >= MakeUBigFloat b -- -- ------------- Simple multiply by an integer ------------- -- -- bn_times multiplies Bigfloat by an Integer, done by repeated addition. -- The *first* parameter is the Integer, which seems more convenient for -- partial parameterisation purposes. -- bn_add normalises the result for us. bn_times n b | n > 0 = bn_times' n b | n < 0 = bn_neg (bn_times' (-n) b) | n == 0 = bn_0 -- reduced case: n always >= 1 -- Optimisation: (2n) x b = (n x b) + (n x b) -- Optimisation to do: multiply all digits by n, then do a super-carry. bn_times' n b | even n = bn_twice (bn_times' (n `div` 2) b) -- n is odd: the (n div 2) here is really ((n-1) div 2). | n > 2 = b `bn_add` (bn_twice (bn_times' (n `div` 2) b)) | n == 1 = b -- Same stuff for mantissae. Caller must guarantee no overflow, and -- maximum n is base-1. m_times :: Int -> Mantissa -> Mantissa m_times n m | sanity_checks && (n < 0 || n >= base) = error("m_times " ++ show n) | otherwise = m_carry2 (map (* n) m) -- -- ------------- Long multiplication ------------- -- -- bn_mul: Multiply two Bigfloats. bn_sqr is a modified version of this. -- -- As an alternative, cleaner and maybe faster, it may be worth using -- Karatsuba's divide-and-conquer algorithm: -- -- Split U in two pieces, U1 and U0, such that -- U = U0 + U1*(B^n), -- and V in V1 and V0, such that -- V = V0 + V1*(B^n). -- -- UV is then computed recursively using the identity -- -- UV = ( B^2n + B^n ) U1 V1 + B^n (U1-U0) (V0-V1) + (B^n + 1) U0 V0 bn_mul (as,ae,am) (bs,be,bm) = bn_normalise (as*bs, e, m) where (e,m) = ubn_mul (ae,am) (be,bm) ubn_mul :: UBigfloat -> UBigfloat -> UBigfloat ubn_mul (ae,am) (be,bm) = (ae+be, m_mul am bm) -- m_mul multiplies two mantissas. -- -- Strategy: -- sum the columns of the skewed cross product and then do an intelligent -- super-carry on the result. The super-carry works by applying a one-digit -- carry to the list of column sums until we can be sure that the maximum -- value of each column is (base * 2). Then it applies m_carry. -- [Alternate version to try: apply one-digit carry until maximum column value -- is (base - 1) ^ 2, then apply m_carry2.] -- -- We do this by giving two lists to the carry-preparation function: -- - the list of column sums -- - the list of maximum values for the column sums. -- and applying the one-digit carry to both lists until the maximum value -- of the leading digits is within the required maximum. m_mul :: [Int]->[Int]->[Int] m_mul a b = supercarry (0:(partial_products a b)) max_mul_vals supercarry list maxvals = m_carry1 (reduce_for_carry list maxvals) partial_products :: [Int]->[Int]->[Int] partial_products a b = add_skew (cross_product a b) -- cross_product: -- given [A,B,C] [a,b,c] -- calculate partial result -- [ [Aa,Ab,Ac], -- [Ba,Bb,Bc], -- [Ca,Cb,Cc] ] cross_product :: [Int]->[Int]->[[Int]] cross_product a b = map (ma_times b) a where ma_times l n = map (* n) l -- add_skew: Return the sum of a potentially infinite number of mantissas, -- offsetting the rows by one place: -- given [ [Aa,Ab,Ac], -- [Ba,Bb,Bc], -- [Ca,Cb,Cc] ] -- return [ [Aa,Ab,Ac], -- + [Ba,Bb,Bc], -- + [Ca,Cb,Cc] ] -- ================ -- [r1,r2,r3,r4,r5] add_skew :: [[Int]] -> [Int] add_skew ((a:x):y) = a : m_add x (add_skew y) -- add_skew ([]:y) = 0 : add_skew y -- never happens (?) -- Sum of a single row is just that row. add_skew [x] = x -- We get handed an empty list when they multiply by bn_0. add_skew [] = [] -- max_mul_vals gives the maximum possible value of a partial product term. -- max_mul_vals = 0 : partial_products nines nines where nines = [9,9..] -- Turns out to be equal to: max_mul_vals = [0, base_1 ^ 2 .. ] -- reduce_for_carry takes a list of column sums and a list of the maximum -- possible values of each column sum, and applies a one-digit carry until -- the terms it returns are sure to be <= (base-1) * 2 -- The second list is always infinite. -- In practice, repeated applications of carry1digit to max_mul_vals do this: -- [0,81,162,243,324,405,486,567,648,729,810,891,972,1053,1134,1215,1296,1377.. -- [8,17,26,35,44,53,62,71,80,90,89,98,107,116,125,134,143,152,161,171,170,179.. -- [9,9,9,9,9,9,9,9,9,8,18,18,18,18,18,18,18,18,18,18,17,27,27,27,27,27,27,27.. -- [9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9.. -- and the ..9,9,9,8,18,18,18.. pattern always repeats itself further down the -- list in successive iterations. -- Thus we can tell when the current digit is gonna be ok by looking for when -- the following digit's maximum possible value exceeds 9 (the 8,17 or 8,18 -- point); at this point a new invocation of carry1digit is necessary from the -- "8" digit onwards. -- Note! Although this reduces all terms in the 2nd parameter to the range -- 0-base-1, it *does not* necessarily reduce all terms of the first parameter -- to this range. For pathological cases like -- bn_sqr (bn "999999999999999999999999999999999999999999"), -- there is a "..9,9,9,9,9,9,10.." near the end that would require many -- iterations of carry1digit to propagate the carry back to its proper place. -- Only the final application of carry1 can resolve this (see m_mul above). reduce_for_carry :: [Int] -> [Int] -> [Int] reduce_for_carry (a:b:x) (c:d:y) | d < base = a : reduce_for_carry (b:x) (d:y) | otherwise = reduce_for_carry (carry1digit (a:b:x)) (carry1digit (c:d:y)) --reduce_for_carry [a] y = [a] --reduce_for_carry [] y = [] reduce_for_carry x y = x -- synthesis of the above two cases -- carry1digit just performs the carry from each digit into the previous -- one without progating further. carry1digit (a:b:x) = (a + b `div` base) : carry1digit (b `mod` base : x) --carry1digit [a] = [a] --carry1digit [] = [] carry1digit x = x -- Perform carry throughout number. Each term must be <= (base-1)^2. -- This ensures that the maximum carry from any digit to the next will be -- (base-2) because: -- max_digit = (base - 1) * (base - 1) = (base^2 - 2 * base + 1) -- max_carry = (max_digit + max_carry) `div` base -- If max_carry is (base - 2) then (max_digit + max_carry) `div` base = -- ( (base^2 - 2 * base + 1) + (base - 2) ) `div` base = -- ( (base^2 - 2 * base + (1 + (base - 2) ) `div` base = -- ( (base^2 - 2 * base + (base - 1) ) `div` base = -- ( (base - 2) + 0 ) = base - 2 -- Example of highest carry: -- [0 81 81 81] -> [0 81 89 1] -> [0 89 9 1] -> [8 9 9 1] m_carry2 :: [Int] -> [Int] m_carry2 (a:b:c:x) -- Max carry from b is (base - 2), so if b `mod` base < 2, the carry -- from x cannot affect a. This turns out to be a win so rarely that -- the extra checking slows us down more than we gain. -- = (a + b `div` base) : m_carry2 ((b `mod` base) : c : x), -- if b `mod` base < 2 -- Carry from x into c can only affect carry from c into b by 1, -- so if (b + (c `div` 10)) `mod` 10 is less than 9, -- x won't affect the carry from b into a. | bc `mod` base < base_1 = (a + bc `div` base) : m_carry2 (bc `mod` base : c `mod` base : x) -- general case | otherwise = (a + carry) : m_carry2 ((b - carry * base) : c : x) where bc = b + c `div` base carry = m_carry2from (b:c:x) m_carry2 [a,b] | b `mod` base > 0 = [a + b `div` base, b `mod` base] | otherwise = [a + b `div` base] m_carry2 [a] | a > 0 = [a] | otherwise = [] m_carry2 [] = [] -- Return the carry from the list of numbers where each number <= (base-1)^2 -- Maximum carry from any digit into the previous one is (base-2). -- To find the carry from A:B:C... in base 10: -- If A `mod` base < 2 then carry = a `div` base -- (carry from B cannot affect carry from A) -- The carry from x into b can only affect the carry from b into a by one. -- So if (a + (b `div` 10)) `mod` 10 is less than 9 -- then x won't affect the carry from a. m_carry2from :: [Int] -> Int m_carry2from (a:b:x) -- This first case gains so infrequently that it ends up a net loss. -- | if a `mod` base < 2 = a `div` base | ab `mod` base < base_1 = ab `div` base | otherwise = (a + m_carry2from(b:x)) `div` base where ab = a + (b `div` base) m_carry2from [a] = a `div` base m_carry2from [] = 0 -- ---------- end of long multiplication stuff ---------- -- bn_raise - simple power function by repeated multiplication -- a is the Bigfloat to raise to a power; -- x is the Integer power to raise it to. -- Optimisation: bn_raise x (2*n) = sqr (bn_raise x n) bn_raise a x | x == 0 && bn_is_zero a = error "Raising 0 to the 0th power is undefined" | bn_is_zero a = bn_0 | x == 0 = bn_1 | x < 0 = bn_div bn_1 (bn_raise' a (-x)) | otherwise = bn_raise' a x bn_raise' a x | even x = bn_sqr (bn_raise' a (x `div` 2)) | (x > 2) = bn_mul a (bn_raise' a (x-1)) | x == 1 = a -- bn_sqr - square a Bigfloat, used frequently in the scientific functions, -- hence optimised here to avoid pointlessly doing certain things twice. -- was: bn_sqr x = x `bn_mul` x bn_sqr (as,ae,am) = bn_normalise (1, e, m) where (e,m) = ubn_sqr (ae,am) ubn_sqr (ae,am) = (ae+ae, m_sqr am) m_sqr :: [Int]->[Int] m_sqr a = supercarry (0:(add_skew2 (square_product a))) max_mul_vals -- square_product is like cross_product except that it folds the lower left -- triangle of partial results onto the upper right triangle (since the -- cross product of a single number is symmetrical about its main diagonal). -- Given [a,b,c] -- calculate partial result -- [ [aa,2ab,2ac,..], -- [bb,2bc,2bd,..], -- [cc,2cd,2ce,..] ] square_product :: [Int]->[[Int]] square_product (a:x) = (a*a : map (* (a+a)) x) : square_product x square_product [] = [] -- add_skew2: Return the sum of a potentially infinite number of mantissas, -- offsetting the rows by two places: -- given [ [aa,2ab,2ac,..], -- [bb,2bc,2bd,..], -- [cc,2cd,2ce,..] ] -- return [ [aa,2ab,2ac,..], -- + [bb,2bc,2bd,..], -- + [cc,2cd,2ce,..] ] -- ================ -- [r1,r2,r3,r4,r5] add_skew2 :: [[Int]] -> [Int] add_skew2 ((a:b:x):y) = a : b : m_add x (add_skew2 y) -- Sum of a single row is just that row. add_skew2 [x] = x -- This happens when squaring numbers with an odd number of digits add_skew2 ([a]:y) = a : 0 : add_skew2 y --add_skew2 ([]:y) = 0 : 0 : add_skew2 y -- never happens (?) -- We get handed an empty list when they square bn_0. add_skew2 [] = [] -- bn_quot: Long division of a Bigfloat by an integer. -- Uses the built-in infinitely-long integers to do all the hard work. -- The time it takes is pretty much independent of the value of q. -- -- bn_quot is the generic divide-by-Integer version. -- bn_quotI is a faster version for Haskell, dividing by an Int, -- used in atan_series. The maximum remainder, stored in an Int, -- is less than (10 * dividend), so the maximum Int you can use it -- to divide by is MaxInt/10. -- It came about so that I could calculate pi with nhc98-1.16 -- even though its divide-by-Integer code gives wrong answers and -- segmentation faults. -- However, it is a win in any case, and also helps with some -- other series-based calculations. bn_quot x y = if integer_division_is_broken then bn_div x (i2bn y) else bn_quot' m_quot x y bn_quotI = bn_quot' m_quotI -- Perform initial sanity checks and conversions, -- then call the appropriate real-work function. bn_quot' :: Integral a => ([Int] -> a -> [Int]) -> Bigfloat -> a -> Bigfloat bn_quot' m_quot_version (s,e,m) q | q == 0 = error "bn_quot cannot divide by 0" | m_is_zero m = bn_0 | q == 1 = (s,e,m) | q == -1 = (-s,e,m) | q > 1 = bn_normalise (s, e, m_quot_version m q) | q < -1 = bn_normalise (-s, e, m_quot_version m (-q)) -- Original Miranda version, now using machine-word Ints m_quotI :: [Int] -> Int -> [Int] m_quotI m q | m_is_zero m = [] | otherwise = (a `div` q) : m_quotI (shiftin ((a `mod` q):x)) q where (a:x) = m shiftin (a:b:x) = (a * base + b) : x shiftin [a] = [a * base] -- Divide a mantissa by an Integer m_quot :: [Int] -> Integer -> [Int] m_quot m q | sanity_checks && q == 0 = error "Internal error: BigFloat.m_quot was told to divide by zero" | sanity_checks && q < 0 = error "Internal error: BigFloat.m_quot was told to divide by a negative number" | m_is_zero m = [] | q == 1 = m | integer_division_is_broken = m_quot2 m q | otherwise -- glue for the Integer-and-rest-of-mantissa version. = m_quot1 (toInteger (head m), tail m) q -- m_quot1 keeps the head of the remainder still to be divided away -- in a separate Integer rather than in the list. m_quot1 :: (Integer,[Int]) -> Integer -> [Int] m_quot1 (a,x) q | a == 0 && m_is_zero x = [] -- remainder==0 => division is exact | otherwise = fromInteger (a `div` q) : m_quot1 (shiftin ((a `mod` q),x)) q where shiftin (b,(c:y)) = ((b * toInteger base + toInteger c), y) shiftin (b,[]) = (b * toInteger base, []) -- Avoid Integer division by using bn_div to get 1/q, -- convert that to a pure mantissa and then multiply by it. -- (q will be > 1) m_quot2 n q = n `m_mul` m_recip_q where m_recip_q = prefix0s (-e) m (s,e,m) = bn_1 `bn_div` (i2bn q) -- -- ---------------- Long division --------------------- -- bn_div (as,ae,am) (bs,be,bm) = bn_normalise (as*bs, ae-be+1, m_div am bm) ubn_div (ae,am) (be,bm) = (ae-be+1, m_div am bm) -- m_div takes two mantissae and divides the first by the second, returning a -- mantissa. -- div_loop requires divisor <= dividend -- - zero divided by anything is 0. -- - div_loop returns the first digit of the result and the remainder of the -- dividend after the divisor has been subtracted from it "dig" times. -- "dend" is the dividend; "sor" is the divisor. -- -- mdiv' seems to be necessary in Haskell cos I don't know how to get a where -- to apply to all 4 cases (nhc98 complains about needing a scope decl) m_div = m_div2 m_div1 dend sor = m_div' dend sor (m_div_loop (0, dend, sor)) m_div' dend sor (dig, dend2, sor2) | m_is_zero dend = [] -- The following line is logically ok but hits so rarely that -- it slows things down more than it speeds them up. -- | if sor `m_gt` dend = 0: m_div dend (0:sor) | dend2 == [] = [dig] -- exact | head dend2 == 0 = dig : m_div (tail dend2) sor | otherwise = dig : m_div dend2 (0:sor) {- was: (works with hugs, not with nhc, which says: ====== Errors after type inference/checking: Context for Prelude.Num needed in left hand pattern. ) m_div dend sor | m_is_zero dend = [] -- The following line is logically ok but hits so rarely that -- it slows things down more than it speeds them up. -- | if sor `m_gt` dend = 0: m_div dend (0:sor) | dend2 == [] = [dig] -- exact | head dend2 == 0 = dig : m_div (tail dend2) sor | otherwise = dig : m_div dend2 (0:sor) where (dig, dend2, sor2) = m_div_loop (0, dend, sor) -} -- Work out how many times the divisor can be subtracted from the dividend -- without the dividend going negative. -- -- "hack" is our conservative first approximation for the first digit, formed -- by dividing the first two digits of dend and sor. -- The slow m_div_loop is then used to arrive at the exact result. -- The >2-2-0-1 order is empirically the fastest using bn_pi as the test piece. m_div_loop (0, dend, sor) | hack > 2 = m_div_loop' (hack, dend `m_subcarry` (m_times hack sor), sor) | hack == 2 = m_div_loop' (2, dend `m_subcarry` (m_twice sor), sor) | hack == 0 = m_div_loop' (0, dend, sor) | hack == 1 = m_div_loop' (1, dend `m_subcarry` sor, sor) where -- The "+ 1" makes sure we don't exceed in our estimate. hack = (take2 dend) `div` (take2 sor + 1) take2 (a:b:x) = a * base + b take2 [a] = a * base take2 [] = 0 -- This is the regular version, done by repeated subtraction. -- Using m_subcarry2 here and above turns out to be >10% slower that m_subcarry -- in this context, probably because the optimisation never comes into effect. m_div_loop' x = until div_loop_done div_loop_body x where div_loop_done (dig, dend, sor) = dend `m_lt` sor div_loop_body (dig, dend, sor) = (dig+1, m_subcarry dend sor, sor) -- WORK IN PROGRESS ON NEW DIVIDE ROUTINE, AVOIDING RELIANCE ON BUILT-IN BIGNUMS m_div2 n d | sanity_checks && m_is_zero d = error "Division by zero" | m_is_zero n = [] | otherwise = m_div_considering_1 n d m_div_considering_1 :: Mantissa -> Mantissa -> Mantissa {- If we can be sure of the first result digit by examining just the first digit of each parameter, then do so. If not, examine more digits. Here's the table of what the first result digit will be, given the first one of the numerator and denominator (n/d) for base 10. '?' means we cannot be sure. n:0 1 2 3 4 5 6 7 8 9 d:0 ? ? ? ? ? ? ? ? ? ? 1 0 ? ? ? ? ? ? ? ? ? 2 0 0 ? 1 ? ? ? ? ? ? 3 0 0 0 ? 1 1 ? ? 2 ? 4 0 0 0 0 ? 1 1 1 ? ? 5 0 0 0 0 0 ? 1 1 1 1 6 0 0 0 0 0 0 ? 1 1 1 7 0 0 0 0 0 0 0 ? 1 1 8 0 0 0 0 0 0 0 0 ? 1 9 0 0 0 0 0 0 0 0 0 ? The '0' and '1' cases are easy enough to isolate; the '2' at 8/3 is more interesting. In higher number bases, more of these cases appear as thin wedges of the same value, of which the '2' here is really just the tip. We can spot these in general by seeing whether (n.999999 / d.000000) == (d.999999 / n.000000) or, in integer maths, whether n / d == n / (d+1). -} -- if the numerator is 0, then the result has come out exact. m_div_considering_1 [] d | sanity_checks && m_is_zero d = error "Internal error in BigFloat.m_div_considering_1 []: Divide by zero" | otherwise = [] m_div_considering_1 n@(na:nx) d@(da:dx) | sanity_checks && m_is_zero d = error "Internal error in BigFloat.m_div_considering_1: Divide by zero" -- if the numerator is 0, then the result has come out exact. | m_is_zero n = [] -- reduce the case of both n & d starting with a 0 digit | na == 0 && da == 0 = m_div_considering_1 nx dx -- The 0 zone of the above figure. Move on to next digit, -- optimising to avoid passing 0:x 0:y. | na < da = 0 : if na == 0 then m_div nx d else m_div n (0:d) -- The diagonal and the top row (avoiding na `div` 0) -- These *could* be included in the "otherwise" case -- but are quick so what the heck. | na == da || da == 0 = m_div_considering 2 n d -- The 1 zone. In fact we deal with all '1' cases regardless of -- parameter length, since m_lt is optimally lazy anyway. | m_half n `m_lt` d = 1 : m_div (n `m_subcarry` d) (0:d) -- Sure cases in the ? ? ? zone | na_div_da == na `div` (da+1) = na_div_da : m_div (n `m_subcarry` (m_times na_div_da d)) (0:d) -- Unsure cases in the ? ? ? zone | otherwise = m_div_considering 2 n d where na_div_da = na `div` da -- No, I don't trust compilers! -- More general case tries to work out the first result digit from the -- first "ndigits" digits of each operand. Much like the above, except that -- terminal and pathological cases have been dealt with. m_div_considering :: Int -> Mantissa -> Mantissa -> Mantissa m_div_considering ndigits n@(na:nx) d@(da:dx) -- If n is less than d but we limit ourselves here to an n-digit -- comparison, we will only recurse with ndigits+1 anyway if it fails -- here for lack of digits, so go straight for full-length comparison. | n `m_lt` d = 0 : if na == 0 then m_div nx d else m_div n (0:d) -- If all the digits of d are comprised in "ndigits" we can always be -- sure of our first digit. | drop ndigits d == [] = sure -- n == d || d == 0 | m_eq (take ndigits n) (take ndigits d) || m_is_zero (take ndigits d) = unsure -- The 1 zone is dealt with completely above. -- Sure cases in the ? ? ? zone | n_div_d == m_divn ndigits n (m_addcarry d (replicate (ndigits-1) 0 ++ [1])) = sure -- Unsure cases in the ? ? ? zone | otherwise = unsure where -- Yes, we are sure of the first digit sure = n_div_d : m_div (n `m_subcarry` (m_times n_div_d d)) (0:d) -- No, we are not sure of the first digit unsure = m_div_considering (ndigits+1) n d n_div_d = m_divn ndigits n d -- Return the first n digits of x divided by the first n digits of y. -- The result will be a single digit from 0 to base-1. -- For now we are lazy and use built-in div function m_divn :: Int -> Mantissa -> Mantissa -> Int m_divn ndigits x y -- Integers go up to 2^31 = 2,147,XXX,XXX | ndigits > 9 = error "Internal error: Integer overflow in m_divn" | otherwise = div (m2Int (take ndigits (x ++ repeat 0))) (m2Int (take ndigits (y ++ repeat 0))) -- END OF WORK IN PROGRESS ON NEW DIVIDE ROUTINE -- -- ---------------- Square root --------------------- -- -- Square root adapted from Mr C. Woo's algorithm for abacus: -- convert number in base 10 to base 100 by pairing digits, -- then each application of sqrt_loop generates a digit of the result. -- For details of the algorithm, see under http://freaknet.org/martin/tape bn_sqrt (s,e,m) | s_is_negative s = error "bn_sqrt of negative number" | otherwise = bn_normalise $ ubn2bn $ ubn_sqrt (e,m) ubn_sqrt (e,m) | even e = (e `div` 2, m_sqrt m) | odd e = ((e+1) `div` 2, m_sqrt (0:m)) -- m_sqrt gives the square root of a mantissa, which must be displaced by an -- even number of digits from the decimal point to give the correct answer. m_sqrt :: [Int] -> [Int] m_sqrt = if integer_division_is_broken then m_sqrt2 else m_sqrt1 m_sqrt1 :: [Int] -> [Int] m_sqrt1 m | m_is_zero m = [] -- sqrt(0) | otherwise = m_sqrt' (0, (toInteger a,x)) where (a:x) = m_sqrbase m -- prime sqrt engine with first digit -- Convert mantissa in base n to mantissa in base n^2. m_sqrbase :: [Int] -> [Int] m_sqrbase (a:b:m) = (a * base + b) : m_sqrbase m m_sqrbase [a] = [a * base] m_sqrbase [] = [] m_sqrt' :: (Integer,(Integer,[Int])) -> [Int] m_sqrt' (r1,(hd_o1,tl_o1)) | hd_o1 == 0 && m_is_zero tl_o1 = [] | otherwise = fromInteger (r2 `mod` ibase) : m_sqrt' (r2 * ibase, m_sqrt_shiftin (hd_o2,tl_o1)) where (r2,hd_o2) = m_sqrt_loop (r1,hd_o1) ibase = toInteger base -- Shift another term of the operand into the calculation. m_sqrt_shiftin :: (Integer, [Int]) -> (Integer, [Int]) m_sqrt_shiftin (a,b:x) = (a * sqr_base + toInteger b, x) m_sqrt_shiftin (a,[]) = (a * sqr_base, []) sqr_base = toInteger (base * base) m_sqrt_loop :: (Integer,Integer) -> (Integer,Integer) m_sqrt_loop ro = until sqrt_loop_done sqrt_loop_body ro where sqrt_loop_done (r, o) = 2*r >= o -- was: 2*r+1 > o sqrt_loop_body (r, o) = (r+1, o - (2*r+1)) -- Second version avoids using Integer mod function -- (which is broken in nhc98-1.16) -- TODO: This is slow cos "+ one" twice in the loop body -- addcarries loads of zeroes followed by a 1. -- Speed it up by working in ubignums instead of mantissae? m_sqrt2 :: [Int] -> [Int] m_sqrt2 m | m_is_zero m = [] | otherwise = m_sqrt2' 1 [0,1] m_0 (m_sqrt2_shiftin ([], m)) m_sqrt2' :: Int -> [Int] -> [Int] -> ([Int],[Int]) -> [Int] m_sqrt2' places one r1 (hd_o1,tl_o1) | m_is_zero hd_o1 && m_is_zero tl_o1 = [] | otherwise = (r2 ++ repeat 0) !! places : m_sqrt2' (places+2) (0:0:one) (0:r2) (m_sqrt2_shiftin (hd_o2,tl_o1)) where (r2,hd_o2) = m_sqrt2_loop one (r1, hd_o1) m_sqrt2_shiftin :: ([Int],[Int]) -> ([Int],[Int]) -- m_sqrt2_shiftin (a,b) = (a ++ take 2 (b++[0,0]), drop 2 b) m_sqrt2_shiftin (a,b) = (a ++ take 2 b, drop 2 b) m_sqrt2_loop :: [Int] -> ([Int],[Int]) -> ([Int],[Int]) m_sqrt2_loop one ro = until sqrt_loop_done sqrt_loop_body ro where sqrt_loop_done (r, o) = m_twice r `m_ge` o -- was: 2*r+1 > o sqrt_loop_body (r, o) = (m_addcarry r one, m_subcarry o (m_twice r `m_addcarry` one)) -- -- ------------ Number generators --------------------- -- -- -- bn_e: A minimally modified version of DAT's edigits -- bn_e = bn_make (1, 1, e_digits) e_digits :: [Int] e_digits = (2 : e_convert(repeat 1)) e_convert :: [Int] -> [Int] e_convert x = (head x') : e_convert (tail x') where x' = e_norm 2 (0:map (10*) x) e_norm :: Int -> [Int] -> [Int] e_norm c (d:e:x) | emodc + base_1 < c = d + edivc: e' `mod` c : x' | otherwise = d + e'divc : e'modc : x' where (e':x') = e_norm (c+1) (e:x) (edivc,emodc) = e `divMod` c (e'divc,e'modc) = e' `divMod` c -- -- ----------------- Logarithms and trigonometric functions ----------------- -- -- From libmath.b, part of the GNU "bc": -- -- exp x = 1 + x + x^2/2! + x^3/3! ..., if x <= 1 -- = (exp (x/2))^2, otherwise -- -- ln x = 2 * (a + a^3/3 + a^5/5 ...), if 0.5 < x < 2 -- where a = (x - 1) / (x + 1) -- = 2 * ln (sqrt x), otherwise -- -- sin x = x - x^3/3! + x^5/5! - x^7/7! ... -- -- cos x = sin (x + pi/2) where pi/2 = (atan 1) * 2 -- -- atan x = x - x^3/3 + x^5/5 - x^7/7 ..., if x < c -- = atan(c) + atan ( (x - c) / (1 + x * c) ), otherwise -- where c = 0.2 -- bn_exp -- -- exp x = 1 + x + x^2/2! + x^3/3! ..., if x <= 1 -- = (exp (x/2))^2, otherwise -- -- CONVERGANCE: -- x <= 1, so we just need to be sure that the divisor gets multiplied by -- at least 10 for each iteration. This is true for the term x^10/10! onward. -- If x <= 0.5 instead, it's true from the term x^5/5! onward. -- For x <= 0.1, it's true for the whole series. -- -- We use sum_series from the third term (x^2/2!) onwards. For the 4th term -- to be 1/10th of the 3rd one, x/3 must be <= 0.1 -- This is a shame because once the bottom reaches 10 * x, the series converges -- ever faster. A better strategy might be to sum the first N terms by hand -- (where N is determined by the value of x) and then apply sum_series to the -- rest, avoiding the call to bn_sqr. -- -- An alternative strategy for negative values of x would be to reduce x -- to -1 by recursion and then pair the terms before summing the series. -- -- We optimise for integral powers by using fast bn_e ^ x bn_exp x | bn_is_zero x = bn_1 | bn_is_neg x = bn_1 `bn_div` (bn_exp (bn_neg x)) | bn_is_integer x = bn_e `bn_raise` bn2Integer x | otherwise = bn_exp' x -- bn_exp' knows that its argument is > 0 bn_exp' x | x `bn_gt` (bn "0.3") = bn_sqr (bn_exp' (bn_half x)) | otherwise = (bn_1 `bn_add` x) `bn_add` (bn_sum_series (zipWith bn_quot (exp_top x) exp_bot)) -- exp_top x = [ x^2, x^3, x^4 .. ] exp_top x = iterate (bn_mul x) (bn_sqr x) -- exp_bot = [2!, 3!, 4! .. ] exp_bot = drop 2 factorials -- factorial function factorial :: Int -> Integer factorial n = factorials !! n -- Result cache for factorial function: -- the list of all factorials: 0! 1! 2! ... factorials :: [Integer] factorials = [a | (a,n) <- iterate (\(a,n)->(a*n,n+1)) (1,1)] ---- WORK IN PROGRESS ON BETTER VERSION OF EXP FUNCTION ---- -- bn_exp bottoms out on some simple values, such as 0.2 -- because an intermediate term in the series tries to sum -- one term ending in 33333.. and one ending in 66666.. -- I tried summing pairs of terms in bn_exp' before calling sum_series -- which fixed some cases while introducing others, including 0.1 -- This is an attempt at a different solution, by rewriting -- exp x = 1 + x + x^2/2! + x^3/3! + x^4/4! ... -- as -- exp x = 1 + x + 1/2 (x^2 + 1/3 (x^3 + 1/4 (x^4 + ... -- For immediately-evaluatable values, this is twice as fast as the -- above version, but as yet gives wrong answers and still bottoms -- out on the same cases as the above. bn_exp2 x | bn_is_zero x = bn_1 | bn_is_neg x = bn_1 `bn_div` (bn_exp2 (bn_neg x)) -- | bn_is_integer x = bn_e `bn_raise` bn2Integer x -- diked out for testing | otherwise = bn_exp2' x -- bn_exp2' knows that its argument is > 0 bn_exp2' x | x `bn_gt` (bn "0.1") = bn_sqr (bn_exp2' (bn_half x)) | otherwise = (bn_1 `bn_add` x) `bn_add` bn_half (bn_sumdiv_series (exp_top x) [3,4..]) ---- END OF WORK IN PROGRESS ON EXP -- bn_ln -- -- From the math library of bc: -- ln x = 2 * (a + a^3/3 + a^5/5 ...), if 0.5 < x < 2 -- where a = (x - 1) / (x + 1) -- = 2 * ln (sqrt x), otherwise -- -- CONVERGANCE: -- 0.5 < x < 2, so (x-1)/(x+1) ranges from (-.5/1.5) to (1/3), ie -1/3 to +1/3. -- Each term is slightly less than a^2 times the previous one... -- unfortunately 1/9th is not quite 1/10th! -- For the series to converge, we need -- a^2 <= 1/10 which is to say -1/sqrt(10) <= a <= 1/sqrt(10) -- Lower limit: -- (x-1)/(x+1) = -1/sqrt(10) -- x-1 = -1/sqrt(10) * (x+1) -- x-1 = (-1/sqrt(10))x + (-1/sqrt(10)) -- x (1 + 1/sqrt(10)) = 1 - 1/sqrt(10) -- x = (1 - 1/sqrt(10)) / (1 + 1/sqrt(10)) -- Upper limit: -- (x-1)/(x+1) = 1/sqrt(10) -- x-1 = (1/sqrt(10)) x + 1/sqrt(10) -- x (1 - 1/sqrt(10)) = 1 + 1/sqrt(10) -- x = (1 + 1/sqrt(10)) / (1 - 1/sqrt(10)) -- -- one_root10 = bn_1 `bn_div` (bn_sqrt (bn "10")) -- one_plus_root10 = bn_1 `bn_add` one_root10 -- one_minus_root10 = bn_1 `bn_sub` one_root10 -- lower_limit = one_minus_root10 `bn_div` one_plus_root10 -- upper_limit = one_plus_root10 `bn_div` one_minus_root10 -- -- Results: -- .51949385329... (we use ...533 to be sure) -- 1.9249505911... -- -- Fortunately, numbers very close to these limits calculate faster -- than other numbers, since the resulting series converges faster. -- We therefore keep these limits as wide as we can. -- -- If we were to sum pairs of terms before calling sum_series, we would get a -- convergance of order a^4 instead of a^2, which allows us to process a wider -- range of values of x without having to recurse (.28013 - 3.56977). This -- proves to be three or four times faster in cases where x has few digits and -- is in the range .281-.519 or 1.925-3.569, because the time taken by bn_ln -- depends on the number of digits in the operand more than anything else, -- and a square root almost always gives an infinite number of digits. -- In all other cases, this strategy turns out to be three or four times slower. -- Limits within which the summation of the series will work ln_lower_limit = bn "0.5194938533" ln_upper_limit = bn "1.9249505911" bn_ln x | bn_is_zero x = error "Can't take log of zero" | bn_is_neg x = error "Can't take log of negative numbers" | x `bn_eq` bn_1 = bn_0 | x `bn_lt` ln_lower_limit || x `bn_gt` ln_upper_limit = bn_twice (bn_ln (bn_sqrt x)) | otherwise = bn_ln' x -- Sum the series bn_ln' x = bn_twice (bn_sum_series (zipWith bn_quotI (ln_top x) [1,3..])) ln_top x = iterate (bn_mul asq) a where asq = bn_sqr a a = (x `bn_sub` bn_1) `bn_div` (x `bn_add` bn_1) -- Now we can do a power function as we should. Mind you, it's incredibly slow! -- We do not optimise integral powers into bn_raise because it makes things like -- (sqrt 2) ^ 3 bottom out needlessly. bn_pow a x = bn_exp (bn_ln a `bn_mul` x) -- bn_sin -- -- sin x = x - x^3/3! + x^5/5! - x^7/7! ... -- -- CONVERGANCE: -- In each term, the top increases by a maximum factor of pi/2^2 (1.570796..) -- and the bottom increases by 2*3 then 4*5 then 6*7, so the minimum -- convergence factor is 1.57 / 6 = 0.1309. For the series formed by summing -- pairs of adjacent terms, it is always less than 0.1. -- -- NOTE! Our value of bn_pi is calculated from the sum of a series, -- so bn_sin will be at its fastest from -pi/2 to +pi/2 -- (Though you can't take the sin of exactly pi/2 because bn_gt in bn_sin' -- goes off to infinity.) -- -- bn_sum_series normalises the result for us. -- A) reduce range of argument to first quadrant: 0..pi/2 -- 1) reduce to positive numbers: sin(n) = -sin(-n) bn_sin n | bn_is_zero n = bn_0 | bn_is_neg n = bn_neg (bn_sin' (bn_neg n)) | otherwise = bn_sin' n -- 2) reduce to 0..2*pi: for n > 2 * pi, sin(n) = sin(n - 2*pi) bn_sin' n | n `bn_gt` bn_2_pi = bn_sin' (bn_sub n bn_2_pi) | otherwise = bn_sin'' n -- 3) reduce to 0..pi: for pi < n <= 2*pi, sin(n) = -sin(n - pi) bn_sin'' n | n `bn_gt` bn_pi = bn_neg (bn_sin''' (bn_sub n bn_pi)) | otherwise = bn_sin''' n -- 4) reduce to 0..pi/2: for pi/2 < n <= pi, sin(n) = sin(pi - n) bn_sin''' n | n `bn_gt` bn_pi_2 = bn_sin'''' (bn_sub bn_pi n) | otherwise = bn_sin'''' n -- B) sum the series: sin x = x - x^3/3! + x^5/5! - x^7/7! bn_sin'''' n = bn_sum_series (sin_series n) -- Create the series for sin -- -- invert the sign of the even terms of the series -- sin_top x = [ x, -x^3, x^5, -x^7 .. ] -- NB: sin_top is re-used below in atan_series -- sin_top x = [ a | a <- x, (bn_mul minus_x_squared) a .. ] sin_top x = iterate (bn_mul minus_x_squared) x where minus_x_squared = bn_neg (bn_sqr x) -- The divisors of the terms of the series: -- sin_bot = [ factorial n | n <- [1,3..] ] sin_bot = everyother factorials where everyother (_:b:x) = b : everyother x sin_series x = addpairs (zipWith bn_quot (sin_top x) sin_bot) -- Add pairs of terms in a series (also used in bn_atan, below) addpairs :: [Bigfloat] -> [Bigfloat] addpairs (a:b:x) = bn_add a b : addpairs x addpairs x = x -- for [a] and [] -- cosine - obvious enough, since we have sin: cos(x) = sin(x + pi/2) -- We add special test for cos(0) since sin(pi/2) never returns. bn_cos x | bn_is_zero x = bn_1 -- avoid bottom-out case | otherwise = bn_sin (x `bn_add` bn_pi_2) -- atan x = x - x^3/3 + x^5/5 - x^7/7 ..., if x < c -- = atan(c) + atan ( (x - c) / (1 + x * c) ), otherwise -- where c = 0.2 -- -- By inspection with a test function, the series converges faster than -- one digit per term (paired terms, that is) with c = 0.5; even faster -- with c = 0.2 (which is the value used inside "GNU bc" whence we robbed the -- algorithm). Since, for us, the cost of doing the mul - `div` etc involved -- in recursing is higher than that of summing the series, we use 0.5. -- Experimentally (with bc), .562341325 converges ok (by a factor of >10) -- while .56241326 doesn't quite. c_lim = bn "0.562341325" -- bn_atan2 takes this two steps further: It lets you specify c at runtime, -- and applies two recursive calls in one go if appropriate. This has three -- advantages: it saves a humungous amount of multiplication and division, -- lets us calculate atan(c) once instead of twice and avoids indefinite -- recursion in intermediate results in many cases when the top and bottom -- of the division attain unfortunately-related infinitely-long values. -- In practice, with c = 0.5: -- 0 <= x <= .5 atan(x) -- .5 < x < 1.3XXX atan(.5) + atan(d) -- 1.3XXX < x < 5.5 atan(.5) + atan(.5) + atan(e) -- 5.5 <= x atan(.5) + atan(.5) + atan(.5) + atan(f) -- -- 0.5 seemed a good general purpose value for c - other values run slower -- except for some special cases - until I found out that it never gives an -- answer to atan(7), atan(9), atan(12), atan(15) and probably others too. -- -- In the case of atan(1), used for pi, c=0.4 runs staggeringly faster than 0.5: -- atan2 1 0.5 => atan 0.5 + atan 0.33333333333.. -- atan2 1 0.4 => (2 * atan 0.4) + atan 0.02439.. which converges much faster -- (requiring 345,003 reductions instead of 1,907,882 for 20 decimal places). -- Your mileage may vary taking the atan of other specific values. bn_atan2 :: Bigfloat -> Bigfloat -> Bigfloat bn_atan2 x c | c `bn_gt` c_lim = error "Internal error: bn_atan2 limit check" | bn_is_zero x = bn_0 | bn_is_neg x = bn_neg (bn_atan2' (bn_neg x) c) | otherwise = bn_atan2' x c -- bn_atan -- Specialised version of bn_atan2, using constant c_use for value of c -- giving speedup for the constant calculations that drop out. -- The speedup is large (more than twice) when several atans are performed -- in one Miranda execution since atan(c) and derived calculation are only -- done once. -- For the generic constant-c version, the value of c to use. -- No single digit version of c works for small positive integers: -- they all (.2, .3, .4, .5) go into infinite regress in the summation -- of the series, so we use a 2-digit value. -- -- Experimentally, using values from .41 to .56 for c_use -- .45 and .48 bottom out in atan(6) -- .56 bottoms out on atan(37) -- .49 was the fastest in tests from 1 to 101 c_use = bn "0.49" bn_atan x | bn_is_zero x = bn_0 | bn_is_neg x = bn_neg (bn_atan' (bn_neg x)) | otherwise = bn_atan' x -- atan2 x c -- = atan x , if x <= c -- = atan c + atan d , if x > c && d <= c -- = atan c + atan c + atan2 e c , otherwise -- where -- d = (x - c) / (1 + x * c) -- e = (d - c) / (1 + d * c) -- -- e = (x - 2c - x c^2) / (1 + 2xc - c^2) -- or (x(1-c^2) - 2c) / (1 + 2xc - c^2) -- -- bn_atan2' knows that its parameter x is not negative. bn_atan2' :: Bigfloat -> Bigfloat -> Bigfloat bn_atan2' x c | x `bn_le` c = bn_atan'' x | d `bn_le` c = bn_atan'' c `bn_add` (bn_atan2' d c) | otherwise = bn_twice (bn_atan'' c) `bn_add` (bn_atan2' e c) where d = (x `bn_sub` c) `bn_div` (bn_1 `bn_add` x_mul_c) e = etop `bn_div` ebot etop = (x `bn_mul` (bn_1 `bn_sub` csq)) `bn_sub` (bn_twice c) ebot = bn_1 `bn_add` (bn_twice x_mul_c) `bn_sub` csq csq = bn_sqr c x_mul_c = x `bn_mul` c -- Now that we have arctangent, we can calculate pi as 4 * atan(1) -- -- "Now I, even I would celebrate -- In rhymes unapt the great -- Immortal Syracusan, rivaled nevermore -- Who, in his wondrous lore -- Passed on before -- Left men his guidance -- How to circles mensurate." bn_pi30 = bn "3.141592653589793238462643383279" bn_2_pi = bn_twice bn_pi -- 2 * pi (used in bn_sin) bn_pi = bn_twice bn_pi_2 -- pi bn_pi_2 = bn_twice bn_pi_4 -- pi / 2 (used in bn_sin) -- bn_pi_4 = bn_atan2 bn_1 (bn "0.4") -- pi / 4 -- for pi/4, calculated with atan(1) using c = 0.4, -- we expand the first iteration of the atan series with constant values: -- d = (1 - .4) / (1 + .4) = .6/1.4 = 3/7 = 0.428571..., which is not < .4 -- so... -- etop = (1 * (1 - .4^2)) - (2 * .4) = 1-.16 - .8 = 0.04 -- ebot = 1 + (2 * 1 * .4) - (.4^2) = 1.8 - .16 = 1.64 -- e = 0.04 / 1.64 = 4 / 164 = 2 / 82 = 1 / 41 -- pi/4 = (2 * atan(0.4)) + atan(1/41) -- this is 4.5% better in space and speed than the bn_pi_4 above. bn_pi_4 = bn_twice (bn_atan'' (bn "0.4")) `bn_add` (bn_atan'' (bn_quotI bn_1 41)) -- Optimised version using constant c_use as the value of c bn_atan' :: Bigfloat -> Bigfloat bn_atan' x | x `bn_le` c_use = bn_atan'' x | d `bn_le` c_use = atan_c `bn_add` (bn_atan' d) | otherwise = twice_atan_c `bn_add` (bn_atan' e) where d = (x `bn_sub` c_use) `bn_div` (bn_1 `bn_add` x_mul_c) e = etop `bn_div` ebot etop = (x `bn_mul` one_sub_sqr_c) `bn_sub` twice_c ebot = one_sub_sqr_c `bn_add` (bn_twice x_mul_c) x_mul_c = x `bn_mul` c_use -- In Miranda, global constant functions seem to fare better than derived -- constant calculations in where clauses. atan_c = bn_atan'' c_use twice_atan_c = bn_twice atan_c sqr_c = bn_sqr c_use twice_c = bn_twice c_use one_sub_sqr_c = bn_1 `bn_sub` sqr_c -- START OF WORK IN PROGRESS ON OPTIMISED ATAN WITH FIXED C=0.5 -- Even more specialised version uses constant 0.5 for value of c -- This lets us use bn_times, bn_quot, bn_half instead of mul and `div` -- -- This works and is 3% faster in space and speed than bn_atan, -- but bottoms out on certain values (7, 9, 12, 15 ...) -- probably because unfortunate pairs occur in the summing series: -- infinitely-long terms that, when subtracted, would give a -- finite-length sum. bn_atan5 x | bn_is_zero x = bn_0 | bn_is_neg x = bn_neg (bn_atan5' (bn_neg x)) | otherwise = bn_atan5' x -- atan5 x -- = atan x , if x <= .5 -- = atan .5 + atan d , if x > .5 && d <= .5 -- = atan .5 + atan .5 + atan5 e , otherwise -- where -- d = (x - .5) / (1 + x * .5) = (x - .5) / (1 + x/2) -- e = (d - .5) / (1 + d * .5) -- -- e = (x - 2*.5 - x .5^2) / (1 + 2x*.5 - .5^2) -- = (x(1-.5^2) - 2*.5) / (1 + x - .25) -- = (.75*x - 1) / (.75 + x) bn_atan5' :: Bigfloat -> Bigfloat bn_atan5' x | x `bn_le` bn_p5 = bn_atan'' x | d `bn_le` bn_p5 = atan_p5 `bn_add` (bn_atan5' d) | otherwise = twice_atan_p5 `bn_add` (bn_atan5' e) where d = (x `bn_sub` bn_p5) `bn_div` (bn_1 `bn_add` (bn_half x)) e = ((bn_half.bn_half)(bn_times 3 x) `bn_sub` bn_1) `bn_div` (bn_p75 `bn_add` x) bn_p5 = bn(".5") bn_p75 = bn(".75") atan_p5 = bn_atan'' bn_p5 twice_atan_p5 = bn_twice atan_p5 -- END OF WORK IN PROGRESS ON OPTIMISED ATAN -- This version always sums the series. x must be <= c_lim -- so that the series converges ok. bn_atan'' x = bn_sum_series (atan_series x) -- The terms of the series for atan(x). -- The dividends are the same as for sin(), so we reuse them. atan_series x = addpairs (zipWith bn_quotI (sin_top x) [1,3..]) bn_tan x = (bn_sin x) `bn_div` (bn_cos x) bn_sinh x = bn_half ( bn_exp x `bn_sub` (bn_exp (bn_neg x)) ) bn_cosh x = bn_half ( bn_exp x `bn_add` (bn_exp (bn_neg x)) ) bn_tanh x = (bn_sinh x) `bn_div` (bn_cosh x) bn_asinh x = bn_ln ( x `bn_add` (bn_sqrt (bn_sqr x `bn_add` bn_1)) ) bn_acosh x = bn_ln ( x `bn_add` (bn_sqrt (bn_sqr x `bn_sub` bn_1)) ) bn_atanh x = bn_half ( bn_ln (bn_add bn_1 x) `bn_sub` (bn_ln (bn_sub bn_1 x)) ) -- bn_sum_series sums a converging series of Bigfloat. It requires that every -- term in the series be less that 1/10th (ok, 1/baseth) of the previous one -- and that all elements in the series be of the same sign (all +ve or all -ve). -- ubn2bn normalises the result for us. -- This copes with series in which the terms are all positive or all negative. -- The check for bn_is_zero (head s) is necessary because if the first element -- is zero, we can't tell whether the elements are all positive or all negative. -- If we are passed a series of mixed positive and negative terms by mistake, -- the calculation will throw a fatal error in bn2ubn. -- Every term in the series has already been normalised. bn_sum_series :: [Bigfloat] -> Bigfloat -- first, eliminate two simple cases -- (These never happen because we are always handed infinite lists) -- bn_sum_series [] = bn_0 -- bn_sum_series [a] = a -- now reduce the job to summing a positive series bn_sum_series s | bn_is_zero (head s) = bn_sum_series (tail s) | bn_is_neg (head s) = bn_neg (bn_sum_series (map bn_neg s)) | otherwise = ubn2bn (ubn_sum_series (map bn2ubn s)) -- unsigned version - tack a leading zero onto all terms so that we can be -- sure that the sum of any term and all the following ones will have the -- same exponent as the term itself. The fact that the first digit of each -- term will now be 0 is taken account of in the unequal-exponent lazy -- lookahead code, so is not as much of a lossage as it might seem. -- Without it, we produce "10" digits from time to time. -- However it is much faster (2/3 of the reductions). Can we use this -- by letting it produce 10 digits and then applying m_carry1? ubn_sum_series :: [UBigfloat] -> UBigfloat ubn_sum_series s = ubn_sum_series2 (map prefix_0 s) where prefix_0 (e,m) = (e+1, 0:m) -- ubn_sum_series2 is the second version of the series summer that uses the same -- kind of logic as ubn_add2 to return the initial digit from the first term -- through inspection only of the first term and the exponent of the second. -- -- When we can no longer be sure that the digits of the result will be the -- same as the digits of the first term, we add the first and second terms -- and repeat the process with the sum as the first term of the list. -- -- We know that zeroes have been prefixed to all terms, guaranteeing that the -- exponent of the sum from a particular term onwards will not grow as a result -- of adding in the following terms, and that the maximum addition from the -- first digit of the second term and all that follow it will be 1. This -- corresponds to the logic of seeing whether there are any non-9 digits in -- ubn_add2. ubn_sum_series2 :: [UBigfloat] -> UBigfloat -- From now on we don't normalise or adjust exponents, so return -- the exponent and work out the mantissa. ubn_sum_series2 ((ae,am):x) = (ae, ubn_sum_series2' ((ae,am):x)) -- These never happen cos all the series we use are infinite... -- ubn_sum_series2 [a] = a -- ubn_sum_series2 [] = ubn_0 -- ubn_sum_series2' just returns the mantissa of the sum, which has the same -- exponent as the first term of the series. ubn_sum_series2' :: [UBigfloat] -> Mantissa -- Most common case first... ubn_sum_series2' ((ae,(a:ax)):(be,bm):x) -- We can be sure of the first digit of a if ae-be >= 2 and there is a -- digit between the first digit of a and the digit above the first of -- b whose value is less that (base-1). Remember that the first digit -- of b will be 0, and so the most that the sum of b onwards could add -- to the final result will be a 1 in the first position of b. -- Note: the any_lt test *includes* the first digit of a that overlaps b. -- I think this is right... Here's the worst case: -- a: 0 8 9 9 9 9 9 9 9 9 9 9 -- b: 0 9 9 9 9 9 9 9 9 9 -- c: 0 9 9 9 9 9 9 9 9 -- b+c: 1 0 9 9 9 9 9 9 9 9 -- abc: 9 -- .... so given the 8, the first digit cannot be affected by adding in -- b,c... ok! In fact, this is worse than the worst case because -- .0099 is not < 1/10th of .0899 | ae-be >= 2 && any_lt_base_1 ax (ae-be) = a : ( ubn_sum_series2' ((ae-1,ax):(be,bm):x) ) -- Otherwise add the first two terms of the series together and repeat the -- process with this as the new first term of the series, which will give -- us more distance between the exponents. | otherwise = ubn_sum_series2' ((ae, (m_addcarry (a:ax) (prefix0s (ae-be) bm))) : x) -- If am is exhausted, eliminate the term and proceed with the rest ubn_sum_series2' ((ae,[]):(be,bm):x) = prefix0s (ae-be) (ubn_sum_series2' ((be,bm):x)) -- Deal with terminal cases ubn_sum_series2' [(e,m)] = m -- One item ubn_sum_series2' [] = [] -- No items ---- WORK IN PROGRESS ON BETTER VERSION FOR EXP FUNCTION ---- -- This specialised version is used to evaluate formulae of the form: -- a0 + 1/d0 * (a1 + 1/d1 * (a2 + 1/d2 * (... -- for example exp(), which is 1 + x + x^2/2 + x^3/3! + x^4/4! -- rewritten as 1 + x + 1/2 (x^2 + 1/3 (x^3 + 1/4 (x^4 + ... -- Convergance: a(n+1) <= a(n)/base, d(x) > 1 and d(n+1) > d(n) -- I hope this means that the RHS of every '+' is always < 1/base of the LHS, -- allowing us to reason as above about the digits of a0 that we can -- safely return by examining only a0 and the exponent of a1. -- The actual division is then applied in the recursive call. -- This is twice as fast as bn_exp for values < 0.1, gives wrong answers, -- and locks solid after 5 digits of exp "0.1" the same as the sum_series -- version. bn_sumdiv_series :: [Bigfloat] -> [Int] -> Bigfloat bn_sumdiv_series s t | bn_is_zero (head s) = bn_sumdiv_series (tail s) (tail t) `bn_quotI` (head t) | bn_is_neg (head s) = bn_neg (bn_sumdiv_series (map bn_neg s) t) | otherwise = ubn2bn (ubn_sumdiv_series (map bn2ubn s) t) ubn_sumdiv_series :: [UBigfloat] -> [Int] -> UBigfloat ubn_sumdiv_series s t = ubn_sumdiv_series2 (map prefix_0 s) t where prefix_0 (e,m) = (e+1, 0:m) ubn_sumdiv_series2 :: [UBigfloat] -> [Int] -> UBigfloat ubn_sumdiv_series2 ((ae,am):x) t = (ae, ubn_sumdiv_series2' ((ae,am):x) t) ubn_sumdiv_series2' :: [UBigfloat] -> [Int] -> Mantissa -- First element exhausted ubn_sumdiv_series2' ((ae,[]):(be,bm):x) (d:y) = prefix0s (ae-be) ( (ubn_sumdiv_series2' ((be,bm):x) y) `m_quotI` d ) ubn_sumdiv_series2' ((a0e,a0m):(a1e,a1m):ax) (d0:d1:dx) | a0e-a1e >= 2 && any_lt_base_1 (tail a0m) (a0e-a1e) = head a0m : ( ubn_sumdiv_series2' ((a0e-1,tail a0m):(a1e,a1m):ax) (d0:d1:dx) ) | d0 > 200000000 = error "Ooverflue" | otherwise -- sumdiv [ a0 + a1/d0, a2, a3 .. ] [ d0 * d1, d2, d3 .. ] = ubn_sumdiv_series2' ( (a0e,m_addcarry a0m (prefix0s (a0e-a1e) a1m `m_quotI` d0)) : ax ) ( d0*d1 : dx ) ubn_sumdiv_series2' [(e,m)] (d:_) = m_quotI m d -- One item ubn_sumdiv_series2' [] _ = [] -- No items ---- END OF WORK IN PROGRESS FOR EXP ---- ---- ...HERE ENDETH ---- ---- End of BigFloat library ----