1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
|
(*************************************************************************)
(* *)
(* OCaml *)
(* *)
(* Pierre Weis, projet Estime, INRIA Rocquencourt *)
(* *)
(* Copyright 2008 Institut National de Recherche en Informatique et *)
(* en Automatique. All rights reserved. This file is distributed *)
(* under the terms of the Q Public License version 1.0. *)
(* *)
(*************************************************************************)
(* 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 ()
;;
|
