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 /****************************************************************************/ 191