[Shootout-list] Improved ocaml pidigits... maybe.
Will M. Farr
farr@MIT.EDU
Tue, 11 Jan 2005 23:05:34 -0500
--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--