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 154