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
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
|
(*s Test program for [Creal]. *)
open Printf
open Gmp
open Creal
(*s Options *)
let prec = ref 50
let display = ref true
let sanity_check = ref true
let _ =
Arg.parse
["-p", Arg.Int ((:=) prec), "-p n set the precision";
"-silent", Arg.Clear display, "-silent no display";
"-check", Arg.Set sanity_check, "-check only sanity checks" ]
(fun s -> raise (Arg.Bad ("unknown option " ^ s)))
"test [-p prec] [silent]"
(*s Sanity checks. Compare two numbers up to the precision. *)
let _ =
if !sanity_check then begin
printf "*** Sanity checks ***\n\n"; flush stdout
end
let check msg x y =
if !sanity_check then begin
printf "%s... " msg; flush stdout;
let delta = Z.sub (approx x !prec) (approx y !prec) in
if Z.cmp_si (Z.abs delta) 1 <= 0 then
printf "ok\n\n"
else begin
printf "FAILED!\n\n"; exit 1
end;
flush stdout
end
let sqrt_2 = sqrt two
let _ = check "sqrt(2)^2 = 2" (sqrt_2 *! sqrt_2) two
let _ = check "1/sqrt(2) = sqrt(2)/2" (inv sqrt_2) (sqrt_2 /! two)
let sqrt_3 = sqrt (of_int 3)
let _ = check "1 = (sqrt(3) + sqrt(2)) * (sqrt(3) - sqrt(2))"
one ((sqrt_3 +! sqrt_2) *! (sqrt_3 -! sqrt_2))
let _ = check "(sqrt(2) ^ sqrt(2)) ^ sqrt(2) = 2"
(pow (pow sqrt_2 sqrt_2) sqrt_2) two
let one_third = of_q (Q.from_ints 1 3)
let root3 x = pow x one_third
let _ = check "54^1/3 - 2^1/3 = 16^1/3"
(root3 (of_int 54) -! root3 two) (root3 (of_int 16))
let _ = check "cos(0)=1" (cos zero) one
let _ = check "cos(pi/2)=0" (cos pi_over_2) zero
let _ = check "sin(0)=0" (sin zero) zero
let _ = check "sin(pi/2)=1" (sin pi_over_2) one
let pi_over_4 = pi /! (of_int 4)
let square x = x *! x
let _ = check "cos^2(pi/4) + sin^2(pi/4) = 1"
(square (cos pi_over_4) +! square (sin pi_over_4)) one
let _ = check "tan(pi/4) = 1" (tan pi_over_4) one
let _ = check "pi/4 = 4arctan(1/5) - arctan(1/239)" pi_over_4
(of_int 4 *! arctan_reciproqual 5 -! arctan_reciproqual 239)
let _ = check "ln(1) = 0" (ln one) zero
let _ = check "ln(e) = 1" (ln e) one
let _ = check "ln(pi*pi) = 2ln(pi)" (ln (square pi)) (two *! ln pi)
let _ = check "exp(-pi) = exp(-pi/2) * exp(-pi/2)"
(exp (neg pi)) (let y = exp (neg pi_over_2) in y *! y)
let _ = if !sanity_check then exit 0
(*s Benchmark. *)
(* Test function: display the real number, if not [silent] ; otherwise,
just compute the approximation (for timings). *)
let _ = printf "\n*** Benchmarks ***\n\n"; flush stdout
let test msg beautiful x =
if !display then begin
printf "%s = " msg; flush stdout;
printf "%s\n\n"
(if beautiful then to_beautiful_string x !prec else to_string x !prec);
flush stdout
end else begin
printf "%s\n" msg;
flush stdout;
ignore (approx x !prec)
end
(*s golden ratio *)
let phi = (one +! sqrt (of_int 5)) /! (of_int 2)
let _ = test "golden ratio" true phi
(* e (predefined in [Creal]) *)
let _ = test "e" true e
(* pi (predefined in [Creal]) *)
let _ = test "pi" true pi
(*s The Exact Arithmetic Competition: Level 0 Tests
http://www.cs.man.ac.uk/arch/dlester/arithmetic/level0t.html *)
(* sqrt(pi) *)
let _ = test "sqrt(pi)" false (sqrt pi)
(* sin(exp(1)) *)
let _ = test "sin(e)" false (sin e)
(* cos(exp(1)) *)
let _ = test "cos(e)" false (cos e)
(* sin(sin(sin(1))) *)
let x = sin (sin (sin one))
let _ = test "sin(sin(sin(1)))" false x
(* cos(cos(cos(1))) *)
let x = cos (cos (cos one))
let _ = test "cos(cos(cos(1)))" false x
(* exp(exp(exp(1))) *)
let x = exp (exp (exp one))
let _ = test "exp(exp(exp(1)))" false x
(* log(pi) *)
let _ = test "ln(pi)" false (ln pi)
(* log(1+log(1+log(1+pi))) *)
let ln_ln_ln_pi = ln (one +! ln (one +! ln (one +! pi)))
let _ = test "ln(1+ln(1+ln(1+pi)))" false ln_ln_ln_pi
(* log(1+log(1+log(1+exp(1)))) *)
let ln_ln_ln_e = ln (one +! ln (one +! ln (one +! e)))
let _ = test "ln(1+ln(1+ln(1+e)))" false ln_ln_ln_e
(*i
(* log(1+log(1+log(1+log(1+log(1+log(1+pi)))))) *)
let x = ln (one +! ln (one +! ln (one +! ln_ln_ln_pi)))
let _ = test "ln(1+ln(1+ln(1+ln(1+ln(1+ln(1+pi))))))" false x
(* log(1+log(1+log(1+log(1+log(1+log(1+exp(1))))))) *)
let x = ln (one +! ln (one +! ln (one +! ln_ln_ln_e)))
let _ = test "ln(1+ln(1+ln(1+ln(1+ln(1+ln(1+e))))))" false x
i*)
(* sin(1e50) *)
let ten_to_50 = pow_int (of_int 10) 50
let x = sin ten_to_50
let _ = test "sin(1e50)" false x
(* cos(1e50) *)
let x = cos ten_to_50
let _ = test "cos(1e50)" false x
(* arctan(1) *)
let _ = test "arctan(1)" false (arctan one)
(*i
(* BUG GMP 2 *)
let q =
Q.from_zs (Z.from_int 1) (Z.from_string "19807040628566084398385987584" 10)
let _ = Q.add q (Q.from_ints 1 2)
(* BUG GMP 3 *)
let q = Q.from_zs (Z.from_string "112803124130337404998606757686274889113032882986303222429756948481" 10) (Z.from_string "5192296858534827628530496329220096" 10)
let q' = Q.add q (Q.from_ints 1 2)
let _ = Z.fdiv_q (Q.get_num q') (Q.get_den q')
let time f x =
let old = Sys.time () in
let y = f x in
Printf.printf "%f\n" (Sys.time () -. old);
y
;;
i*)
|
