1
2/*
3 *  M_APM  -  mapm_lg2.c
4 *
5 *  Copyright (C) 2003 - 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: mapm_lg2.c,v 1.9 2007/12/03 01:42:06 mike Exp $
23 *
24 *      This file contains the iterative function to compute the LOG
25 *	This is an internal function to the library and is not intended
26 *	to be called directly by the user.
27 *
28 *      $Log: mapm_lg2.c,v $
29 *      Revision 1.9  2007/12/03 01:42:06  mike
30 *      Update license
31 *
32 *      Revision 1.8  2003/05/12 17:52:15  mike
33 *      rearrange some logic
34 *
35 *      Revision 1.7  2003/05/03 11:24:51  mike
36 *      optimize decimal places
37 *
38 *      Revision 1.6  2003/05/01 21:58:27  mike
39 *      remove math.h
40 *
41 *      Revision 1.5  2003/05/01 20:53:38  mike
42 *      implement a new algorithm for log
43 *
44 *      Revision 1.4  2003/04/09 20:21:29  mike
45 *      fix rare corner condition by intentionally inducing a
46 *      10 ^ -5 error in the initial guess.
47 *
48 *      Revision 1.3  2003/03/31 22:13:15  mike
49 *      call generic error handling function
50 *
51 *      Revision 1.2  2003/03/30 21:27:22  mike
52 *      add comments
53 *
54 *      Revision 1.1  2003/03/30 21:18:07  mike
55 *      Initial revision
56 */
57
58#include "m_apm_lc.h"
59
60/****************************************************************************/
61
62/*
63 *      compute rr = log(nn)
64 *
65 *	input is assumed to not exceed the exponent range of a normal
66 *	'C' double ( |exponent| must be < 308)
67 */
68
69/****************************************************************************/
70void	M_log_solve_cubic(M_APM rr, int places, M_APM nn)
71{
72M_APM   tmp0, tmp1, tmp2, tmp3, guess;
73int	ii, maxp, tolerance, local_precision;
74
75guess = M_get_stack_var();
76tmp0  = M_get_stack_var();
77tmp1  = M_get_stack_var();
78tmp2  = M_get_stack_var();
79tmp3  = M_get_stack_var();
80
81M_get_log_guess(guess, nn);
82
83tolerance       = -(places + 4);
84maxp            = places + 16;
85local_precision = 18;
86
87/*    Use the following iteration to solve for log :
88
89                        exp(X) - N
90      X     =  X - 2 * ------------
91       n+1              exp(X) + N
92
93
94      this is a cubically convergent algorithm
95      (each iteration yields 3X more digits)
96*/
97
98ii = 0;
99
100while (TRUE)
101  {
102   m_apm_exp(tmp1, local_precision, guess);
103
104   m_apm_subtract(tmp3, tmp1, nn);
105   m_apm_add(tmp2, tmp1, nn);
106
107   m_apm_divide(tmp1, local_precision, tmp3, tmp2);
108   m_apm_multiply(tmp0, MM_Two, tmp1);
109   m_apm_subtract(tmp3, guess, tmp0);
110
111   if (ii != 0)
112     {
113      if (((3 * tmp0->m_apm_exponent) < tolerance) || (tmp0->m_apm_sign == 0))
114        break;
115     }
116
117   m_apm_round(guess, local_precision, tmp3);
118
119   local_precision *= 3;
120
121   if (local_precision > maxp)
122     local_precision = maxp;
123
124   ii = 1;
125  }
126
127m_apm_round(rr, places, tmp3);
128M_restore_stack(5);
129}
130/****************************************************************************/
131/*
132 *      find log(N)
133 *
134 *      if places < 360
135 *         solve with cubically convergent algorithm above
136 *
137 *      else
138 *
139 *      let 'X' be 'close' to the solution   (we use ~110 decimal places)
140 *
141 *      let Y = N * exp(-X) - 1
142 *
143 *	then
144 *
145 *	log(N) = X + log(1 + Y)
146 *
147 *      since 'Y' will be small, we can use the efficient log_near_1 algorithm.
148 *
149 */
150void	M_log_basic_iteration(M_APM rr, int places, M_APM nn)
151{
152M_APM   tmp0, tmp1, tmp2, tmpX;
153
154if (places < 360)
155  {
156   M_log_solve_cubic(rr, places, nn);
157  }
158else
159  {
160   tmp0 = M_get_stack_var();
161   tmp1 = M_get_stack_var();
162   tmp2 = M_get_stack_var();
163   tmpX = M_get_stack_var();
164
165   M_log_solve_cubic(tmpX, 110, nn);
166
167   m_apm_negate(tmp0, tmpX);
168   m_apm_exp(tmp1, (places + 8), tmp0);
169   m_apm_multiply(tmp2, tmp1, nn);
170   m_apm_subtract(tmp1, tmp2, MM_One);
171
172   M_log_near_1(tmp0, (places - 104), tmp1);
173
174   m_apm_add(tmp1, tmpX, tmp0);
175   m_apm_round(rr, places, tmp1);
176
177   M_restore_stack(4);
178  }
179}
180/****************************************************************************/
181