"Fossies" - the Fresh Open Source Software Archive 
Member "mapm_4.9.5a/mapm_lg3.c" (21 Feb 2010, 6190 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_lg3.c" see the
Fossies "Dox" file reference documentation.
1
2 /*
3 * M_APM - mapm_lg3.c
4 *
5 * Copyright (C) 2003 - 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_lg3.c,v 1.7 2007/12/03 01:42:59 mike Exp $
23 *
24 * This file contains the function to compute log(2), log(10),
25 * and 1/log(10) to the desired precision using an AGM algorithm.
26 *
27 * $Log: mapm_lg3.c,v $
28 * Revision 1.7 2007/12/03 01:42:59 mike
29 * Update license
30 *
31 * Revision 1.6 2003/12/09 01:25:06 mike
32 * actually compute the first term of the AGM iteration instead
33 * of assuming the inputs a=1 and b=10^-N.
34 *
35 * Revision 1.5 2003/12/04 03:19:16 mike
36 * rearrange logic in AGM to be more straight-forward
37 *
38 * Revision 1.4 2003/05/01 22:04:37 mike
39 * rearrange some code
40 *
41 * Revision 1.3 2003/05/01 21:58:31 mike
42 * remove math.h
43 *
44 * Revision 1.2 2003/03/30 22:14:58 mike
45 * add comments
46 *
47 * Revision 1.1 2003/03/30 21:18:04 mike
48 * Initial revision
49 */
50
51 #include "m_apm_lc.h"
52
53 /*
54 * using the 'R' function (defined below) for 'N' decimal places :
55 *
56 *
57 * -N -N
58 * log(2) = R(1, 0.5 * 10 ) - R(1, 10 )
59 *
60 *
61 * -N -N
62 * log(10) = R(1, 0.1 * 10 ) - R(1, 10 )
63 *
64 *
65 * In general:
66 *
67 * -N -N
68 * log(x) = R(1, 10 / x) - R(1, 10 )
69 *
70 *
71 * I found this on a web site which went into considerable detail
72 * on the history of log(2). This formula is algebraically identical
73 * to the formula specified in J. Borwein and P. Borwein's book
74 * "PI and the AGM". (reference algorithm 7.2)
75 */
76
77 /****************************************************************************/
78 /*
79 * check if our local copy of log(2) & log(10) is precise
80 * enough for our purpose. if not, calculate them so it's
81 * as precise as desired, accurate to at least 'places'.
82 */
83 void M_check_log_places(int places)
84 {
85 M_APM tmp6, tmp7, tmp8, tmp9;
86 int dplaces;
87
88 dplaces = places + 4;
89
90 if (dplaces > MM_lc_log_digits)
91 {
92 MM_lc_log_digits = dplaces + 4;
93
94 tmp6 = M_get_stack_var();
95 tmp7 = M_get_stack_var();
96 tmp8 = M_get_stack_var();
97 tmp9 = M_get_stack_var();
98
99 dplaces += 6 + (int)log10((double)places);
100
101 m_apm_copy(tmp7, MM_One);
102 tmp7->m_apm_exponent = -places;
103
104 M_log_AGM_R_func(tmp8, dplaces, MM_One, tmp7);
105
106 m_apm_multiply(tmp6, tmp7, MM_0_5);
107
108 M_log_AGM_R_func(tmp9, dplaces, MM_One, tmp6);
109
110 m_apm_subtract(MM_lc_log2, tmp9, tmp8); /* log(2) */
111
112 tmp7->m_apm_exponent -= 1; /* divide by 10 */
113
114 M_log_AGM_R_func(tmp9, dplaces, MM_One, tmp7);
115
116 m_apm_subtract(MM_lc_log10, tmp9, tmp8); /* log(10) */
117 m_apm_reciprocal(MM_lc_log10R, dplaces, MM_lc_log10); /* 1 / log(10) */
118
119 M_restore_stack(4);
120 }
121 }
122 /****************************************************************************/
123
124 /*
125 * define a notation for a function 'R' :
126 *
127 *
128 *
129 * 1
130 * R (a0, b0) = ------------------------------
131 *
132 * ----
133 * \
134 * \ n-1 2 2
135 * 1 - | 2 * (a - b )
136 * / n n
137 * /
138 * ----
139 * n >= 0
140 *
141 *
142 * where a, b are the classic AGM iteration :
143 *
144 *
145 * a = 0.5 * (a + b )
146 * n+1 n n
147 *
148 *
149 * b = sqrt(a * b )
150 * n+1 n n
151 *
152 *
153 *
154 * define a variable 'c' for more efficient computation :
155 *
156 * 2 2 2
157 * c = 0.5 * (a - b ) , c = a - b
158 * n+1 n n n n n
159 *
160 */
161
162 /****************************************************************************/
163 void M_log_AGM_R_func(M_APM rr, int places, M_APM aa, M_APM bb)
164 {
165 M_APM tmp1, tmp2, tmp3, tmp4, tmpC2, sum, pow_2, tmpA0, tmpB0;
166 int tolerance, dplaces;
167
168 tmpA0 = M_get_stack_var();
169 tmpB0 = M_get_stack_var();
170 tmpC2 = M_get_stack_var();
171 tmp1 = M_get_stack_var();
172 tmp2 = M_get_stack_var();
173 tmp3 = M_get_stack_var();
174 tmp4 = M_get_stack_var();
175 sum = M_get_stack_var();
176 pow_2 = M_get_stack_var();
177
178 tolerance = places + 8;
179 dplaces = places + 16;
180
181 m_apm_copy(tmpA0, aa);
182 m_apm_copy(tmpB0, bb);
183 m_apm_copy(pow_2, MM_0_5);
184
185 m_apm_multiply(tmp1, aa, aa); /* 0.5 * [ a ^ 2 - b ^ 2 ] */
186 m_apm_multiply(tmp2, bb, bb);
187 m_apm_subtract(tmp3, tmp1, tmp2);
188 m_apm_multiply(sum, MM_0_5, tmp3);
189
190 while (TRUE)
191 {
192 m_apm_subtract(tmp1, tmpA0, tmpB0); /* C n+1 = 0.5 * [ An - Bn ] */
193 m_apm_multiply(tmp4, MM_0_5, tmp1); /* C n+1 */
194 m_apm_multiply(tmpC2, tmp4, tmp4); /* C n+1 ^ 2 */
195
196 /* do the AGM */
197
198 m_apm_add(tmp1, tmpA0, tmpB0);
199 m_apm_multiply(tmp3, MM_0_5, tmp1);
200
201 m_apm_multiply(tmp2, tmpA0, tmpB0);
202 m_apm_sqrt(tmpB0, dplaces, tmp2);
203
204 m_apm_round(tmpA0, dplaces, tmp3);
205
206 /* end AGM */
207
208 m_apm_multiply(tmp2, MM_Two, pow_2);
209 m_apm_copy(pow_2, tmp2);
210
211 m_apm_multiply(tmp1, tmpC2, pow_2);
212 m_apm_add(tmp3, sum, tmp1);
213
214 if ((tmp1->m_apm_sign == 0) ||
215 ((-2 * tmp1->m_apm_exponent) > tolerance))
216 break;
217
218 m_apm_round(sum, dplaces, tmp3);
219 }
220
221 m_apm_subtract(tmp4, MM_One, tmp3);
222 m_apm_reciprocal(rr, places, tmp4);
223
224 M_restore_stack(9);
225 }
226 /****************************************************************************/