1
2/*
3 *  M_APM  -  mapm_lg3.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_lg3.c,v 1.7 2007/12/03 01:42:59 mike Exp $
23 *
24 *      This file contains the function to compute log(2), log(10),
25 *	and 1/log(10) to the desired precision using an AGM algorithm.
26 *
27 *      $Log: mapm_lg3.c,v $
28 *      Revision 1.7  2007/12/03 01:42:59  mike
29 *      Update license
30 *
31 *      Revision 1.6  2003/12/09 01:25:06  mike
32 *      actually compute the first term of the AGM iteration instead
33 *      of assuming the inputs a=1 and b=10^-N.
34 *
35 *      Revision 1.5  2003/12/04 03:19:16  mike
36 *      rearrange logic in AGM to be more straight-forward
37 *
38 *      Revision 1.4  2003/05/01 22:04:37  mike
39 *      rearrange some code
40 *
41 *      Revision 1.3  2003/05/01 21:58:31  mike
42 *      remove math.h
43 *
44 *      Revision 1.2  2003/03/30 22:14:58  mike
45 *      add comments
46 *
47 *      Revision 1.1  2003/03/30 21:18:04  mike
48 *      Initial revision
49 */
50
51#include "m_apm_lc.h"
52
53/*
54 *  using the 'R' function (defined below) for 'N' decimal places :
55 *
56 *
57 *                          -N             -N
58 *  log(2)  =  R(1, 0.5 * 10  )  -  R(1, 10  )
59 *
60 *
61 *                          -N             -N
62 *  log(10) =  R(1, 0.1 * 10  )  -  R(1, 10  )
63 *
64 *
65 *  In general:
66 *
67 *                    -N                -N
68 *  log(x)  =  R(1, 10  / x)  -  R(1, 10  )
69 *
70 *
71 *  I found this on a web site which went into considerable detail
72 *  on the history of log(2). This formula is algebraically identical
73 *  to the formula specified in J. Borwein and P. Borwein's book
74 *  "PI and the AGM". (reference algorithm 7.2)
75 */
76
77/****************************************************************************/
78/*
79 *	check if our local copy of log(2) & log(10) is precise
80 *      enough for our purpose. if not, calculate them so it's
81 *	as precise as desired, accurate to at least 'places'.
82 */
83void	M_check_log_places(int places)
84{
85M_APM   tmp6, tmp7, tmp8, tmp9;
86int     dplaces;
87
88dplaces = places + 4;
89
90if (dplaces > MM_lc_log_digits)
91  {
92   MM_lc_log_digits = dplaces + 4;
93
94   tmp6 = M_get_stack_var();
95   tmp7 = M_get_stack_var();
96   tmp8 = M_get_stack_var();
97   tmp9 = M_get_stack_var();
98
99   dplaces += 6 + (int)log10((double)places);
100
101   m_apm_copy(tmp7, MM_One);
102   tmp7->m_apm_exponent = -places;
103
104   M_log_AGM_R_func(tmp8, dplaces, MM_One, tmp7);
105
106   m_apm_multiply(tmp6, tmp7, MM_0_5);
107
108   M_log_AGM_R_func(tmp9, dplaces, MM_One, tmp6);
109
110   m_apm_subtract(MM_lc_log2, tmp9, tmp8);               /* log(2) */
111
112   tmp7->m_apm_exponent -= 1;                            /* divide by 10 */
113
114   M_log_AGM_R_func(tmp9, dplaces, MM_One, tmp7);
115
116   m_apm_subtract(MM_lc_log10, tmp9, tmp8);              /* log(10) */
117   m_apm_reciprocal(MM_lc_log10R, dplaces, MM_lc_log10); /* 1 / log(10) */
118
119   M_restore_stack(4);
120  }
121}
122/****************************************************************************/
123
124/*
125 *	define a notation for a function 'R' :
126 *
127 *
128 *
129 *                                    1
130 *      R (a0, b0)  =  ------------------------------
131 *
132 *                          ----
133 *                           \
134 *                            \     n-1      2    2
135 *                      1  -   |   2    *  (a  - b )
136 *                            /              n    n
137 *                           /
138 *                          ----
139 *                         n >= 0
140 *
141 *
142 *      where a, b are the classic AGM iteration :
143 *
144 *
145 *      a    =  0.5 * (a  + b )
146 *       n+1            n    n
147 *
148 *
149 *      b    =  sqrt(a  * b )
150 *       n+1          n    n
151 *
152 *
153 *
154 *      define a variable 'c' for more efficient computation :
155 *
156 *                                      2     2     2
157 *      c    =  0.5 * (a  - b )    ,   c  =  a  -  b
158 *       n+1            n    n          n     n     n
159 *
160 */
161
162/****************************************************************************/
163void	M_log_AGM_R_func(M_APM rr, int places, M_APM aa, M_APM bb)
164{
165M_APM   tmp1, tmp2, tmp3, tmp4, tmpC2, sum, pow_2, tmpA0, tmpB0;
166int	tolerance, dplaces;
167
168tmpA0 = M_get_stack_var();
169tmpB0 = M_get_stack_var();
170tmpC2 = M_get_stack_var();
171tmp1  = M_get_stack_var();
172tmp2  = M_get_stack_var();
173tmp3  = M_get_stack_var();
174tmp4  = M_get_stack_var();
175sum   = M_get_stack_var();
176pow_2 = M_get_stack_var();
177
178tolerance = places + 8;
179dplaces   = places + 16;
180
181m_apm_copy(tmpA0, aa);
182m_apm_copy(tmpB0, bb);
183m_apm_copy(pow_2, MM_0_5);
184
185m_apm_multiply(tmp1, aa, aa);		    /* 0.5 * [ a ^ 2 - b ^ 2 ] */
186m_apm_multiply(tmp2, bb, bb);
187m_apm_subtract(tmp3, tmp1, tmp2);
188m_apm_multiply(sum, MM_0_5, tmp3);
189
190while (TRUE)
191  {
192   m_apm_subtract(tmp1, tmpA0, tmpB0);      /* C n+1 = 0.5 * [ An - Bn ] */
193   m_apm_multiply(tmp4, MM_0_5, tmp1);      /* C n+1 */
194   m_apm_multiply(tmpC2, tmp4, tmp4);       /* C n+1 ^ 2 */
195
196   /* do the AGM */
197
198   m_apm_add(tmp1, tmpA0, tmpB0);
199   m_apm_multiply(tmp3, MM_0_5, tmp1);
200
201   m_apm_multiply(tmp2, tmpA0, tmpB0);
202   m_apm_sqrt(tmpB0, dplaces, tmp2);
203
204   m_apm_round(tmpA0, dplaces, tmp3);
205
206   /* end AGM */
207
208   m_apm_multiply(tmp2, MM_Two, pow_2);
209   m_apm_copy(pow_2, tmp2);
210
211   m_apm_multiply(tmp1, tmpC2, pow_2);
212   m_apm_add(tmp3, sum, tmp1);
213
214   if ((tmp1->m_apm_sign == 0) ||
215      ((-2 * tmp1->m_apm_exponent) > tolerance))
216     break;
217
218   m_apm_round(sum, dplaces, tmp3);
219  }
220
221m_apm_subtract(tmp4, MM_One, tmp3);
222m_apm_reciprocal(rr, places, tmp4);
223
224M_restore_stack(9);
225}
226/****************************************************************************/
227