BigNums in Caml Light 0.5

BigNums in Caml Light 0.5

[HEDDEN]

Dear Caml Light users,

I would like to see some form of "arbitrary precision" math or
"big nums" added to Caml Light.  Below is a very preliminary
attempt at such an implementation written in Caml Light.  At this
point, I am not so much concerned with speed as I am with
obtaining the capabilities. 

I am using vectors and a "base 10000" implementation to represent
floating point numbers.  So far I have written routines for +, -
and *.  (They were easy.)  I would like to expand it to include
division, mod/remainder, gcd, sqrt, exp, log, trig functions,
etc.

If anyone would like to help me in this effort, I'd appreciate 
it.  If anyone can supply/recommend algorithms, that would be 
appreciated, too.  Alternate implementations and improvements 
will be accepted, as well.  In short, I am completely open to 
suggestion.  Please help.

I will probably need to stick with a Caml Light implementation 
because the 'C' interface is only available on Unix, and I am 
working on an Amiga.  Besides, this would make the code useable 
by all.

Jerry D. Hedden
hedden@esdsdf.dnet.ge.com

====================== Cut Here ===============================


(* Bignums module *)
(* Functions implemented so far:

        renormalize : bignum -> bignum = <fun>
          Internal -- used to "fixup" the result of a calculation

        add_big : bignum -> bignum -> bignum = <fun>
        sub_big : bignum -> bignum -> bignum = <fun>
        mult_big : bignum -> bignum -> bignum = <fun>
          Basic operations

        string_of_bignum : bignum -> string = <fun>
        bignum_of_string : string -> bignum = <fun>
          Conversion functions

        big_fib : int -> bignum = <fun>
        big_fact : int -> bignum = <fun>
          Example functions that use bignums
*)

type bignum = {Neg:bool; Point:int; Value:int vect};;

let rec renormalize {Neg=neg;Point=pt;Value=v} =
  let vl = vect_length v in
  begin
    for i = (pred vl) downto 1 do
      let x = v.(i) in
      if x >= 10000
         then (v.(i)      <- x mod 10000;
               v.(pred i) <- v.(pred i) + x/10000)
         else
      if x < 0
         then (v.(i)      <- 10000 + (x mod 10000);
               v.(pred i) <- v.(pred i) - (succ (x/10000)))
    done;
    if v.(0) < 0
       then renormalize {Neg=not neg; Point=pt;
                         Value=(map_vect (function a -> -a) v)}
       else let first =
              let ff =
                let rec find_first n =
                  if n < vl
                     then if (eq_int v.(n) 0)
                             then find_first (succ n)
                             else n
                     else n
                in find_first 0
              in min ff pt
            and last =
              let ll =
                let rec find_last n =
                  if n >= 0
                     then if (eq_int v.(n) 0)
                             then find_last (pred n)
                             else n
                     else n
                in find_last (pred vl)
              in max ll pt
            in if first >= vl
                  then {Neg=false;Point=0;Value=[|0|]}
                  else if (eq_int first 0) & (eq_int last (pred vl))
                          then {Neg=neg;Point=pt;Value=v}
                          else let vvl = succ (last-first) in
                               let vv = make_vect vvl 0 in
                               begin
                                 blit_vect v first vv 0 vvl;
                                 {Neg=neg;Point=(pt-first);Value=vv}
                               end
  end
;;


let add_big {Neg=nx;Point=px;Value=vx} {Neg=ny;Point=py;Value=vy} =
  let xl = vect_length vx
  and yl = vect_length vy
  and dpt = succ (max px py) in
  let len = dpt + (max (xl-px) (yl-py)) in
  let xx = make_vect len 0
  and yy = make_vect len 0
  in begin
    blit_vect (if nx then (map_vect (function a -> -a) vx)
                        else vx)
              0 xx (dpt-px) xl;
    blit_vect (if ny then (map_vect (function a -> -a) vy)
                        else vy)
              0 yy (dpt-py) yl;
    renormalize {Neg=false; Point=dpt;
                 Value=vect_of_list (map2 (prefix +) (list_of_vect xx)
                                                     (list_of_vect yy))}
  end
;;


let sub_big {Neg=nx;Point=px;Value=vx} {Neg=ny;Point=py;Value=vy} =
  let xl = vect_length vx
  and yl = vect_length vy
  and dpt = succ (max px py) in
  let len = dpt + (max (xl-px) (yl-py)) in
  let xx = make_vect len 0
  and yy = make_vect len 0
  in begin
    blit_vect (if nx then (map_vect (function a -> -a) vx)
                        else vx)
              0 xx (dpt-px) xl;
    blit_vect (if ny then (map_vect (function a -> -a) vy)
                        else vy)
              0 yy (dpt-py) yl;
    renormalize {Neg=false; Point=dpt;
                  Value=vect_of_list (map2 (prefix -) (list_of_vect xx)
                                                      (list_of_vect yy))}
  end
;;


let mult_big {Neg=nx;Point=px;Value=vx} {Neg=ny;Point=py;Value=vy} =
  let xl = vect_length vx
  and yl = vect_length vy
  in
  let z = make_vect (xl+yl) 0
  in begin
    for j = (yl-1) downto 0 do
      for i = (xl-1) downto 0 do
        z.(i+j+1) <- z.(i+j+1) + (vx.(i) * vy.(j))
      done;
      for i = (xl+j) downto (j+1) do
        let q = z.(i) in
        if q >= 10000
           then (z.(i)      <- q mod 10000;
                 z.(pred i) <- z.(pred i) + q/10000)
      done
    done;
    renormalize {Neg=not(nx=ny);Point=(px+py+1);Value=z}
  end
;;


let string_of_bignum {Neg=neg;Point=pt;Value=v} =
  let vl = vect_length v
  and n4 n =
    let s = string_of_int n in
    (make_string (4-(string_length s)) `0`) ^ s
  in
  let rec no_trail s =
    let ln = string_length s in
    if (eq_int ln 0) then ""
                     else if (nth_char s (pred ln)) = `0`
                             then no_trail (sub_string s 0 (pred ln))
                             else s
  and nstr a b =
    if a<=b then (n4 v.(a)) ^ (nstr (succ a) b)
            else ""
  in
  let sign = if neg then "-" else ""
  and point = if (succ pt) < vl then "." else ""
  and last = if (succ pt) < vl
                then no_trail (n4 v.(pred vl))
                else ""
  in sign ^ (string_of_int v.(0)) ^ (nstr 1 pt) ^
              point ^ (nstr (succ pt) (vl-2)) ^ last
;;


let bignum_of_string s =
  let last = pred (string_length s)
  and neg,first =
    match nth_char s 0 with
        `-` -> true,1
      | `+` -> false,1
      | _   -> false,0
  in
  let dp = find_dp 0
    where rec find_dp n =
      if n > last
         then n
         else if (nth_char s n) = `.`
                 then n
                 else find_dp (succ n)
  in
  let asize = let size = (3+dp-first)/4 in
              if (eq_int size 0) then 1 else size
  and bsize = (3+last-dp)/4
  in
  let x = make_vect (asize+bsize) 0 in
  let rec set_num n y =
    if (n-4) <= first
       then  x.(y) <- int_of_string (sub_string s first (n-first))
       else (x.(y) <- int_of_string (sub_string s (n-4) 4);
             set_num (n-4) (pred y))
  and set_dec n y =
    if (n+3) >= last
       then  x.(y) <- int_of_string (sub_string (s^"000") n 4)
       else (x.(y) <- int_of_string (sub_string s n 4);
             set_dec (n+4) (succ y))
  in
  begin
    (if dp > first then set_num dp (pred asize));
    (if dp < last then set_dec (succ dp) asize);
    {Neg=neg;Point=(pred asize);Value=x}
  end
;;


let big_fib = fib_rec (bignum_of_string "0") (bignum_of_string "1")
  where rec fib_rec a b n = if n<=0 then a else fib_rec b (add_big a b) (pred n)
;;


let big_fact = fact_rec (bignum_of_string "1")
  where rec fact_rec ans n =
    if n<=0
       then ans
       else fact_rec (mult_big ans (bignum_of_string (string_of_int n))) (pred n)
;;

BigNums in Caml Light 0.5

[Vale'rie Me'nissier-Morain]

You want to implement bigfloats but you seems to mix big integers
features (Fibonacci numbers, fact, mod and gcd) with really floating
point/real features (sqrt, exp, log, trig functions).

It is classical to work in a 2^b-basis (generally b=32). It costs more
to read or write a number but pure computation is much faster.

Here is some bibliography references which may interest you even if
they refer essentially to classical limited floating point numbers:

  David Goldberg "What every computer scientist should know about
  floating point arithmetic", Computing Surveys, 23 (1), March 1991.

  In proceedings of SIGPLAN '90 PLDI, there are two papers by Clinger,
  Steele and another one about how to read and write floating point
  numbers accurately

Otherwise, there exist arbitrarily large rational numbers in Caml V3.1
(the strong version on Unix) (numerator and denominator with less than
79803 decimal digits) with exact computation. In the other hand, we
have only built-in floating point numbers. It will be intesting to
have arbitrary large floating point numbers arithmetic . Furthermore,
it will be interesting to extend rational arithmetic to transcendantal
functions with the required precision as additional argument. 


------------------- Vale'rie ME'NISSIER - MORAIN -----------------------------
INRIA Projet Formel			        Tel: (33 1) (16 1) 39 63 55 98
Domaine de Voluceau		                Fax: (33 1) (16 1) 39 63 53 30
78153 Rocquencourt			     e-mail: menissie@margaux.inria.fr
BP 105
------------------------------------------------------------------------------