|| bignum - Infinite precision calculator, Martin Guy, Catania, October 2001.
|| All original, except for edigits, modified from DAT's example script.
||
|| Numbers are represented as as infinite series of numbers which in a decimal
|| system would be referred to as decimal places.
|| The number represented by the bignum (4, [5, 2, 7]) is 0.527 * 10^4.
|| Support is not guaranteed for bases other than 10, but no base is
|| assumed by almost any of the functions.
||
|| bignum is a tuple containing sign, exponent and mantissa in that order.
||
|| in the code we have three principal sorts of functions:
|| bn_* treat bignums (sign, exp, mantissa)
|| ubn_* treat unsigned bignums (exp, mantissa)
|| m_* just treat mantissae [num]
|| and normally the bn_ functions use ubn functions to achieve their result,
|| which in turn call m_ functions to perform the detailed calculation.
%export bignum bn showbignum bn_0 bn_1 bn_2 bn_is_neg bn_is_zero bn_cmp bn_eq bn_ne bn_gt bn_lt bn_ge bn_le bn_neg bn_add bn_sub bn_twice bn_half bn_times bn_mul bn_raise bn_quot bn_div bn_sqr bn_sqrt bn_rnd bn_e bn_pi bn_ln bn_exp bn_pow bn_sin bn_cos bn_atan
|| 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.
base = 10
|| Constants derived from base. Literal values are notably faster than values
|| derived from simple formulae of constants (which is strange - shouldn't
|| Miranda reduce them to a singularity on first evaluation?)
|| base_1 = base - 1
base_1 = 9
|| base_2 = base ^2
base_2 = 100
|| base_div_2 = base div 2
base_div_2 = 5
|| The number of digits that bn_show prints. if scale = 0, infinite.
scale = 0
|| The truncation function used in the printing functions below.
scaletrunc x
= take scale x, if scale > 0
= x, otherwise
|| Exported abstract type bignum represents an infinite-precision number.
abstype bignum
with
|| Internal functions, not exported
|| signed <-> unsigned bignum converters
bn2ubn :: bignum -> ubignum
ubn2bn :: ubignum -> bignum
|| normaliser strips leading zeroes from mantissa, adjusting
|| exponent accordingly. It also corrects the "-0" case.
bn_normalise :: bignum -> bignum
|| bn_trim strips trailing zeroes from the mastissa.
bn_trim :: bignum -> bignum
|| Exported functions dealing directly with bignums:
|| bn_make and bn_unmake construct and destruct a bignum.
bn_make :: (num,num,[num]) -> bignum
bn_unmake :: bignum -> (num,num,[num])
|| signof returns the exponent of a bignum
sign_of :: bignum -> sign
|| exponent_of returns the exponent of a bignum
exponent_of :: bignum -> exponent
|| mantissa_of returns the mantissa of a bignum
mantissa_of :: bignum -> mantissa
|| bn converts a decimal string representation of a number
|| into a bignum. It's the easiest way of entering constants into
|| an equation, and accepts all the formats that "showbignum" outputs.
bn :: [char] -> bignum
|| showbignum converts a bignum into a printable form: an initial
|| '-' if it's negative, followed by one of three basic formats:
|| - a string of decimal digits if the bignum has no fractional part
|| (ie represents an integer)
|| - a string of digits + a period + a string of digits if the bignum
|| is >= 1 and is not integral
|| - "0." ++ a string of digits if the bignum is < 1.
|| showbignum is picked up by the Miranda system as the standard
|| function for outputting bignums.
showbignum :: bignum -> [char]
|| bn_is_neg tells you whether a bignum is a negative number
|| (0 is not negative!)
bn_is_neg :: bignum -> bool
bn_is_zero :: bignum -> bool
|| bn_cmp compares two bignums and returns:
|| 0 if the two numbers have the same value,
|| 1 if the first number is greater than the second
|| -1 if the second number is greater than the first.
bn_cmp :: bignum -> bignum -> num
|| The six usual comparison functions:
|| == ~= > < >= <=
bn_eq :: bignum -> bignum -> bool
bn_ne :: bignum -> bignum -> bool
bn_gt :: bignum -> bignum -> bool
bn_lt :: bignum -> bignum -> bool
bn_ge :: bignum -> bignum -> bool
bn_le :: bignum -> bignum -> bool
|| bn_rnd takes a positive integer seed parameter and returns a
|| pseudo-random bignum in the range 0 <= bn_rnd n < 1
bn_rnd :: num -> bignum
|| bn_neg negates a bignum
bn_neg :: bignum -> bignum
|| bn_add adds two bignums together, returning their sum.
bn_add :: bignum -> bignum -> bignum
|| bn_sub subtracts its second bignum parameter from its first,
|| returning the difference. Since negative numbers are not
|| representable (yet), a negative result is a fatal error.
bn_sub :: bignum -> bignum -> bignum
|| bn_twice multiplies by two.
bn_twice :: bignum -> bignum
|| bn_half divides by two.
bn_half :: bignum -> bignum
|| bn_times multiplies an integer by a bignum.
|| Done by repeated addition.
bn_times :: num -> bignum -> bignum
|| bn_mul multiplies two bignums.
bn_mul :: bignum -> bignum -> bignum
|| bn_raise raises a bignum to an integral power
bn_raise :: bignum -> num -> bignum
|| bn_quot divides a bignum by a positive integer
bn_quot :: bignum -> num -> bignum
|| bn_div performs division on two bignums
bn_div :: bignum -> bignum -> bignum
|| bn_sqr returns the square of its argument.
bn_sqr :: bignum -> bignum
|| bn_sqrt returns the square root of its argument.
bn_sqrt :: bignum -> bignum
|| Logarithms and trigonometric functions
bn_ln :: bignum -> bignum
bn_exp :: bignum -> bignum
bn_pow :: bignum -> bignum -> bignum
bn_sin :: bignum -> bignum
bn_cos :: bignum -> bignum
bn_atan :: bignum -> bignum
|| Here beginneth...
bignum == ( sign, exponent, mantissa )
ubignum == ( exponent, mantissa )
sign == num
exponent == num
mantissa == [num]
|| 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 bignum constants
ubn_0 :: ubignum
ubn_0 = (0, [])
||
|| ---------- simple conversion functions ----------
||
|| make a bignum 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 bignum
sign_of (s,e,m) = s
exponent_of (s,e,m) = e
mantissa_of (s,e,m) = m
|| functions to select parts of a ubignum
uexponent_of (e,m) = e
umantissa_of (e,m) = m
|| Convert bignum to unsigned bignum and vice versa
bn2ubn (1,e,m) = (e,m)
bn2ubn other = error "Internal error: bn2ubn of negative number"
ubn2bn (e,m) = bn_normalise (1,e,m)
|| bn_normalise converts a bignum 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
|| m_trim strips trailing zero digits from a mantissa. It's tempting to use
|| it liberally because trailing zeroes cause pointless calculation.
|| In practice, it's hardly ever worth doing. Instead, we test for a few cases
|| here and there where it costs us next to nothing to check for a generated
|| trailing zero (e.g. in bn_twice and bn_carry1), and then just apply m_trim
|| in showbignum to suppress trailing zeroes in the output.
m_trim :: mantissa -> mantissa
m_trim (0:x)
= [], if m_is_zero x
= 0 : m_trim x, otherwise
m_trim (a:x) = a : m_trim x
m_trim [] = []
|| ubignum and bignum versions:
ubn_trim :: ubignum -> ubignum
ubn_trim (e,m) = (e, m_trim m)
bn_trim (s,e,m) = (s, e, m_trim m)
||
|| ---------- text <-> bignum conversion functions ----------
||
|| Convert bignum to text.
showbignum (s,e,m)
= sign ++ ubn_show (e, m_trim m) || suppress trailing zeroes
where
sign = "-", if s < 0
= "", otherwise
|| We produce output in three forms: .456, 123 and 123.456.
|| ubn_show deals with the first form, ubn_show' with the other two forms.
ubn_show :: (num,[num]) -> [char]
ubn_show (e,m)
= "0", if m_is_zero m
= ubn_show' (e,m), if e > 0
= "." ++ scaletrunc ((rep (-e) '0') ++ m_show m), if e <= 0
|| Deal with bignums whose exponent > 0 (i.e. >= 1.0)
ubn_show' :: (num,[num]) -> [char]
ubn_show' (e, []) = rep e '0' || Do final zeroes (if any) of integral bn.
ubn_show' (0, m) = '.' : scaletrunc (m_show m)
ubn_show' (e, (a:x)) = mkdigit a : ubn_show' (e-1, x)
m_show :: [num] -> [char]
m_show m = map mkdigit m
mkdigit :: num -> char
mkdigit n = decode(n + code '0'), if 0 <= n < 10
= decode(n - 10 + code 'a'), if 10 <= n < 36
|| bn: Convert text to bignum
|| Syntax: ['-'] [digit] [digit]* \/ ['-'] [digit]* ++ '.' ++ [digit]*
|| Once we've found a '.', we can accept an infinite number of decimal places.
|| Until we find a '.' (or the end of the string), we don't know what the
|| exponent will be.
|| For simplicity, we also allow "", "-", "." and "-." as synonyms for 0.
bn a = bn_normalise (bn_trim (bn_text a))
where
|| We don't make bn_text public because it returns unnormalised bignums
|| (thus removing a trap for the unwary).
bn_text ('-':x) = (-1,e,m) where (e,m) = ubn_text x
bn_text x = (1,e,m) where (e,m) = ubn_text x
|| Convert text to unsigned bignum
ubn_text :: [char] -> ubignum
ubn_text (a:x)
= (0, m_text x), if a = '.'
= (xe+1, (code a - code '0'):xm), if isdigit a
= error "Non-digit in bignum text", otherwise
where (xe,xm) = ubn_text x
ubn_text [] = (0,[])
m_text :: [char] -> mantissa
m_text (a:x)
= (code a - code '0'): m_text x, if isdigit a
= error "Non-digit in bignum text", otherwise
m_text [] = []
isdigit :: char -> bool
isdigit c = code '0' <= code c <= code '9'
|| ubn_same_exp takes a pair of bignums 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 bignum with the smaller exponent.
ubn_same_exp :: (ubignum, ubignum) -> (ubignum, ubignum)
ubn_same_exp (a, b)
= ((ae,am), (be,bm)) , if ae = be
= (ubn_abnormalise be a, b) , if ae < be
= (a, ubn_abnormalise ae b) , if ae > be
where (ae,am) = a
(be,bm) = b
|| ubn_abnormalise sets the exponent of its second parameter to its first
|| parameter, giving as a result an un-normalised bignum. 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 -> ubignum -> ubignum
ubn_abnormalise new_e (old_e, m)
= (new_e, prefix0s (new_e - old_e) m) , if new_e >= old_e
= error "ubn_abnormalise called with duff parameters", otherwise
|| Negate a bignum, avoiding the "-0" syndrome.
bn_neg (s,e,m)
= bn_0, if m_is_zero m
= (-s,e,m), otherwise
||
|| ----------- Stuff for addition/subtraction of bignums ------------
||
|| bn_add gives the arithmetic sum of two bignums.
|| 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, ubn_add (ae,am) (be,bm)), if as = bs
= (as, ubn_sub (ae,am) (be,bm)), if (ae,am) $ubn_ge (be,bm)
= (bs, ubn_sub (be,bm) (ae,am)), otherwise
|| bn_sub performs subtraction of 2 bignums.
bn_sub a b = bn_add a (bn_neg b)
|| ubn_add gives the arithmetic sum of two unsigned bignums.
|| 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.
|| To do: take advantage of unequal exponents and do the sum-carry only on the
|| parts of the numbers that overlap.
ubn_add :: ubignum -> ubignum -> ubignum
ubn_add = ubn_add2
ubn_add1 :: ubignum -> ubignum -> ubignum
ubn_add1 a b
= (ae+1, m_addcarry0 am bm)
where
((ae,am),(be,bm)) = ubn_same_exp (a, b)
|| ubn_add2 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.
|| It would be better if we returned the first sure digit without having to
|| inspect *all* the sure digits (and the digit following the last one).
ubn_add2 :: ubignum -> ubignum -> ubignum
ubn_add2 (ae,am) (be,bm)
|| First, make sure ae >= be
= (ae+1, m_addcarry0 am bm), if ae = be
= ubn_add' (ae,am) (be,bm), if ae > be
= ubn_add' (be,bm) (ae,am), otherwise
|| We can only exploit a difference >= 2 in the exponents
ubn_add' :: ubignum -> ubignum -> ubignum
ubn_add' (ae,am) (be,bm)
= (ae+1, m_addcarry0 am bm), if ae = be
= (ae+1, m_addcarry0 am (0:bm)), if ae = be + 1
|| OK! the difference in the exponents is >= 2:
= (ae, (take nsd am) ++ m_addcarry
(drop nsd am)
(prefix0s (ae-be-nsd) bm)), otherwise
where
nsd = n_sure_digits (ae,am) be
|| END OF WORK IN PROGRESS on new add routine
|| ubn_sub performs subtraction of 2 unsigned bignums.
|| 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.
|| To do: take advantage of ae > be and return the initial part of am
|| unaltered rather than subtracting a string of zeroes from it.
ubn_sub :: ubignum -> ubignum -> ubignum
ubn_sub a b
= (ae, m_subcarry2 am bm)
where
((ae,am),(be,bm)) = ubn_same_exp (a, b)
|| 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 prefizes 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 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 neg 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 neg 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 :: [num] -> [num]
m_carry1 (a:b:x)
= a : m_carry1 (b : x) , if b < base_1
= (a + 1) : m_carry1 ((b - base) : x) , if b > base_1
= (a + carry): m_carry1 ((b-carry*base):x), otherwise
where carry = m_carry1from (b:x)
m_carry1 [a]
= [a], if a > 0
= [], otherwise
m_carry1 [] = []
|| mcarry1a is like mcarry1 except that it never removes trailing zeroes
m_carry1a :: [num] -> [num]
m_carry1a (a:b:x)
= a : m_carry1a (b : x) , if b < base_1
= (a + 1) : m_carry1a ((b - base) : x) , if b > base_1
= (a + carry): m_carry1a ((b-carry*base):x), otherwise
where carry = m_carry1from (b:x)
m_carry1a [a] = [a]
m_carry1a [] = []
|| 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 :: [num] -> num
m_carry1from (a:x)
= 0 , if a < base_1
= 1 , if a > base_1
= m_carry1from x, otherwise
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 :: [num] -> [num]
m_borrow1 (a:b:x)
= a : m_borrow1 (b:x) , if b > 0
= (a-1) : m_borrow1 ((b+base):x) , if b < 0
= (a-borrow) : m_borrow1 ((b + borrow * base) : x), otherwise
where borrow = m_borrow1from (b:x)
m_borrow1 [a]
= [a], if a > 0
= [], if a = 0
= error "Negative result in m_borrow", otherwise || Never happens
m_borrow1 [] = []
|| m_borrow1a is like m_borrow1 except that it doesn't strip trailing zeroes.
m_borrow1a :: [num] -> [num]
m_borrow1a (a:b:x)
= a : m_borrow1a (b:x) , if b > 0
= (a-1) : m_borrow1a ((b+base):x) , if b < 0
= (a-borrow) : m_borrow1a ((b + borrow * base) : x), otherwise
where borrow = m_borrow1from (b:x)
m_borrow1a [a] = [a]
m_borrow1a [] = []
m_borrow1from (a:x)
= 0 , if a > 0
= 1 , if a < 0
= m_borrow1from x, otherwise
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.
bn_twice (s,e,m)
= bn_0, if m_is_zero m
= bn_add (s,e,m) (s,e,m), if base mod 2 ~= 0
= (s, e, m_twice m), if hd m < base_div_2
= (s, e+1, m_twice (0:m)), otherwise
|| m_twice is always applied to a list whose first element is < base/2
m_twice :: [num] -> [num]
m_twice (a:b:x)
= (a+a) : m_twice (b:x), if b < base_div_2
= (a+a+1) : m_twice ((b - base_div_2) : x), otherwise
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)
|| m_half only works for even number bases
= bn_normalise (s, e, m_half m), if base mod 2 = 0
= bn_quot (s,e,m) 2, otherwise
m_half (a:x)
= (a div 2) : m_half x, if a mod 2 = 0
= (a div 2) : (m_add [base_div_2] (m_half x)), otherwise
m_half [] = []
||
|| --------- Simple comparison functions ----------
||
|| Is a bignum negative?
bn_is_neg (s,e,m) = (s < 0)
|| Is it zero?
bn_is_zero (s,e,m) = m_is_zero m
|| Return 0 if the two numbers have the same value,
|| return 1 if the first number is greater than the second
|| return -1 if the second number is greater than the first.
|| Assumes that the parameters are normalised.
bn_cmp (as,ae,am) (bs,be,bm)
= as * ubn_cmp (ae,am) (be,bm), if as = bs
= as, otherwise
|| tests functions for equality, inequality, less then, less that or equal,
|| greater than, greater than or equal.
bn_eq a b = bn_cmp a b = 0
bn_ne a b = bn_cmp a b ~= 0
bn_lt a b = bn_cmp a b < 0
bn_le a b = bn_cmp a b <= 0
bn_gt a b = bn_cmp a b > 0
bn_ge a b = bn_cmp a b >= 0
|| Unsigned bignum comparison; result as per bn_cmp
ubn_cmp :: ubignum -> ubignum -> num
ubn_cmp (ae,am) (be,bm)
= 1 , if ae > be
= -1 , if ae < be
= m_cmp am bm , if ae = be
|| tests functions for equality, inequality, less then, less that or equal,
|| greater than, greater than or equal.
|| Type declarations are inherited from ubn_cmp.
ubn_eq a b = ubn_cmp a b = 0
ubn_ne a b = ubn_cmp a b ~= 0
ubn_lt a b = ubn_cmp a b < 0
ubn_le a b = ubn_cmp a b <= 0
ubn_gt a b = ubn_cmp a b > 0
ubn_ge a b = ubn_cmp a b >= 0
m_cmp :: mantissa -> mantissa -> num
|| Compare two mantissae.
m_cmp (a:x) (b:y) = 1 , if a > b
= -1 , if a < b
= m_cmp x y , if a = b
m_cmp a [] = 0 , if m_is_zero a
= 1 , otherwise
m_cmp [] b = 0 , if m_is_zero b
= -1 , otherwise
|| equal, not equal, less then [or equal to], greater than [or equal to]
|| for mantissae.
|| Type declarations are inherited from m_cmp.
m_eq a b = m_cmp a b = 0
m_ne a b = m_cmp a b ~= 0
|| m_lt a b = m_cmp a b < 0
|| m_le a b = m_cmp a b <= 0
|| m_gt a b = m_cmp a b > 0
|| m_ge a b = m_cmp a b >= 0
|| optimised m_lt, used in m_div
m_lt (a:x) (b:y)
= (a < b), if a ~= b
= m_lt x y, otherwise
m_lt a [] = False
m_lt [] b = ~ (m_is_zero b)
|| reuse code
m_gt a b = m_lt b a
m_ge a b = ~(m_lt a b)
m_le a b = ~(m_lt b a)
m_is_zero :: mantissa -> bool
|| Are all elements of a mantissa equal to zero?
|| It'd be nice to use a "fold" function here, but I've seen no guarantee that
|| the "&" and "\/" operators are lazy on the right hand parameter when the left
|| hand one determines the result (or vice versa), even though it seems so
|| experimentally.
|| Was:
|| m_is_zero [] = True
|| m_is_zero (a:x)
|| = False, if a > 0
|| = m_is_zero x, if a = 0
|| = error "negative number in mantissa!", otherwise
|| Note: The following version applied to [0,0..] provokes
|| "impossibile tag (-9) in reduce" rather than "Out of heap space" from the
|| Linux version 2.025 (22 mar 97) of the Miranda system but is much faster.
m_is_zero x = and (map (= 0) x)
||
|| ------------- Simple multiply by an integer -------------
||
|| bn_times multiplies bignum 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
= error "bn_times can only multiply by integers.", if ~ (integer n)
= bn_times' n b, if n > 0
= bn_neg (bn_times' (-n) b), if n < 0
= bn_0, if n = 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
= bn_twice (bn_times' (n div 2) b), if n mod 2 = 0
= b $bn_add (bn_times' (n-1) b), if n > 2
= b, if n = 1
|| Same stuff for mantissae. Caller must guarantee no overflow, and
|| maximum n is base-1.
m_times :: num -> mantissa -> mantissa
m_times n m
= m_carry2 (map (* n) m), if 0 <= n < base
= error ("m_times " ++ show n), otherwise
||
|| ------------- Long multiplication -------------
||
|| bn_mul: Multiply two bignums. bn_sqr is a modified version of this.
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 :: ubignum -> ubignum -> ubignum
|| Add a leading 0 so that there's somewhere for the carry to propagate into.
|| If we were to simply m_multiply the mantissas without prefixing a 0, the
|| exponent of the result would be (ae + be - 1) (think: 0.1 * 0.1 = 0.01).
|| The extra digit that m_mul prefixes to the result cancels the "- 1".
ubn_mul (ae,am) (be,bm)
= (ae+be, m_mul am bm)
|| m_mul multiplies two mantissas.
||
|| NOTE! m_mul prefixes a digit to the result so that there's somewhere for
|| the carry from the first result digit to propagate into,
|| so callers must adjust result exponents accordingly.
||
|| 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 maximul 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 :: [num]->[num]->[num]
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 :: [num]->[num]->[num]
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 :: [num]->[num]->[[num]]
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 :: [[num]] -> [num]
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 :: [num] -> [num] -> [num]
reduce_for_carry (a:b:x) (c:d:y)
= a : reduce_for_carry (b:x) (d:y), if d < base
= reduce_for_carry (carry1digit (a:b:x)) (carry1digit (c:d:y)), otherwise
||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 :: [num] -> [num]
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.
= (a + bc div base) : m_carry2 (bc mod base : c mod base : x),
if bc mod base < base_1
|| general case
= (a + carry) : m_carry2 ((b - carry * base) : c : x),
otherwise
where
bc = b + c div base
carry = m_carry2from (b:c:x)
m_carry2 [a,b]
= [a + b div base, b mod base], if b mod base > 0
= [a + b div base], otherwise
m_carry2 [a]
= [a], if a > 0
= [], 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, 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,
|| x won't affect the carry from a.
m_carry2from :: [num] -> num
m_carry2from (a:b:x)
|| This first case gains so infrequently that it ends up a net loss.
|| = a div base , if a mod base < 2
= ab div base , if ab mod base < base_1
= (a + m_carry2from(b:x)) div base, otherwise
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 bignum to raise to a power;
|| x is the integral power to raise it to.
|| Optimisation: bn_raise x (2*n) = sqr (bn_raise x n)
bn_raise a x
= error "bn_raise only does integer powers", if ~ (integer x)
= bn_div bn_1 (bn_raise' a (-x)), if x < 0
= bn_1, if x = 0
= bn_raise' a x, otherwise
bn_raise' a x
= bn_sqr (bn_raise' a (x div 2)), if x mod 2 = 0
= bn_mul a (bn_raise' a (x-1)), if x >= 2
= a, if x = 1
|| bn_sqr - square a bignum, 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 :: [num]->[num]
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 :: [num]->[[num]]
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 :: [[num]] -> [num]
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 multiply by bn_0.
add_skew2 [] = []
||
|| bn_quot: Long division of a bignum by an integer.
|| Uses Miranda's marvellous infinitely-long integers to do all the hard work.
|| The time it takes is pretty much independent of the value of q.
||
bn_quot (s,e,m) q
= error "bn_quot can only divide by integers", if ~(integer q)
= error "bn_quot cannot divide by 0", if q = 0
= bn_0, if m_is_zero m
= (s,e,m), if q = 1
= (-s,e,m), if q = -1
= bn_normalise (s, e, m_quot (m,q)), if q > 1
= bn_normalise (-s, e, m_quot (m,-q)), if q < -1
m_quot (m,q)
= [], if m_is_zero m
= (a div q) : m_quot (shiftin ((a mod q):x), q), otherwise
where
(a:x) = m
shiftin (a:b:x) = (a * base + b) : x
shiftin [a] = [a * base]
||
|| bn_div: 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.
m_div dend sor
= [], if 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.
|| = 0: m_div dend (0:sor), if sor $m_gt dend
= [dig], if dend2 = [] || exact
= dig : m_div (tl dend2) sor, if hd dend2 = 0
= dig : m_div dend2 (0:sor), otherwise
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)
= m_div_loop' (hack, dend $m_subcarry (m_times hack sor), sor)
, if hack > 2
= m_div_loop' (2, dend $m_subcarry (m_twice sor), sor)
, if hack = 2
= m_div_loop' (0, dend, sor)
, if hack = 0
= m_div_loop' (1, dend $m_subcarry sor, sor)
, if hack = 1
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 turns out to be >10% slower that plain 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)
||
|| ---------------- 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.
bn_sqrt (s,e,m)
= error "bn_sqrt of negative number", if s ~= 1
= bn_normalise (1, e div 2, m_sqrt m) , if e mod 2 = 0
= bn_normalise (1, (e+1) div 2, m_sqrt (0:m)) , if e mod 2 = 1
m_sqrt [] = [] || sqrt(0)
m_sqrt m = m_sqrt' (0, m_sqrbase m)
|| Convert mantissa in base b to mantissa in base b^2.
m_sqrbase (a:b:m) = (a * base + b) : m_sqrbase m
m_sqrbase [a] = [a * base]
m_sqrbase [] = []
||m_sqrt' :: (num,[num]) -> [num]
m_sqrt' (r1,(hd_o1:tl_o1))
= [] , if hd_o1 = 0 & m_is_zero tl_o1
= (r2 mod base) : m_sqrt' (r2 * base, m_sqrt_shiftin (hd_o2:tl_o1))
, otherwise
where
(r2,hd_o2) = m_sqrt_loop (r1,hd_o1)
|| Shift another term of the operand into the calculation.
m_sqrt_shiftin (a:b:x) = (a * base_2 + b) : x
m_sqrt_shiftin [a] = [a * base_2]
m_sqrt_loop :: (num,num) -> (num,num)
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))
||
|| ------------ Number generators ---------------------
||
|| bn_rnd gives a pseudo-random bignum in the range 0 <= bn_rnd x < 1
|| based on the seed, an integer >= 0.
bn_rnd seed
= bn_normalise (1, 0, rnd seed)
where rnd seed
= newseed mod base : rnd newseed
where newseed = (seed * 9197 + 7473) mod 32768
|| A minimally modified version of DAT's edigits
||
|| A better function may be the classic SUM (n <- 0..) (1 / n!)
|| Since running edigits for thousands of digits generates some enormous
|| Miranda integers which make the heap grow indefinitely.
bn_e = bn_make (1, 1, e_digits)
e_digits = (2: e_convert(repeat 1))
e_convert x
= (hd x'):e_convert (tl x')
where x' = e_norm 2 (0:map (10*) x)
e_norm c (d:e:x)
= d + e div c: e' mod c : x', if e mod c + 9 < c
= d + e' div c : e' mod c : x', otherwise
where
(e':x') = e_norm (c+1) (e:x)
|| "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_pi = bn "3.141592653589793238462643383279"
||
|| pi = 4 * atan(1)
|| See commentary under bn_atan2.
bn_pi = bn_times 4 bn_pi_4
bn_pi_4 = bn_atan2 bn_1 (bn "0.4") || pi / 4
|| Old literal constant value, used in benchmarks
bn_pi60 = bn_times 4 (bn ".785398163397448309615660845819875721049292349843776455243736")
||
|| ----------------- Logarithms and trigonometric functions -----------------
||
|| Da 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.
bn_exp x
= bn_1, if bn_is_zero x
= bn_1 $bn_div (bn_exp (bn_neg x)), if bn_is_neg x
= bn_sqr (bn_exp (bn_half x)), if x $bn_gt (bn "0.3")
= (bn_1 $bn_add x) $bn_add (bn_sum_series (map2 bn_quot (exp_top x) exp_bot)), otherwise
exp_top x = [ a | a <- bn_sqr x, bn_mul x a .. ]
exp_bot = [ a | (a,n) <- (2,3), (a*n,n+1) .. ]
|| 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:
|| .5194938532
|| 1.9249505911
||
|| 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 if digits.
|| In all other cases, this strategy turns out to be three or four times slower.
bn_ln x
= error "Can't take log of zero",
if x = bn_0
= error "Can't take log of negative numbers",
if bn_is_neg x
= bn_0,
if x $bn_eq bn_1
= bn_twice (bn_ln (bn_sqrt x)),
if x $bn_lt lower_limit \/ x $bn_gt upper_limit
= bn_ln' x,
otherwise
where
|| Limits within which the summation of the series will work
lower_limit = bn ".5194938533"
upper_limit = bn "1.9249505911"
|| Sum the series
bn_ln' x
= bn_twice (bn_sum_series (map2 bn_quot (ln_top x) [1,3..]))
ln_top x
= [n | n <- a, n $bn_mul asq ..]
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!
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_0, if n = bn_0
= bn_neg (bn_sin' (bn_neg n)), if n $bn_lt bn_0
= bn_sin' n, otherwise
|| 2) reduce to 0..2*pi: for n > 2 * pi, sin(n) = sin(n - 2*pi)
bn_sin' n = bn_sin' (bn_sub n bn_2pi), if n $bn_gt bn_2pi
= bn_sin'' n, otherwise
where bn_2pi = bn_twice bn_pi
|| 3) reduce to 0..pi: for pi < n <= 2*pi, sin(n) = -sin(n - pi)
bn_sin'' n = bn_neg (bn_sin''' (bn_sub n bn_pi)), if n $bn_gt bn_pi
= bn_sin''' n, otherwise
|| 4) reduce to 0..pi/2: for pi/2 < n <= pi, sin(n) = sin(pi - n)
bn_sin''' n = bn_sin'''' (bn_sub bn_pi n), if n $bn_gt bn_pi_2
= bn_sin'''' n, otherwise
where bn_pi_2 = bn_twice bn_pi_4
|| 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 .. ]
where minus_x_squared = bn_neg (bn_mul x x)
|| The divisors of the terms of the series:
|| sin_bot = [ factorial n | n <- [1,3..] ]
sin_bot = [ a | (a,n) <- (1,2), (a*n*(n+1),n+2) .. ]
sin_series x = addpairs (map2 bn_quot (sin_top x) sin_bot)
|| 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_1, if x = bn_0 || avoid bottom-out case
= bn_sin (x $bn_add (bn_half bn_pi)), otherwise
|| 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.
atan_c_lim = bn ".562341325"
bn_atan x = bn_atan2 x (bn ".5")
|| 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 seems a good general purpose value for c - other values run slower
|| except for some special cases.
|| In the case of atan(1), used for pi, 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 arctan of other specific values.
bn_atan2 x c
= error "Internal error: bn_atan2 limit check", if c $bn_gt atan_c_lim
= bn_0, if x = bn_0
= bn_neg (bn_atan2' (bn_neg x) c), if bn_is_neg x
= bn_atan2' x c, otherwise
|| 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)
||
|| We tried using the knowledge that c = 0.5 here to optimise the expressions
|| but it gives a negligable speedup - the time is mostly spent in division
|| and in generating and summing the series - so we keep the general case.
||
|| bn_atan2' knows that its parameter x is not negative.
bn_atan2' :: bignum -> bignum -> bignum
bn_atan2' x c
= bn_atan'' x, if x $bn_le c
= bn_atan'' c $bn_add (bn_atan2' d c) , if d $bn_le c
= bn_twice (bn_atan'' c) $bn_add (bn_atan2' (e x) c), otherwise
where
d = (x $bn_sub c) $bn_div (bn_1 $bn_add (c $bn_mul x))
e x = (etop x) $bn_div (ebot x)
etop x = (x $bn_sub (bn_twice c)) $bn_sub (x $bn_mul csq)
ebot x = bn_1 $bn_add (bn_twice (bn_mul x c)) $bn_sub csq
csq = bn_sqr c
|| This version always sums the series. x must be <= atan_c_lim
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 (map2 bn_quot (sin_top x) [1,3..])
addpairs :: [bignum] -> [bignum]
addpairs (a:b:x) = bn_add a b : addpairs x
addpairs [a] = [a]
addpairs [] = []
|| bn_sum_series sums a converging series of bignum. 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 (hd 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, the
|| calculation will cause a fatal error in bn2ubn.
bn_sum_series :: [bignum] -> bignum
bn_sum_series s
= bn_sum_series (tl s), if bn_is_zero (hd s)
= bn_neg (bn_sum_series (map bn_neg s)), if bn_is_neg (hd s)
= ubn2bn (ubn_sum_series (map bn2ubn s)), otherwise
|| unsigned version - tack a leading zero on so that carry from first digit
|| (if any) has somewhere to go.
ubn_sum_series :: [ubignum] -> ubignum
ubn_sum_series x
= sum_series (map prefix_0 x)
where
prefix_0 (e,m) = (e+1, 0:m)
sum_series :: [ubignum] -> ubignum
sum_series ((ae,am):(be,bm):x)
= (ae, (take nsd am) ++ umantissa_of (sum_series ((ubn_add_no_overflow (ae-nsd, drop nsd am) (be,bm)) : x)))
where
nsd = n_sure_digits (ae,am) be
sum_series [a] = a || Sum of one number is that number!
sum_series [] = ubn_0
|| ubn_add_no_overflow knows that there will be no overflow from the result,
|| and so we don't prefix a zero. It also knows that ae >= be.
ubn_add_no_overflow :: ubignum -> ubignum -> ubignum
ubn_add_no_overflow (ae,am) (be,bm)
|| We don't use m_addcarry here because the args almost certainly
|| have an infinite number of digits.
= (ae, m_carry1 (m_add am (prefix0s (ae-be) bm)))
prefix0s n x = take n [0,0..] ++ x
|| Given an ubignum and the exponent of the next term in the series, work out
|| how many digits of the result depend on just the first term.
|| We can be sure of some digits if the exponent of the first term is greater
|| than that of the second:
|| - None if there is only a difference of one in the exponents (1st digit of
|| 2nd term can cause carry into 1st digit of 1st term)
|| - if the first term has so few digits that there is no overlap with the
|| mantissa of the second AND there's a clear digit between the end of the
|| first and the start of the 2nd:
|| a = 0.123
|| b = 0.0000987
|| then the whole of the first term is sure (any overflow from the sum of b
|| and the rest of the terms of the series would fall into the gap)
|| - If the last digit of a before it hits the mantissa of b is < 9, then all
|| digits before that digit are sure.
|| - If that digit *is* 9 then repeat the test for the previous digit etc etc
n_sure_digits :: ubignum -> num -> num
n_sure_digits (ae,am) be
= 0, if ae-be <= 1
= len_am, if len_am < ae-be
= ae - be - 1, if am ! (ae-be-1) < base_1
= n_sure_digits (ae,am) (be+1), otherwise
where
len_am = len_n am (ae-be)
|| len_n is a length-limited length-of-list function, necessary since we
|| cannot take the length of an infinite list.
len_n :: [*] -> num -> num
len_n (a:x) n
= 1 + (len_n x (n-1)), if n > 0
= 0, otherwise
len_n [] n = 0
||
|| ----------------- Test functions -----------------
||
bn_third = bn_make (1,0,[3,3..])
testaddsub = bn_sub (bn_add bn_e bn_third) bn_e
testmuldiv = bn_mul (bn_div bn_third bn_e) bn_e
bnshow a = "(" ++ show (sign_of a) ++
"," ++ show (exponent_of a) ++
"," ++ show (mantissa_of a) ++ ")"