double-math.c
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-2024 J. Marcel van der Veer [algol68g@xs4all.nl].
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 [http://www.gnu.org/licenses/].
21
22 //! @section Synopsis
23 //!
24 //! LONG REAL, LONG COMPLEX routines.
25
26 // References:
27 // M. Abramowitz and I. Stegun, Handbook of Mathematical Functions,
28 // Dover Publications, New York [1970]
29 // https://en.wikipedia.org/wiki/Abramowitz_and_Stegun
30
31 #include "a68g.h"
32 #include "a68g-double.h"
33 #include "a68g-numbers.h"
34
35 #if (A68_LEVEL >= 3)
36
37 // Infinity.
38
39 DOUBLE_T a68_div_double (DOUBLE_T x, DOUBLE_T y)
40 {
41 return x / y;
42 }
43
44 DOUBLE_T a68_posinf_double (void)
45 {
46 return a68_div_double (+1.0, 0.0);
47 }
48
49 DOUBLE_T a68_neginf_double (void)
50 {
51 return a68_div_double (-1.0, 0.0);
52 }
53
54 //! @brief PROC (LONG REAL, LONG REAL) LONG REAL long hypot
55
56 DOUBLE_T a68_hypot_double (DOUBLE_T x, DOUBLE_T y)
57 {
58 // Sqrt (x^2 + y^2) that does not needlessly overflow.
59 DOUBLE_T xabs = ABSQ (x), yabs = ABSQ (y), min, max;
60 if (xabs < yabs) {
61 min = xabs;
62 max = yabs;
63 } else {
64 min = yabs;
65 max = xabs;
66 }
67 if (min == 0.0q) {
68 return max;
69 } else {
70 DOUBLE_T u = min / max;
71 return max * sqrt_double (1.0q + u * u);
72 }
73 }
74
75 //! @brief PROC (LONG REAL, LONG REAL, LONG REAL) LONG REAL long beta inc
76
77 DOUBLE_T a68_beta_inc_double (DOUBLE_T s, DOUBLE_T t, DOUBLE_T x)
78 {
79 // Incomplete beta function I{x}(s, t) by continued fraction,
80 // see dlmf.nist.gov/8.17; Lentz's algorithm.
81 if (x < 0.0q || x > 1.0q) {
82 errno = ERANGE;
83 return -1.0q;
84 } else {
85 const INT_T lim = 16 * sizeof (DOUBLE_T);
86 // Rapid convergence when x <= (s+1)/(s+t+2) or else recursion.
87 if (x > (s + 1.0q) / (s + t + 2.0q)) {
88 // B{x}(s, t) = 1 - B{1-x}(t, s)
89 return 1.0q - a68_beta_inc_double (s, t, 1.0q - x);
90 }
91 // Lentz's algorithm for continued fraction.
92 DOUBLE_T W = 1.0q, F = 1.0q, c = 1.0q, d = 0.0q;
93 BOOL_T cont = A68_TRUE;
94 for (int N = 0, m = 0; cont && N < lim; N++) {
95 DOUBLE_T T;
96 if (N == 0) {
97 T = 1.0q;
98 } else if (N % 2 == 0) {
99 // d{2m} := x m(t-m)/((s+2m-1)(s+2m))
100 T = x * m * (t - m) / (s + 2.0q * m - 1.0q) / (s + 2.0q * m);
101 } else {
102 // d{2m+1} := -x (s+m)(s+t+m)/((s+2m+1)(s+2m))
103 T = -x * (s + m) * (s + t + m) / (s + 2.0q * m + 1.0q) / (s + 2.0q * m);
104 m++;
105 }
106 d = 1.0q / (T * d + 1.0q);
107 c = T / c + 1.0q;
108 F *= c * d;
109 if (F == W) {
110 cont = A68_FALSE;
111 } else {
112 W = F;
113 }
114 }
115 // I{x}(s,t)=x^s(1-x)^t / s / B(s,t) F
116 DOUBLE_T beta = exp_double (lgamma_double (s) + lgamma_double (t) - lgamma_double (s + t));
117 return pow_double (x, s) * pow_double (1.0q - x, t) / s / beta * (F - 1.0q);
118 }
119 }
120
121 //! @brief PROC (LONG REAL) LONG REAL long csc
122
123 DOUBLE_T a68_csc_double (DOUBLE_T x)
124 {
125 DOUBLE_T z = sin_double (x);
126 A68_OVERFLOW (z == 0.0q);
127 return 1.0q / z;
128 }
129
130 //! @brief PROC (LONG REAL) LONG REAL long acsc
131
132 DOUBLE_T a68_acsc_double (DOUBLE_T x)
133 {
134 A68_OVERFLOW (x == 0.0q);
135 return asin_double (1.0q / x);
136 }
137
138 //! @brief PROC (LONG REAL) LONG REAL long sec
139
140 DOUBLE_T a68_sec_double (DOUBLE_T x)
141 {
142 DOUBLE_T z = cos_double (x);
143 A68_OVERFLOW (z == 0.0q);
144 return 1.0q / z;
145 }
146
147 //! @brief PROC (LONG REAL) LONG REAL long asec
148
149 DOUBLE_T a68_asec_double (DOUBLE_T x)
150 {
151 A68_OVERFLOW (x == 0.0q);
152 return acos_double (1.0q / x);
153 }
154
155 //! @brief PROC (LONG REAL) LONG REAL long cot
156
157 DOUBLE_T a68_cot_double (DOUBLE_T x)
158 {
159 DOUBLE_T z = sin_double (x);
160 A68_OVERFLOW (z == 0.0q);
161 return cos_double (x) / z;
162 }
163
164 //! @brief PROC (LONG REAL) LONG REAL long acot
165
166 DOUBLE_T a68_acot_double (DOUBLE_T x)
167 {
168 A68_OVERFLOW (x == 0.0q);
169 return atan_double (1 / x);
170 }
171
172 //! brief PROC (LONG REAL) LONG REAL long cas
173
174 DOUBLE_T a68_cas_double (DOUBLE_T x)
175 {
176 // Hartley kernel, which Hartley named cas (cosine and sine).
177 return cos_double (x) + sin_double (x);
178 }
179
180 //! brief PROC (LONG REAL) LONG REAL long sindg
181
182 DOUBLE_T a68_sindg_double (DOUBLE_T x)
183 {
184 return sin_double (x * CONST_PI_OVER_180_Q);
185 }
186
187 //! brief PROC (LONG REAL) LONG REAL long cosdg
188
189 DOUBLE_T a68_cosdg_double (DOUBLE_T x)
190 {
191 return cos_double (x * CONST_PI_OVER_180_Q);
192 }
193
194 //! brief PROC (LONG REAL) LONG REAL long tandg
195
196 DOUBLE_T a68_tandg_double (DOUBLE_T x)
197 {
198 return tan_double (x * CONST_PI_OVER_180_Q);
199 }
200
201 //! brief PROC (LONG REAL) LONG REAL long asindg
202
203 DOUBLE_T a68_asindg_double (DOUBLE_T x)
204 {
205 return asin_double (x) * CONST_180_OVER_PI_Q;
206 }
207
208 //! brief PROC (LONG REAL) LONG REAL long acosdg
209
210 DOUBLE_T a68_acosdg_double (DOUBLE_T x)
211 {
212 return acos_double (x) * CONST_180_OVER_PI_Q;
213 }
214
215 //! brief PROC (LONG REAL) LONG REAL long atandg
216
217 DOUBLE_T a68_atandg_double (DOUBLE_T x)
218 {
219 return atan_double (x) * CONST_180_OVER_PI_Q;
220 }
221 //! @brief PROC (LONG REAL) LONG REAL long cscdg
222
223 DOUBLE_T a68_cscdg_double (DOUBLE_T x)
224 {
225 DOUBLE_T z = a68_sindg_double (x);
226 A68_OVERFLOW (z == 0.0q);
227 return 1.0q / z;
228 }
229
230 //! @brief PROC (LONG REAL) LONG REAL long acscdg
231
232 DOUBLE_T a68_acscdg_double (DOUBLE_T x)
233 {
234 A68_OVERFLOW (x == 0.0q);
235 return a68_asindg_double (1.0q / x);
236 }
237
238 //! @brief PROC (LONG REAL) LONG REAL long secdg
239
240 DOUBLE_T a68_secdg_double (DOUBLE_T x)
241 {
242 DOUBLE_T z = a68_cosdg_double (x);
243 A68_OVERFLOW (z == 0.0q);
244 return 1.0q / z;
245 }
246
247 //! @brief PROC (LONG REAL) LONG REAL long asecdg
248
249 DOUBLE_T a68_asecdg_double (DOUBLE_T x)
250 {
251 A68_OVERFLOW (x == 0.0q);
252 return a68_acosdg_double (1.0q / x);
253 }
254
255 // PROC (LONG REAL) LONG REAL long cotdg
256
257 DOUBLE_T a68_cotdg_double (DOUBLE_T x)
258 {
259 DOUBLE_T z = a68_sindg_double (x);
260 A68_OVERFLOW (z == 0);
261 return a68_cosdg_double (x) / z;
262 }
263
264 // PROC (LONG REAL) LONG REAL long acotdg
265
266 DOUBLE_T a68_acotdg_double (DOUBLE_T z)
267 {
268 A68_OVERFLOW (z == 0);
269 return a68_atandg_double (1 / z);
270 }
271
272 // @brief PROC (LONG REAL) LONG REAL long sinpi
273
274 DOUBLE_T a68_sinpi_double (DOUBLE_T x)
275 {
276 x = fmod_double (x, 2.0q);
277 if (x <= -1.0q) {
278 x += 2.0q;
279 } else if (x > 1.0q) {
280 x -= 2.0q;
281 }
282 // x in <-1, 1].
283 if (x == 0.0q || x == 1.0q) {
284 return 0.0q;
285 } else if (x == 0.5q) {
286 return 1.0q;
287 }
288 if (x == -0.5q) {
289 return -1.0q;
290 } else {
291 return sin_double (CONST_PI_Q * x);
292 }
293 }
294
295 // @brief PROC (LONG REAL) LONG REAL long cospi
296
297 DOUBLE_T a68_cospi_double (DOUBLE_T x)
298 {
299 x = fmod_double (fabs_double (x), 2.0q);
300 // x in [0, 2>.
301 if (x == 0.5q || x == 1.5q) {
302 return 0.0q;
303 } else if (x == 0.0q) {
304 return 1.0q;
305 } else if (x == 1.0q) {
306 return -1.0q;
307 } else {
308 return cos_double (CONST_PI_Q * x);
309 }
310 }
311
312 // @brief PROC (LONG REAL) LONG REAL long tanpi
313
314 DOUBLE_T a68_tanpi_double (DOUBLE_T x)
315 {
316 x = fmod_double (x, 1.0q);
317 if (x <= -0.5q) {
318 x += 1.0q;
319 } else if (x > 0.5q) {
320 x -= 1.0q;
321 }
322 // x in <-1/2, 1/2].
323 A68_OVERFLOW (x == 0.5q);
324 if (x == -0.25q) {
325 return -1.0q;
326 } else if (x == 0) {
327 return 0.0q;
328 } else if (x == 0.25q) {
329 return 1.0q;
330 } else {
331 return a68_sinpi_double (x) / a68_cospi_double (x);
332 }
333 }
334
335 // @brief PROC (LONG REAL) LONG REAL long cotpi
336
337 DOUBLE_T a68_cotpi_double (DOUBLE_T x)
338 {
339 x = fmod_double (x, 1.0q);
340 if (x <= -0.5q) {
341 x += 1.0q;
342 } else if (x > 0.5q) {
343 x -= 1.0q;
344 }
345 // x in <-1/2, 1/2].
346 A68_OVERFLOW (x == 0.0q);
347 if (x == -0.25q) {
348 return -1.0q;
349 } else if (x == 0.25q) {
350 return 1.0q;
351 } else if (x == 0.5q) {
352 return 0.0q;
353 } else {
354 return a68_cospi_double (x) / a68_sinpi_double (x);
355 }
356 }
357
358 #endif
