/* [<][>][^][v][top][bottom][index][help] */
DEFINITIONS
This source file includes following definitions.
- Spigot_calculate_pi
- Spigot_calculate_e
- Scm_Init_spigot
1 /*
2 * spigot.c - calculate pi and e by spigot algorithm
3 *
4 * Written by Shiro Kawai (shiro@acm.org)
5 * I put this program in public domain. Use it as you like.
6 *
7 * $Id: spigot.c,v 1.3 2003/06/29 08:21:07 shirok Exp $
8 */
9
10 #include <math.h>
11 #include <gauche.h>
12 #include <gauche/extend.h>
13 #include "spigot.h"
14
15 /* The spigot algorithm for pi:
16
17 pi = 2 + 1/3*(2 + 2/5*(2 + 3/7*(2 + ... (2 + k/(2k+1)*(2 + ...)) ...))) :[1]
18
19 The last term can be approximated by:
20
21 2 + k/(2k+1)*(4)
22
23 where
24
25 k = n * log2 (10)
26
27 for a desired precision of n decimal digits.
28
29 Now, pi is 3.14159 ... in decimal. That means:
30
31 pi = 3 + 1/10*(1 + 1/10*(4 + 1/10*(1 + 1/10*(5 + 1/10*(9 + ....))))) :[2]
32
33 Note the similality of [1] and [2]. If we think a special variable
34 base system that the k-th digit after point has the weight of k/(2k+1)
35 of the previous digit, [1] means pi can be writen as [2.222222...]
36 in such a base system.
37
38 Base: [1/3, 2/5, 3/7, ... , k/(2k+1), ...]
39 Number: [2 . 222222222 ...]
40
41 Base: [1/10, 1/10, 1/10, ....]
42 Number: [3 . 141592653 ...]
43
44 Now, calculating pi is effectively a base conversion problem.
45
46 The actual calculation process is explained in detail at:
47 http://membres.lycos.fr/bgourevitch/mathematiciens/goutte/goutte.html
48
49 It uses a table, but the rows "A" and "B" are mechanically calculated,
50 and the rows "retenue", "somme", and "reste" can share the same array
51 for those columns are calculated left to right for each iteration.
52 Thus we need to keep only one array to hold intermediate calculation.
53
54 E can be calculated in the same way, starting from the following number:
55
56 Base: [1, 1/2, 1/3, 1/4, ... , 1/(k-1), ...]
57 Number: [1 . 11111111 ...]
58
59 The following code is a very naive implementation of the algorithm.
60 You would get less memory footprint and faster calculation by using
61 big digit, like 10000 instead of 10.
62 */
63
64 ScmObj Spigot_calculate_pi(int digits)
65 {
66 int k, i, j, l, b, q, r, *array;
67 ScmObj rvec, *relts;
68
69 if (digits <= 0) Scm_Error("digits must be a positive integer");
70
71 /* Scheme vector to keep the result */
72 rvec = Scm_MakeVector(digits, SCM_MAKE_INT(0));
73 relts = SCM_VECTOR_ELEMENTS(rvec);
74
75 /* Prepare the array for variable base system */
76 k = (int)floor(digits * 3.3219280948873626);
77 array = SCM_NEW_ATOMIC2(int *, (k+1)*sizeof(int));
78 for (i=0; i<k; i++) array[i] = 2;
79 array[k] = 4;
80
81 for (i=0, b=2; i<digits; i++) {
82 q = 0;
83 for (j=k; j>0; j--) {
84 q += array[j] * 10;
85 array[j] = q % (2*j+1);
86 q /= 2*j+1;
87 q *= j;
88 }
89 r = b + q/10;
90 b = q % 10;
91 /* Here, we have the i-th digit in r.
92 In rare occasions, r becomes more than 10, and we need to back-up
93 to increment the previous digit(s). (It's rarely the case that
94 this back-up cascades for more than one digit). */
95 if (r < 10) {
96 relts[i] = SCM_MAKE_INT(r);
97 } else {
98 relts[i] = SCM_MAKE_INT(r%10);
99 for (l=i-1, r/=10; r && l>=0; l--, r/=10) {
100 r += SCM_INT_VALUE(relts[l]);
101 relts[l] = SCM_MAKE_INT(r%10);
102 }
103 }
104 }
105 return rvec;
106 }
107
108 ScmObj Spigot_calculate_e(int digits)
109 {
110 int k, i, j, l, b, q, r, *array;
111 ScmObj rvec, *relts;
112
113 if (digits <= 0) Scm_Error("digits must be a positive integer");
114
115 /* Scheme vector to keep the result */
116 rvec = Scm_MakeVector(digits, SCM_MAKE_INT(0));
117 relts = SCM_VECTOR_ELEMENTS(rvec);
118
119 /* Prepare the array for variable base system */
120 k = (int)floor(digits * 3.3219280948873626);
121 array = SCM_NEW_ATOMIC2(int *, (k+1)*sizeof(int));
122 for (i=0; i<k; i++) array[i] = 1;
123 array[k] = 2;
124
125 for (i=0, b=1; i<digits; i++) {
126 q = 0;
127 for (j=k; j>0; j--) {
128 q += array[j] * 10;
129 array[j] = q % j;
130 q /= j;
131 }
132 r = b + q/10;
133 b = q % 10;
134 /* Here, we have the i-th digit in r.
135 In rare occasions, r becomes more than 10, and we need to back-up
136 to increment the previous digit(s). (It's rarely the case that
137 this back-up cascades for more than one digit). */
138 if (r < 10) {
139 relts[i] = SCM_MAKE_INT(r);
140 } else {
141 relts[i] = SCM_MAKE_INT(r%10);
142 for (l=i-1, r/=10; r && l>=0; l--, r/=10) {
143 r += SCM_INT_VALUE(relts[l]);
144 relts[l] = SCM_MAKE_INT(r%10);
145 }
146 }
147 }
148 return rvec;
149 }
150
151 /*
152 * Module initialization function.
153 * This is called when spigot.so is dynamically loaded into gosh.
154 */
155 ScmObj Scm_Init_spigot(void)
156 {
157 ScmModule *mod;
158
159 /* Register this DSO to Gauche */
160 SCM_INIT_EXTENSION(spigot);
161
162 /* Create "spigot" module if it doesn't exist yet. */
163 mod = SCM_MODULE(SCM_FIND_MODULE("spigot", TRUE));
164
165 /* Register stub-generated procedures inside "spigot" module */
166 Scm_Init_spigotlib(mod);
167 }
168