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 /****************************************************************************/
m_apm_cbrt(M_APM rr,int places,M_APM aa)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
154