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 /****************************************************************************/
M_get_sqrt_guess(M_APM r,M_APM a)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 /****************************************************************************/
M_get_cbrt_guess(M_APM r,M_APM a)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 /****************************************************************************/
M_get_log_guess(M_APM r,M_APM a)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 */
M_get_asin_guess(M_APM r,M_APM a)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 /****************************************************************************/
M_get_acos_guess(M_APM r,M_APM a)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 */
m_apm_set_double(M_APM atmp,double dd)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