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 /****************************************************************************/ 70 void M_log_solve_cubic(M_APM rr, int places, M_APM nn) 71 { 72 M_APM tmp0, tmp1, tmp2, tmp3, guess; 73 int ii, maxp, tolerance, local_precision; 74 75 guess = M_get_stack_var(); 76 tmp0 = M_get_stack_var(); 77 tmp1 = M_get_stack_var(); 78 tmp2 = M_get_stack_var(); 79 tmp3 = M_get_stack_var(); 80 81 M_get_log_guess(guess, nn); 82 83 tolerance = -(places + 4); 84 maxp = places + 16; 85 local_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 98 ii = 0; 99 100 while (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 127 m_apm_round(rr, places, tmp3); 128 M_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 */ 150 void M_log_basic_iteration(M_APM rr, int places, M_APM nn) 151 { 152 M_APM tmp0, tmp1, tmp2, tmpX; 153 154 if (places < 360) 155 { 156 M_log_solve_cubic(rr, places, nn); 157 } 158 else 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