File indexing completed on 2024-05-12 17:22:33
0001 /* floatexp.c: exponential function 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 "floatconst.h" 0032 #include "floatcommon.h" 0033 #include "floatseries.h" 0034 #include "floatexp.h" 0035 0036 /* uses the addition theorem 0037 cosh(2x)-1 == 2*(cosh x - 1)*(cosh x + 1) 0038 to reduce the argument to the range |x| < 0.01. 0039 Starting with x == 1, you need 7 reduction steps 0040 to achieve the desired magnitude. 0041 The relative error is < 8e-100 for a 100 digit result. 0042 The return value is 0, if the result underflows. 0043 |x| < 1, otherwise the final process, where 0044 the reductions are unwinded, becomes too 0045 unstable */ 0046 char 0047 _coshminus1lt1( 0048 floatnum x, 0049 int digits) 0050 { 0051 floatstruct tmp; 0052 int reductions; 0053 0054 if (float_iszero(x)) 0055 return 1; 0056 float_abs(x); 0057 reductions = 0; 0058 while(float_getexponent(x) >= -2) 0059 { 0060 float_mul(x, x, &c1Div2, digits+1); 0061 ++reductions; 0062 } 0063 if (!coshminus1near0(x, digits) && reductions == 0) 0064 return !float_iszero(x); 0065 float_create(&tmp); 0066 for(; reductions-- > 0;) 0067 { 0068 float_mul(&tmp, x, x, digits); 0069 float_add(x, x, x, digits+2); 0070 float_add(x, x, &tmp, digits+2); 0071 float_add(x, x, x, digits+2); 0072 } 0073 float_free(&tmp); 0074 return 1; 0075 } 0076 0077 /* sinh x == sqrt((cosh x - 1) * (cosh x + 1)) */ 0078 static void 0079 _sinhfromcoshminus1( 0080 floatnum x, 0081 int digits) 0082 { 0083 floatstruct tmp; 0084 0085 float_create(&tmp); 0086 float_add(&tmp, x, &c2, digits); 0087 float_mul(x, &tmp, x, digits+1); 0088 float_sqrt(x, digits+1); 0089 float_free(&tmp); 0090 } 0091 0092 /* sinh x for |x| < 1. Derived from cosh x - 1. 0093 The relative error is < 8e-100 for a 100 digit result */ 0094 void 0095 _sinhlt1( 0096 floatnum x, 0097 int digits) 0098 { 0099 signed char sgn; 0100 if (float_getexponent(x) < -digits) 0101 /* for very small x: sinh(x) approx. == x. */ 0102 return; 0103 sgn = float_getsign(x); 0104 _coshminus1lt1(x, digits); 0105 _sinhfromcoshminus1(x, digits); 0106 float_setsign(x, sgn); 0107 } 0108 0109 /* evaluates exp(x) - 1. This value can be obtained 0110 by exp(x) - 1 == sinh(x) + cosh(x) - 1 0111 relative error < 8e-100 for a 100 digit result */ 0112 void 0113 _expminus1lt1( 0114 floatnum x, 0115 int digits) 0116 { 0117 floatstruct tmp; 0118 signed char sgn; 0119 0120 if (float_getexponent(x) < -digits || float_iszero(x)) 0121 /* for very small x: exp(x)-1 approx.== x */ 0122 return; 0123 float_create(&tmp); 0124 sgn = float_getsign(x); 0125 _coshminus1lt1(x, digits); 0126 float_copy(&tmp, x, EXACT); 0127 _sinhfromcoshminus1(x, digits); 0128 float_setsign(x, sgn); 0129 float_add(x, x, &tmp, digits+1); 0130 float_free(&tmp); 0131 } 0132 0133 /* exp(x) for 0 <= x < ln 10 0134 relative error < 5e-100 */ 0135 void 0136 _expltln10( 0137 floatnum x, 0138 int digits) 0139 { 0140 int expx; 0141 int factor; 0142 char sgnf; 0143 0144 expx = float_getexponent(x); 0145 factor = 1; 0146 if (expx >= -1) 0147 { 0148 sgnf = leadingdigits(x, 2 + expx); 0149 if (sgnf > 4) 0150 { 0151 if (sgnf < 9) 0152 { 0153 factor = 2; 0154 float_sub(x, x, &cLn2, digits+1); 0155 } 0156 else if (sgnf < 14) 0157 { 0158 factor = 3; 0159 float_sub(x, x, &cLn3, digits+1); 0160 } 0161 else if (sgnf < 21) 0162 { 0163 factor = 7; 0164 float_sub(x, x, &cLn7, digits+1); 0165 } 0166 else 0167 { 0168 factor = 10; 0169 float_sub(x, x, &cLn10, digits+1); 0170 } 0171 } 0172 } 0173 _expminus1lt1(x, digits); 0174 float_add(x, x, &c1, digits+1); 0175 if (factor != 1) 0176 float_muli(x, x, factor, digits+1); 0177 } 0178 0179 /* exp(x) for all x. Underflow or overflow is 0180 indicated by the return value (0, if error) 0181 relative error for 100 digit results is 5e-100 */ 0182 char 0183 _exp( 0184 floatnum x, 0185 int digits) 0186 { 0187 floatstruct exp, tmp; 0188 int expx, extra; 0189 char ok; 0190 0191 if (float_iszero(x)) 0192 { 0193 float_copy(x, &c1, EXACT); 0194 return 1; 0195 } 0196 expx = float_getexponent(x); 0197 if (expx >= (int)(BITS_IN_EXP >> 1)) 0198 /* obvious overflow or underflow */ 0199 return 0; 0200 0201 float_create(&exp); 0202 float_create(&tmp); 0203 float_setzero(&exp); 0204 0205 if (expx >= 0) 0206 { 0207 float_div(&exp, x, &cLn10, expx+1); 0208 float_int(&exp); 0209 extra = float_getexponent(&exp)+1; 0210 float_mul(&tmp, &exp, &cLn10, digits+extra); 0211 float_sub(x, x, &tmp, digits+extra); 0212 if (float_cmp(x, &cLn10) >= 0) 0213 { 0214 float_add(&exp, &exp, &c1, EXACT); 0215 float_sub(x, x, &cLn10, digits); 0216 } 0217 } 0218 if (float_getsign(x) < 0) 0219 { 0220 float_sub(&exp, &exp, &c1, EXACT); 0221 float_add(x, x, &cLn10, digits); 0222 } 0223 /* when we get here 0 <= x < ln 10 */ 0224 _expltln10(x, digits); 0225 /* just in case rounding leads to a value >= 10 */ 0226 expx = float_getexponent(x); 0227 if (expx != 0) 0228 float_addi(&exp, &exp, expx, EXACT); 0229 0230 ok = 1; 0231 if (!float_iszero(&exp)) 0232 { 0233 expx = float_asinteger(&exp); 0234 ok = expx != 0; 0235 float_setexponent(x, expx); 0236 } 0237 float_free(&exp); 0238 float_free(&tmp); 0239 return ok && !float_isnan(x); 0240 } 0241 0242 static char 0243 _0_5exp( 0244 floatnum x, 0245 int digits) 0246 { 0247 float_sub(x, x, &cLn2, digits + (3*logexp(x)/10)+1); 0248 return _exp(x, digits); 0249 } 0250 0251 /* exp(x)-1 for all x. Overflow is 0252 indicated by the return value (0, if error) 0253 relative error for 100 digit results is 8e-100 */ 0254 char 0255 _expminus1( 0256 floatnum x, 0257 int digits) 0258 { 0259 int expr; 0260 0261 if (float_abscmp(x, &c1Div2) < 0) 0262 { 0263 _expminus1lt1(x, digits); 0264 return 1; 0265 } 0266 if (float_getsign(x) < 0) 0267 { 0268 expr = (2*float_getexponent(x)/5); 0269 if (expr >= digits) 0270 float_setinteger(x, -1); 0271 else 0272 { 0273 _exp(x, digits-expr); 0274 float_sub(x, x, &c1, digits); 0275 } 0276 return 1; 0277 } 0278 if (!_exp(x, digits)) 0279 return 0; 0280 float_sub(x, x, &c1, digits); 0281 return 1; 0282 } 0283 0284 static void 0285 _addreciproc( 0286 floatnum x, 0287 int digits, 0288 signed char sgn) 0289 { 0290 floatstruct tmp; 0291 int expx; 0292 0293 expx = float_getexponent(x); 0294 if (2*expx < digits) 0295 { 0296 float_create(&tmp); 0297 float_muli(&tmp, x, 4, digits-2*expx); 0298 float_reciprocal(&tmp, digits-2*expx); 0299 float_setsign(&tmp, sgn); 0300 float_add(x, x, &tmp, digits+1); 0301 float_free(&tmp); 0302 } 0303 } 0304 0305 /* cosh(x)-1 for all x. Underflow or overflow is 0306 indicated by the return value (0, if error) 0307 relative error for 100 digit results is 6e-100 */ 0308 0309 char 0310 _coshminus1( 0311 floatnum x, 0312 int digits) 0313 { 0314 if (float_getexponent(x) < 0 || float_iszero(x)) 0315 return _coshminus1lt1(x, digits); 0316 if(!_0_5exp(x, digits)) 0317 return 0; 0318 _addreciproc(x, digits, 1); 0319 float_sub(x, x, &c1, digits); 0320 return 1; 0321 } 0322 0323 /* sinh(x) for all x. Overflow is 0324 indicated by the return value (0, if error) 0325 relative error for 100 digit results is < 8e-100 */ 0326 char 0327 _sinh( 0328 floatnum x, 0329 int digits) 0330 { 0331 if (float_getexponent(x) < 0 || float_iszero(x)) 0332 _sinhlt1(x, digits); 0333 else 0334 { 0335 if(!_0_5exp(x, digits)) 0336 return 0; 0337 _addreciproc(x, digits, -1); 0338 } 0339 return 1; 0340 } 0341 0342 /* tanh(x) for |x| <= 0.5. 0343 relative error for 100 digit results is < 7e-100 */ 0344 void 0345 _tanhlt0_5( 0346 floatnum x, 0347 int digits) 0348 { 0349 floatstruct tmp; 0350 signed char sgn; 0351 0352 float_create(&tmp); 0353 sgn = float_getsign(x); 0354 float_abs(x); 0355 float_add(x, x, x, digits+1); 0356 _expminus1lt1(x, digits); 0357 float_add(&tmp, x, &c2, digits); 0358 float_div(x, x, &tmp, digits); 0359 float_setsign(x, sgn); 0360 float_free(&tmp); 0361 } 0362 0363 /* tanh(x)-1 for x > 0. 0364 relative error for 100 digit results is < 9e-100 */ 0365 char 0366 _tanhminus1gt0( 0367 floatnum x, 0368 int digits) 0369 { 0370 if (float_add(x, x, x, digits+1) && _0_5exp(x, digits)) 0371 { 0372 float_add(x, x, &c1Div2, digits+1); 0373 float_reciprocal(x, digits+1); 0374 float_setsign(x, -1); 0375 return 1; 0376 } 0377 return 0; 0378 } 0379 0380 void 0381 _tanhgt0_5( 0382 floatnum x, 0383 int digits) 0384 { 0385 int expx; 0386 0387 expx = float_getexponent(x); 0388 if (5*expx >= digits) 0389 float_copy(x, &c1, EXACT); 0390 else 0391 { 0392 _tanhminus1gt0(x, digits - 5*expx); 0393 float_add(x, x, &c1, digits); 0394 } 0395 } 0396 0397 char 0398 _power10( 0399 floatnum x, 0400 int digits) 0401 { 0402 int exp; 0403 0404 if (float_isinteger(x)) 0405 { 0406 exp = float_asinteger(x); 0407 if (exp == 0 && !float_iszero(x)) 0408 return 0; 0409 float_copy(x, &c1, EXACT); 0410 float_setexponent(x, exp); 0411 return !float_isnan(x); 0412 } 0413 return float_mul(x, x, &cLn10, digits+2) && _exp(x, digits); 0414 }