|| 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) ++ ")"