File indexing completed on 2024-05-12 05:55:15
0001 /* floatseries.c: basic series, based on floatnum. */ 0002 /* 0003 Copyright (C) 2007 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 "floatseries.h" 0032 #include "floatconst.h" 0033 #include "floatcommon.h" 0034 #include "floatexp.h" 0035 0036 /* Though all these serieses can be used with arguments |x| < 1 or 0037 even more, the computation time increases rapidly with x. 0038 Tests show, that for 100 digit results, it is best to limit x 0039 to |x| < 0.01..0.02, and use reduction formulas for greater x */ 0040 0041 /* the Taylor series of arctan/arctanh x at x == 0. For small 0042 |x| < 0.01 this series converges very fast, yielding 4 or 0043 more digits of the result with every summand. The working 0044 precision is adjusted, so that the relative error for 0045 100-digit arguments is around 5.0e-100. This means, the error 0046 is 1 in the 100-th place (or less) */ 0047 void 0048 arctanseries( 0049 floatnum x, 0050 int digits, 0051 char alternating) 0052 { 0053 int expx; 0054 int expsqrx; 0055 int pwrsz; 0056 int addsz; 0057 int i; 0058 floatstruct xsqr; 0059 floatstruct pwr; 0060 floatstruct smd; 0061 floatstruct sum; 0062 0063 /* upper limit of log(x) and log(result) */ 0064 expx = float_getexponent(x)+1; 0065 0066 /* the summands of the series from the second on are 0067 bounded by x^(2*i-1)/3. So the summation yields a 0068 result bounded by (x^3/(1-x*x))/3. 0069 For x < sqrt(1/3) approx.= 0.5, this is less than 0.5*x^3. 0070 We need to sum up only, if the first <digits> places 0071 of the result (roughly x) are touched. Ignoring the effect of 0072 a possile carry, this is only the case, if 0073 x*x >= 2*10^(-digits) > 10^(-digits) 0074 Example: for x = 9e-51, a 100-digits result covers 0075 the decimal places from 1e-51 to 1e-150. x^3/3 0076 is roughly 3e-151, and so is the sum of the series. 0077 So we can ignore the sum, but we couldn't for x = 9e-50 */ 0078 if (float_iszero(x) || 2*expx < -digits) 0079 /* for very tiny arguments arctan/arctanh x is approx.== x */ 0080 return; 0081 float_create(&xsqr); 0082 float_create(&pwr); 0083 float_create(&smd); 0084 float_create(&sum); 0085 0086 /* we adapt the working precision to the decreasing 0087 summands, saving time when multiplying. Unfortunately, 0088 there is no error bound given for the operations of 0089 bc_num. Tests show, that the last digit in an incomplete 0090 multiplication is usually not correct up to 5 ULP's. */ 0091 pwrsz = digits + 2*expx + 1; 0092 /* the precision of the addition must not decrease, of course */ 0093 addsz = pwrsz; 0094 i = 3; 0095 float_mul(&xsqr, x, x, pwrsz); 0096 float_setsign(&xsqr, alternating? -1 : 1); 0097 expsqrx = float_getexponent(&xsqr); 0098 float_copy(&pwr, x, pwrsz); 0099 float_setzero(&sum); 0100 0101 for(; pwrsz > 0; ) 0102 { 0103 /* x^i */ 0104 float_mul(&pwr, &pwr, &xsqr, pwrsz+1); 0105 /* x^i/i */ 0106 float_divi(&smd, &pwr, i, pwrsz); 0107 /* The addition virtually does not introduce errors */ 0108 float_add(&sum, &sum, &smd, addsz); 0109 /* reduce the working precision according to the decreasing powers */ 0110 pwrsz = digits - expx + float_getexponent(&smd) + expsqrx + 3; 0111 i += 2; 0112 } 0113 /* add the first summand */ 0114 float_add(x, x, &sum, digits+1); 0115 0116 float_free(&xsqr); 0117 float_free(&pwr); 0118 float_free(&smd); 0119 float_free(&sum); 0120 } 0121 0122 /* series expansion of cos/cosh - 1 used for small x, 0123 |x| <= 0.01. 0124 The function returns 0, if an underflow occurs. 0125 The relative error seems to be less than 5e-100 for 0126 a 100-digit calculation with |x| < 0.01 */ 0127 char 0128 cosminus1series( 0129 floatnum x, 0130 int digits, 0131 char alternating) 0132 { 0133 floatstruct sum, smd; 0134 int expsqrx, pwrsz, addsz, i; 0135 0136 expsqrx = 2 * float_getexponent(x); 0137 float_setexponent(x, 0); 0138 float_mul(x, x, x, digits+1); 0139 float_mul(x, x, &c1Div2, digits+1); 0140 float_setsign(x, alternating? -1 : 1); 0141 expsqrx += float_getexponent(x); 0142 if (float_iszero(x) || expsqrx < EXPMIN) 0143 { 0144 /* underflow */ 0145 float_setzero(x); 0146 return expsqrx == 0; 0147 } 0148 float_setexponent(x, expsqrx); 0149 pwrsz = digits + expsqrx + 2; 0150 if (pwrsz <= 0) 0151 /* for very small x, cos/cosh(x) - 1 = (-/+)0.5*x*x */ 0152 return 1; 0153 addsz = pwrsz; 0154 float_create(&sum); 0155 float_create(&smd); 0156 float_copy(&smd, x, pwrsz); 0157 float_setzero(&sum); 0158 i = 2; 0159 while (pwrsz > 0) 0160 { 0161 float_mul(&smd, &smd, x, pwrsz+1); 0162 float_divi(&smd, &smd, i*(2*i-1), pwrsz); 0163 float_add(&sum, &sum, &smd, addsz); 0164 ++i; 0165 pwrsz = digits + float_getexponent(&smd); 0166 } 0167 float_add(x, x, &sum, digits+1); 0168 float_free(&sum); 0169 float_free(&smd); 0170 return 1; 0171 }