root/examples/spigot/spigot.c

/* [<][>][^][v][top][bottom][index][help] */

DEFINITIONS

This source file includes following definitions.
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 + ...)) ...)))  :
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 + ....)))))    :
32
33   Note the similality of  and .  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,  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

/* [<][>][^][v][top][bottom][index][help] */