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/****************************************************************************/
93void	M_get_sqrt_guess(M_APM r, M_APM a)
94{
95char	buf[48];
96double  dd;
97
98m_apm_to_string(buf, 15, a);
99dd = atof(buf);                     /* sqrt algorithm actually finds 1/sqrt */
100m_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/****************************************************************************/
112void	M_get_cbrt_guess(M_APM r, M_APM a)
113{
114char	buf[48];
115double  dd;
116
117m_apm_to_string(buf, 15, a);
118dd = atof(buf);
119dd = log(dd) / 3.0;                 /* cbrt algorithm actually finds 1/cbrt */
120m_apm_set_double(r, (1.00001 / exp(dd)));
121}
122/****************************************************************************/
123void	M_get_log_guess(M_APM r, M_APM a)
124{
125char	buf[48];
126double  dd;
127
128m_apm_to_string(buf, 15, a);
129dd = atof(buf);
130m_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 */
138void	M_get_asin_guess(M_APM r, M_APM a)
139{
140char	buf[48];
141double  dd;
142
143m_apm_to_string(buf, 15, a);
144dd = atof(buf);
145m_apm_set_double(r, (1.00001 * asin(dd)));       /* induce error of 10 ^ -5 */
146}
147/****************************************************************************/
148void	M_get_acos_guess(M_APM r, M_APM a)
149{
150char	buf[48];
151double  dd;
152
153m_apm_to_string(buf, 15, a);
154dd = atof(buf);
155m_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*/
161void	m_apm_set_double(M_APM atmp, double dd)
162{
163char	*cp, *p, *ps, buf[64];
164
165if (dd == 0.0)                     /* special case for 0 exactly */
166   M_set_to_zero(atmp);
167else
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