"Fossies" - the Fresh Open Source Software Archive 
Member "mapm_4.9.5a/mapm_cpi.c" (21 Feb 2010, 4278 Bytes) of package /linux/misc/old/mapm-4.9.5a.tar.gz:
As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style:
standard) with prefixed line numbers and
code folding option.
Alternatively you can here
view or
download the uninterpreted source code file.
For more information about "mapm_cpi.c" see the
Fossies "Dox" file reference documentation.
1
2 /*
3 * M_APM - mapm_cpi.c
4 *
5 * Copyright (C) 1999 - 2007 Michael C. Ring
6 *
7 * Permission to use, copy, and distribute this software and its
8 * documentation for any purpose with or without fee is hereby granted,
9 * provided that the above copyright notice appear in all copies and
10 * that both that copyright notice and this permission notice appear
11 * in supporting documentation.
12 *
13 * Permission to modify the software is granted. Permission to distribute
14 * the modified code is granted. Modifications are to be distributed by
15 * using the file 'license.txt' as a template to modify the file header.
16 * 'license.txt' is available in the official MAPM distribution.
17 *
18 * This software is provided "as is" without express or implied warranty.
19 */
20
21 /*
22 * $Id: mapm_cpi.c,v 1.4 2007/12/03 01:34:29 mike Exp $
23 *
24 * This file contains the PI related functions.
25 *
26 * $Log: mapm_cpi.c,v $
27 * Revision 1.4 2007/12/03 01:34:29 mike
28 * Update license
29 *
30 * Revision 1.3 2002/11/05 23:10:14 mike
31 * streamline the PI AGM algorithm
32 *
33 * Revision 1.2 2002/11/03 21:56:21 mike
34 * Updated function parameters to use the modern style
35 *
36 * Revision 1.1 2001/03/25 21:01:53 mike
37 * Initial revision
38 */
39
40 #include "m_apm_lc.h"
41
42 /****************************************************************************/
43 /*
44 * check if our local copy of PI is precise enough
45 * for our purpose. if not, calculate PI so it's
46 * as precise as desired, accurate to 'places' decimal
47 * places.
48 */
49 void M_check_PI_places(int places)
50 {
51 int dplaces;
52
53 dplaces = places + 2;
54
55 if (dplaces > MM_lc_PI_digits)
56 {
57 MM_lc_PI_digits = dplaces + 2;
58
59 /* compute PI using the AGM (see right below) */
60
61 M_calculate_PI_AGM(MM_lc_PI, (dplaces + 5));
62
63 m_apm_multiply(MM_lc_HALF_PI, MM_0_5, MM_lc_PI);
64 m_apm_multiply(MM_lc_2_PI, MM_Two, MM_lc_PI);
65 }
66 }
67 /****************************************************************************/
68 /*
69 * Calculate PI using the AGM (Arithmetic-Geometric Mean)
70 *
71 * Init : A0 = 1
72 * B0 = 1 / sqrt(2)
73 * Sum = 1
74 *
75 * Iterate: n = 1...
76 *
77 *
78 * A = 0.5 * [ A + B ]
79 * n n-1 n-1
80 *
81 *
82 * B = sqrt [ A * B ]
83 * n n-1 n-1
84 *
85 *
86 *
87 * C = 0.5 * [ A - B ]
88 * n n-1 n-1
89 *
90 *
91 * 2 n+1
92 * Sum = Sum - C * 2
93 * n
94 *
95 *
96 * At the end when C is 'small enough' :
97 * n
98 *
99 * 2
100 * PI = 4 * A / Sum
101 * n+1
102 *
103 * -OR-
104 *
105 * 2
106 * PI = ( A + B ) / Sum
107 * n n
108 *
109 */
110 void M_calculate_PI_AGM(M_APM outv, int places)
111 {
112 M_APM tmp1, tmp2, a0, b0, c0, a1, b1, sum, pow_2;
113 int dplaces, nn;
114
115 tmp1 = M_get_stack_var();
116 tmp2 = M_get_stack_var();
117 a0 = M_get_stack_var();
118 b0 = M_get_stack_var();
119 c0 = M_get_stack_var();
120 a1 = M_get_stack_var();
121 b1 = M_get_stack_var();
122 sum = M_get_stack_var();
123 pow_2 = M_get_stack_var();
124
125 dplaces = places + 16;
126
127 m_apm_copy(a0, MM_One);
128 m_apm_copy(sum, MM_One);
129 m_apm_copy(pow_2, MM_Four);
130 m_apm_sqrt(b0, dplaces, MM_0_5); /* sqrt(0.5) */
131
132 while (TRUE)
133 {
134 m_apm_add(tmp1, a0, b0);
135 m_apm_multiply(a1, MM_0_5, tmp1);
136
137 m_apm_multiply(tmp1, a0, b0);
138 m_apm_sqrt(b1, dplaces, tmp1);
139
140 m_apm_subtract(tmp1, a0, b0);
141 m_apm_multiply(c0, MM_0_5, tmp1);
142
143 /*
144 * the net 'PI' calculated from this iteration will
145 * be accurate to ~4 X the value of (c0)'s exponent.
146 * this was determined experimentally.
147 */
148
149 nn = -4 * c0->m_apm_exponent;
150
151 m_apm_multiply(tmp1, c0, c0);
152 m_apm_multiply(tmp2, tmp1, pow_2);
153 m_apm_subtract(tmp1, sum, tmp2);
154 m_apm_round(sum, dplaces, tmp1);
155
156 if (nn >= dplaces)
157 break;
158
159 m_apm_copy(a0, a1);
160 m_apm_copy(b0, b1);
161
162 m_apm_multiply(tmp1, pow_2, MM_Two);
163 m_apm_copy(pow_2, tmp1);
164 }
165
166 m_apm_add(tmp1, a1, b1);
167 m_apm_multiply(tmp2, tmp1, tmp1);
168 m_apm_divide(tmp1, dplaces, tmp2, sum);
169 m_apm_round(outv, places, tmp1);
170
171 M_restore_stack(9);
172 }
173 /****************************************************************************/