"Fossies" - the Fresh Open Source Software Archive 
Member "mapm_4.9.5a/mapmcbrt.c" (21 Feb 2010, 3874 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 "mapmcbrt.c" see the
Fossies "Dox" file reference documentation.
1
2 /*
3 * M_APM - mapmcbrt.c
4 *
5 * Copyright (C) 2000 - 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: mapmcbrt.c,v 1.8 2007/12/03 01:50:32 mike Exp $
23 *
24 * This file contains the CBRT (cube root) function.
25 *
26 * $Log: mapmcbrt.c,v $
27 * Revision 1.8 2007/12/03 01:50:32 mike
28 * Update license
29 *
30 * Revision 1.7 2003/05/05 18:17:38 mike
31 * simplify some logic
32 *
33 * Revision 1.6 2003/04/16 16:55:58 mike
34 * use new faster algorithm which finds 1 / cbrt(N)
35 *
36 * Revision 1.5 2002/11/03 21:34:34 mike
37 * Updated function parameters to use the modern style
38 *
39 * Revision 1.4 2000/10/30 16:42:22 mike
40 * minor speed optimization
41 *
42 * Revision 1.3 2000/07/11 18:03:39 mike
43 * make better estimate for initial precision
44 *
45 * Revision 1.2 2000/04/08 18:34:35 mike
46 * added some more comments
47 *
48 * Revision 1.1 2000/04/03 17:58:04 mike
49 * Initial revision
50 */
51
52 #include "m_apm_lc.h"
53
54 /****************************************************************************/
55 void m_apm_cbrt(M_APM rr, int places, M_APM aa)
56 {
57 M_APM last_x, guess, tmpN, tmp7, tmp8, tmp9;
58 int ii, nexp, bflag, tolerance, maxp, local_precision;
59
60 /* result is 0 if input is 0 */
61
62 if (aa->m_apm_sign == 0)
63 {
64 M_set_to_zero(rr);
65 return;
66 }
67
68 last_x = M_get_stack_var();
69 guess = M_get_stack_var();
70 tmpN = M_get_stack_var();
71 tmp7 = M_get_stack_var();
72 tmp8 = M_get_stack_var();
73 tmp9 = M_get_stack_var();
74
75 /* compute the cube root of the positive number, we'll fix the sign later */
76
77 m_apm_absolute_value(tmpN, aa);
78
79 /*
80 normalize the input number (make the exponent near 0) so
81 the 'guess' function will not over/under flow on large
82 magnitude exponents.
83 */
84
85 nexp = aa->m_apm_exponent / 3;
86 tmpN->m_apm_exponent -= 3 * nexp;
87
88 M_get_cbrt_guess(guess, tmpN);
89
90 tolerance = places + 4;
91 maxp = places + 16;
92 bflag = FALSE;
93 local_precision = 14;
94
95 m_apm_negate(last_x, MM_Ten);
96
97 /* Use the following iteration to calculate 1 / cbrt(N) :
98
99 4
100 X = [ 4 * X - N * X ] / 3
101 n+1
102 */
103
104 ii = 0;
105
106 while (TRUE)
107 {
108 m_apm_multiply(tmp8, guess, guess);
109 m_apm_multiply(tmp7, tmp8, tmp8);
110 m_apm_round(tmp8, local_precision, tmp7);
111 m_apm_multiply(tmp9, tmpN, tmp8);
112
113 m_apm_multiply(tmp8, MM_Four, guess);
114 m_apm_subtract(tmp7, tmp8, tmp9);
115 m_apm_divide(guess, local_precision, tmp7, MM_Three);
116
117 if (bflag)
118 break;
119
120 /* force at least 2 iterations so 'last_x' has valid data */
121
122 if (ii != 0)
123 {
124 m_apm_subtract(tmp8, guess, last_x);
125
126 if (tmp8->m_apm_sign == 0)
127 break;
128
129 if ((-4 * tmp8->m_apm_exponent) > tolerance)
130 bflag = TRUE;
131 }
132
133 local_precision *= 2;
134
135 if (local_precision > maxp)
136 local_precision = maxp;
137
138 m_apm_copy(last_x, guess);
139 ii = 1;
140 }
141
142 /* final cbrt = N * guess ^ 2 */
143
144 m_apm_multiply(tmp9, guess, guess);
145 m_apm_multiply(tmp8, tmp9, tmpN);
146 m_apm_round(rr, places, tmp8);
147
148 rr->m_apm_exponent += nexp;
149 rr->m_apm_sign = aa->m_apm_sign;
150 M_restore_stack(6);
151 }
152 /****************************************************************************/
153