mp-math.c
1 //! @file mp-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 [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] LONG REAL math functions.
25
26 #include "a68g.h"
27
28 #include "a68g-double.h"
29 #include "a68g-mp.h"
30 #include "a68g-prelude.h"
31
32 //! @brief Test on |"z"| > 0.001 for argument reduction in "sin" and "exp".
33
34 static inline BOOL_T eps_mp (MP_T * z, int digs)
35 {
36 if (MP_DIGIT (z, 1) == 0) {
37 return A68_FALSE;
38 } else if (MP_EXPONENT (z) > -1) {
39 return A68_TRUE;
40 } else if (MP_EXPONENT (z) < -1) {
41 return A68_FALSE;
42 } else {
43 // More or less optimised for LONG and default LONG LONG precisions.
44 return (BOOL_T) (digs <= 10 ? ABS (MP_DIGIT (z, 1)) > 0.01 * MP_RADIX : ABS (MP_DIGIT (z, 1)) > 0.001 * MP_RADIX);
45 }
46 }
47
48 //! @brief PROC (LONG REAL) LONG REAL sqrt
49
50 MP_T *sqrt_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
51 {
52 ADDR_T pop_sp = A68_SP;
53 if (MP_DIGIT (x, 1) == 0) {
54 A68_SP = pop_sp;
55 SET_MP_ZERO (z, digs);
56 return z;
57 }
58 if (MP_DIGIT (x, 1) < 0) {
59 A68_SP = pop_sp;
60 errno = EDOM;
61 return NaN_MP;
62 }
63 int gdigs = FUN_DIGITS (digs), hdigs;
64 BOOL_T reciprocal = A68_FALSE;
65 MP_T *z_g = nil_mp (p, gdigs), *x_g = len_mp (p, x, digs, gdigs);
66 MP_T *tmp = nil_mp (p, gdigs);
67 // Scaling for small x; sqrt (x) = 1 / sqrt (1 / x).
68 if ((reciprocal = (BOOL_T) (MP_EXPONENT (x_g) < 0)) == A68_TRUE) {
69 (void) rec_mp (p, x_g, x_g, gdigs);
70 }
71 if (ABS (MP_EXPONENT (x_g)) >= 2) {
72 // For extreme arguments we want accurate results as well.
73 int expo = (int) MP_EXPONENT (x_g);
74 MP_EXPONENT (x_g) = (MP_T) (expo % 2);
75 (void) sqrt_mp (p, z_g, x_g, gdigs);
76 MP_EXPONENT (z_g) += (MP_T) (expo / 2);
77 } else {
78 // Argument is in range. Estimate the root as double.
79 #if (A68_LEVEL >= 3)
80 DOUBLE_T x_d = mp_to_double (p, x_g, gdigs);
81 (void) double_to_mp (p, z_g, sqrt_double (x_d), gdigs);
82 #else
83 REAL_T x_d = mp_to_real (p, x_g, gdigs);
84 (void) real_to_mp (p, z_g, sqrt (x_d), gdigs);
85 #endif
86 // Newton's method: x<n+1> = (x<n> + a / x<n>) / 2.
87 int decimals = DOUBLE_ACCURACY;
88 do {
89 decimals <<= 1;
90 hdigs = MINIMUM (1 + decimals / LOG_MP_RADIX, gdigs);
91 (void) div_mp (p, tmp, x_g, z_g, hdigs);
92 (void) add_mp (p, tmp, z_g, tmp, hdigs);
93 (void) half_mp (p, z_g, tmp, hdigs);
94 } while (decimals < 2 * gdigs * LOG_MP_RADIX);
95 }
96 if (reciprocal) {
97 (void) rec_mp (p, z_g, z_g, digs);
98 }
99 (void) shorten_mp (p, z, digs, z_g, gdigs);
100 // Exit.
101 A68_SP = pop_sp;
102 return z;
103 }
104
105 //! @brief PROC (LONG REAL) LONG REAL curt
106
107 MP_T *curt_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
108 {
109 ADDR_T pop_sp = A68_SP;
110 if (MP_DIGIT (x, 1) == 0) {
111 A68_SP = pop_sp;
112 SET_MP_ZERO (z, digs);
113 return z;
114 }
115 BOOL_T change_sign = A68_FALSE;
116 if (MP_DIGIT (x, 1) < 0) {
117 change_sign = A68_TRUE;
118 MP_DIGIT (x, 1) = -MP_DIGIT (x, 1);
119 }
120 int gdigs = FUN_DIGITS (digs), hdigs;
121 BOOL_T reciprocal = A68_FALSE;
122 MP_T *z_g = nil_mp (p, gdigs), *x_g = len_mp (p, x, digs, gdigs);
123 MP_T *tmp = nil_mp (p, gdigs);
124 // Scaling for small x; curt (x) = 1 / curt (1 / x).
125 if ((reciprocal = (BOOL_T) (MP_EXPONENT (x_g) < 0)) == A68_TRUE) {
126 (void) rec_mp (p, x_g, x_g, gdigs);
127 }
128 if (ABS (MP_EXPONENT (x_g)) >= 3) {
129 // For extreme arguments we want accurate results as well.
130 int expo = (int) MP_EXPONENT (x_g);
131 MP_EXPONENT (x_g) = (MP_T) (expo % 3);
132 (void) curt_mp (p, z_g, x_g, gdigs);
133 MP_EXPONENT (z_g) += (MP_T) (expo / 3);
134 } else {
135 // Argument is in range. Estimate the root as double.
136 #if (A68_LEVEL >= 3)
137 DOUBLE_T x_d = mp_to_double (p, x_g, gdigs);
138 (void) double_to_mp (p, z_g, cbrt_double (x_d), gdigs);
139 #else
140 REAL_T x_d = mp_to_real (p, x_g, gdigs);
141 (void) real_to_mp (p, z_g, cbrt (x_d), gdigs);
142 #endif
143 // Newton's method: x<n+1> = (2 x<n> + a / x<n> ^ 2) / 3.
144 int decimals = DOUBLE_ACCURACY;
145 do {
146 decimals <<= 1;
147 hdigs = MINIMUM (1 + decimals / LOG_MP_RADIX, gdigs);
148 (void) mul_mp (p, tmp, z_g, z_g, hdigs);
149 (void) div_mp (p, tmp, x_g, tmp, hdigs);
150 (void) add_mp (p, tmp, z_g, tmp, hdigs);
151 (void) add_mp (p, tmp, z_g, tmp, hdigs);
152 (void) div_mp_digit (p, z_g, tmp, (MP_T) 3, hdigs);
153 } while (decimals < gdigs * LOG_MP_RADIX);
154 }
155 if (reciprocal) {
156 (void) rec_mp (p, z_g, z_g, digs);
157 }
158 (void) shorten_mp (p, z, digs, z_g, gdigs);
159 // Exit.
160 A68_SP = pop_sp;
161 if (change_sign) {
162 MP_DIGIT (z, 1) = -MP_DIGIT (z, 1);
163 }
164 return z;
165 }
166
167 //! @brief PROC (LONG REAL) LONG REAL hypot
168
169 MP_T *hypot_mp (NODE_T * p, MP_T * z, MP_T * x, MP_T * y, int digs)
170 {
171 // sqrt (x^2 + y^2).
172 ADDR_T pop_sp = A68_SP;
173 MP_T *t = nil_mp (p, digs), *u = nil_mp (p, digs), *v = nil_mp (p, digs);
174 (void) move_mp (u, x, digs);
175 (void) move_mp (v, y, digs);
176 MP_DIGIT (u, 1) = ABS (MP_DIGIT (u, 1));
177 MP_DIGIT (v, 1) = ABS (MP_DIGIT (v, 1));
178 if (IS_ZERO_MP (u)) {
179 (void) move_mp (z, v, digs);
180 } else if (IS_ZERO_MP (v)) {
181 (void) move_mp (z, u, digs);
182 } else {
183 SET_MP_ONE (t, digs);
184 (void) sub_mp (p, z, u, v, digs);
185 if (MP_DIGIT (z, 1) > 0) {
186 (void) div_mp (p, z, v, u, digs);
187 (void) mul_mp (p, z, z, z, digs);
188 (void) add_mp (p, z, t, z, digs);
189 (void) sqrt_mp (p, z, z, digs);
190 (void) mul_mp (p, z, u, z, digs);
191 } else {
192 (void) div_mp (p, z, u, v, digs);
193 (void) mul_mp (p, z, z, z, digs);
194 (void) add_mp (p, z, t, z, digs);
195 (void) sqrt_mp (p, z, z, digs);
196 (void) mul_mp (p, z, v, z, digs);
197 }
198 }
199 A68_SP = pop_sp;
200 return z;
201 }
202
203 //! @brief PROC (LONG REAL) LONG REAL exp
204
205 MP_T *exp_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
206 {
207 // Argument is reduced by using exp (z / (2 ** n)) ** (2 ** n) = exp(z).
208 int gdigs = FUN_DIGITS (digs);
209 ADDR_T pop_sp = A68_SP;
210 if (MP_DIGIT (x, 1) == 0) {
211 SET_MP_ONE (z, digs);
212 return z;
213 }
214 MP_T *x_g = len_mp (p, x, digs, gdigs);
215 MP_T *pwr = nil_mp (p, gdigs), *fac = nil_mp (p, gdigs), *sum = nil_mp (p, gdigs);
216 MP_T *tmp = nil_mp (p, gdigs);
217 int m = 0;
218 // Scale x down.
219 while (eps_mp (x_g, gdigs)) {
220 m++;
221 (void) half_mp (p, x_g, x_g, gdigs);
222 }
223 // Calculate Taylor sum exp (z) = 1 + z / 1 ! + z ** 2 / 2 ! + ..
224 SET_MP_ONE (sum, gdigs);
225 (void) add_mp (p, sum, sum, x_g, gdigs);
226 (void) mul_mp (p, pwr, x_g, x_g, gdigs);
227 (void) half_mp (p, tmp, pwr, gdigs);
228 (void) add_mp (p, sum, sum, tmp, gdigs);
229 (void) mul_mp (p, pwr, pwr, x_g, gdigs);
230 (void) div_mp_digit (p, tmp, pwr, (MP_T) 6, gdigs);
231 (void) add_mp (p, sum, sum, tmp, gdigs);
232 (void) mul_mp (p, pwr, pwr, x_g, gdigs);
233 (void) div_mp_digit (p, tmp, pwr, (MP_T) 24, gdigs);
234 (void) add_mp (p, sum, sum, tmp, gdigs);
235 (void) mul_mp (p, pwr, pwr, x_g, gdigs);
236 (void) div_mp_digit (p, tmp, pwr, (MP_T) 120, gdigs);
237 (void) add_mp (p, sum, sum, tmp, gdigs);
238 (void) mul_mp (p, pwr, pwr, x_g, gdigs);
239 (void) div_mp_digit (p, tmp, pwr, (MP_T) 720, gdigs);
240 (void) add_mp (p, sum, sum, tmp, gdigs);
241 (void) mul_mp (p, pwr, pwr, x_g, gdigs);
242 (void) div_mp_digit (p, tmp, pwr, (MP_T) 5040, gdigs);
243 (void) add_mp (p, sum, sum, tmp, gdigs);
244 (void) mul_mp (p, pwr, pwr, x_g, gdigs);
245 (void) div_mp_digit (p, tmp, pwr, (MP_T) 40320, gdigs);
246 (void) add_mp (p, sum, sum, tmp, gdigs);
247 (void) mul_mp (p, pwr, pwr, x_g, gdigs);
248 (void) div_mp_digit (p, tmp, pwr, (MP_T) 362880, gdigs);
249 (void) add_mp (p, sum, sum, tmp, gdigs);
250 (void) mul_mp (p, pwr, pwr, x_g, gdigs);
251 (void) set_mp (fac, (MP_T) (MP_T) 3628800, 0, gdigs);
252 int n = 10;
253 BOOL_T iter = (BOOL_T) (MP_DIGIT (pwr, 1) != 0);
254 while (iter) {
255 (void) div_mp (p, tmp, pwr, fac, gdigs);
256 if (MP_EXPONENT (tmp) <= (MP_EXPONENT (sum) - gdigs)) {
257 iter = A68_FALSE;
258 } else {
259 (void) add_mp (p, sum, sum, tmp, gdigs);
260 (void) mul_mp (p, pwr, pwr, x_g, gdigs);
261 n++;
262 (void) mul_mp_digit (p, fac, fac, (MP_T) n, gdigs);
263 }
264 }
265 // Square exp (x) up.
266 while (m--) {
267 (void) mul_mp (p, sum, sum, sum, gdigs);
268 }
269 (void) shorten_mp (p, z, digs, sum, gdigs);
270 // Exit.
271 A68_SP = pop_sp;
272 return z;
273 }
274
275 //! @brief PROC (LONG REAL) LONG REAL exp (x) - 1
276 // assuming "x" to be close to 0.
277
278 MP_T *expm1_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
279 {
280 int gdigs = FUN_DIGITS (digs);
281 ADDR_T pop_sp = A68_SP;
282 if (MP_DIGIT (x, 1) == 0) {
283 SET_MP_ONE (z, digs);
284 return z;
285 }
286 MP_T *x_g = len_mp (p, x, digs, gdigs);
287 MP_T *pwr = nil_mp (p, gdigs), *fac = nil_mp (p, gdigs), *sum = nil_mp (p, gdigs);
288 MP_T *tmp = nil_mp (p, gdigs);
289 // Calculate Taylor sum expm1 (z) = z / 1 ! + z ** 2 / 2 ! + ...
290 (void) add_mp (p, sum, sum, x_g, gdigs);
291 (void) mul_mp (p, pwr, x_g, x_g, gdigs);
292 (void) half_mp (p, tmp, pwr, gdigs);
293 (void) add_mp (p, sum, sum, tmp, gdigs);
294 (void) mul_mp (p, pwr, pwr, x_g, gdigs);
295 (void) div_mp_digit (p, tmp, pwr, (MP_T) 6, gdigs);
296 (void) add_mp (p, sum, sum, tmp, gdigs);
297 (void) mul_mp (p, pwr, pwr, x_g, gdigs);
298 (void) div_mp_digit (p, tmp, pwr, (MP_T) 24, gdigs);
299 (void) add_mp (p, sum, sum, tmp, gdigs);
300 (void) mul_mp (p, pwr, pwr, x_g, gdigs);
301 (void) div_mp_digit (p, tmp, pwr, (MP_T) 120, gdigs);
302 (void) add_mp (p, sum, sum, tmp, gdigs);
303 (void) mul_mp (p, pwr, pwr, x_g, gdigs);
304 (void) div_mp_digit (p, tmp, pwr, (MP_T) 720, gdigs);
305 (void) add_mp (p, sum, sum, tmp, gdigs);
306 (void) mul_mp (p, pwr, pwr, x_g, gdigs);
307 (void) div_mp_digit (p, tmp, pwr, (MP_T) 5040, gdigs);
308 (void) add_mp (p, sum, sum, tmp, gdigs);
309 (void) mul_mp (p, pwr, pwr, x_g, gdigs);
310 (void) div_mp_digit (p, tmp, pwr, (MP_T) 40320, gdigs);
311 (void) add_mp (p, sum, sum, tmp, gdigs);
312 (void) mul_mp (p, pwr, pwr, x_g, gdigs);
313 (void) div_mp_digit (p, tmp, pwr, (MP_T) 362880, gdigs);
314 (void) add_mp (p, sum, sum, tmp, gdigs);
315 (void) mul_mp (p, pwr, pwr, x_g, gdigs);
316 (void) set_mp (fac, (MP_T) (MP_T) 3628800, 0, gdigs);
317 int n = 10;
318 BOOL_T iter = (BOOL_T) (MP_DIGIT (pwr, 1) != 0);
319 while (iter) {
320 (void) div_mp (p, tmp, pwr, fac, gdigs);
321 if (MP_EXPONENT (tmp) <= (MP_EXPONENT (sum) - gdigs)) {
322 iter = A68_FALSE;
323 } else {
324 (void) add_mp (p, sum, sum, tmp, gdigs);
325 (void) mul_mp (p, pwr, pwr, x_g, gdigs);
326 n++;
327 (void) mul_mp_digit (p, fac, fac, (MP_T) n, gdigs);
328 }
329 }
330 (void) shorten_mp (p, z, digs, sum, gdigs);
331 // Exit.
332 A68_SP = pop_sp;
333 return z;
334 }
335
336 //! @brief ln scale
337
338 MP_T *mp_ln_scale (NODE_T * p, MP_T * z, int digs)
339 {
340 ADDR_T pop_sp = A68_SP;
341 int gdigs = FUN_DIGITS (digs);
342 MP_T *z_g = nil_mp (p, gdigs);
343 // First see if we can restore a previous calculation.
344 if (gdigs <= A68_MP (mp_ln_scale_size)) {
345 (void) move_mp (z_g, A68_MP (mp_ln_scale), gdigs);
346 } else {
347 // No luck with the kept value, we generate a longer one.
348 (void) set_mp (z_g, (MP_T) 1, 1, gdigs);
349 (void) ln_mp (p, z_g, z_g, gdigs);
350 A68_MP (mp_ln_scale) = (MP_T *) get_heap_space ((unt) SIZE_MP (gdigs));
351 (void) move_mp (A68_MP (mp_ln_scale), z_g, gdigs);
352 A68_MP (mp_ln_scale_size) = gdigs;
353 }
354 (void) shorten_mp (p, z, digs, z_g, gdigs);
355 A68_SP = pop_sp;
356 return z;
357 }
358
359 //! @brief ln 10
360
361 MP_T *mp_ln_10 (NODE_T * p, MP_T * z, int digs)
362 {
363 ADDR_T pop_sp = A68_SP;
364 int gdigs = FUN_DIGITS (digs);
365 MP_T *z_g = nil_mp (p, gdigs);
366 // First see if we can restore a previous calculation.
367 if (gdigs <= A68_MP (mp_ln_10_size)) {
368 (void) move_mp (z_g, A68_MP (mp_ln_10), gdigs);
369 } else {
370 // No luck with the kept value, we generate a longer one.
371 (void) set_mp (z_g, (MP_T) 10, 0, gdigs);
372 (void) ln_mp (p, z_g, z_g, gdigs);
373 A68_MP (mp_ln_10) = (MP_T *) get_heap_space ((unt) SIZE_MP (gdigs));
374 (void) move_mp (A68_MP (mp_ln_10), z_g, gdigs);
375 A68_MP (mp_ln_10_size) = gdigs;
376 }
377 (void) shorten_mp (p, z, digs, z_g, gdigs);
378 A68_SP = pop_sp;
379 return z;
380 }
381
382 //! @brief PROC (LONG REAL) LONG REAL ln
383
384 MP_T *ln_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
385 {
386 // Depending on the argument we choose either Taylor or Newton.
387 ADDR_T pop_sp = A68_SP;
388 int gdigs = FUN_DIGITS (digs);
389 BOOL_T neg, scale;
390 MP_T expo = 0;
391 if (MP_DIGIT (x, 1) <= 0) {
392 errno = EDOM;
393 return NaN_MP;
394 }
395 MP_T *x_g = len_mp (p, x, digs, gdigs), *z_g = nil_mp (p, gdigs);
396 // We use ln (1 / x) = - ln (x).
397 neg = (BOOL_T) (MP_EXPONENT (x_g) < 0);
398 if (neg) {
399 (void) rec_mp (p, x_g, x_g, digs);
400 }
401 // We want correct results for extreme arguments. We scale when "x_g" exceeds
402 // "MP_RADIX ** +- 2", using ln (x * MP_RADIX ** n) = ln (x) + n * ln (MP_RADIX).
403 scale = (BOOL_T) (ABS (MP_EXPONENT (x_g)) >= 2);
404 if (scale) {
405 expo = MP_EXPONENT (x_g);
406 MP_EXPONENT (x_g) = (MP_T) 0;
407 }
408 if (MP_EXPONENT (x_g) == 0 && MP_DIGIT (x_g, 1) == 1 && MP_DIGIT (x_g, 2) == 0) {
409 // Taylor sum for x close to unity.
410 // ln (x) = (x - 1) - (x - 1) ** 2 / 2 + (x - 1) ** 3 / 3 - ...
411 // This is faster for small x and avoids cancellation.
412 MP_T *pwr = nil_mp (p, gdigs), *tmp = nil_mp (p, gdigs);
413 (void) minus_one_mp (p, x_g, x_g, gdigs);
414 (void) mul_mp (p, pwr, x_g, x_g, gdigs);
415 (void) move_mp (z_g, x_g, gdigs);
416 int n = 2;
417 BOOL_T iter = (BOOL_T) (MP_DIGIT (pwr, 1) != 0);
418 while (iter) {
419 (void) div_mp_digit (p, tmp, pwr, (MP_T) n, gdigs);
420 if (MP_EXPONENT (tmp) <= (MP_EXPONENT (z_g) - gdigs)) {
421 iter = A68_FALSE;
422 } else {
423 MP_DIGIT (tmp, 1) = (EVEN (n) ? -MP_DIGIT (tmp, 1) : MP_DIGIT (tmp, 1));
424 (void) add_mp (p, z_g, z_g, tmp, gdigs);
425 (void) mul_mp (p, pwr, pwr, x_g, gdigs);
426 n++;
427 }
428 }
429 } else {
430 // Newton's method: x<n+1> = x<n> - 1 + a / exp (x<n>).
431 MP_T *tmp = nil_mp (p, gdigs);
432 // Construct an estimate.
433 #if (A68_LEVEL >= 3)
434 (void) double_to_mp (p, z_g, log_double (mp_to_double (p, x_g, gdigs)), gdigs);
435 #else
436 (void) real_to_mp (p, z_g, log (mp_to_real (p, x_g, gdigs)), gdigs);
437 #endif
438 int decimals = DOUBLE_ACCURACY;
439 do {
440 decimals <<= 1;
441 int hdigs = MINIMUM (1 + decimals / LOG_MP_RADIX, gdigs);
442 (void) exp_mp (p, tmp, z_g, hdigs);
443 (void) div_mp (p, tmp, x_g, tmp, hdigs);
444 (void) minus_one_mp (p, z_g, z_g, hdigs);
445 (void) add_mp (p, z_g, z_g, tmp, hdigs);
446 } while (decimals < gdigs * LOG_MP_RADIX);
447 }
448 // Inverse scaling.
449 if (scale) {
450 // ln (x * MP_RADIX ** n) = ln (x) + n * ln (MP_RADIX).
451 MP_T *ln_base = nil_mp (p, gdigs);
452 (void) mp_ln_scale (p, ln_base, gdigs);
453 (void) mul_mp_digit (p, ln_base, ln_base, (MP_T) expo, gdigs);
454 (void) add_mp (p, z_g, z_g, ln_base, gdigs);
455 }
456 if (neg) {
457 MP_DIGIT (z_g, 1) = -MP_DIGIT (z_g, 1);
458 }
459 (void) shorten_mp (p, z, digs, z_g, gdigs);
460 // Exit.
461 A68_SP = pop_sp;
462 return z;
463 }
464
465 //! @brief PROC (LONG REAL) LONG REAL log
466
467 MP_T *log_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
468 {
469 ADDR_T pop_sp = A68_SP;
470 MP_T *ln_10 = nil_mp (p, digs);
471 if (ln_mp (p, z, x, digs) == NaN_MP) {
472 errno = EDOM;
473 return NaN_MP;
474 }
475 (void) mp_ln_10 (p, ln_10, digs);
476 (void) div_mp (p, z, z, ln_10, digs);
477 A68_SP = pop_sp;
478 return z;
479 }
480
481 //! @brief sinh ("z") and cosh ("z")
482
483 MP_T *hyp_mp (NODE_T * p, MP_T * sh, MP_T * ch, MP_T * z, int digs)
484 {
485 ADDR_T pop_sp = A68_SP;
486 MP_T *x_g = nil_mp (p, digs), *y_g = nil_mp (p, digs), *z_g = nil_mp (p, digs);
487 (void) move_mp (z_g, z, digs);
488 (void) exp_mp (p, x_g, z_g, digs);
489 (void) rec_mp (p, y_g, x_g, digs);
490 (void) add_mp (p, ch, x_g, y_g, digs);
491 // Avoid cancellation for sinh.
492 if ((MP_DIGIT (x_g, 1) == 1 && MP_DIGIT (x_g, 2) == 0) || (MP_DIGIT (y_g, 1) == 1 && MP_DIGIT (y_g, 2) == 0)) {
493 (void) expm1_mp (p, x_g, z_g, digs);
494 MP_DIGIT (z_g, 1) = -MP_DIGIT (z_g, 1);
495 (void) expm1_mp (p, y_g, z_g, digs);
496 }
497 (void) sub_mp (p, sh, x_g, y_g, digs);
498 (void) half_mp (p, sh, sh, digs);
499 (void) half_mp (p, ch, ch, digs);
500 A68_SP = pop_sp;
501 return z;
502 }
503
504 //! @brief PROC (LONG REAL) LONG REAL sinh
505
506 MP_T *sinh_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
507 {
508 ADDR_T pop_sp = A68_SP;
509 int gdigs = FUN_DIGITS (digs);
510 MP_T *x_g = len_mp (p, x, digs, gdigs), *y_g = nil_mp (p, gdigs), *z_g = nil_mp (p, gdigs);
511 (void) hyp_mp (p, z_g, y_g, x_g, gdigs);
512 (void) shorten_mp (p, z, digs, z_g, gdigs);
513 A68_SP = pop_sp;
514 return z;
515 }
516
517 //! @brief PROC (LONG REAL) LONG REAL asinh
518
519 MP_T *asinh_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
520 {
521 if (IS_ZERO_MP (x)) {
522 SET_MP_ZERO (z, digs);
523 return z;
524 } else {
525 ADDR_T pop_sp = A68_SP;
526 int gdigs;
527 if (MP_EXPONENT (x) >= -1) {
528 gdigs = FUN_DIGITS (digs);
529 } else {
530 // Extra precision when x^2+1 gets close to 1.
531 gdigs = 2 * FUN_DIGITS (digs);
532 }
533 MP_T *x_g = len_mp (p, x, digs, gdigs), *y_g = nil_mp (p, gdigs), *z_g = nil_mp (p, gdigs);
534 (void) mul_mp (p, z_g, x_g, x_g, gdigs);
535 SET_MP_ONE (y_g, gdigs);
536 (void) add_mp (p, y_g, z_g, y_g, gdigs);
537 (void) sqrt_mp (p, y_g, y_g, gdigs);
538 (void) add_mp (p, y_g, y_g, x_g, gdigs);
539 (void) ln_mp (p, z_g, y_g, gdigs);
540 if (IS_ZERO_MP (z_g)) {
541 (void) move_mp (z, x, digs);
542 } else {
543 (void) shorten_mp (p, z, digs, z_g, gdigs);
544 }
545 A68_SP = pop_sp;
546 return z;
547 }
548 }
549
550 //! @brief PROC (LONG REAL) LONG REAL cosh
551
552 MP_T *cosh_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
553 {
554 ADDR_T pop_sp = A68_SP;
555 int gdigs = FUN_DIGITS (digs);
556 MP_T *x_g = len_mp (p, x, digs, gdigs), *y_g = nil_mp (p, gdigs), *z_g = nil_mp (p, gdigs);
557 (void) hyp_mp (p, y_g, z_g, x_g, gdigs);
558 (void) shorten_mp (p, z, digs, z_g, gdigs);
559 A68_SP = pop_sp;
560 return z;
561 }
562
563 //! @brief PROC (LONG REAL) LONG REAL acosh
564
565 MP_T *acosh_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
566 {
567 ADDR_T pop_sp = A68_SP;
568 int gdigs;
569 if (MP_DIGIT (x, 1) == 1 && MP_DIGIT (x, 2) == 0) {
570 // Extra precision when x^2-1 gets close to 0.
571 gdigs = 2 * FUN_DIGITS (digs);
572 } else {
573 gdigs = FUN_DIGITS (digs);
574 }
575 MP_T *x_g = len_mp (p, x, digs, gdigs), *y_g = nil_mp (p, gdigs), *z_g = nil_mp (p, gdigs);
576 (void) mul_mp (p, z_g, x_g, x_g, gdigs);
577 SET_MP_ONE (y_g, gdigs);
578 (void) sub_mp (p, y_g, z_g, y_g, gdigs);
579 (void) sqrt_mp (p, y_g, y_g, gdigs);
580 (void) add_mp (p, y_g, y_g, x_g, gdigs);
581 (void) ln_mp (p, z_g, y_g, gdigs);
582 (void) shorten_mp (p, z, digs, z_g, gdigs);
583 A68_SP = pop_sp;
584 return z;
585 }
586
587 //! @brief PROC (LONG REAL) LONG REAL tanh
588
589 MP_T *tanh_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
590 {
591 ADDR_T pop_sp = A68_SP;
592 int gdigs = FUN_DIGITS (digs);
593 MP_T *x_g = len_mp (p, x, digs, gdigs), *y_g = nil_mp (p, gdigs), *z_g = nil_mp (p, gdigs);
594 (void) hyp_mp (p, y_g, z_g, x_g, gdigs);
595 (void) div_mp (p, z_g, y_g, z_g, gdigs);
596 (void) shorten_mp (p, z, digs, z_g, gdigs);
597 A68_SP = pop_sp;
598 return z;
599 }
600
601 //! @brief PROC (LONG REAL) LONG REAL atanh
602
603 MP_T *atanh_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
604 {
605 ADDR_T pop_sp = A68_SP;
606 int gdigs = FUN_DIGITS (digs);
607 MP_T *x_g = len_mp (p, x, digs, gdigs), *y_g = nil_mp (p, gdigs), *z_g = nil_mp (p, gdigs);
608 SET_MP_ONE (y_g, gdigs);
609 (void) add_mp (p, z_g, y_g, x_g, gdigs);
610 (void) sub_mp (p, y_g, y_g, x_g, gdigs);
611 (void) div_mp (p, y_g, z_g, y_g, gdigs);
612 (void) ln_mp (p, z_g, y_g, gdigs);
613 (void) half_mp (p, z_g, z_g, gdigs);
614 (void) shorten_mp (p, z, digs, z_g, gdigs);
615 A68_SP = pop_sp;
616 return z;
617 }
618
619 //! @brief PROC (LONG REAL) LONG REAL sin
620
621 MP_T *sin_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
622 {
623 // Use triple-angle relation to reduce argument.
624 ADDR_T pop_sp = A68_SP;
625 int gdigs = FUN_DIGITS (digs);
626 // We will use "pi".
627 MP_T *pi = nil_mp (p, gdigs), *tpi = nil_mp (p, gdigs), *hpi = nil_mp (p, gdigs);
628 (void) mp_pi (p, pi, MP_PI, gdigs);
629 (void) mp_pi (p, tpi, MP_TWO_PI, gdigs);
630 (void) mp_pi (p, hpi, MP_HALF_PI, gdigs);
631 // Argument reduction (1): sin (x) = sin (x mod 2 pi).
632 MP_T *x_g = len_mp (p, x, digs, gdigs);
633 (void) mod_mp (p, x_g, x_g, tpi, gdigs);
634 // Argument reduction (2): sin (-x) = sin (x)
635 // sin (x) = - sin (x - pi); pi < x <= 2 pi
636 // sin (x) = sin (pi - x); pi / 2 < x <= pi
637 BOOL_T neg = (BOOL_T) (MP_DIGIT (x_g, 1) < 0);
638 if (neg) {
639 MP_DIGIT (x_g, 1) = -MP_DIGIT (x_g, 1);
640 }
641 MP_T *tmp = nil_mp (p, gdigs);
642 (void) sub_mp (p, tmp, x_g, pi, gdigs);
643 BOOL_T flip = (BOOL_T) (MP_DIGIT (tmp, 1) > 0);
644 if (flip) { // x > pi
645 (void) sub_mp (p, x_g, x_g, pi, gdigs);
646 }
647 (void) sub_mp (p, tmp, x_g, hpi, gdigs);
648 if (MP_DIGIT (tmp, 1) > 0) { // x > pi / 2
649 (void) sub_mp (p, x_g, pi, x_g, gdigs);
650 }
651 // Argument reduction (3): (follows from De Moivre's theorem)
652 // sin (3x) = sin (x) * (3 - 4 sin ^ 2 (x))
653 int m = 0;
654 while (eps_mp (x_g, gdigs)) {
655 m++;
656 (void) div_mp_digit (p, x_g, x_g, (MP_T) 3, gdigs);
657 }
658 // Taylor sum.
659 MP_T *sqr = nil_mp (p, gdigs), *z_g = nil_mp (p, gdigs);
660 MP_T *pwr = nil_mp (p, gdigs), *fac = nil_mp (p, gdigs);
661 (void) mul_mp (p, sqr, x_g, x_g, gdigs);
662 (void) mul_mp (p, pwr, sqr, x_g, gdigs);
663 (void) move_mp (z_g, x_g, gdigs);
664 (void) div_mp_digit (p, tmp, pwr, (MP_T) 6, gdigs);
665 (void) sub_mp (p, z_g, z_g, tmp, gdigs);
666 (void) mul_mp (p, pwr, pwr, sqr, gdigs);
667 (void) div_mp_digit (p, tmp, pwr, (MP_T) 120, gdigs);
668 (void) add_mp (p, z_g, z_g, tmp, gdigs);
669 (void) mul_mp (p, pwr, pwr, sqr, gdigs);
670 (void) div_mp_digit (p, tmp, pwr, (MP_T) 5040, gdigs);
671 (void) sub_mp (p, z_g, z_g, tmp, gdigs);
672 (void) mul_mp (p, pwr, pwr, sqr, gdigs);
673 (void) set_mp (fac, (MP_T) 362880, 0, gdigs);
674 int n = 9;
675 BOOL_T even = A68_TRUE, iter = (BOOL_T) (MP_DIGIT (pwr, 1) != 0);
676 while (iter) {
677 (void) div_mp (p, tmp, pwr, fac, gdigs);
678 if (MP_EXPONENT (tmp) <= (MP_EXPONENT (z_g) - gdigs)) {
679 iter = A68_FALSE;
680 } else {
681 if (even) {
682 (void) add_mp (p, z_g, z_g, tmp, gdigs);
683 even = A68_FALSE;
684 } else {
685 (void) sub_mp (p, z_g, z_g, tmp, gdigs);
686 even = A68_TRUE;
687 }
688 (void) mul_mp (p, pwr, pwr, sqr, gdigs);
689 (void) mul_mp_digit (p, fac, fac, (MP_T) (++n), gdigs);
690 (void) mul_mp_digit (p, fac, fac, (MP_T) (++n), gdigs);
691 }
692 }
693 // Inverse scaling using sin (3x) = sin (x) * (3 - 4 sin ** 2 (x)).
694 // Use existing mp's for intermediates.
695 (void) set_mp (fac, (MP_T) 3, 0, gdigs);
696 while (m--) {
697 (void) mul_mp (p, pwr, z_g, z_g, gdigs);
698 (void) mul_mp_digit (p, pwr, pwr, (MP_T) 4, gdigs);
699 (void) sub_mp (p, pwr, fac, pwr, gdigs);
700 (void) mul_mp (p, z_g, pwr, z_g, gdigs);
701 }
702 (void) shorten_mp (p, z, digs, z_g, gdigs);
703 if (neg ^ flip) {
704 MP_DIGIT (z, 1) = -MP_DIGIT (z, 1);
705 }
706 // Exit.
707 A68_SP = pop_sp;
708 return z;
709 }
710
711 //! @brief PROC (LONG REAL) LONG REAL cas
712
713 MP_T *cas_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
714 {
715 // Hartley kernel, which Hartley named cas (cosine and sine).
716 ADDR_T pop_sp = A68_SP;
717 MP_T *sinx = nil_mp (p, digs), *cosx = nil_mp (p, digs);
718 cos_mp (p, cosx, x, digs);
719 sin_mp (p, sinx, x, digs);
720 (void) add_mp (p, z, cosx, sinx, digs);
721 A68_SP = pop_sp;
722 return z;
723 }
724
725 //! @brief PROC (LONG REAL) LONG REAL cos
726
727 MP_T *cos_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
728 {
729 // Use cos (x) = sin (pi / 2 - x).
730 // Compute x mod 2 pi before subtracting to avoid cancellation.
731 ADDR_T pop_sp = A68_SP;
732 int gdigs = FUN_DIGITS (digs);
733 MP_T *hpi = nil_mp (p, gdigs), *tpi = nil_mp (p, gdigs);
734 MP_T *x_g = len_mp (p, x, digs, gdigs), *y = nil_mp (p, digs);
735 (void) mp_pi (p, hpi, MP_HALF_PI, gdigs);
736 (void) mp_pi (p, tpi, MP_TWO_PI, gdigs);
737 (void) mod_mp (p, x_g, x_g, tpi, gdigs);
738 (void) sub_mp (p, x_g, hpi, x_g, gdigs);
739 (void) shorten_mp (p, y, digs, x_g, gdigs);
740 (void) sin_mp (p, z, y, digs);
741 A68_SP = pop_sp;
742 return z;
743 }
744
745 //! @brief PROC (LONG REAL) LONG REAL tan
746
747 MP_T *tan_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
748 {
749 // Use tan (x) = sin (x) / sqrt (1 - sin ^ 2 (x)).
750 ADDR_T pop_sp = A68_SP;
751 int gdigs = FUN_DIGITS (digs);
752 MP_T *pi = nil_mp (p, gdigs), *hpi = nil_mp (p, gdigs);
753 (void) mp_pi (p, pi, MP_PI, gdigs);
754 (void) mp_pi (p, hpi, MP_HALF_PI, gdigs);
755 MP_T *x_g = len_mp (p, x, digs, gdigs), *y_g = nil_mp (p, gdigs);
756 MP_T *sns = nil_mp (p, digs), *cns = nil_mp (p, digs);
757 // Argument mod pi.
758 BOOL_T neg;
759 (void) mod_mp (p, x_g, x_g, pi, gdigs);
760 if (MP_DIGIT (x_g, 1) >= 0) {
761 (void) sub_mp (p, y_g, x_g, hpi, gdigs);
762 neg = (BOOL_T) (MP_DIGIT (y_g, 1) > 0);
763 } else {
764 (void) add_mp (p, y_g, x_g, hpi, gdigs);
765 neg = (BOOL_T) (MP_DIGIT (y_g, 1) < 0);
766 }
767 (void) shorten_mp (p, x, digs, x_g, gdigs);
768 // tan(x) = sin(x) / sqrt (1 - sin ** 2 (x)).
769 (void) sin_mp (p, sns, x, digs);
770 (void) mul_mp (p, cns, sns, sns, digs);
771 (void) one_minus_mp (p, cns, cns, digs);
772 (void) sqrt_mp (p, cns, cns, digs);
773 if (div_mp (p, z, sns, cns, digs) == NaN_MP) {
774 errno = EDOM;
775 A68_SP = pop_sp;
776 return NaN_MP;
777 }
778 A68_SP = pop_sp;
779 if (neg) {
780 MP_DIGIT (z, 1) = -MP_DIGIT (z, 1);
781 }
782 return z;
783 }
784
785 //! @brief PROC (LONG REAL) LONG REAL arcsin
786
787 MP_T *asin_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
788 {
789 ADDR_T pop_sp = A68_SP;
790 int gdigs = FUN_DIGITS (digs);
791 MP_T *x_g = len_mp (p, x, digs, gdigs), *z_g = nil_mp (p, gdigs);
792 MP_T *y = nil_mp (p, digs);
793 (void) mul_mp (p, z_g, x_g, x_g, gdigs);
794 (void) one_minus_mp (p, z_g, z_g, gdigs);
795 if (sqrt_mp (p, z_g, z_g, digs) == NaN_MP) {
796 errno = EDOM;
797 A68_SP = pop_sp;
798 return NaN_MP;
799 }
800 if (MP_DIGIT (z_g, 1) == 0) {
801 (void) mp_pi (p, z, MP_HALF_PI, digs);
802 MP_DIGIT (z, 1) = (MP_DIGIT (x_g, 1) >= 0 ? MP_DIGIT (z, 1) : -MP_DIGIT (z, 1));
803 A68_SP = pop_sp;
804 return z;
805 }
806 if (div_mp (p, x_g, x_g, z_g, gdigs) == NaN_MP) {
807 errno = EDOM;
808 A68_SP = pop_sp;
809 return NaN_MP;
810 }
811 (void) shorten_mp (p, y, digs, x_g, gdigs);
812 (void) atan_mp (p, z, y, digs);
813 A68_SP = pop_sp;
814 return z;
815 }
816
817 //! @brief PROC (LONG REAL) LONG REAL arccos
818
819 MP_T *acos_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
820 {
821 ADDR_T pop_sp = A68_SP;
822 int gdigs = FUN_DIGITS (digs);
823 BOOL_T neg = (BOOL_T) (MP_DIGIT (x, 1) < 0);
824 if (MP_DIGIT (x, 1) == 0) {
825 (void) mp_pi (p, z, MP_HALF_PI, digs);
826 A68_SP = pop_sp;
827 return z;
828 }
829 MP_T *x_g = len_mp (p, x, digs, gdigs), *z_g = nil_mp (p, gdigs);
830 MP_T *y = nil_mp (p, digs);
831 (void) mul_mp (p, z_g, x_g, x_g, gdigs);
832 (void) one_minus_mp (p, z_g, z_g, gdigs);
833 if (sqrt_mp (p, z_g, z_g, digs) == NaN_MP) {
834 errno = EDOM;
835 A68_SP = pop_sp;
836 return NaN_MP;
837 }
838 if (div_mp (p, x_g, z_g, x_g, gdigs) == NaN_MP) {
839 errno = EDOM;
840 A68_SP = pop_sp;
841 return NaN_MP;
842 }
843 (void) shorten_mp (p, y, digs, x_g, gdigs);
844 (void) atan_mp (p, z, y, digs);
845 if (neg) {
846 (void) mp_pi (p, y, MP_PI, digs);
847 (void) add_mp (p, z, z, y, digs);
848 }
849 A68_SP = pop_sp;
850 return z;
851 }
852
853 //! @brief PROC (LONG REAL) LONG REAL arctan
854
855 MP_T *atan_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
856 {
857 // Depending on the argument we choose either Taylor or Newton.
858 ADDR_T pop_sp = A68_SP;
859 if (MP_DIGIT (x, 1) == 0) {
860 A68_SP = pop_sp;
861 SET_MP_ZERO (z, digs);
862 return z;
863 }
864 int gdigs = FUN_DIGITS (digs);
865 MP_T *x_g = len_mp (p, x, digs, gdigs), *z_g = nil_mp (p, gdigs);
866 BOOL_T neg = (BOOL_T) (MP_DIGIT (x_g, 1) < 0);
867 if (neg) {
868 MP_DIGIT (x_g, 1) = -MP_DIGIT (x_g, 1);
869 }
870 // For larger arguments we use atan(x) = pi/2 - atan(1/x).
871 BOOL_T flip = (BOOL_T) (((MP_EXPONENT (x_g) > 0) || (MP_EXPONENT (x_g) == 0 && MP_DIGIT (x_g, 1) > 1)) && (MP_DIGIT (x_g, 1) != 0));
872 if (flip) {
873 (void) rec_mp (p, x_g, x_g, gdigs);
874 }
875 if (MP_EXPONENT (x_g) < -1 || (MP_EXPONENT (x_g) == -1 && MP_DIGIT (x_g, 1) < MP_RADIX / 100)) {
876 // Taylor sum for x close to zero.
877 // arctan (x) = x - x ** 3 / 3 + x ** 5 / 5 - x ** 7 / 7 + ..
878 // This is faster for small x and avoids cancellation
879 MP_T *pwr = nil_mp (p, gdigs), *sqr = nil_mp (p, gdigs), *tmp = nil_mp (p, gdigs);
880 (void) mul_mp (p, sqr, x_g, x_g, gdigs);
881 (void) mul_mp (p, pwr, sqr, x_g, gdigs);
882 (void) move_mp (z_g, x_g, gdigs);
883 int n = 3;
884 BOOL_T even = A68_FALSE, iter = (BOOL_T) (MP_DIGIT (pwr, 1) != 0);
885 while (iter) {
886 (void) div_mp_digit (p, tmp, pwr, (MP_T) n, gdigs);
887 if (MP_EXPONENT (tmp) <= (MP_EXPONENT (z_g) - gdigs)) {
888 iter = A68_FALSE;
889 } else {
890 if (even) {
891 (void) add_mp (p, z_g, z_g, tmp, gdigs);
892 even = A68_FALSE;
893 } else {
894 (void) sub_mp (p, z_g, z_g, tmp, gdigs);
895 even = A68_TRUE;
896 }
897 (void) mul_mp (p, pwr, pwr, sqr, gdigs);
898 n += 2;
899 }
900 }
901 } else {
902 // Newton's method: x<n+1> = x<n> - cos (x<n>) * (sin (x<n>) - a cos (x<n>)).
903 MP_T *sns = nil_mp (p, gdigs), *cns = nil_mp (p, gdigs), *tmp = nil_mp (p, gdigs);
904 // Construct an estimate.
905 #if (A68_LEVEL >= 3)
906 (void) double_to_mp (p, z_g, atan_double (mp_to_double (p, x_g, gdigs)), gdigs);
907 #else
908 (void) real_to_mp (p, z_g, atan (mp_to_real (p, x_g, gdigs)), gdigs);
909 #endif
910 int decimals = DOUBLE_ACCURACY;
911 do {
912 decimals <<= 1;
913 int hdigs = MINIMUM (1 + decimals / LOG_MP_RADIX, gdigs);
914 (void) sin_mp (p, sns, z_g, hdigs);
915 (void) mul_mp (p, tmp, sns, sns, hdigs);
916 (void) one_minus_mp (p, tmp, tmp, hdigs);
917 (void) sqrt_mp (p, cns, tmp, hdigs);
918 (void) mul_mp (p, tmp, x_g, cns, hdigs);
919 (void) sub_mp (p, tmp, sns, tmp, hdigs);
920 (void) mul_mp (p, tmp, tmp, cns, hdigs);
921 (void) sub_mp (p, z_g, z_g, tmp, hdigs);
922 } while (decimals < gdigs * LOG_MP_RADIX);
923 }
924 if (flip) {
925 MP_T *hpi = nil_mp (p, gdigs);
926 (void) sub_mp (p, z_g, mp_pi (p, hpi, MP_HALF_PI, gdigs), z_g, gdigs);
927 }
928 (void) shorten_mp (p, z, digs, z_g, gdigs);
929 MP_DIGIT (z, 1) = (neg ? -MP_DIGIT (z, 1) : MP_DIGIT (z, 1));
930 // Exit.
931 A68_SP = pop_sp;
932 return z;
933 }
934
935 //! @brief PROC (LONG REAL, LONG REAL) LONG REAL atan2
936
937 MP_T *atan2_mp (NODE_T * p, MP_T * z, MP_T * x, MP_T * y, int digs)
938 {
939 ADDR_T pop_sp = A68_SP;
940 MP_T *t = nil_mp (p, digs);
941 if (MP_DIGIT (x, 1) == 0 && MP_DIGIT (y, 1) == 0) {
942 errno = EDOM;
943 A68_SP = pop_sp;
944 return NaN_MP;
945 } else {
946 BOOL_T flip = (BOOL_T) (MP_DIGIT (y, 1) < 0);
947 MP_DIGIT (y, 1) = ABS (MP_DIGIT (y, 1));
948 if (IS_ZERO_MP (x)) {
949 (void) mp_pi (p, z, MP_HALF_PI, digs);
950 } else {
951 BOOL_T flop = (BOOL_T) (MP_DIGIT (x, 1) <= 0);
952 MP_DIGIT (x, 1) = ABS (MP_DIGIT (x, 1));
953 (void) div_mp (p, z, y, x, digs);
954 (void) atan_mp (p, z, z, digs);
955 if (flop) {
956 (void) mp_pi (p, t, MP_PI, digs);
957 (void) sub_mp (p, z, t, z, digs);
958 }
959 }
960 if (flip) {
961 MP_DIGIT (z, 1) = -MP_DIGIT (z, 1);
962 }
963 }
964 A68_SP = pop_sp;
965 return z;
966 }
967
968 //! @brief PROC (LONG REAL) LONG REAL csc
969
970 MP_T *csc_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
971 {
972 sin_mp (p, z, x, digs);
973 if (rec_mp (p, z, z, digs) == NaN_MP) {
974 return NaN_MP;
975 }
976 return z;
977 }
978
979 //! @brief PROC (LONG REAL) LONG REAL cscdg
980
981 MP_T *cscdg_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
982 {
983 sindg_mp (p, z, x, digs);
984 if (rec_mp (p, z, z, digs) == NaN_MP) {
985 return NaN_MP;
986 }
987 return z;
988 }
989
990 //! @brief PROC (LONG REAL) LONG REAL acsc
991
992 MP_T *acsc_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
993 {
994 if (rec_mp (p, z, x, digs) == NaN_MP) {
995 return NaN_MP;
996 }
997 return asin_mp (p, z, z, digs);
998 }
999
1000 //! @brief PROC (LONG REAL) LONG REAL acscdg
1001
1002 MP_T *acscdg_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1003 {
1004 if (rec_mp (p, z, x, digs) == NaN_MP) {
1005 return NaN_MP;
1006 }
1007 return asindg_mp (p, z, z, digs);
1008 }
1009
1010 //! @brief PROC (LONG REAL) LONG REAL sec
1011
1012 MP_T *sec_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1013 {
1014 cos_mp (p, z, x, digs);
1015 if (rec_mp (p, z, z, digs) == NaN_MP) {
1016 return NaN_MP;
1017 }
1018 return z;
1019 }
1020
1021 //! @brief PROC (LONG REAL) LONG REAL secdg
1022
1023 MP_T *secdg_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1024 {
1025 cosdg_mp (p, z, x, digs);
1026 if (rec_mp (p, z, z, digs) == NaN_MP) {
1027 return NaN_MP;
1028 }
1029 return z;
1030 }
1031
1032 //! @brief PROC (LONG REAL) LONG REAL asec
1033
1034 MP_T *asec_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1035 {
1036 if (rec_mp (p, z, x, digs) == NaN_MP) {
1037 return NaN_MP;
1038 }
1039 return acos_mp (p, z, z, digs);
1040 }
1041
1042 //! @brief PROC (LONG REAL) LONG REAL asec
1043
1044 MP_T *asecdg_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1045 {
1046 if (rec_mp (p, z, x, digs) == NaN_MP) {
1047 return NaN_MP;
1048 }
1049 return acosdg_mp (p, z, z, digs);
1050 }
1051
1052 //! @brief PROC (LONG REAL) LONG REAL cot
1053
1054 MP_T *cot_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1055 {
1056 // Use tan (x) = sin (x) / sqrt (1 - sin ^ 2 (x)).
1057 ADDR_T pop_sp = A68_SP;
1058 int gdigs = FUN_DIGITS (digs);
1059 MP_T *pi = nil_mp (p, gdigs), *hpi = nil_mp (p, gdigs);
1060 (void) mp_pi (p, pi, MP_PI, gdigs);
1061 (void) mp_pi (p, hpi, MP_HALF_PI, gdigs);
1062 MP_T *x_g = len_mp (p, x, digs, gdigs), *y_g = nil_mp (p, gdigs);
1063 MP_T *sns = nil_mp (p, digs), *cns = nil_mp (p, digs);
1064 // Argument mod pi.
1065 BOOL_T neg;
1066 (void) mod_mp (p, x_g, x_g, pi, gdigs);
1067 if (MP_DIGIT (x_g, 1) >= 0) {
1068 (void) sub_mp (p, y_g, x_g, hpi, gdigs);
1069 neg = (BOOL_T) (MP_DIGIT (y_g, 1) > 0);
1070 } else {
1071 (void) add_mp (p, y_g, x_g, hpi, gdigs);
1072 neg = (BOOL_T) (MP_DIGIT (y_g, 1) < 0);
1073 }
1074 (void) shorten_mp (p, x, digs, x_g, gdigs);
1075 // tan(x) = sin(x) / sqrt (1 - sin ** 2 (x)).
1076 (void) sin_mp (p, sns, x, digs);
1077 (void) mul_mp (p, cns, sns, sns, digs);
1078 (void) one_minus_mp (p, cns, cns, digs);
1079 (void) sqrt_mp (p, cns, cns, digs);
1080 if (div_mp (p, z, cns, sns, digs) == NaN_MP) {
1081 errno = EDOM;
1082 A68_SP = pop_sp;
1083 return NaN_MP;
1084 }
1085 A68_SP = pop_sp;
1086 if (neg) {
1087 MP_DIGIT (z, 1) = -MP_DIGIT (z, 1);
1088 }
1089 return z;
1090 }
1091
1092 //! @brief PROC (LONG REAL) LONG REAL arccot
1093
1094 MP_T *acot_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1095 {
1096 ADDR_T pop_sp = A68_SP;
1097 MP_T *f = nil_mp (p, digs);
1098 if (rec_mp (p, f, x, digs) == NaN_MP) {
1099 errno = EDOM;
1100 A68_SP = pop_sp;
1101 return NaN_MP;
1102 } else {
1103 (void) atan_mp (p, z, f, digs);
1104 A68_SP = pop_sp;
1105 return z;
1106 }
1107 }
1108
1109 //! @brief PROC (LONG REAL) LONG REAL sindg
1110
1111 MP_T *sindg_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1112 {
1113 ADDR_T pop_sp = A68_SP;
1114 MP_T *f = nil_mp (p, digs), *g = nil_mp (p, digs);
1115 (void) mp_pi (p, f, MP_PI_OVER_180, digs);
1116 (void) mul_mp (p, g, x, f, digs);
1117 (void) sin_mp (p, z, g, digs);
1118 A68_SP = pop_sp;
1119 return z;
1120 }
1121
1122 //! @brief PROC (LONG REAL) LONG REAL cosdg
1123
1124 MP_T *cosdg_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1125 {
1126 ADDR_T pop_sp = A68_SP;
1127 MP_T *f = nil_mp (p, digs), *g = nil_mp (p, digs);
1128 (void) mp_pi (p, f, MP_PI_OVER_180, digs);
1129 (void) mul_mp (p, g, x, f, digs);
1130 (void) cos_mp (p, z, g, digs);
1131 A68_SP = pop_sp;
1132 return z;
1133 }
1134
1135 //! @brief PROC (LONG REAL) LONG REAL tandg
1136
1137 MP_T *tandg_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1138 {
1139 ADDR_T pop_sp = A68_SP;
1140 MP_T *f = nil_mp (p, digs), *g = nil_mp (p, digs);
1141 (void) mp_pi (p, f, MP_PI_OVER_180, digs);
1142 (void) mul_mp (p, g, x, f, digs);
1143 if (tan_mp (p, z, g, digs) == NaN_MP) {
1144 errno = EDOM;
1145 A68_SP = pop_sp;
1146 return NaN_MP;
1147 } else {
1148 A68_SP = pop_sp;
1149 return z;
1150 }
1151 }
1152
1153 //! @brief PROC (LONG REAL) LONG REAL cotdg
1154
1155 MP_T *cotdg_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1156 {
1157 ADDR_T pop_sp = A68_SP;
1158 MP_T *f = nil_mp (p, digs), *g = nil_mp (p, digs);
1159 (void) mp_pi (p, f, MP_PI_OVER_180, digs);
1160 (void) mul_mp (p, g, x, f, digs);
1161 if (cot_mp (p, z, g, digs) == NaN_MP) {
1162 errno = EDOM;
1163 A68_SP = pop_sp;
1164 return NaN_MP;
1165 } else {
1166 A68_SP = pop_sp;
1167 return z;
1168 }
1169 }
1170
1171 //! @brief PROC (LONG REAL) LONG REAL arcsindg
1172
1173 MP_T *asindg_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1174 {
1175 ADDR_T pop_sp = A68_SP;
1176 MP_T *f = nil_mp (p, digs), *g = nil_mp (p, digs);
1177 (void) asin_mp (p, f, x, digs);
1178 (void) mp_pi (p, g, MP_180_OVER_PI, digs);
1179 (void) mul_mp (p, z, f, g, digs);
1180 A68_SP = pop_sp;
1181 return z;
1182 }
1183
1184 //! @brief PROC (LONG REAL) LONG REAL arccosdg
1185
1186 MP_T *acosdg_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1187 {
1188 ADDR_T pop_sp = A68_SP;
1189 MP_T *f = nil_mp (p, digs), *g = nil_mp (p, digs);
1190 (void) acos_mp (p, f, x, digs);
1191 (void) mp_pi (p, g, MP_180_OVER_PI, digs);
1192 (void) mul_mp (p, z, f, g, digs);
1193 A68_SP = pop_sp;
1194 return z;
1195 }
1196
1197 //! @brief PROC (LONG REAL) LONG REAL arctandg
1198
1199 MP_T *atandg_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1200 {
1201 ADDR_T pop_sp = A68_SP;
1202 MP_T *f = nil_mp (p, digs), *g = nil_mp (p, digs);
1203 (void) atan_mp (p, f, x, digs);
1204 (void) mp_pi (p, g, MP_180_OVER_PI, digs);
1205 (void) mul_mp (p, z, f, g, digs);
1206 A68_SP = pop_sp;
1207 return z;
1208 }
1209
1210 //! @brief PROC (LONG REAL) LONG REAL arccotdg
1211
1212 MP_T *acotdg_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1213 {
1214 ADDR_T pop_sp = A68_SP;
1215 MP_T *f = nil_mp (p, digs), *g = nil_mp (p, digs);
1216 if (acot_mp (p, f, x, digs) == NaN_MP) {
1217 errno = EDOM;
1218 A68_SP = pop_sp;
1219 return NaN_MP;
1220 } else {
1221 (void) mp_pi (p, g, MP_180_OVER_PI, digs);
1222 (void) mul_mp (p, z, f, g, digs);
1223 A68_SP = pop_sp;
1224 return z;
1225 }
1226 }
1227
1228 //! @brief PROC (LONG REAL, LONG REAL) LONG REAL atan2dg
1229
1230 MP_T *atan2dg_mp (NODE_T * p, MP_T * z, MP_T * x, MP_T * y, int digs)
1231 {
1232 ADDR_T pop_sp = A68_SP;
1233 MP_T *f = nil_mp (p, digs), *g = nil_mp (p, digs);
1234 if (atan2_mp (p, f, x, y, digs) == NaN_MP) {
1235 errno = EDOM;
1236 A68_SP = pop_sp;
1237 return NaN_MP;
1238 } else {
1239 (void) mp_pi (p, g, MP_180_OVER_PI, digs);
1240 (void) mul_mp (p, z, f, g, digs);
1241 A68_SP = pop_sp;
1242 return z;
1243 }
1244 }
1245
1246 //! @brief PROC (LONG REAL) LONG REAL sinpi
1247
1248 MP_T *sinpi_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1249 {
1250 ADDR_T pop_sp = A68_SP;
1251 MP_T *f = nil_mp (p, digs), *g = lit_mp (p, 1, 0, digs);
1252 (void) mod_mp (p, f, x, g, digs);
1253 if (IS_ZERO_MP (f)) {
1254 SET_MP_ZERO (z, digs);
1255 A68_SP = pop_sp;
1256 return z;
1257 }
1258 (void) mp_pi (p, f, MP_PI, digs);
1259 (void) mul_mp (p, g, x, f, digs);
1260 (void) sin_mp (p, z, g, digs);
1261 A68_SP = pop_sp;
1262 return z;
1263 }
1264
1265 //! @brief PROC (LONG REAL) LONG REAL acospi
1266
1267 MP_T *cospi_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1268 {
1269 ADDR_T pop_sp = A68_SP;
1270 MP_T *f = nil_mp (p, digs), *g = lit_mp (p, 1, 0, digs);
1271 (void) mod_mp (p, f, x, g, digs);
1272 abs_mp (p, f, f, digs);
1273 SET_MP_HALF (g, digs);
1274 (void) sub_mp (p, g, f, g, digs);
1275 if (IS_ZERO_MP (g)) {
1276 SET_MP_ZERO (z, digs);
1277 A68_SP = pop_sp;
1278 return z;
1279 }
1280 (void) mp_pi (p, f, MP_PI, digs);
1281 (void) mul_mp (p, g, x, f, digs);
1282 (void) cos_mp (p, z, g, digs);
1283 A68_SP = pop_sp;
1284 return z;
1285 }
1286
1287 //! @brief PROC (LONG REAL) LONG REAL tanpi
1288
1289 MP_T *tanpi_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1290 {
1291 ADDR_T pop_sp = A68_SP;
1292 MP_T *f = nil_mp (p, digs), *g = lit_mp (p, 1, 0, digs);
1293 MP_T *h = nil_mp (p, digs), *half = nil_mp (p, digs);
1294 SET_MP_ONE (g, digs);
1295 (void) mod_mp (p, f, x, g, digs);
1296 if (IS_ZERO_MP (f)) {
1297 SET_MP_ZERO (z, digs);
1298 A68_SP = pop_sp;
1299 return z;
1300 }
1301 SET_MP_MINUS_HALF (half, digs);
1302 (void) sub_mp (p, h, f, half, digs);
1303 if (MP_DIGIT (h, 1) < 0) {
1304 (void) plus_one_mp (p, f, f, digs);
1305 } else {
1306 SET_MP_HALF (half, digs);
1307 (void) sub_mp (p, h, f, half, digs);
1308 if (MP_DIGIT (h, 1) < 0) {
1309 (void) minus_one_mp (p, f, f, digs);
1310 }
1311 }
1312 BOOL_T neg = MP_DIGIT (f, 1) < 0;
1313 (void) abs_mp (p, f, f, digs);
1314 SET_MP_QUART (g, digs);
1315 (void) sub_mp (p, g, f, g, digs);
1316 if (IS_ZERO_MP (g)) {
1317 if (neg) {
1318 SET_MP_MINUS_ONE (z, digs);
1319 } else {
1320 SET_MP_ONE (z, digs);
1321 }
1322 A68_SP = pop_sp;
1323 return z;
1324 }
1325 (void) mp_pi (p, f, MP_PI, digs);
1326 (void) mul_mp (p, g, x, f, digs);
1327 if (tan_mp (p, z, g, digs) == NaN_MP) {
1328 errno = EDOM;
1329 A68_SP = pop_sp;
1330 return NaN_MP;
1331 } else {
1332 A68_SP = pop_sp;
1333 return z;
1334 }
1335 }
1336
1337 //! @brief PROC (LONG REAL) LONG REAL cotpi
1338
1339 MP_T *cotpi_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1340 {
1341 ADDR_T pop_sp = A68_SP;
1342 MP_T *f = nil_mp (p, digs), *g = lit_mp (p, 1, 0, digs);
1343 MP_T *h = nil_mp (p, digs), *half = nil_mp (p, digs);
1344 (void) mod_mp (p, f, x, g, digs);
1345 if (IS_ZERO_MP (f)) {
1346 errno = EDOM;
1347 A68_SP = pop_sp;
1348 return NaN_MP;
1349 }
1350 SET_MP_MINUS_HALF (half, digs);
1351 (void) sub_mp (p, h, f, half, digs);
1352 if (MP_DIGIT (h, 1) < 0) {
1353 (void) plus_one_mp (p, f, f, digs);
1354 } else {
1355 SET_MP_HALF (half, digs);
1356 (void) sub_mp (p, h, f, half, digs);
1357 if (MP_DIGIT (h, 1) < 0) {
1358 (void) minus_one_mp (p, f, f, digs);
1359 }
1360 }
1361 BOOL_T neg = MP_DIGIT (f, 1) < 0;
1362 (void) abs_mp (p, f, f, digs);
1363 SET_MP_QUART (g, digs);
1364 (void) sub_mp (p, g, f, g, digs);
1365 if (IS_ZERO_MP (g)) {
1366 if (neg) {
1367 SET_MP_MINUS_ONE (z, digs);
1368 } else {
1369 SET_MP_ONE (z, digs);
1370 }
1371 A68_SP = pop_sp;
1372 return z;
1373 }
1374 (void) mp_pi (p, f, MP_PI, digs);
1375 (void) mul_mp (p, g, x, f, digs);
1376 if (cot_mp (p, z, g, digs) == NaN_MP) {
1377 errno = EDOM;
1378 A68_SP = pop_sp;
1379 return NaN_MP;
1380 } else {
1381 A68_SP = pop_sp;
1382 return z;
1383 }
1384 }
