74 lines
1.6 KiB
OCaml
74 lines
1.6 KiB
OCaml
|
|
(* 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 usage () =
|
|
prerr_endline "Usage: pi_num <number of digits to compute for pi>";
|
|
exit 2
|
|
;;
|
|
|
|
let main () =
|
|
let args = Sys.argv in
|
|
if Array.length args <> 2 then usage () else
|
|
digits (int_of_string Sys.argv.(1))
|
|
;;
|
|
|
|
main ()
|
|
;;
|