"Fossies" - the Fresh Open Source Software Archive 
Member "mapm_4.9.5a/mapmsqrt.c" (21 Feb 2010, 5225 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 "mapmsqrt.c" see the
Fossies "Dox" file reference documentation.
1
2 /*
3 * M_APM - mapmsqrt.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: mapmsqrt.c,v 1.19 2007/12/03 01:57:31 mike Exp $
23 *
24 * This file contains the SQRT function.
25 *
26 * $Log: mapmsqrt.c,v $
27 * Revision 1.19 2007/12/03 01:57:31 mike
28 * Update license
29 *
30 * Revision 1.18 2003/07/21 20:39:00 mike
31 * Modify error messages to be in a consistent format.
32 *
33 * Revision 1.17 2003/05/07 16:36:04 mike
34 * simplify 'nexp' logic
35 *
36 * Revision 1.16 2003/03/31 21:50:14 mike
37 * call generic error handling function
38 *
39 * Revision 1.15 2003/03/11 21:29:00 mike
40 * round an intermediate result for faster runtime.
41 *
42 * Revision 1.14 2002/11/03 22:00:46 mike
43 * Updated function parameters to use the modern style
44 *
45 * Revision 1.13 2001/07/10 22:50:31 mike
46 * minor optimization
47 *
48 * Revision 1.12 2000/09/26 18:32:04 mike
49 * use new algorithm which only uses multiply and subtract
50 * (avoids the slower version which used division)
51 *
52 * Revision 1.11 2000/07/11 17:56:22 mike
53 * make better estimate for initial precision
54 *
55 * Revision 1.10 1999/07/21 02:48:45 mike
56 * added some comments
57 *
58 * Revision 1.9 1999/07/19 00:25:44 mike
59 * adjust local precision again
60 *
61 * Revision 1.8 1999/07/19 00:09:41 mike
62 * adjust local precision during loop
63 *
64 * Revision 1.7 1999/07/18 22:57:08 mike
65 * change to dynamically changing local precision and
66 * change tolerance checks using integers
67 *
68 * Revision 1.6 1999/06/19 21:18:00 mike
69 * changed local static variables to MAPM stack variables
70 *
71 * Revision 1.5 1999/05/31 01:40:39 mike
72 * minor update to normalizing the exponent
73 *
74 * Revision 1.4 1999/05/31 01:21:41 mike
75 * optimize for large exponents
76 *
77 * Revision 1.3 1999/05/12 20:59:35 mike
78 * use a better 'guess' function
79 *
80 * Revision 1.2 1999/05/10 21:15:26 mike
81 * added some comments
82 *
83 * Revision 1.1 1999/05/10 20:56:31 mike
84 * Initial revision
85 */
86
87 #include "m_apm_lc.h"
88
89 /****************************************************************************/
90 void m_apm_sqrt(M_APM rr, int places, M_APM aa)
91 {
92 M_APM last_x, guess, tmpN, tmp7, tmp8, tmp9;
93 int ii, bflag, nexp, tolerance, dplaces;
94
95 if (aa->m_apm_sign <= 0)
96 {
97 if (aa->m_apm_sign == -1)
98 {
99 M_apm_log_error_msg(M_APM_RETURN, "\'m_apm_sqrt\', Negative argument");
100 }
101
102 M_set_to_zero(rr);
103 return;
104 }
105
106 last_x = M_get_stack_var();
107 guess = M_get_stack_var();
108 tmpN = M_get_stack_var();
109 tmp7 = M_get_stack_var();
110 tmp8 = M_get_stack_var();
111 tmp9 = M_get_stack_var();
112
113 m_apm_copy(tmpN, aa);
114
115 /*
116 normalize the input number (make the exponent near 0) so
117 the 'guess' function will not over/under flow on large
118 magnitude exponents.
119 */
120
121 nexp = aa->m_apm_exponent / 2;
122 tmpN->m_apm_exponent -= 2 * nexp;
123
124 M_get_sqrt_guess(guess, tmpN); /* actually gets 1/sqrt guess */
125
126 tolerance = places + 4;
127 dplaces = places + 16;
128 bflag = FALSE;
129
130 m_apm_negate(last_x, MM_Ten);
131
132 /* Use the following iteration to calculate 1 / sqrt(N) :
133
134 X = 0.5 * X * [ 3 - N * X^2 ]
135 n+1
136 */
137
138 ii = 0;
139
140 while (TRUE)
141 {
142 m_apm_multiply(tmp9, tmpN, guess);
143 m_apm_multiply(tmp8, tmp9, guess);
144 m_apm_round(tmp7, dplaces, tmp8);
145 m_apm_subtract(tmp9, MM_Three, tmp7);
146 m_apm_multiply(tmp8, tmp9, guess);
147 m_apm_multiply(tmp9, tmp8, MM_0_5);
148
149 if (bflag)
150 break;
151
152 m_apm_round(guess, dplaces, tmp9);
153
154 /* force at least 2 iterations so 'last_x' has valid data */
155
156 if (ii != 0)
157 {
158 m_apm_subtract(tmp7, guess, last_x);
159
160 if (tmp7->m_apm_sign == 0)
161 break;
162
163 /*
164 * if we are within a factor of 4 on the error term,
165 * we will be accurate enough after the *next* iteration
166 * is complete. (note that the sign of the exponent on
167 * the error term will be a negative number).
168 */
169
170 if ((-4 * tmp7->m_apm_exponent) > tolerance)
171 bflag = TRUE;
172 }
173
174 m_apm_copy(last_x, guess);
175 ii++;
176 }
177
178 /*
179 * multiply by the starting number to get the final
180 * sqrt and then adjust the exponent since we found
181 * the sqrt of the normalized number.
182 */
183
184 m_apm_multiply(tmp8, tmp9, tmpN);
185 m_apm_round(rr, places, tmp8);
186 rr->m_apm_exponent += nexp;
187
188 M_restore_stack(6);
189 }
190 /****************************************************************************/