1 /* floatseries.c: basic series, based on floatnum. */
2 /*
3 Copyright (C) 2007 Wolf Lammen.
4
5 This program is free software; you can redistribute it and/or modify
6 it under the terms of the GNU General Public License as published by
7 the Free Software Foundation; either version 2 of the License , or
8 (at your option) any later version.
9
10 This program is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU General Public License for more details.
14
15 You should have received a copy of the GNU General Public License
16 along with this program; see the file COPYING. If not, write to:
17
18 The Free Software Foundation, Inc.
19 59 Temple Place, Suite 330
20 Boston, MA 02111-1307 USA.
21
22
23 You may contact the author by:
24 e-mail: ookami1 <at> gmx <dot> de
25 mail: Wolf Lammen
26 Oertzweg 45
27 22307 Hamburg
28 Germany
29 *************************************************************************/
30
31 #include "floatseries.h"
32 #include "floatconst.h"
33 #include "floatcommon.h"
34 #include "floatexp.h"
35
36 /* Though all these serieses can be used with arguments |x| < 1 or
37 even more, the computation time increases rapidly with x.
38 Tests show, that for 100 digit results, it is best to limit x
39 to |x| < 0.01..0.02, and use reduction formulas for greater x */
40
41 /* the Taylor series of arctan/arctanh x at x == 0. For small
42 |x| < 0.01 this series converges very fast, yielding 4 or
43 more digits of the result with every summand. The working
44 precision is adjusted, so that the relative error for
45 100-digit arguments is around 5.0e-100. This means, the error
46 is 1 in the 100-th place (or less) */
47 void
arctanseries(floatnum x,int digits,char alternating)48 arctanseries(
49 floatnum x,
50 int digits,
51 char alternating)
52 {
53 int expx;
54 int expsqrx;
55 int pwrsz;
56 int addsz;
57 int i;
58 floatstruct xsqr;
59 floatstruct pwr;
60 floatstruct smd;
61 floatstruct sum;
62
63 /* upper limit of log(x) and log(result) */
64 expx = float_getexponent(x)+1;
65
66 /* the summands of the series from the second on are
67 bounded by x^(2*i-1)/3. So the summation yields a
68 result bounded by (x^3/(1-x*x))/3.
69 For x < sqrt(1/3) approx.= 0.5, this is less than 0.5*x^3.
70 We need to sum up only, if the first <digits> places
71 of the result (roughly x) are touched. Ignoring the effect of
72 a possile carry, this is only the case, if
73 x*x >= 2*10^(-digits) > 10^(-digits)
74 Example: for x = 9e-51, a 100-digits result covers
75 the decimal places from 1e-51 to 1e-150. x^3/3
76 is roughly 3e-151, and so is the sum of the series.
77 So we can ignore the sum, but we couldn't for x = 9e-50 */
78 if (float_iszero(x) || 2*expx < -digits)
79 /* for very tiny arguments arctan/arctanh x is approx.== x */
80 return;
81 float_create(&xsqr);
82 float_create(&pwr);
83 float_create(&smd);
84 float_create(&sum);
85
86 /* we adapt the working precision to the decreasing
87 summands, saving time when multiplying. Unfortunately,
88 there is no error bound given for the operations of
89 bc_num. Tests show, that the last digit in an incomplete
90 multiplication is usually not correct up to 5 ULP's. */
91 pwrsz = digits + 2*expx + 1;
92 /* the precision of the addition must not decrease, of course */
93 addsz = pwrsz;
94 i = 3;
95 float_mul(&xsqr, x, x, pwrsz);
96 float_setsign(&xsqr, alternating? -1 : 1);
97 expsqrx = float_getexponent(&xsqr);
98 float_copy(&pwr, x, pwrsz);
99 float_setzero(&sum);
100
101 for(; pwrsz > 0; )
102 {
103 /* x^i */
104 float_mul(&pwr, &pwr, &xsqr, pwrsz+1);
105 /* x^i/i */
106 float_divi(&smd, &pwr, i, pwrsz);
107 /* The addition virtually does not introduce errors */
108 float_add(&sum, &sum, &smd, addsz);
109 /* reduce the working precision according to the decreasing powers */
110 pwrsz = digits - expx + float_getexponent(&smd) + expsqrx + 3;
111 i += 2;
112 }
113 /* add the first summand */
114 float_add(x, x, &sum, digits+1);
115
116 float_free(&xsqr);
117 float_free(&pwr);
118 float_free(&smd);
119 float_free(&sum);
120 }
121
122 /* series expansion of cos/cosh - 1 used for small x,
123 |x| <= 0.01.
124 The function returns 0, if an underflow occurs.
125 The relative error seems to be less than 5e-100 for
126 a 100-digit calculation with |x| < 0.01 */
127 char
cosminus1series(floatnum x,int digits,char alternating)128 cosminus1series(
129 floatnum x,
130 int digits,
131 char alternating)
132 {
133 floatstruct sum, smd;
134 int expsqrx, pwrsz, addsz, i;
135
136 expsqrx = 2 * float_getexponent(x);
137 float_setexponent(x, 0);
138 float_mul(x, x, x, digits+1);
139 float_mul(x, x, &c1Div2, digits+1);
140 float_setsign(x, alternating? -1 : 1);
141 expsqrx += float_getexponent(x);
142 if (float_iszero(x) || expsqrx < EXPMIN)
143 {
144 /* underflow */
145 float_setzero(x);
146 return expsqrx == 0;
147 }
148 float_setexponent(x, expsqrx);
149 pwrsz = digits + expsqrx + 2;
150 if (pwrsz <= 0)
151 /* for very small x, cos/cosh(x) - 1 = (-/+)0.5*x*x */
152 return 1;
153 addsz = pwrsz;
154 float_create(&sum);
155 float_create(&smd);
156 float_copy(&smd, x, pwrsz);
157 float_setzero(&sum);
158 i = 2;
159 while (pwrsz > 0)
160 {
161 float_mul(&smd, &smd, x, pwrsz+1);
162 float_divi(&smd, &smd, i*(2*i-1), pwrsz);
163 float_add(&sum, &sum, &smd, addsz);
164 ++i;
165 pwrsz = digits + float_getexponent(&smd);
166 }
167 float_add(x, x, &sum, digits+1);
168 float_free(&sum);
169 float_free(&smd);
170 return 1;
171 }
172