[Shootout-list] Improved ocaml pidigits... maybe.

[Will M. Farr]

--Apple-Mail-1-493999647
Content-Transfer-Encoding: 7bit
Content-Type: text/plain; charset=US-ASCII; delsp=yes; format=flowed

Hello,

I borrowed a few ideas from William Douglas Neumann (see  
http://lists.alioth.debian.org/pipermail/shootout-list/2004-December/ 
000926.html ) for making my personal version of pidigits in ocaml  
faster.  Using the create_ratio and floor_ratio sped my code up by  
factors of *thousands* (gcd is very slow, and apparently the // routine  
from Num uses gcd to reduce the fraction, while make_ratio doesn't).   
I'm glad WDN noticed this, or I would have quit while taking hours to  
compute thousands of digits!

The implementation below runs roughly two times faster on my computer  
than the one posted from WDN, BUT it uses a suggestion for optimization  
found in the text, but not the pseudocode, of the referenced paper.   
During the computation, I reduce the linear fractional transform by  
eliminating common factors from the components (I found that the  
optimum seemed to be to do this about 4 times through the course of the  
computation -- this seems strange to me, but it really did work best on  
my computer).  Technically, this is not the same algorithm as in the  
pseudocode from the paper, but I'm submitting it anyway because the  
optimization is mentioned in the paper.

Will

let ( */ ) = Num.( */ );;
let ( +/ ) = Num.( +/ );;
let ( // ) = Num.( // );;
let ( -/ ) = Num.( -/ );;

type lft = {
   q : Big_int.big_int;
   r : Big_int.big_int;
   (* s = 0 always! *)
   t : Big_int.big_int
}

let unit = {q = Big_int.unit_big_int;
             r = Big_int.zero_big_int;
             t = Big_int.unit_big_int}

let comp = function {q = q1; r = r1; t = t1} ->
   function {q = q2; r = r2; t = t2} ->
     {q = Big_int.mult_big_int q1 q2;
      r = Big_int.add_big_int
            (Big_int.mult_big_int q1 r2)
            (Big_int.mult_big_int r1 t2);
      t = Big_int.mult_big_int t1 t2};;

let apply = function {q = qq; r = rr; t = tt} ->
   fun x ->
     let numer = Big_int.add_big_int
                   (Big_int.mult_int_big_int x qq) rr in
       Ratio.create_ratio numer tt;;

let pi_series_term i =
   Some {q = Big_int.big_int_of_int (i+1);
         r = Big_int.big_int_of_int (4*(i+1)+2);
         t = Big_int.big_int_of_int (2*(i+1)+1)};;

let extract_int l x =
   Big_int.int_of_big_int (Ratio.floor_ratio (apply l x));;

let remove_digit l n =
   comp {q = Big_int.big_int_of_int 10;
         r = Big_int.big_int_of_int (~-n*10);
         t = Big_int.unit_big_int}
     l;;

let lcd a b =
   let g = Big_int.gcd_big_int a b in
     Big_int.mult_big_int g (Big_int.mult_big_int
                               (Big_int.div_big_int a g)
                               (Big_int.div_big_int b g));;

(** Uses the fact that s := 0 always to reduce the lft. *)
let reduce = function {q = qq; r = rr; t = tt} ->
   let gcd1 = Big_int.gcd_big_int qq tt and
     gcd2 = Big_int.gcd_big_int rr tt in
   let q1 = Big_int.div_big_int qq gcd1 and
     t1 = Big_int.div_big_int tt gcd1 and
     r2 = Big_int.div_big_int rr gcd2 and
     t2 = Big_int.div_big_int tt gcd2 in
   let denom = lcd t1 t2 in
     {q = Big_int.mult_big_int q1 (Big_int.div_big_int denom t1);
      r = Big_int.mult_big_int r2 (Big_int.div_big_int denom t2);
      t = denom};;

let make_pi_stream n =
   let reduce_mod = n / 4 in
   let series_state = ref unit in
   let current_digit = ref ~-1 in
   let pi_series_stream = Stream.from pi_series_term in
   let next_state () =
     series_state := comp !series_state (Stream.next pi_series_stream)
   in
     fun i ->
       if i = !current_digit + 1 then
         begin
           current_digit := i;
           next_state ();
           if i mod reduce_mod = 0 then
             series_state := reduce !series_state;
           let rec loop three_value =
             if three_value = (extract_int !series_state 4) then
               begin
                 series_state := remove_digit !series_state three_value;
                 Some three_value
               end
             else
               begin
                 next_state ();
                 loop (extract_int !series_state 3)
               end
           in
             loop (extract_int !series_state 3)
         end
       else
         None;;


let pi_digits n =
   let pi_stream = Stream.from (make_pi_stream n) in
   let n_mod = n / 10 in
     for i = 1 to n_mod do
       for j = 1 to 10 do
         Printf.printf "%1d" (Stream.next pi_stream)
       done;
       Printf.printf "\t:%d\n" (10*i)
     done;
     let n_rem = n mod 10 in
       if n_rem <> 0 then
         begin
           for i = 1 to n_rem do
             Printf.printf "%1d" (Stream.next pi_stream)
           done;
           Printf.printf "\t:%d\n" n
         end

let _ = pi_digits (int_of_string Sys.argv.(1));;

--Apple-Mail-1-493999647
content-type: application/pgp-signature; x-mac-type=70674453;
	name=PGP.sig
content-description: This is a digitally signed message part
content-disposition: inline; filename=PGP.sig
content-transfer-encoding: 7bit

-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.4.0 (Darwin)

iD8DBQFB5KIVjFCrhUweU3MRAlteAKCetD4RAcVlqqyckkQGmu9m1ZyNTACdHsZa
xGtHVxM6VdxU2l2BamyYrao=
=EnMA
-----END PGP SIGNATURE-----

--Apple-Mail-1-493999647--