"Fossies" - the Fresh Open Source Software Archive 
Member "mapm_4.9.5a/mapm_lg2.c" (21 Feb 2010, 4577 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_lg2.c" see the
Fossies "Dox" file reference documentation.
1
2 /*
3 * M_APM - mapm_lg2.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_lg2.c,v 1.9 2007/12/03 01:42:06 mike Exp $
23 *
24 * This file contains the iterative function to compute the LOG
25 * This is an internal function to the library and is not intended
26 * to be called directly by the user.
27 *
28 * $Log: mapm_lg2.c,v $
29 * Revision 1.9 2007/12/03 01:42:06 mike
30 * Update license
31 *
32 * Revision 1.8 2003/05/12 17:52:15 mike
33 * rearrange some logic
34 *
35 * Revision 1.7 2003/05/03 11:24:51 mike
36 * optimize decimal places
37 *
38 * Revision 1.6 2003/05/01 21:58:27 mike
39 * remove math.h
40 *
41 * Revision 1.5 2003/05/01 20:53:38 mike
42 * implement a new algorithm for log
43 *
44 * Revision 1.4 2003/04/09 20:21:29 mike
45 * fix rare corner condition by intentionally inducing a
46 * 10 ^ -5 error in the initial guess.
47 *
48 * Revision 1.3 2003/03/31 22:13:15 mike
49 * call generic error handling function
50 *
51 * Revision 1.2 2003/03/30 21:27:22 mike
52 * add comments
53 *
54 * Revision 1.1 2003/03/30 21:18:07 mike
55 * Initial revision
56 */
57
58 #include "m_apm_lc.h"
59
60 /****************************************************************************/
61
62 /*
63 * compute rr = log(nn)
64 *
65 * input is assumed to not exceed the exponent range of a normal
66 * 'C' double ( |exponent| must be < 308)
67 */
68
69 /****************************************************************************/
70 void M_log_solve_cubic(M_APM rr, int places, M_APM nn)
71 {
72 M_APM tmp0, tmp1, tmp2, tmp3, guess;
73 int ii, maxp, tolerance, local_precision;
74
75 guess = M_get_stack_var();
76 tmp0 = M_get_stack_var();
77 tmp1 = M_get_stack_var();
78 tmp2 = M_get_stack_var();
79 tmp3 = M_get_stack_var();
80
81 M_get_log_guess(guess, nn);
82
83 tolerance = -(places + 4);
84 maxp = places + 16;
85 local_precision = 18;
86
87 /* Use the following iteration to solve for log :
88
89 exp(X) - N
90 X = X - 2 * ------------
91 n+1 exp(X) + N
92
93
94 this is a cubically convergent algorithm
95 (each iteration yields 3X more digits)
96 */
97
98 ii = 0;
99
100 while (TRUE)
101 {
102 m_apm_exp(tmp1, local_precision, guess);
103
104 m_apm_subtract(tmp3, tmp1, nn);
105 m_apm_add(tmp2, tmp1, nn);
106
107 m_apm_divide(tmp1, local_precision, tmp3, tmp2);
108 m_apm_multiply(tmp0, MM_Two, tmp1);
109 m_apm_subtract(tmp3, guess, tmp0);
110
111 if (ii != 0)
112 {
113 if (((3 * tmp0->m_apm_exponent) < tolerance) || (tmp0->m_apm_sign == 0))
114 break;
115 }
116
117 m_apm_round(guess, local_precision, tmp3);
118
119 local_precision *= 3;
120
121 if (local_precision > maxp)
122 local_precision = maxp;
123
124 ii = 1;
125 }
126
127 m_apm_round(rr, places, tmp3);
128 M_restore_stack(5);
129 }
130 /****************************************************************************/
131 /*
132 * find log(N)
133 *
134 * if places < 360
135 * solve with cubically convergent algorithm above
136 *
137 * else
138 *
139 * let 'X' be 'close' to the solution (we use ~110 decimal places)
140 *
141 * let Y = N * exp(-X) - 1
142 *
143 * then
144 *
145 * log(N) = X + log(1 + Y)
146 *
147 * since 'Y' will be small, we can use the efficient log_near_1 algorithm.
148 *
149 */
150 void M_log_basic_iteration(M_APM rr, int places, M_APM nn)
151 {
152 M_APM tmp0, tmp1, tmp2, tmpX;
153
154 if (places < 360)
155 {
156 M_log_solve_cubic(rr, places, nn);
157 }
158 else
159 {
160 tmp0 = M_get_stack_var();
161 tmp1 = M_get_stack_var();
162 tmp2 = M_get_stack_var();
163 tmpX = M_get_stack_var();
164
165 M_log_solve_cubic(tmpX, 110, nn);
166
167 m_apm_negate(tmp0, tmpX);
168 m_apm_exp(tmp1, (places + 8), tmp0);
169 m_apm_multiply(tmp2, tmp1, nn);
170 m_apm_subtract(tmp1, tmp2, MM_One);
171
172 M_log_near_1(tmp0, (places - 104), tmp1);
173
174 m_apm_add(tmp1, tmpX, tmp0);
175 m_apm_round(rr, places, tmp1);
176
177 M_restore_stack(4);
178 }
179 }
180 /****************************************************************************/