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/****************************************************************************/
55void	m_apm_cbrt(M_APM rr, int places, M_APM aa)
56{
57M_APM   last_x, guess, tmpN, tmp7, tmp8, tmp9;
58int	ii, nexp, bflag, tolerance, maxp, local_precision;
59
60/* result is 0 if input is 0 */
61
62if (aa->m_apm_sign == 0)
63  {
64   M_set_to_zero(rr);
65   return;
66  }
67
68last_x = M_get_stack_var();
69guess  = M_get_stack_var();
70tmpN   = M_get_stack_var();
71tmp7   = M_get_stack_var();
72tmp8   = M_get_stack_var();
73tmp9   = M_get_stack_var();
74
75/* compute the cube root of the positive number, we'll fix the sign later */
76
77m_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
85nexp = aa->m_apm_exponent / 3;
86tmpN->m_apm_exponent -= 3 * nexp;
87
88M_get_cbrt_guess(guess, tmpN);
89
90tolerance       = places + 4;
91maxp            = places + 16;
92bflag           = FALSE;
93local_precision = 14;
94
95m_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
104ii = 0;
105
106while (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
144m_apm_multiply(tmp9, guess, guess);
145m_apm_multiply(tmp8, tmp9, tmpN);
146m_apm_round(rr, places, tmp8);
147
148rr->m_apm_exponent += nexp;
149rr->m_apm_sign = aa->m_apm_sign;
150M_restore_stack(6);
151}
152/****************************************************************************/
153
154