File indexing completed on 2024-05-12 17:22:35
0001 /* floatlog.c: logarithm and friends, based on floatnum. */ 0002 /* 0003 Copyright (C) 2007, 2008 Wolf Lammen. 0004 0005 This program is free software; you can redistribute it and/or modify 0006 it under the terms of the GNU General Public License as published by 0007 the Free Software Foundation; either version 2 of the License , or 0008 (at your option) any later version. 0009 0010 This program is distributed in the hope that it will be useful, 0011 but WITHOUT ANY WARRANTY; without even the implied warranty of 0012 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0013 GNU General Public License for more details. 0014 0015 You should have received a copy of the GNU General Public License 0016 along with this program; see the file COPYING. If not, write to: 0017 0018 The Free Software Foundation, Inc. 0019 59 Temple Place, Suite 330 0020 Boston, MA 02111-1307 USA. 0021 0022 0023 You may contact the author by: 0024 e-mail: ookami1 <at> gmx <dot> de 0025 mail: Wolf Lammen 0026 Oertzweg 45 0027 22307 Hamburg 0028 Germany 0029 *************************************************************************/ 0030 0031 #include "floatlog.h" 0032 #include "floatconst.h" 0033 #include "floatcommon.h" 0034 #include "floatseries.h" 0035 0036 typedef struct 0037 { 0038 char c2; 0039 char c3; 0040 char c5; 0041 char c7; 0042 char c10; 0043 } _lincomb; 0044 0045 static _lincomb _lincombtbl[] = 0046 { 0047 {0,3,3,2,5},{12,4,0,2,7},{4,0,0,0,1},{0,2,2,1,3},{13,3,0,1,6}, /* -.40 */ 0048 {7,5,0,2,6},{10,1,0,2,5},{4,3,0,3,5},{1,6,0,0,3},{4,2,0,0,2}, /* -.35 */ 0049 {0,4,2,1,4},{1,0,0,1,1},{0,9,0,1,5},{3,5,0,1,4},{6,1,0,1,3}, /* -.30 */ 0050 {0,3,0,2,3},{8,6,0,1,6},{11,2,0,1,5},{0,6,2,1,5},{22,1,0,0,7}, /* -.25 */ 0051 {9,5,0,0,5},{0,0,2,2,3},{0,5,1,0,3},{2,1,0,0,1},{0,3,4,1,5}, /* -.20 */ 0052 {0,1,8,0,6},{0,3,3,3,6},{5,6,0,2,6},{8,2,0,2,5},{0,13,0,1,7}, /* -.15 */ 0053 {12,3,0,0,5},{0,7,1,0,4},{0,2,1,4,5},{0,7,0,2,5},{3,3,0,2,4}, /* -.10 */ 0054 {20,0,0,0,6},{7,4,0,0,4},{10,0,0,0,3},{8,4,0,2,6},{0,0,0,0,0}, /* -.05 */ 0055 {0,0,0,0,0},{0,9,1,0,5},{13,5,0,2,8},{9,3,0,1,5},{0,7,4,1,7}, /* 0.00 */ 0056 {10,3,0,3,7},{0,1,5,0,4},{0,12,2,1,8},{10,2,0,0,4},{8,6,0,2,7},/* 0.05 */ 0057 {8,1,0,6,8},{7,0,0,1,3},{5,4,0,3,6},{22,1,0,1,8},{2,7,0,0,4}, /* 0.10 */ 0058 {2,2,0,4,5},{0,0,2,3,4},{0,5,1,1,4},{0,3,5,0,5},{14,6,0,1,8}, /* 0.15 */ 0059 {10,4,0,0,5},{3,1,0,3,4},{4,6,0,1,5},{0,4,0,0,2},{14,0,0,2,6}, /* 0.20 */ 0060 {0,13,1,0,7},{8,2,0,3,6},{4,0,0,2,3},{2,4,0,4,6},{0,2,2,3,5}, /* 0.25 */ 0061 {0,0,6,2,6},{0,5,5,0,6},{23,2,0,0,8},{0,1,2,0,2},{17,4,0,1,8}, /* 0.30 */ 0062 {0,10,3,0,7},{11,6,0,2,8},{4,3,0,5,7},{14,2,0,2,7},{10,0,0,1,4},/* 0.35 */ 0063 {8,4,0,3,7},{1,1,0,6,6},{11,0,0,3,6},{12,5,0,1,7},{2,1,0,8,8}, /* 0.40 */ 0064 {15,1,0,1,6},{0,7,5,0,7},{6,2,0,6,8},{2,0,0,5,5},{1,14,0,1,8}, /* 0.45 */ 0065 {0,12,3,0,8},{20,2,0,1,8},{0,8,0,0,4},{14,4,0,2,8},{3,4,0,0,3},/* 0.50 */ 0066 {0,1,4,3,6},{6,0,0,0,2},{4,4,0,2,5},{0,2,0,1,2},{8,0,1,2,5}, /* 0.55 */ 0067 {8,5,0,0,5},{15,3,0,1,7},{0,9,5,0,8},{0,0,3,2,4},{0,5,2,0,4}, /* 0.60 */ 0068 {5,3,0,1,4},{0,0,2,4,5},{13,6,0,0,7},{3,2,0,7,8},{16,2,0,0,6}, /* 0.65 */ 0069 {0,1,9,0,7},{3,6,0,0,4},{10,4,0,1,6},{0,0,0,8,7},{13,0,0,1,5}, /* 0.70 */ 0070 {0,6,7,0,8},{0,4,0,1,3},{0,2,4,0,4},{3,0,0,1,2},{15,5,0,1,8}, /* 0.75 */ 0071 {8,2,0,4,7},{2,9,0,1,6},{4,0,0,3,4},{2,4,0,5,7},{12,3,0,2,7}, /* 0.80 */ 0072 {1,3,0,0,2},{0,0,6,3,7},{6,5,0,3,7},{16,4,0,0,7},{0,3,9,0,8}, /* 0.85 */ 0073 {0,1,2,1,3},{10,6,0,1,7},{0,2,0,8,8},{0,10,3,1,8},{0,1,1,3,4}, /* 0.90 */ 0074 {9,0,0,0,3},{4,3,0,6,8},{0,4,4,0,5},{8,9,0,0,7},{10,0,0,2,5} /* 0.95 */ 0075 }; 0076 0077 /* artanh x and ln(x+1) are closely related. Evaluate the logarithm 0078 using the artanh series, which converges twice as fast as the 0079 logarithm series. 0080 With 100 digits, if -0.0198 <= x <= 0.02, the relative error is 0081 less than 5.0e-100, or at most 1 unit in the 100th place */ 0082 void 0083 _lnxplus1near0( 0084 floatnum x, 0085 int digits) 0086 { 0087 floatstruct tmp; 0088 0089 float_create(&tmp); 0090 float_add(&tmp, &c2, x, digits); 0091 float_div(x, x, &tmp, digits); 0092 artanhnear0(x, digits); 0093 float_add(x, x, x, EXACT); 0094 float_free(&tmp); 0095 } 0096 0097 /* helpers */ 0098 static void 0099 _addcoef( 0100 floatnum x, 0101 int coef, 0102 floatnum cnst, 0103 int digits) 0104 { 0105 floatstruct tmp; 0106 0107 if (coef != 0) 0108 { 0109 float_create(&tmp); 0110 float_muli(&tmp, cnst, coef, digits); 0111 float_add(x, x, &tmp, digits); 0112 float_free(&tmp); 0113 } 0114 } 0115 0116 static int 0117 _factor( 0118 int idx) 0119 { 0120 int factor; 0121 int i; 0122 0123 factor = 1; 0124 for (i = _lincombtbl[idx].c2; --i >= 0;) 0125 factor *= 2; 0126 for (i = _lincombtbl[idx].c3; --i >= 0;) 0127 factor *= 3; 0128 for (i = _lincombtbl[idx].c5; --i >= 0;) 0129 factor *= 5; 0130 for (i = _lincombtbl[idx].c7; --i >= 0;) 0131 factor *= 7; 0132 return factor; 0133 } 0134 0135 static void 0136 _lnguess( 0137 floatnum dest, 0138 int digits, 0139 int idx) 0140 { 0141 _addcoef(dest, _lincombtbl[idx].c2-_lincombtbl[idx].c5, &cLn2, digits); 0142 _addcoef(dest, _lincombtbl[idx].c3, &cLn3, digits); 0143 _addcoef(dest, _lincombtbl[idx].c7, &cLn7, digits); 0144 _addcoef(dest, _lincombtbl[idx].c5-_lincombtbl[idx].c10, &cLn10, digits); 0145 } 0146 0147 /* reduces x using a special factor whose prime factors 0148 are 2, 3, 5 and 7 only. x is multiplied by this 0149 factor, yielding a value near a power of ten. 0150 Then x is divided by this power of ten. 0151 The logarithm of this factor is returned in 0152 lnguess. Valid for -0.4 <= x < 1. Relative error < 5e-100 for 0153 100 digit result. 0154 This algorithm reduces x to a value < 0.01, which 0155 is appropriately small for submitting to a series evaluation */ 0156 void 0157 _lnreduce( 0158 floatnum x, 0159 floatnum lnguess, 0160 int digits) 0161 { 0162 floatstruct tmp1, tmp2; 0163 int idx; 0164 int expx; 0165 signed char pos; 0166 0167 float_setzero(lnguess); 0168 expx = float_getexponent(x); 0169 if (expx < -2) 0170 return; 0171 float_create(&tmp1); 0172 float_create(&tmp2); 0173 0174 idx = leadingdigits(x, 3 + expx) + 39; 0175 if (idx < 0) 0176 idx = 0; 0177 pos = idx >= 40? 1 : 0; 0178 idx += pos; 0179 float_setinteger(&tmp1, _factor(idx)); 0180 float_setexponent(&tmp1, -pos); 0181 float_sub(&tmp1, &tmp1, &c1, EXACT); 0182 float_mul(&tmp2, x, &tmp1, digits); 0183 float_add(&tmp1, &tmp1, &tmp2, digits+1); 0184 float_add(x, x, &tmp1, digits); 0185 _lnguess(lnguess, digits+4, idx); 0186 float_free(&tmp1); 0187 float_free(&tmp2); 0188 return; 0189 } 0190 0191 /* for -0.4 <= x < 1.0 */ 0192 void 0193 _lnxplus1lt1( 0194 floatnum x, 0195 int digits) 0196 { 0197 floatstruct lnfactor; 0198 0199 float_create(&lnfactor); 0200 _lnreduce(x, &lnfactor, digits); 0201 _lnxplus1near0(x, digits); 0202 float_sub(x, x, &lnfactor, digits+1); 0203 float_free(&lnfactor); 0204 } 0205 0206 /* the general purpose routine evaluating ln(x) for all 0207 positive arguments. It uses multiplication and division to 0208 reduce the argument to a value near 1. The factors are 0209 always small, so these operations do not take much time. 0210 -- 0211 One often sees an argument reduction by taking the 0212 square root several times in succession. But the square root is 0213 slow in number.c, so we avoid it best. 0214 -- 0215 Super-fast algorithm like the AGM based ones, do not 0216 pay off for limited lengths because of the same reason: They 0217 require taking square roots several times. A test showed, 0218 a hundred digits calculation using AGM is 8 times slower than 0219 the algorithm here. For extreme precision, of course, AGM will be 0220 superior. 0221 The relative error seems to be less than 7e-101 for 100-digits 0222 computations */ 0223 void 0224 _ln( 0225 floatnum x, 0226 int digits) 0227 { 0228 floatstruct tmp; 0229 int coef10; 0230 char coef3; 0231 char dgt; 0232 0233 float_create(&tmp); 0234 coef10 = float_getexponent(x); 0235 0236 /* reducing the significand to 0.6 <= x < 2 0237 by simple multiplication */ 0238 dgt = leadingdigits(x, 1); 0239 coef3 = 0; 0240 if (dgt == 1) 0241 float_setexponent(x, 0); 0242 else 0243 { 0244 ++coef10; 0245 float_setexponent(x, - 1); 0246 coef3 = (dgt < 6); 0247 if (coef3) 0248 float_mul(x, x, &c3, digits+1); 0249 } 0250 float_sub(x, x, &c1, digits+1); 0251 _lnxplus1lt1(x, digits); 0252 if (coef10 != 0) 0253 { 0254 float_muli(&tmp, &cLn10, coef10, digits+1); 0255 float_add(x, x, &tmp, digits+1); 0256 } 0257 if (coef3) 0258 float_sub(x, x, &cLn3, digits); 0259 float_free(&tmp); 0260 } 0261 0262 void 0263 _lnxplus1( 0264 floatnum x, 0265 int digits) 0266 { 0267 if (float_cmp(x, &cMinus0_4) >= 0 && float_cmp(x, &c1) < 0) 0268 _lnxplus1lt1(x, digits); 0269 else 0270 { 0271 float_add(x, x, &c1, digits+1); 0272 _ln(x, digits); 0273 } 0274 } 0275 0276 /* The artanh function has a pole at x == 1. 0277 For values close to 1 the general purpose 0278 evaluation is too unstable. We use a dedicated 0279 algorithm here to get around the problems. 0280 Avoid using this function for values of x > 0.5, 0281 _artanh is the better choice then. 0282 Values near the other pole of artanh at -1 can 0283 be derived from this function using 0284 artanh(-1+x) = -artanh(1-x). 0285 For x < 0.5 and a 100-digit computation, the maximum 0286 relative error is in the order of 1e-99 */ 0287 void 0288 _artanh1minusx( 0289 floatnum x, 0290 int digits) 0291 { 0292 floatstruct tmp; 0293 0294 float_create(&tmp); 0295 float_sub(&tmp, &c2, x, digits); 0296 float_div(x, &tmp, x, digits); 0297 _ln(x, digits); 0298 float_mul(x, x, &c1Div2, digits); 0299 float_free(&tmp); 0300 } 0301 0302 /* designed for |x| <= 0.5. The evaluation is not 0303 symmetric with respect to 0, i.e. it might be 0304 -artanh x != artanh -x. The difference 0305 is in the order of the last digit, of course. 0306 Evaluation of positive values yield a slightly 0307 better relative error than that of negative 0308 values. 0309 The maximum relative error is the maximum of 0310 that of artanhnear0 and 1e-99 */ 0311 0312 void 0313 _artanhlt0_5( 0314 floatnum x, 0315 int digits) 0316 { 0317 floatstruct tmp; 0318 0319 if (float_getexponent(x) < -2 || float_iszero(x)) 0320 artanhnear0(x, digits); 0321 else 0322 { 0323 float_create(&tmp); 0324 float_sub(&tmp, &c1, x, digits+1); 0325 float_add(x, x, x, digits); 0326 float_div(x, x, &tmp, digits); 0327 _lnxplus1(x, digits); 0328 float_mul(x, x, &c1Div2, digits); 0329 float_free(&tmp); 0330 } 0331 } 0332 0333 /* the general purpose routine for evaluating artanh. 0334 The evaluation is symmetric with respect to 0, i.e. 0335 it is always -artanh(x) == artanh -x. 0336 Valid for |x| < 1, but unstable for |x| approx. == 1 */ 0337 void 0338 _artanh( 0339 floatnum x, 0340 int digits) 0341 { 0342 signed char sgn; 0343 0344 sgn = float_getsign(x); 0345 float_abs(x); 0346 if (float_cmp(x, &c1Div2) <= 0) 0347 _artanhlt0_5(x, digits); 0348 else 0349 { 0350 float_sub(x, &c1, x, digits+1); 0351 _artanh1minusx(x, digits); 0352 } 0353 float_setsign(x, sgn); 0354 } 0355 0356 /* evaluates ln(2*x), the asymptotic function of arsinh and 0357 arcosh for large x */ 0358 static void 0359 _ln2x( 0360 floatnum x, 0361 int digits) 0362 { 0363 _ln(x, digits); 0364 float_add(x, &cLn2, x, digits); 0365 } 0366 0367 /* valid for all x. The relative error is less than 1e-99. 0368 The evaluation is symmetric with respect to 0: 0369 arsinh -x == - arsinh x */ 0370 void 0371 _arsinh( 0372 floatnum x, 0373 int digits) 0374 { 0375 floatstruct tmp; 0376 int expxsqr; 0377 signed char sgn; 0378 0379 expxsqr = 2*float_getexponent(x)+2; 0380 /* for extreme small x, sqrt(1+x*x) is approx == 1 + x*x/2 - x^4/8, 0381 so arsinh x == ln(1 + x + x*x/2 - x^4/8 +...) 0382 approx == x - x*x*x/6 + ... 0383 If we approximate arsinh x by x, the relative error is 0384 roughly |x*x/6|. This is less than acceptable 10^(-<digits>), if 0385 log(x*x) < -<digits> */ 0386 if (expxsqr < -digits || float_iszero(x)) 0387 return; 0388 float_create(&tmp); 0389 0390 /* arsinh -x = -arsinh x, so we can derive the 0391 arsinh for negative arguments from that 0392 of positive arguments */ 0393 sgn = float_getsign(x); 0394 float_abs(x); 0395 if (expxsqr-2 > digits) 0396 /* for very large x, use the asymptotic formula ln(2*x) */ 0397 _ln2x(x, digits); 0398 else 0399 { 0400 float_mul(&tmp, x, x, digits + (expxsqr>=0? 0 : expxsqr)); 0401 float_add(&tmp, &tmp, &c1, digits+1); 0402 float_sqrt(&tmp, digits); 0403 if (float_getexponent(x) < 0) 0404 { 0405 /* for small x, use arsinh x == artanh (x/sqrt(1+x*x)). 0406 Stable for x < 1, but not for large x*/ 0407 float_div(x, x, &tmp, digits); 0408 _artanh(x, digits+1); 0409 } 0410 else 0411 { 0412 /* arsinh x = ln(x+sqrt(1+x*x)), stable for x >= 1, but 0413 not for small x */ 0414 float_add(x, x, &tmp, digits); 0415 _ln(x, digits); 0416 } 0417 } 0418 float_setsign(x, sgn); 0419 float_free(&tmp); 0420 } 0421 0422 /* arcosh(x+1), x >= 0, is the stable variant of 0423 arcosh(x), x >= 1. 0424 The relative error is less than 1e-99. */ 0425 void 0426 _arcoshxplus1( 0427 floatnum x, 0428 int digits) 0429 { 0430 floatstruct tmp; 0431 0432 float_create(&tmp); 0433 if (2*float_getexponent(x) > digits) 0434 { 0435 /* for very large x, use the asymptotic formula ln(2*(x+1)) */ 0436 float_add(x, x, &c1, digits+1); 0437 _ln2x(x, digits); 0438 } 0439 else 0440 { 0441 /* arcosh(x+1) = ln(1+(x+sqrt(x*(x+2)))), stable 0442 for all positive x, except for extreme large x, where 0443 x*(x+2) might overflow */ 0444 0445 /* get sinh(arcosh (1+x)) = sqrt(x*(x+2)) */ 0446 float_add(&tmp, x, &c2, digits); 0447 float_mul(&tmp, x, &tmp, digits); 0448 float_sqrt(&tmp, digits); 0449 0450 float_add(x, x, &tmp, digits); 0451 _lnxplus1(x, digits); 0452 } 0453 float_free(&tmp); 0454 }