1 2 /* 3 * M_APM - mapmgues.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: mapmgues.c,v 1.20 2007/12/03 01:52:55 mike Exp $ 23 * 24 * This file contains the functions that generate the initial 25 * 'guesses' for the sqrt, cbrt, log, arcsin, and arccos functions. 26 * 27 * $Log: mapmgues.c,v $ 28 * Revision 1.20 2007/12/03 01:52:55 mike 29 * Update license 30 * 31 * Revision 1.19 2003/07/21 20:03:49 mike 32 * check for invalid inputs to _set_double 33 * 34 * Revision 1.18 2003/05/01 21:58:45 mike 35 * remove math.h 36 * 37 * Revision 1.17 2003/04/16 16:52:47 mike 38 * change cbrt guess to use reciprocal value for new cbrt algorithm 39 * 40 * Revision 1.16 2003/04/11 14:18:13 mike 41 * add comments 42 * 43 * Revision 1.15 2003/04/10 22:28:35 mike 44 * . 45 * 46 * Revision 1.14 2003/04/09 21:33:19 mike 47 * induce same error for asin and acos 48 * 49 * Revision 1.13 2003/04/09 20:11:38 mike 50 * induce error of 10 ^ -5 in log guess for known starting 51 * point in the iterative algorithm 52 * 53 * Revision 1.12 2003/03/27 19:32:59 mike 54 * simplify log guess since caller guarantee's limited input magnitude 55 * 56 * Revision 1.11 2002/11/03 21:45:53 mike 57 * Updated function parameters to use the modern style 58 * 59 * Revision 1.10 2001/03/20 22:08:27 mike 60 * delete unneeded logic in asin guess 61 * 62 * Revision 1.9 2000/09/26 17:05:11 mike 63 * guess 1/sqrt instead of sqrt due to new sqrt algorithm 64 * 65 * Revision 1.8 2000/04/10 21:13:13 mike 66 * minor tweaks to log_guess 67 * 68 * Revision 1.7 2000/04/03 17:25:45 mike 69 * added function to estimate the cube root (cbrt) 70 * 71 * Revision 1.6 1999/07/18 01:57:35 mike 72 * adjust arc-sin guess for small exponents 73 * 74 * Revision 1.5 1999/07/09 22:32:50 mike 75 * optimize some functions 76 * 77 * Revision 1.4 1999/05/12 21:22:27 mike 78 * add more comments 79 * 80 * Revision 1.3 1999/05/12 21:00:51 mike 81 * added new sqrt guess function 82 * 83 * Revision 1.2 1999/05/11 02:10:12 mike 84 * added some comments 85 * 86 * Revision 1.1 1999/05/10 20:56:31 mike 87 * Initial revision 88 */ 89 90 #include "m_apm_lc.h" 91 92 /****************************************************************************/ 93 void M_get_sqrt_guess(M_APM r, M_APM a) 94 { 95 char buf[48]; 96 double dd; 97 98 m_apm_to_string(buf, 15, a); 99 dd = atof(buf); /* sqrt algorithm actually finds 1/sqrt */ 100 m_apm_set_double(r, (1.0 / sqrt(dd))); 101 } 102 /****************************************************************************/ 103 /* 104 * for cbrt, log, asin, and acos we induce an error of 10 ^ -5. 105 * this enables the iterative routine to be more efficient 106 * by knowing exactly how accurate the initial guess is. 107 * 108 * this also prevents some corner conditions where the iterative 109 * functions may terminate too soon. 110 */ 111 /****************************************************************************/ 112 void M_get_cbrt_guess(M_APM r, M_APM a) 113 { 114 char buf[48]; 115 double dd; 116 117 m_apm_to_string(buf, 15, a); 118 dd = atof(buf); 119 dd = log(dd) / 3.0; /* cbrt algorithm actually finds 1/cbrt */ 120 m_apm_set_double(r, (1.00001 / exp(dd))); 121 } 122 /****************************************************************************/ 123 void M_get_log_guess(M_APM r, M_APM a) 124 { 125 char buf[48]; 126 double dd; 127 128 m_apm_to_string(buf, 15, a); 129 dd = atof(buf); 130 m_apm_set_double(r, (1.00001 * log(dd))); /* induce error of 10 ^ -5 */ 131 } 132 /****************************************************************************/ 133 /* 134 * the implementation of the asin & acos functions 135 * guarantee that 'a' is always < 0.85, so it is 136 * safe to multiply by a number > 1 137 */ 138 void M_get_asin_guess(M_APM r, M_APM a) 139 { 140 char buf[48]; 141 double dd; 142 143 m_apm_to_string(buf, 15, a); 144 dd = atof(buf); 145 m_apm_set_double(r, (1.00001 * asin(dd))); /* induce error of 10 ^ -5 */ 146 } 147 /****************************************************************************/ 148 void M_get_acos_guess(M_APM r, M_APM a) 149 { 150 char buf[48]; 151 double dd; 152 153 m_apm_to_string(buf, 15, a); 154 dd = atof(buf); 155 m_apm_set_double(r, (1.00001 * acos(dd))); /* induce error of 10 ^ -5 */ 156 } 157 /****************************************************************************/ 158 /* 159 convert a C 'double' into an M_APM value. 160 */ 161 void m_apm_set_double(M_APM atmp, double dd) 162 { 163 char *cp, *p, *ps, buf[64]; 164 165 if (dd == 0.0) /* special case for 0 exactly */ 166 M_set_to_zero(atmp); 167 else 168 { 169 sprintf(buf, "%.14E", dd); 170 171 if ((cp = strstr(buf, "E")) == NULL) 172 { 173 M_apm_log_error_msg(M_APM_RETURN, 174 "\'m_apm_set_double\', Invalid double input (likely a NAN or +/- INF)"); 175 176 M_set_to_zero(atmp); 177 return; 178 } 179 180 if (atoi(cp + sizeof(char)) == 0) 181 *cp = '\0'; 182 183 p = cp; 184 185 while (TRUE) 186 { 187 p--; 188 if (*p == '0' || *p == '.') 189 *p = ' '; 190 else 191 break; 192 } 193 194 ps = buf; 195 p = buf; 196 197 while (TRUE) 198 { 199 if ((*p = *ps) == '\0') 200 break; 201 202 if (*ps++ != ' ') 203 p++; 204 } 205 206 m_apm_set_string(atmp, buf); 207 } 208 } 209 /****************************************************************************/ 210