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/****************************************************************************/
90void	m_apm_sqrt(M_APM rr, int places, M_APM aa)
91{
92M_APM   last_x, guess, tmpN, tmp7, tmp8, tmp9;
93int	ii, bflag, nexp, tolerance, dplaces;
94
95if (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
106last_x = M_get_stack_var();
107guess  = M_get_stack_var();
108tmpN   = M_get_stack_var();
109tmp7   = M_get_stack_var();
110tmp8   = M_get_stack_var();
111tmp9   = M_get_stack_var();
112
113m_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
121nexp = aa->m_apm_exponent / 2;
122tmpN->m_apm_exponent -= 2 * nexp;
123
124M_get_sqrt_guess(guess, tmpN);    /* actually gets 1/sqrt guess */
125
126tolerance = places + 4;
127dplaces   = places + 16;
128bflag     = FALSE;
129
130m_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
138ii = 0;
139
140while (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
184m_apm_multiply(tmp8, tmp9, tmpN);
185m_apm_round(rr, places, tmp8);
186rr->m_apm_exponent += nexp;
187
188M_restore_stack(6);
189}
190/****************************************************************************/
191