double-math.c
You can download the current version of Algol 68 Genie and its documentation here.
1 //! @file double-math.c
2 //! @author J. Marcel van der Veer
3 //!
4 //! @section Copyright
5 //!
6 //! This file is part of Algol68G - an Algol 68 compiler-interpreter.
7 //! Copyright 2001-2023 J. Marcel van der Veer .
8 //!
9 //! @section License
10 //!
11 //! This program is free software; you can redistribute it and/or modify it
12 //! under the terms of the GNU General Public License as published by the
13 //! Free Software Foundation; either version 3 of the License, or
14 //! (at your option) any later version.
15 //!
16 //! This program is distributed in the hope that it will be useful, but
17 //! WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
18 //! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
19 //! more details. You should have received a copy of the GNU General Public
20 //! License along with this program. If not, see .
21
22 //! @section Synopsis
23 //!
24 //! LONG REAL, LONG COMPLEX routines.
25
26 // References:
27 //
28 // Milton Abramowitz and Irene Stegun, Handbook of Mathematical Functions,
29 // Dover Publications, New York [1970]
30 // https://en.wikipedia.org/wiki/Abramowitz_and_Stegun
31
32 #include "a68g.h"
33 #include "a68g-genie.h"
34 #include "a68g-prelude.h"
35 #include "a68g-double.h"
36 #include "a68g-numbers.h"
37 #include "a68g-math.h"
38
39 #if (A68_LEVEL >= 3)
40
41 DOUBLE_T a68_beta_inc_16 (DOUBLE_T s, DOUBLE_T t, DOUBLE_T x)
42 {
43 // Incomplete beta function I{x}(s, t).
44 // Continued fraction, see dlmf.nist.gov/8.17; Lentz's algorithm.
45 if (x < 0.0q || x > 1.0q) {
46 errno = ERANGE;
47 return -1.0q;
48 } else {
49 const INT_T lim = 16 * sizeof (DOUBLE_T);
50 // Rapid convergence when x <= (s+1)/(s+t+2) or else recursion.
51 if (x > (s + 1.0q) / (s + t + 2.0q)) {
52 // B{x}(s, t) = 1 - B{1-x}(t, s)
53 return 1.0q - a68_beta_inc_16 (s, t, 1.0q - x);
54 }
55 // Lentz's algorithm for continued fraction.
56 DOUBLE_T W = 1.0q, F = 1.0q, c = 1.0q, d = 0.0q;
57 BOOL_T cont = A68_TRUE;
58 for (int N = 0, m = 0; cont && N < lim; N++) {
59 DOUBLE_T T;
60 if (N == 0) {
61 T = 1.0q;
62 } else if (N % 2 == 0) {
63 // d{2m} := x m(t-m)/((s+2m-1)(s+2m))
64 T = x * m * (t - m) / (s + 2.0q * m - 1.0q) / (s + 2.0q * m);
65 } else {
66 // d{2m+1} := -x (s+m)(s+t+m)/((s+2m+1)(s+2m))
67 T = -x * (s + m) * (s + t + m) / (s + 2.0q * m + 1.0q) / (s + 2.0q * m);
68 m++;
69 }
70 d = 1.0q / (T * d + 1.0q);
71 c = T / c + 1.0q;
72 F *= c * d;
73 if (F == W) {
74 cont = A68_FALSE;
75 } else {
76 W = F;
77 }
78 }
79 // I{x}(s,t)=x^s(1-x)^t / s / B(s,t) F
80 DOUBLE_T beta = expq (lgammaq (s) + lgammaq (t) - lgammaq (s + t));
81 return powq (x, s) * powq (1.0q - x, t) / s / beta * (F - 1.0q);
82 }
83 }
84
85 //! @brief PROC (LONG REAL) LONG REAL csc
86
87 DOUBLE_T a68_csc_16 (DOUBLE_T x)
88 {
89 DOUBLE_T z = sinq (x);
90 A68_OVERFLOW (z == 0.0q);
91 return 1.0q / z;
92 }
93
94 //! @brief PROC (LONG REAL) LONG REAL acsc
95
96 DOUBLE_T a68_acsc_16 (DOUBLE_T x)
97 {
98 A68_OVERFLOW (x == 0.0q);
99 return asinq (1.0q / x);
100 }
101
102 //! @brief PROC (LONG REAL) LONG REAL sec
103
104 DOUBLE_T a68_sec_16 (DOUBLE_T x)
105 {
106 DOUBLE_T z = cosq (x);
107 A68_OVERFLOW (z == 0.0q);
108 return 1.0q / z;
109 }
110
111 //! @brief PROC (LONG REAL) LONG REAL asec
112
113 DOUBLE_T a68_asec_16 (DOUBLE_T x)
114 {
115 A68_OVERFLOW (x == 0.0q);
116 return acosq (1.0q / x);
117 }
118
119 //! @brief PROC (LONG REAL) LONG REAL cot
120
121 DOUBLE_T a68_cot_16 (DOUBLE_T x)
122 {
123 DOUBLE_T z = sinq (x);
124 A68_OVERFLOW (z == 0.0q);
125 return cosq (x) / z;
126 }
127
128 //! @brief PROC (LONG REAL) LONG REAL acot
129
130 DOUBLE_T a68_acot_16 (DOUBLE_T x)
131 {
132 A68_OVERFLOW (x == 0.0q);
133 return atanq (1 / x);
134 }
135
136 //! brief PROC (LONG REAL) LONG REAL sindg
137
138 DOUBLE_T a68_sindg_16 (DOUBLE_T x)
139 {
140 return sin (x * CONST_PI_OVER_180_Q);
141 }
142
143 //! brief PROC (LONG REAL) LONG REAL cosdg
144
145 DOUBLE_T a68_cosdg_16 (DOUBLE_T x)
146 {
147 return cos (x * CONST_PI_OVER_180_Q);
148 }
149
150 //! brief PROC (LONG REAL) LONG REAL tandg
151
152 DOUBLE_T a68_tandg_16 (DOUBLE_T x)
153 {
154 return tan (x * CONST_PI_OVER_180_Q);
155 }
156
157 //! brief PROC (LONG REAL) LONG REAL asindg
158
159 DOUBLE_T a68_asindg_16 (DOUBLE_T x)
160 {
161 return asin (x) * CONST_180_OVER_PI_Q;
162 }
163
164 //! brief PROC (LONG REAL) LONG REAL acosdg
165
166 DOUBLE_T a68_acosdg_16 (DOUBLE_T x)
167 {
168 return acos (x) * CONST_180_OVER_PI_Q;
169 }
170
171 //! brief PROC (LONG REAL) LONG REAL atandg
172
173 DOUBLE_T a68_atandg_16 (DOUBLE_T x)
174 {
175 return atan (x) * CONST_180_OVER_PI_Q;
176 }
177
178 // PROC (LONG REAL) LONG REAL cotdg
179
180 DOUBLE_T a68_cotdg_16 (DOUBLE_T x)
181 {
182 DOUBLE_T z = a68_sindg (x);
183 A68_OVERFLOW (z == 0);
184 return a68_cosdg (x) / z;
185 }
186
187 // PROC (LONG REAL) LONG REAL acotdg
188
189 DOUBLE_T a68_acotdg_16 (DOUBLE_T z)
190 {
191 A68_OVERFLOW (z == 0);
192 return a68_atandg (1 / z);
193 }
194
195 // @brief PROC (LONG REAL) LONG REAL sinpi
196
197 DOUBLE_T a68_sinpi_16 (DOUBLE_T x)
198 {
199 x = fmodq (x, 2.0q);
200 if (x <= -1.0q) {
201 x += 2.0q;
202 } else if (x > 1.0q) {
203 x -= 2.0q;
204 }
205 // x in <-1, 1].
206 if (x == 0.0q || x == 1.0q) {
207 return 0.0q;
208 } else if (x == 0.5q) {
209 return 1.0q;
210 }
211 if (x == -0.5q) {
212 return -1.0q;
213 } else {
214 return sinq (CONST_PI_Q * x);
215 }
216 }
217
218 // @brief PROC (LONG REAL) LONG REAL cospi
219
220 DOUBLE_T a68_cospi_16 (DOUBLE_T x)
221 {
222 x = fmodq (fabsq (x), 2.0q);
223 // x in [0, 2>.
224 if (x == 0.5q || x == 1.5q) {
225 return 0.0q;
226 } else if (x == 0.0q) {
227 return 1.0q;
228 } else if (x == 1.0q) {
229 return -1.0q;
230 } else {
231 return cosq (CONST_PI_Q * x);
232 }
233 }
234
235 // @brief PROC (LONG REAL) LONG REAL tanpi
236
237 DOUBLE_T a68_tanpi_16 (DOUBLE_T x)
238 {
239 x = fmodq (x, 1.0q);
240 if (x <= -0.5q) {
241 x += 1.0q;
242 } else if (x > 0.5q) {
243 x -= 1.0q;
244 }
245 // x in <-1/2, 1/2].
246 A68_OVERFLOW (x == 0.5q);
247 if (x == -0.25q) {
248 return -1.0q;
249 } else if (x == 0) {
250 return 0.0q;
251 } else if (x == 0.25q) {
252 return 1.0q;
253 } else {
254 return a68_sinpi_16 (x) / a68_cospi_16 (x);
255 }
256 }
257
258 // @brief PROC (LONG REAL) LONG REAL cotpi
259
260 DOUBLE_T a68_cotpi_16 (DOUBLE_T x)
261 {
262 x = fmodq (x, 1.0q);
263 if (x <= -0.5q) {
264 x += 1.0q;
265 } else if (x > 0.5q) {
266 x -= 1.0q;
267 }
268 // x in <-1/2, 1/2].
269 A68_OVERFLOW (x == 0.0q);
270 if (x == -0.25q) {
271 return -1.0q;
272 } else if (x == 0.25q) {
273 return 1.0q;
274 } else if (x == 0.5q) {
275 return 0.0q;
276 } else {
277 return a68_cospi_16 (x) / a68_sinpi_16 (x);
278 }
279 }
280
281 #endif
jmvdveer@xs4all.nl