l4/pkg/ocaml/ocaml/contrib/otherlibs/num/test/pi_big_int.ml

   1 (* Pi digits computed with the sreaming algorithm given on pages 4, 6
   2    & 7 of "Unbounded Spigot Algorithms for the Digits of Pi", Jeremy
   3    Gibbons, August 2004. *)
   4
   5 open Printf;;
   6 open Big_int;;
   7
   8 let ( !$ ) = Big_int.big_int_of_int
   9 and ( +$ ) = Big_int.add_big_int
  10 and ( *$ ) = Big_int.mult_big_int
  11 and ( =$ ) = Big_int.eq_big_int
  12 ;;
  13
  14 let zero = Big_int.zero_big_int
  15 and one = Big_int.unit_big_int
  16 and three = !$ 3
  17 and four = !$ 4
  18 and ten = !$ 10
  19 and neg_ten = !$(-10)
  20 ;;
  21
  22 (* Linear Fractional (aka M=F6bius) Transformations *)
  23 module LFT = struct
  24
  25   let floor_ev (q, r, s, t) x = div_big_int (q *$ x +$ r) (s *$ x +$ t);;
  26
  27   let unit = (one, zero, zero, one);;
  28
  29   let comp (q, r, s, t) (q', r', s', t') =
  30     (q *$ q' +$ r *$ s', q *$ r' +$ r *$ t',
  31      s *$ q' +$ t *$ s', s *$ r' +$ t *$ t')
  32 ;;
  33
  34 end
  35 ;;
  36
  37 let next z = LFT.floor_ev z three
  38 and safe z n = (n =$ LFT.floor_ev z four)
  39 and prod z n = LFT.comp (ten, neg_ten *$ n, zero, one) z
  40 and cons z k =
  41   let den = 2 * k + 1 in
  42   LFT.comp z (!$ k, !$(2 * den), zero, !$ den)
  43 ;;
  44
  45 let rec digit k z n row col =
  46   if n > 0 then
  47     let y = next z in
  48     if safe z y then
  49       if col = 10 then (
  50         let row = row + 10 in
  51         printf "\t:%i\n%s" row (string_of_big_int y);
  52         digit k (prod z y) (n - 1) row 1
  53       )
  54       else (
  55         print_string(string_of_big_int y);
  56         digit k (prod z y) (n - 1) row (col + 1)
  57       )
  58     else digit (k + 1) (cons z k) n row col
  59   else
  60     printf "%*s\t:%i\n" (10 - col) "" (row + col)
  61 ;;
  62
  63 let digits n = digit 1 LFT.unit n 0 0
  64 ;;
  65
  66 let usage () =
  67   prerr_endline "Usage: pi_big_int <number of digits to compute for pi>";
  68   exit 2
  69 ;;
  70
  71 let main () =
  72   let args = Sys.argv in
  73   if Array.length args <> 2 then usage () else
  74   digits (int_of_string Sys.argv.(1))
  75 ;;
  76
  77 main ()
  78 ;;