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 */ 83 void M_check_log_places(int places) 84 { 85 M_APM tmp6, tmp7, tmp8, tmp9; 86 int dplaces; 87 88 dplaces = places + 4; 89 90 if (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 /****************************************************************************/ 163 void M_log_AGM_R_func(M_APM rr, int places, M_APM aa, M_APM bb) 164 { 165 M_APM tmp1, tmp2, tmp3, tmp4, tmpC2, sum, pow_2, tmpA0, tmpB0; 166 int tolerance, dplaces; 167 168 tmpA0 = M_get_stack_var(); 169 tmpB0 = M_get_stack_var(); 170 tmpC2 = M_get_stack_var(); 171 tmp1 = M_get_stack_var(); 172 tmp2 = M_get_stack_var(); 173 tmp3 = M_get_stack_var(); 174 tmp4 = M_get_stack_var(); 175 sum = M_get_stack_var(); 176 pow_2 = M_get_stack_var(); 177 178 tolerance = places + 8; 179 dplaces = places + 16; 180 181 m_apm_copy(tmpA0, aa); 182 m_apm_copy(tmpB0, bb); 183 m_apm_copy(pow_2, MM_0_5); 184 185 m_apm_multiply(tmp1, aa, aa); /* 0.5 * [ a ^ 2 - b ^ 2 ] */ 186 m_apm_multiply(tmp2, bb, bb); 187 m_apm_subtract(tmp3, tmp1, tmp2); 188 m_apm_multiply(sum, MM_0_5, tmp3); 189 190 while (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 221 m_apm_subtract(tmp4, MM_One, tmp3); 222 m_apm_reciprocal(rr, places, tmp4); 223 224 M_restore_stack(9); 225 } 226 /****************************************************************************/ 227