1
+ -- let's compute pi
2
+
3
+ import prelude;
4
+
5
+ -- Stores numbers in a base-100 representation so they are easy to print out
6
+ -- returns whether the result is not zero
7
+ small_div {n} (out dest :: u8[n]) (num :: u8[n]) (div :: u32) :: bool :=
8
+ ( r :: u32 := 0;
9
+ is_zero := True;
10
+ for i in n -> (
11
+ s := r * 100 + num[i];
12
+ dest[i] <- s / div;
13
+ r <- s % div;
14
+ is_zero <- is_zero and dest[i] == 0;
15
+ );
16
+ return not is_zero;
17
+ );
18
+
19
+ bigsum {n} (out dest :: u8[n]) (a :: u8[n]) (b :: u8[n]) :: () :=
20
+ ( carry :: u32 := 0;
21
+ for i in n-1..0 : -1 -> (
22
+ s := a[i] + b[i] + carry;
23
+ dest[i] <- s % 100;
24
+ carry <- s / 100;
25
+ );
26
+ );
27
+ bigdiff {n} (out dest :: u8[n]) (a :: u8[n]) (b :: u8[n]) :: () :=
28
+ ( borrow :: u32 := 1;
29
+ for i in n-1..0 : -1 -> (
30
+ d := a[i] + 99 - b[i] + borrow;
31
+ dest[i] <- d % 100;
32
+ borrow <- d / 100;
33
+ );
34
+ );
35
+
36
+ mkInt {N} (z :: u8) :: u8[N] :=
37
+ vec i in : -> if i == 0 then z else 0;
38
+
39
+ -- Computes 2N digits after the decimal point
40
+ compute (N :: u32) :: () :=
41
+ ( num :: u8[N + 2];
42
+ powers :: u8[N + 2];
43
+ scratch :: u8[N + 2];
44
+
45
+ -- Compute 16arctan(1/5)
46
+ powers <- mkInt 16;
47
+ small_div (out powers) powers 5;
48
+ num <- powers;
49
+
50
+ divisor := 3;
51
+ while (small_div (out powers) powers 25
52
+ and small_div (out scratch) powers divisor) -> (
53
+ if divisor % 4 == 3 then
54
+ bigdiff (out num) num scratch
55
+ else
56
+ bigsum (out num) num scratch;
57
+ divisor <- divisor + 2;
58
+ );
59
+
60
+ -- Compute 4arctan(1/239)
61
+ powers <- mkInt 4;
62
+ small_div (out powers) powers 239;
63
+ bigdiff (out num) num powers;
64
+
65
+ divisor <- 3;
66
+ while (small_div (out powers) powers (239*239)
67
+ and small_div (out scratch) powers divisor) -> (
68
+ if divisor % 4 == 3 then
69
+ bigsum (out num) num scratch
70
+ else
71
+ bigdiff (out num) num scratch;
72
+ divisor <- divisor + 2;
73
+ );
74
+
75
+ printf "%01u." num[0];
76
+ -- the last two digits are probably wrong, so don't print them
77
+ for i in 1:N -> printf "%02u" num[i];
78
+ void $ printf "\n";
79
+ );
80
+
81
+ main () :: int :=
82
+ ( compute 500;
83
+ return 0;
84
+ );
0 commit comments