Num.quo_num

[Christophe TROESTLER <Christophe.Troestler@umh.ac.be>]

[-- Attachment #1: Type: Text/Plain, Size: 958 bytes --]

Hi,

Attached are two implementations of the algorithm for finding the
digits of pi described here:
http://web.comlab.ox.ac.uk/oucl/work/jeremy.gibbons/publications/spigot.pdf
The two implementations differ only by the fact that the first one
uses Big_int while the second one uses Num.  However, the efficiency
figures are rather different (here for 100 digits) :

Big_int
  real    0m0.009s
  user    0m0.007s
  sys     0m0.003s

Num
  real    0m0.487s
  user    0m0.484s
  sys     0m0.000s

Upon investigation, it appears that the culpright is the following:

  let quo_num x y = floor_num (div_num x y)

Indeed "Num.div_num" calls "Num.num_of_ratio" wich in turn calls
"Ratio.normalize_ratio" which I guess accounts for the huge time
difference.  Now, "Num.normalize_ratio" is actually useless for
"quo_num", since it does "num_of_big_int" anyway.  What about writing
a specialized and vastly more efficient implementation for "quo_num"?

Regards,
ChriS

[-- Attachment #2: pidigits.ml --]
[-- Type: Text/Plain, Size: 1498 bytes --]

(* Pi digits computed with the sreaming algorithm given on pages 4, 6
   & 7 of "Unbounded Spigot Algorithms for the Digits of Pi", Jeremy
   Gibbons, August 2004.  *)

open Printf
open Big_int

let ( !$ ) = Big_int.big_int_of_int
let ( +$ ) = Big_int.add_big_int
let ( *$ ) = Big_int.mult_big_int
let ( =$ ) = Big_int.eq_big_int
let zero = Big_int.zero_big_int
and one  = Big_int.unit_big_int
and three = !$ 3
and four = !$ 4
and ten  = !$ 10
and neg_ten = !$(-10)

(* Linear Fractional (aka Möbius) Transformations *)
module LFT =
struct
  let floor_ev (q,r,s,t) x = div_big_int (q *$ x +$ r) (s *$ x +$ t)

  let unit = (one, zero, zero, one)

  let comp (q,r,s,t) (q',r',s',t') =
    (q *$ q' +$ r *$ s',  q *$ r' +$ r *$ t',
     s *$ q' +$ t *$ s',  s *$ r' +$ t *$ t')
end

let next z = LFT.floor_ev z three
and safe z n = (n =$ LFT.floor_ev z four)
and prod z n = LFT.comp (ten, neg_ten *$ n, zero, one) z
and cons z k =
  let den = 2*k+1 in LFT.comp z (!$ k, !$(2*den), zero, !$ den)

let rec digit k z n row col =
  if n > 0 then
    let y = next z in
    if safe z y then
      if col = 10 then (
	let row = row + 10 in
	printf "\t:%i\n%s" row (string_of_big_int y);
	digit k (prod z y) (n-1) row 1
      )
      else (
	print_string(string_of_big_int y);
	digit k (prod z y) (n-1) row (col+1)
      )
    else digit (k+1) (cons z k) n row col
  else
    printf "%*s\t:%i\n" (10 - col) "" (row + col)

let digits n = digit 1 LFT.unit n 0 0

let () =
  digits(int_of_string Sys.argv.(1))

[-- Attachment #3: pidigits-num.ml --]
[-- Type: Text/Plain, Size: 1375 bytes --]

(* Pi digits computed with the sreaming algorithm given on pages 4, 6
   & 7 of "Unbounded Spigot Algorithms for the Digits of Pi", Jeremy
   Gibbons, August 2004.  *)

open Printf
open Num

let zero = num_of_int 0
and one  = num_of_int 1
and three = num_of_int 3
and four = num_of_int 4
and ten  = num_of_int 10
and neg_ten = num_of_int(-10)

(* Linear Fractional Transformation *)
module LFT =
struct
  let floor_ev (q,r,s,t) x = quo_num (q */ x +/ r) (s */ x +/ t)

  let unit = (one, zero, zero, one)

  let comp (q,r,s,t) (q',r',s',t') =
    (q */ q' +/ r */ s',  q */ r' +/ r */ t',
     s */ q' +/ t */ s',  s */ r' +/ t */ t')
end

let next z = LFT.floor_ev z three
and safe z n = (n =/ LFT.floor_ev z four)
and prod z n = LFT.comp (ten, neg_ten */ n, zero, one) z
and cons z k =
  let den = 2*k+1 in
  LFT.comp z (num_of_int k, num_of_int(2*den), zero, num_of_int den)

let rec digit k z n row col =
  if n > 0 then
    let y = next z in
    if safe z y then
      if col = 10 then (
	let row = row + 10 in
	printf "\t:%i\n%s" row (string_of_num y);
	digit k (prod z y) (n-1) row 1
      )
      else (
	print_string(string_of_num y);
	digit k (prod z y) (n-1) row (col+1)
      )
    else digit (k+1) (cons z k) n row col
  else
    printf "%*s\t:%i\n" (10 - col) "" (row + col)

let digits n = digit 1 LFT.unit n 0 0

let () =
  digits(int_of_string Sys.argv.(1))