## mp-math.c

You can download the current version of Algol 68 Genie and its documentation here.

```   1 //! @file mp-math.c
2 //! @author J. Marcel van der Veer
3 //!
5 //!
6 //! This file is part of Algol68G - an Algol 68 compiler-interpreter.
7 //! Copyright 2001-2023 J. Marcel van der Veer .
8 //!
10 //!
11 //! This program is free software; you can redistribute it and/or modify it
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 LONG REAL math routines.
25
26 #include "a68g.h"
27 #include "a68g-genie.h"
28 #include "a68g-prelude.h"
29 #include "a68g-double.h"
30 #include "a68g-mp.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 {
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);
66   MP_T *x_g = len_mp (p, x, digs, gdigs);
67   MP_T *tmp = nil_mp (p, gdigs);
68 // Scaling for small x; sqrt (x) = 1 / sqrt (1 / x).
69   if ((reciprocal = (BOOL_T) (MP_EXPONENT (x_g) < 0)) == A68_TRUE) {
70     (void) rec_mp (p, x_g, x_g, gdigs);
71   }
72   if (ABS (MP_EXPONENT (x_g)) >= 2) {
73 // For extreme arguments we want accurate results as well.
74     int expo = (int) MP_EXPONENT (x_g);
75     MP_EXPONENT (x_g) = (MP_T) (expo % 2);
76     (void) sqrt_mp (p, z_g, x_g, gdigs);
77     MP_EXPONENT (z_g) += (MP_T) (expo / 2);
78   } else {
79 // Argument is in range. Estimate the root as double.
80 #if (A68_LEVEL >= 3)
81     DOUBLE_T x_d = mp_to_real_16 (p, x_g, gdigs);
82     (void) real_16_to_mp (p, z_g, sqrtq (x_d), gdigs);
83 #else
84     REAL_T x_d = mp_to_real (p, x_g, gdigs);
85     (void) real_to_mp (p, z_g, sqrt (x_d), gdigs);
86 #endif
87 // Newton's method: x = (x + a / x) / 2.
88     int decimals = DOUBLE_ACCURACY;
89     do {
90       decimals <<= 1;
91       hdigs = MINIMUM (1 + decimals / LOG_MP_RADIX, gdigs);
92       (void) div_mp (p, tmp, x_g, z_g, hdigs);
93       (void) add_mp (p, tmp, z_g, tmp, hdigs);
94       (void) half_mp (p, z_g, tmp, hdigs);
95     } while (decimals < 2 * gdigs * LOG_MP_RADIX);
96   }
97   if (reciprocal) {
98     (void) rec_mp (p, z_g, z_g, digs);
99   }
100   (void) shorten_mp (p, z, digs, z_g, gdigs);
101 // Exit.
102   A68_SP = pop_sp;
103   return z;
104 }
105
106 //! @brief PROC (LONG REAL) LONG REAL curt
107
108 MP_T *curt_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
109 {
111   if (MP_DIGIT (x, 1) == 0) {
112     A68_SP = pop_sp;
113     SET_MP_ZERO (z, digs);
114     return z;
115   }
116   BOOL_T change_sign = A68_FALSE;
117   if (MP_DIGIT (x, 1) < 0) {
118     change_sign = A68_TRUE;
119     MP_DIGIT (x, 1) = -MP_DIGIT (x, 1);
120   }
121   int gdigs = FUN_DIGITS (digs), hdigs;
122   BOOL_T reciprocal = A68_FALSE;
123   MP_T *z_g = nil_mp (p, gdigs);
124   MP_T *x_g = len_mp (p, x, digs, gdigs);
125   MP_T *tmp = nil_mp (p, gdigs);
126 // Scaling for small x; curt (x) = 1 / curt (1 / x).
127   if ((reciprocal = (BOOL_T) (MP_EXPONENT (x_g) < 0)) == A68_TRUE) {
128     (void) rec_mp (p, x_g, x_g, gdigs);
129   }
130   if (ABS (MP_EXPONENT (x_g)) >= 3) {
131 // For extreme arguments we want accurate results as well.
132     int expo = (int) MP_EXPONENT (x_g);
133     MP_EXPONENT (x_g) = (MP_T) (expo % 3);
134     (void) curt_mp (p, z_g, x_g, gdigs);
135     MP_EXPONENT (z_g) += (MP_T) (expo / 3);
136   } else {
137 // Argument is in range. Estimate the root as double.
138     int decimals;
139 #if (A68_LEVEL >= 3)
140     DOUBLE_T x_d = mp_to_real_16 (p, x_g, gdigs);
141     (void) real_16_to_mp (p, z_g, cbrtq (x_d), gdigs);
142 #else
143     REAL_T x_d = mp_to_real (p, x_g, gdigs);
144     (void) real_to_mp (p, z_g, cbrt (x_d), gdigs);
145 #endif
146 // Newton's method: x = (2 x + a / x ^ 2) / 3.
147     decimals = DOUBLE_ACCURACY;
148     do {
149       decimals <<= 1;
150       hdigs = MINIMUM (1 + decimals / LOG_MP_RADIX, gdigs);
151       (void) mul_mp (p, tmp, z_g, z_g, hdigs);
152       (void) div_mp (p, tmp, x_g, tmp, hdigs);
153       (void) add_mp (p, tmp, z_g, tmp, hdigs);
154       (void) add_mp (p, tmp, z_g, tmp, hdigs);
155       (void) div_mp_digit (p, z_g, tmp, (MP_T) 3, hdigs);
156     } while (decimals < gdigs * LOG_MP_RADIX);
157   }
158   if (reciprocal) {
159     (void) rec_mp (p, z_g, z_g, digs);
160   }
161   (void) shorten_mp (p, z, digs, z_g, gdigs);
162 // Exit.
163   A68_SP = pop_sp;
164   if (change_sign) {
165     MP_DIGIT (z, 1) = -MP_DIGIT (z, 1);
166   }
167   return z;
168 }
169
170 //! @brief PROC (LONG REAL) LONG REAL hypot
171
172 MP_T *hypot_mp (NODE_T * p, MP_T * z, MP_T * x, MP_T * y, int digs)
173 {
174 // sqrt (x^2 + y^2).
176   MP_T *t = nil_mp (p, digs);
177   MP_T *u = nil_mp (p, digs);
178   MP_T *v = nil_mp (p, digs);
179   (void) move_mp (u, x, digs);
180   (void) move_mp (v, y, digs);
181   MP_DIGIT (u, 1) = ABS (MP_DIGIT (u, 1));
182   MP_DIGIT (v, 1) = ABS (MP_DIGIT (v, 1));
183   if (IS_ZERO_MP (u)) {
184     (void) move_mp (z, v, digs);
185   } else if (IS_ZERO_MP (v)) {
186     (void) move_mp (z, u, digs);
187   } else {
188     SET_MP_ONE (t, digs);
189     (void) sub_mp (p, z, u, v, digs);
190     if (MP_DIGIT (z, 1) > 0) {
191       (void) div_mp (p, z, v, u, digs);
192       (void) mul_mp (p, z, z, z, digs);
193       (void) add_mp (p, z, t, z, digs);
194       (void) sqrt_mp (p, z, z, digs);
195       (void) mul_mp (p, z, u, z, digs);
196     } else {
197       (void) div_mp (p, z, u, v, digs);
198       (void) mul_mp (p, z, z, z, digs);
199       (void) add_mp (p, z, t, z, digs);
200       (void) sqrt_mp (p, z, z, digs);
201       (void) mul_mp (p, z, v, z, digs);
202     }
203   }
204   A68_SP = pop_sp;
205   return z;
206 }
207
208 //! @brief PROC (LONG REAL) LONG REAL exp
209
210 MP_T *exp_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
211 {
212 // Argument is reduced by using exp (z / (2 ** n)) ** (2 ** n) = exp(z).
213   int m, n, gdigs = FUN_DIGITS (digs);
215   BOOL_T iterate;
216   if (MP_DIGIT (x, 1) == 0) {
217     SET_MP_ONE (z, digs);
218     return z;
219   }
220   MP_T *x_g = len_mp (p, x, digs, gdigs);
221   MP_T *sum = nil_mp (p, gdigs);
222   MP_T *a68_pow = nil_mp (p, gdigs);
223   MP_T *fac = nil_mp (p, gdigs);
224   MP_T *tmp = nil_mp (p, gdigs);
225   m = 0;
226 // Scale x down.
227   while (eps_mp (x_g, gdigs)) {
228     m++;
229     (void) half_mp (p, x_g, x_g, gdigs);
230   }
231 // Calculate Taylor sum exp (z) = 1 + z / 1 ! + z ** 2 / 2 ! + ..
232   SET_MP_ONE (sum, gdigs);
233   (void) add_mp (p, sum, sum, x_g, gdigs);
234   (void) mul_mp (p, a68_pow, x_g, x_g, gdigs);
235   (void) half_mp (p, tmp, a68_pow, gdigs);
236   (void) add_mp (p, sum, sum, tmp, gdigs);
237   (void) mul_mp (p, a68_pow, a68_pow, x_g, gdigs);
238   (void) div_mp_digit (p, tmp, a68_pow, (MP_T) 6, gdigs);
239   (void) add_mp (p, sum, sum, tmp, gdigs);
240   (void) mul_mp (p, a68_pow, a68_pow, x_g, gdigs);
241   (void) div_mp_digit (p, tmp, a68_pow, (MP_T) 24, gdigs);
242   (void) add_mp (p, sum, sum, tmp, gdigs);
243   (void) mul_mp (p, a68_pow, a68_pow, x_g, gdigs);
244   (void) div_mp_digit (p, tmp, a68_pow, (MP_T) 120, gdigs);
245   (void) add_mp (p, sum, sum, tmp, gdigs);
246   (void) mul_mp (p, a68_pow, a68_pow, x_g, gdigs);
247   (void) div_mp_digit (p, tmp, a68_pow, (MP_T) 720, gdigs);
248   (void) add_mp (p, sum, sum, tmp, gdigs);
249   (void) mul_mp (p, a68_pow, a68_pow, x_g, gdigs);
250   (void) div_mp_digit (p, tmp, a68_pow, (MP_T) 5040, gdigs);
251   (void) add_mp (p, sum, sum, tmp, gdigs);
252   (void) mul_mp (p, a68_pow, a68_pow, x_g, gdigs);
253   (void) div_mp_digit (p, tmp, a68_pow, (MP_T) 40320, gdigs);
254   (void) add_mp (p, sum, sum, tmp, gdigs);
255   (void) mul_mp (p, a68_pow, a68_pow, x_g, gdigs);
256   (void) div_mp_digit (p, tmp, a68_pow, (MP_T) 362880, gdigs);
257   (void) add_mp (p, sum, sum, tmp, gdigs);
258   (void) mul_mp (p, a68_pow, a68_pow, x_g, gdigs);
259   (void) set_mp (fac, (MP_T) (MP_T) 3628800, 0, gdigs);
260   n = 10;
261   iterate = (BOOL_T) (MP_DIGIT (a68_pow, 1) != 0);
262   while (iterate) {
263     (void) div_mp (p, tmp, a68_pow, fac, gdigs);
264     if (MP_EXPONENT (tmp) <= (MP_EXPONENT (sum) - gdigs)) {
265       iterate = A68_FALSE;
266     } else {
267       (void) add_mp (p, sum, sum, tmp, gdigs);
268       (void) mul_mp (p, a68_pow, a68_pow, x_g, gdigs);
269       n++;
270       (void) mul_mp_digit (p, fac, fac, (MP_T) n, gdigs);
271     }
272   }
273 // Square exp (x) up.
274   while (m--) {
275     (void) mul_mp (p, sum, sum, sum, gdigs);
276   }
277   (void) shorten_mp (p, z, digs, sum, gdigs);
278 // Exit.
279   A68_SP = pop_sp;
280   return z;
281 }
282
283 //! @brief PROC (LONG REAL) LONG REAL exp (x) - 1
284 // assuming "x" to be close to 0.
285
286 MP_T *expm1_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
287 {
288   int n, gdigs = FUN_DIGITS (digs);
290   BOOL_T iterate;
291   if (MP_DIGIT (x, 1) == 0) {
292     SET_MP_ONE (z, digs);
293     return z;
294   }
295   MP_T *x_g = len_mp (p, x, digs, gdigs);
296   MP_T *sum = nil_mp (p, gdigs);
297   MP_T *a68_pow = nil_mp (p, gdigs);
298   MP_T *fac = nil_mp (p, gdigs);
299   MP_T *tmp = nil_mp (p, gdigs);
300 // Calculate Taylor sum expm1 (z) = z / 1 ! + z ** 2 / 2 ! + ...
301   (void) add_mp (p, sum, sum, x_g, gdigs);
302   (void) mul_mp (p, a68_pow, x_g, x_g, gdigs);
303   (void) half_mp (p, tmp, a68_pow, gdigs);
304   (void) add_mp (p, sum, sum, tmp, gdigs);
305   (void) mul_mp (p, a68_pow, a68_pow, x_g, gdigs);
306   (void) div_mp_digit (p, tmp, a68_pow, (MP_T) 6, gdigs);
307   (void) add_mp (p, sum, sum, tmp, gdigs);
308   (void) mul_mp (p, a68_pow, a68_pow, x_g, gdigs);
309   (void) div_mp_digit (p, tmp, a68_pow, (MP_T) 24, gdigs);
310   (void) add_mp (p, sum, sum, tmp, gdigs);
311   (void) mul_mp (p, a68_pow, a68_pow, x_g, gdigs);
312   (void) div_mp_digit (p, tmp, a68_pow, (MP_T) 120, gdigs);
313   (void) add_mp (p, sum, sum, tmp, gdigs);
314   (void) mul_mp (p, a68_pow, a68_pow, x_g, gdigs);
315   (void) div_mp_digit (p, tmp, a68_pow, (MP_T) 720, gdigs);
316   (void) add_mp (p, sum, sum, tmp, gdigs);
317   (void) mul_mp (p, a68_pow, a68_pow, x_g, gdigs);
318   (void) div_mp_digit (p, tmp, a68_pow, (MP_T) 5040, gdigs);
319   (void) add_mp (p, sum, sum, tmp, gdigs);
320   (void) mul_mp (p, a68_pow, a68_pow, x_g, gdigs);
321   (void) div_mp_digit (p, tmp, a68_pow, (MP_T) 40320, gdigs);
322   (void) add_mp (p, sum, sum, tmp, gdigs);
323   (void) mul_mp (p, a68_pow, a68_pow, x_g, gdigs);
324   (void) div_mp_digit (p, tmp, a68_pow, (MP_T) 362880, gdigs);
325   (void) add_mp (p, sum, sum, tmp, gdigs);
326   (void) mul_mp (p, a68_pow, a68_pow, x_g, gdigs);
327   (void) set_mp (fac, (MP_T) (MP_T) 3628800, 0, gdigs);
328   n = 10;
329   iterate = (BOOL_T) (MP_DIGIT (a68_pow, 1) != 0);
330   while (iterate) {
331     (void) div_mp (p, tmp, a68_pow, fac, gdigs);
332     if (MP_EXPONENT (tmp) <= (MP_EXPONENT (sum) - gdigs)) {
333       iterate = A68_FALSE;
334     } else {
335       (void) add_mp (p, sum, sum, tmp, gdigs);
336       (void) mul_mp (p, a68_pow, a68_pow, x_g, gdigs);
337       n++;
338       (void) mul_mp_digit (p, fac, fac, (MP_T) n, gdigs);
339     }
340   }
341   (void) shorten_mp (p, z, digs, sum, gdigs);
342 // Exit.
343   A68_SP = pop_sp;
344   return z;
345 }
346
347 //! @brief ln scale
348
349 MP_T *mp_ln_scale (NODE_T * p, MP_T * z, int digs)
350 {
352   int gdigs = FUN_DIGITS (digs);
353   MP_T *z_g = nil_mp (p, gdigs);
354 // First see if we can restore a previous calculation.
355   if (gdigs <= A68_MP (mp_ln_scale_size)) {
356     (void) move_mp (z_g, A68_MP (mp_ln_scale), gdigs);
357   } else {
358 // No luck with the kept value, we generate a longer one.
359     (void) set_mp (z_g, (MP_T) 1, 1, gdigs);
360     (void) ln_mp (p, z_g, z_g, gdigs);
361     A68_MP (mp_ln_scale) = (MP_T *) get_heap_space ((unt) SIZE_MP (gdigs));
362     (void) move_mp (A68_MP (mp_ln_scale), z_g, gdigs);
363     A68_MP (mp_ln_scale_size) = gdigs;
364   }
365   (void) shorten_mp (p, z, digs, z_g, gdigs);
366   A68_SP = pop_sp;
367   return z;
368 }
369
370 //! @brief ln 10
371
372 MP_T *mp_ln_10 (NODE_T * p, MP_T * z, int digs)
373 {
375   int gdigs = FUN_DIGITS (digs);
376   MP_T *z_g = nil_mp (p, gdigs);
377 // First see if we can restore a previous calculation.
378   if (gdigs <= A68_MP (mp_ln_10_size)) {
379     (void) move_mp (z_g, A68_MP (mp_ln_10), gdigs);
380   } else {
381 // No luck with the kept value, we generate a longer one.
382     (void) set_mp (z_g, (MP_T) 10, 0, gdigs);
383     (void) ln_mp (p, z_g, z_g, gdigs);
384     A68_MP (mp_ln_10) = (MP_T *) get_heap_space ((unt) SIZE_MP (gdigs));
385     (void) move_mp (A68_MP (mp_ln_10), z_g, gdigs);
386     A68_MP (mp_ln_10_size) = gdigs;
387   }
388   (void) shorten_mp (p, z, digs, z_g, gdigs);
389   A68_SP = pop_sp;
390   return z;
391 }
392
393 //! @brief PROC (LONG REAL) LONG REAL ln
394
395 MP_T *ln_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
396 {
397 // Depending on the argument we choose either Taylor or Newton.
399   int gdigs = FUN_DIGITS (digs);
400   BOOL_T negative, scale;
401   MP_T expo = 0;
402   if (MP_DIGIT (x, 1) <= 0) {
403     errno = EDOM;
404     return NaN_MP;
405   }
406   MP_T *x_g = len_mp (p, x, digs, gdigs);
407   MP_T *z_g = nil_mp (p, gdigs);
408 // We use ln (1 / x) = - ln (x).
409   negative = (BOOL_T) (MP_EXPONENT (x_g) < 0);
410   if (negative) {
411     (void) rec_mp (p, x_g, x_g, digs);
412   }
413 // We want correct results for extreme arguments. We scale when "x_g" exceeds
414 // "MP_RADIX ** +- 2", using ln (x * MP_RADIX ** n) = ln (x) + n * ln (MP_RADIX).
415   scale = (BOOL_T) (ABS (MP_EXPONENT (x_g)) >= 2);
416   if (scale) {
417     expo = MP_EXPONENT (x_g);
418     MP_EXPONENT (x_g) = (MP_T) 0;
419   }
420   if (MP_EXPONENT (x_g) == 0 && MP_DIGIT (x_g, 1) == 1 && MP_DIGIT (x_g, 2) == 0) {
421 // Taylor sum for x close to unity.
422 // ln (x) = (x - 1) - (x - 1) ** 2 / 2 + (x - 1) ** 3 / 3 - ...
423 // This is faster for small x and avoids cancellation.
424     int n = 2;
425     BOOL_T iterate;
426     MP_T *tmp = nil_mp (p, gdigs);
427     MP_T *a68_pow = nil_mp (p, gdigs);
428     (void) minus_one_mp (p, x_g, x_g, gdigs);
429     (void) mul_mp (p, a68_pow, x_g, x_g, gdigs);
430     (void) move_mp (z_g, x_g, gdigs);
431     iterate = (BOOL_T) (MP_DIGIT (a68_pow, 1) != 0);
432     while (iterate) {
433       (void) div_mp_digit (p, tmp, a68_pow, (MP_T) n, gdigs);
434       if (MP_EXPONENT (tmp) <= (MP_EXPONENT (z_g) - gdigs)) {
435         iterate = A68_FALSE;
436       } else {
437         MP_DIGIT (tmp, 1) = (EVEN (n) ? -MP_DIGIT (tmp, 1) : MP_DIGIT (tmp, 1));
438         (void) add_mp (p, z_g, z_g, tmp, gdigs);
439         (void) mul_mp (p, a68_pow, a68_pow, x_g, gdigs);
440         n++;
441       }
442     }
443   } else {
444 // Newton's method: x = x - 1 + a / exp (x).
445     MP_T *tmp = nil_mp (p, gdigs);
446 // Construct an estimate.
447 #if (A68_LEVEL >= 3)
448     (void) real_16_to_mp (p, z_g, logq (mp_to_real_16 (p, x_g, gdigs)), gdigs);
449 #else
450     (void) real_to_mp (p, z_g, log (mp_to_real (p, x_g, gdigs)), gdigs);
451 #endif
452     int decimals = DOUBLE_ACCURACY;
453     do {
454       int hdigs;
455       decimals <<= 1;
456       hdigs = MINIMUM (1 + decimals / LOG_MP_RADIX, gdigs);
457       (void) exp_mp (p, tmp, z_g, hdigs);
458       (void) div_mp (p, tmp, x_g, tmp, hdigs);
459       (void) minus_one_mp (p, z_g, z_g, hdigs);
460       (void) add_mp (p, z_g, z_g, tmp, hdigs);
461     } while (decimals < gdigs * LOG_MP_RADIX);
462   }
463 // Inverse scaling.
464   if (scale) {
465 // ln (x * MP_RADIX ** n) = ln (x) + n * ln (MP_RADIX).
466     MP_T *ln_base = nil_mp (p, gdigs);
467     (void) mp_ln_scale (p, ln_base, gdigs);
468     (void) mul_mp_digit (p, ln_base, ln_base, (MP_T) expo, gdigs);
469     (void) add_mp (p, z_g, z_g, ln_base, gdigs);
470   }
471   if (negative) {
472     MP_DIGIT (z_g, 1) = -MP_DIGIT (z_g, 1);
473   }
474   (void) shorten_mp (p, z, digs, z_g, gdigs);
475 // Exit.
476   A68_SP = pop_sp;
477   return z;
478 }
479
480 //! @brief PROC (LONG REAL) LONG REAL log
481
482 MP_T *log_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
483 {
485   MP_T *ln_10 = nil_mp (p, digs);
486   if (ln_mp (p, z, x, digs) == NaN_MP) {
487     errno = EDOM;
488     return NaN_MP;
489   }
490   (void) mp_ln_10 (p, ln_10, digs);
491   (void) div_mp (p, z, z, ln_10, digs);
492   A68_SP = pop_sp;
493   return z;
494 }
495
496 //! @brief sinh ("z") and cosh ("z")
497
498 MP_T *hyp_mp (NODE_T * p, MP_T * sh, MP_T * ch, MP_T * z, int digs)
499 {
501   MP_T *x_g = nil_mp (p, digs);
502   MP_T *y_g = nil_mp (p, digs);
503   MP_T *z_g = nil_mp (p, digs);
504   (void) move_mp (z_g, z, digs);
505   (void) exp_mp (p, x_g, z_g, digs);
506   (void) rec_mp (p, y_g, x_g, digs);
507   (void) add_mp (p, ch, x_g, y_g, digs);
508 // Avoid cancellation for sinh.
509   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)) {
510     (void) expm1_mp (p, x_g, z_g, digs);
511     MP_DIGIT (z_g, 1) = -MP_DIGIT (z_g, 1);
512     (void) expm1_mp (p, y_g, z_g, digs);
513   }
514   (void) sub_mp (p, sh, x_g, y_g, digs);
515   (void) half_mp (p, sh, sh, digs);
516   (void) half_mp (p, ch, ch, digs);
517   A68_SP = pop_sp;
518   return z;
519 }
520
521 //! @brief PROC (LONG REAL) LONG REAL sinh
522
523 MP_T *sinh_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
524 {
526   int gdigs = FUN_DIGITS (digs);
527   MP_T *x_g = len_mp (p, x, digs, gdigs);
528   MP_T *y_g = nil_mp (p, gdigs);
529   MP_T *z_g = nil_mp (p, gdigs);
530   (void) hyp_mp (p, z_g, y_g, x_g, gdigs);
531   (void) shorten_mp (p, z, digs, z_g, gdigs);
532   A68_SP = pop_sp;
533   return z;
534 }
535
536 //! @brief PROC (LONG REAL) LONG REAL asinh
537
538 MP_T *asinh_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
539 {
540   if (IS_ZERO_MP (x)) {
541     SET_MP_ZERO (z, digs);
542     return z;
543   } else {
545     int gdigs;
546     if (MP_EXPONENT (x) >= -1) {
547       gdigs = FUN_DIGITS (digs);
548     } else {
549 // Extra precision when x^2+1 gets close to 1.
550       gdigs = 2 * FUN_DIGITS (digs);
551     }
552     MP_T *x_g = len_mp (p, x, digs, gdigs);
553     MP_T *y_g = nil_mp (p, gdigs);
554     MP_T *z_g = nil_mp (p, gdigs);
555     (void) mul_mp (p, z_g, x_g, x_g, gdigs);
556     SET_MP_ONE (y_g, gdigs);
557     (void) add_mp (p, y_g, z_g, y_g, gdigs);
558     (void) sqrt_mp (p, y_g, y_g, gdigs);
559     (void) add_mp (p, y_g, y_g, x_g, gdigs);
560     (void) ln_mp (p, z_g, y_g, gdigs);
561     if (IS_ZERO_MP (z_g)) {
562       (void) move_mp (z, x, digs);
563     } else {
564       (void) shorten_mp (p, z, digs, z_g, gdigs);
565     }
566     A68_SP = pop_sp;
567     return z;
568   }
569 }
570
571 //! @brief PROC (LONG REAL) LONG REAL cosh
572
573 MP_T *cosh_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
574 {
576   int gdigs = FUN_DIGITS (digs);
577   MP_T *x_g = len_mp (p, x, digs, gdigs);
578   MP_T *y_g = nil_mp (p, gdigs);
579   MP_T *z_g = nil_mp (p, gdigs);
580   (void) hyp_mp (p, y_g, z_g, x_g, gdigs);
581   (void) shorten_mp (p, z, digs, z_g, gdigs);
582   A68_SP = pop_sp;
583   return z;
584 }
585
586 //! @brief PROC (LONG REAL) LONG REAL acosh
587
588 MP_T *acosh_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
589 {
591   int gdigs;
592   if (MP_DIGIT (x, 1) == 1 && MP_DIGIT (x, 2) == 0) {
593 // Extra precision when x^2-1 gets close to 0.
594     gdigs = 2 * FUN_DIGITS (digs);
595   } else {
596     gdigs = FUN_DIGITS (digs);
597   }
598   MP_T *x_g = len_mp (p, x, digs, gdigs);
599   MP_T *y_g = nil_mp (p, gdigs);
600   MP_T *z_g = nil_mp (p, gdigs);
601   (void) mul_mp (p, z_g, x_g, x_g, gdigs);
602   SET_MP_ONE (y_g, gdigs);
603   (void) sub_mp (p, y_g, z_g, y_g, gdigs);
604   (void) sqrt_mp (p, y_g, y_g, gdigs);
605   (void) add_mp (p, y_g, y_g, x_g, gdigs);
606   (void) ln_mp (p, z_g, y_g, gdigs);
607   (void) shorten_mp (p, z, digs, z_g, gdigs);
608   A68_SP = pop_sp;
609   return z;
610 }
611
612 //! @brief PROC (LONG REAL) LONG REAL tanh
613
614 MP_T *tanh_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
615 {
617   int gdigs = FUN_DIGITS (digs);
618   MP_T *x_g = len_mp (p, x, digs, gdigs);
619   MP_T *y_g = nil_mp (p, gdigs);
620   MP_T *z_g = nil_mp (p, gdigs);
621   (void) hyp_mp (p, y_g, z_g, x_g, gdigs);
622   (void) div_mp (p, z_g, y_g, z_g, gdigs);
623   (void) shorten_mp (p, z, digs, z_g, gdigs);
624   A68_SP = pop_sp;
625   return z;
626 }
627
628 //! @brief PROC (LONG REAL) LONG REAL atanh
629
630 MP_T *atanh_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
631 {
633   int gdigs = FUN_DIGITS (digs);
634   MP_T *x_g = len_mp (p, x, digs, gdigs);
635   MP_T *y_g = nil_mp (p, gdigs);
636   MP_T *z_g = nil_mp (p, gdigs);
637   SET_MP_ONE (y_g, gdigs);
638   (void) add_mp (p, z_g, y_g, x_g, gdigs);
639   (void) sub_mp (p, y_g, y_g, x_g, gdigs);
640   (void) div_mp (p, y_g, z_g, y_g, gdigs);
641   (void) ln_mp (p, z_g, y_g, gdigs);
642   (void) half_mp (p, z_g, z_g, gdigs);
643   (void) shorten_mp (p, z, digs, z_g, gdigs);
644   A68_SP = pop_sp;
645   return z;
646 }
647
648 //! @brief PROC (LONG REAL) LONG REAL sin
649
650 MP_T *sin_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
651 {
652 // Use triple-angle relation to reduce argument.
654   int m, n, gdigs = FUN_DIGITS (digs);
655   BOOL_T flip, negative, iterate, even;
656 // We will use "pi".
657   MP_T *pi = nil_mp (p, gdigs);
658   MP_T *tpi = nil_mp (p, gdigs);
659   MP_T *hpi = nil_mp (p, gdigs);
660   (void) mp_pi (p, pi, MP_PI, gdigs);
661   (void) mp_pi (p, tpi, MP_TWO_PI, gdigs);
662   (void) mp_pi (p, hpi, MP_HALF_PI, gdigs);
663 // Argument reduction (1): sin (x) = sin (x mod 2 pi).
664   MP_T *x_g = len_mp (p, x, digs, gdigs);
665   (void) mod_mp (p, x_g, x_g, tpi, gdigs);
666 // Argument reduction (2): sin (-x) = sin (x)
667 //                         sin (x) = - sin (x - pi); pi < x <= 2 pi
668 //                         sin (x) = sin (pi - x);   pi / 2 < x <= pi
669   negative = (BOOL_T) (MP_DIGIT (x_g, 1) < 0);
670   if (negative) {
671     MP_DIGIT (x_g, 1) = -MP_DIGIT (x_g, 1);
672   }
673   MP_T *tmp = nil_mp (p, gdigs);
674   (void) sub_mp (p, tmp, x_g, pi, gdigs);
675   flip = (BOOL_T) (MP_DIGIT (tmp, 1) > 0);
676   if (flip) {                   // x > pi
677     (void) sub_mp (p, x_g, x_g, pi, gdigs);
678   }
679   (void) sub_mp (p, tmp, x_g, hpi, gdigs);
680   if (MP_DIGIT (tmp, 1) > 0) {  // x > pi / 2
681     (void) sub_mp (p, x_g, pi, x_g, gdigs);
682   }
683 // Argument reduction (3): (follows from De Moivre's theorem)
684 // sin (3x) = sin (x) * (3 - 4 sin ^ 2 (x))
685   m = 0;
686   while (eps_mp (x_g, gdigs)) {
687     m++;
688     (void) div_mp_digit (p, x_g, x_g, (MP_T) 3, gdigs);
689   }
690 // Taylor sum.
691   MP_T *sqr = nil_mp (p, gdigs);
692   MP_T *a68_pow = nil_mp (p, gdigs);
693   MP_T *fac = nil_mp (p, gdigs);
694   MP_T *z_g = nil_mp (p, gdigs);
695   (void) mul_mp (p, sqr, x_g, x_g, gdigs);
696   (void) mul_mp (p, a68_pow, sqr, x_g, gdigs);
697   (void) move_mp (z_g, x_g, gdigs);
698   (void) div_mp_digit (p, tmp, a68_pow, (MP_T) 6, gdigs);
699   (void) sub_mp (p, z_g, z_g, tmp, gdigs);
700   (void) mul_mp (p, a68_pow, a68_pow, sqr, gdigs);
701   (void) div_mp_digit (p, tmp, a68_pow, (MP_T) 120, gdigs);
702   (void) add_mp (p, z_g, z_g, tmp, gdigs);
703   (void) mul_mp (p, a68_pow, a68_pow, sqr, gdigs);
704   (void) div_mp_digit (p, tmp, a68_pow, (MP_T) 5040, gdigs);
705   (void) sub_mp (p, z_g, z_g, tmp, gdigs);
706   (void) mul_mp (p, a68_pow, a68_pow, sqr, gdigs);
707   (void) set_mp (fac, (MP_T) 362880, 0, gdigs);
708   n = 9;
709   even = A68_TRUE;
710   iterate = (BOOL_T) (MP_DIGIT (a68_pow, 1) != 0);
711   while (iterate) {
712     (void) div_mp (p, tmp, a68_pow, fac, gdigs);
713     if (MP_EXPONENT (tmp) <= (MP_EXPONENT (z_g) - gdigs)) {
714       iterate = A68_FALSE;
715     } else {
716       if (even) {
717         (void) add_mp (p, z_g, z_g, tmp, gdigs);
718         even = A68_FALSE;
719       } else {
720         (void) sub_mp (p, z_g, z_g, tmp, gdigs);
721         even = A68_TRUE;
722       }
723       (void) mul_mp (p, a68_pow, a68_pow, sqr, gdigs);
724       (void) mul_mp_digit (p, fac, fac, (MP_T) (++n), gdigs);
725       (void) mul_mp_digit (p, fac, fac, (MP_T) (++n), gdigs);
726     }
727   }
728 // Inverse scaling using sin (3x) = sin (x) * (3 - 4 sin ** 2 (x)).
729 // Use existing mp's for intermediates.
730   (void) set_mp (fac, (MP_T) 3, 0, gdigs);
731   while (m--) {
732     (void) mul_mp (p, a68_pow, z_g, z_g, gdigs);
733     (void) mul_mp_digit (p, a68_pow, a68_pow, (MP_T) 4, gdigs);
734     (void) sub_mp (p, a68_pow, fac, a68_pow, gdigs);
735     (void) mul_mp (p, z_g, a68_pow, z_g, gdigs);
736   }
737   (void) shorten_mp (p, z, digs, z_g, gdigs);
738   if (negative ^ flip) {
739     MP_DIGIT (z, 1) = -MP_DIGIT (z, 1);
740   }
741 // Exit.
742   A68_SP = pop_sp;
743   return z;
744 }
745
746 //! @brief PROC (LONG REAL) LONG REAL cos
747
748 MP_T *cos_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
749 {
750 // Use cos (x) = sin (pi / 2 - x).
751 // Compute x mod 2 pi before subtracting to avoid cancellation.
753   int gdigs = FUN_DIGITS (digs);
754   MP_T *hpi = nil_mp (p, gdigs);
755   MP_T *tpi = nil_mp (p, gdigs);
756   MP_T *x_g = len_mp (p, x, digs, gdigs);
757   MP_T *y = nil_mp (p, digs);
758   (void) mp_pi (p, hpi, MP_HALF_PI, gdigs);
759   (void) mp_pi (p, tpi, MP_TWO_PI, gdigs);
760   (void) mod_mp (p, x_g, x_g, tpi, gdigs);
761   (void) sub_mp (p, x_g, hpi, x_g, gdigs);
762   (void) shorten_mp (p, y, digs, x_g, gdigs);
763   (void) sin_mp (p, z, y, digs);
764   A68_SP = pop_sp;
765   return z;
766 }
767
768 //! @brief PROC (LONG REAL) LONG REAL tan
769
770 MP_T *tan_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
771 {
772 // Use tan (x) = sin (x) / sqrt (1 - sin ^ 2 (x)).
774   int gdigs = FUN_DIGITS (digs);
775   BOOL_T negate;
776   MP_T *pi = nil_mp (p, gdigs);
777   MP_T *hpi = nil_mp (p, gdigs);
778   (void) mp_pi (p, pi, MP_PI, gdigs);
779   (void) mp_pi (p, hpi, MP_HALF_PI, gdigs);
780   MP_T *x_g = len_mp (p, x, digs, gdigs);
781   MP_T *y_g = nil_mp (p, gdigs);
782   MP_T *sns = nil_mp (p, digs);
783   MP_T *cns = nil_mp (p, digs);
784 // Argument mod pi.
785   (void) mod_mp (p, x_g, x_g, pi, gdigs);
786   if (MP_DIGIT (x_g, 1) >= 0) {
787     (void) sub_mp (p, y_g, x_g, hpi, gdigs);
788     negate = (BOOL_T) (MP_DIGIT (y_g, 1) > 0);
789   } else {
790     (void) add_mp (p, y_g, x_g, hpi, gdigs);
791     negate = (BOOL_T) (MP_DIGIT (y_g, 1) < 0);
792   }
793   (void) shorten_mp (p, x, digs, x_g, gdigs);
794 // tan(x) = sin(x) / sqrt (1 - sin ** 2 (x)).
795   (void) sin_mp (p, sns, x, digs);
796   (void) mul_mp (p, cns, sns, sns, digs);
797   (void) one_minus_mp (p, cns, cns, digs);
798   (void) sqrt_mp (p, cns, cns, digs);
799   if (div_mp (p, z, sns, cns, digs) == NaN_MP) {
800     errno = EDOM;
801     A68_SP = pop_sp;
802     return NaN_MP;
803   }
804   A68_SP = pop_sp;
805   if (negate) {
806     MP_DIGIT (z, 1) = -MP_DIGIT (z, 1);
807   }
808   return z;
809 }
810
811 //! @brief PROC (LONG REAL) LONG REAL arcsin
812
813 MP_T *asin_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
814 {
816   int gdigs = FUN_DIGITS (digs);
817   MP_T *y = nil_mp (p, digs);
818   MP_T *x_g = len_mp (p, x, digs, gdigs);
819   MP_T *z_g = nil_mp (p, gdigs);
820   (void) mul_mp (p, z_g, x_g, x_g, gdigs);
821   (void) one_minus_mp (p, z_g, z_g, gdigs);
822   if (sqrt_mp (p, z_g, z_g, digs) == NaN_MP) {
823     errno = EDOM;
824     A68_SP = pop_sp;
825     return NaN_MP;
826   }
827   if (MP_DIGIT (z_g, 1) == 0) {
828     (void) mp_pi (p, z, MP_HALF_PI, digs);
829     MP_DIGIT (z, 1) = (MP_DIGIT (x_g, 1) >= 0 ? MP_DIGIT (z, 1) : -MP_DIGIT (z, 1));
830     A68_SP = pop_sp;
831     return z;
832   }
833   if (div_mp (p, x_g, x_g, z_g, gdigs) == NaN_MP) {
834     errno = EDOM;
835     A68_SP = pop_sp;
836     return NaN_MP;
837   }
838   (void) shorten_mp (p, y, digs, x_g, gdigs);
839   (void) atan_mp (p, z, y, digs);
840   A68_SP = pop_sp;
841   return z;
842 }
843
844 //! @brief PROC (LONG REAL) LONG REAL arccos
845
846 MP_T *acos_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
847 {
849   int gdigs = FUN_DIGITS (digs);
850   BOOL_T negative = (BOOL_T) (MP_DIGIT (x, 1) < 0);
851   if (MP_DIGIT (x, 1) == 0) {
852     (void) mp_pi (p, z, MP_HALF_PI, digs);
853     A68_SP = pop_sp;
854     return z;
855   }
856   MP_T *y = nil_mp (p, digs);
857   MP_T *x_g = len_mp (p, x, digs, gdigs);
858   MP_T *z_g = nil_mp (p, gdigs);
859   (void) mul_mp (p, z_g, x_g, x_g, gdigs);
860   (void) one_minus_mp (p, z_g, z_g, gdigs);
861   if (sqrt_mp (p, z_g, z_g, digs) == NaN_MP) {
862     errno = EDOM;
863     A68_SP = pop_sp;
864     return NaN_MP;
865   }
866   if (div_mp (p, x_g, z_g, x_g, gdigs) == NaN_MP) {
867     errno = EDOM;
868     A68_SP = pop_sp;
869     return NaN_MP;
870   }
871   (void) shorten_mp (p, y, digs, x_g, gdigs);
872   (void) atan_mp (p, z, y, digs);
873   if (negative) {
874     (void) mp_pi (p, y, MP_PI, digs);
875     (void) add_mp (p, z, z, y, digs);
876   }
877   A68_SP = pop_sp;
878   return z;
879 }
880
881 //! @brief PROC (LONG REAL) LONG REAL arctan
882
883 MP_T *atan_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
884 {
885 // Depending on the argument we choose either Taylor or Newton.
887   if (MP_DIGIT (x, 1) == 0) {
888     A68_SP = pop_sp;
889     SET_MP_ZERO (z, digs);
890     return z;
891   }
892   int gdigs = FUN_DIGITS (digs);
893   MP_T *x_g = len_mp (p, x, digs, gdigs);
894   MP_T *z_g = nil_mp (p, gdigs);
895   BOOL_T negative = (BOOL_T) (MP_DIGIT (x_g, 1) < 0);
896   if (negative) {
897     MP_DIGIT (x_g, 1) = -MP_DIGIT (x_g, 1);
898   }
899 // For larger arguments we use atan(x) = pi/2 - atan(1/x).
900   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));
901   if (flip) {
902     (void) rec_mp (p, x_g, x_g, gdigs);
903   }
904   if (MP_EXPONENT (x_g) < -1 || (MP_EXPONENT (x_g) == -1 && MP_DIGIT (x_g, 1) < MP_RADIX / 100)) {
905 // Taylor sum for x close to zero.
906 // arctan (x) = x - x ** 3 / 3 + x ** 5 / 5 - x ** 7 / 7 + ..
907 // This is faster for small x and avoids cancellation
908     int n = 3;
909     BOOL_T iterate, even;
910     MP_T *tmp = nil_mp (p, gdigs);
911     MP_T *a68_pow = nil_mp (p, gdigs);
912     MP_T *sqr = nil_mp (p, gdigs);
913     (void) mul_mp (p, sqr, x_g, x_g, gdigs);
914     (void) mul_mp (p, a68_pow, sqr, x_g, gdigs);
915     (void) move_mp (z_g, x_g, gdigs);
916     even = A68_FALSE;
917     iterate = (BOOL_T) (MP_DIGIT (a68_pow, 1) != 0);
918     while (iterate) {
919       (void) div_mp_digit (p, tmp, a68_pow, (MP_T) n, gdigs);
920       if (MP_EXPONENT (tmp) <= (MP_EXPONENT (z_g) - gdigs)) {
921         iterate = A68_FALSE;
922       } else {
923         if (even) {
924           (void) add_mp (p, z_g, z_g, tmp, gdigs);
925           even = A68_FALSE;
926         } else {
927           (void) sub_mp (p, z_g, z_g, tmp, gdigs);
928           even = A68_TRUE;
929         }
930         (void) mul_mp (p, a68_pow, a68_pow, sqr, gdigs);
931         n += 2;
932       }
933     }
934   } else {
935 // Newton's method: x = x - cos (x) * (sin (x) - a cos (x)).
936     int decimals, hdigs;
937     MP_T *tmp = nil_mp (p, gdigs);
938     MP_T *sns = nil_mp (p, gdigs);
939     MP_T *cns = nil_mp (p, gdigs);
940 // Construct an estimate.
941 #if (A68_LEVEL >= 3)
942     (void) real_16_to_mp (p, z_g, atanq (mp_to_real_16 (p, x_g, gdigs)), gdigs);
943 #else
944     (void) real_to_mp (p, z_g, atan (mp_to_real (p, x_g, gdigs)), gdigs);
945 #endif
946     decimals = DOUBLE_ACCURACY;
947     do {
948       decimals <<= 1;
949       hdigs = MINIMUM (1 + decimals / LOG_MP_RADIX, gdigs);
950       (void) sin_mp (p, sns, z_g, hdigs);
951       (void) mul_mp (p, tmp, sns, sns, hdigs);
952       (void) one_minus_mp (p, tmp, tmp, hdigs);
953       (void) sqrt_mp (p, cns, tmp, hdigs);
954       (void) mul_mp (p, tmp, x_g, cns, hdigs);
955       (void) sub_mp (p, tmp, sns, tmp, hdigs);
956       (void) mul_mp (p, tmp, tmp, cns, hdigs);
957       (void) sub_mp (p, z_g, z_g, tmp, hdigs);
958     } while (decimals < gdigs * LOG_MP_RADIX);
959   }
960   if (flip) {
961     MP_T *hpi = nil_mp (p, gdigs);
962     (void) sub_mp (p, z_g, mp_pi (p, hpi, MP_HALF_PI, gdigs), z_g, gdigs);
963   }
964   (void) shorten_mp (p, z, digs, z_g, gdigs);
965   MP_DIGIT (z, 1) = (negative ? -MP_DIGIT (z, 1) : MP_DIGIT (z, 1));
966 // Exit.
967   A68_SP = pop_sp;
968   return z;
969 }
970
971 //! @brief PROC (LONG REAL, LONG REAL) LONG REAL atan2
972
973 MP_T *atan2_mp (NODE_T * p, MP_T * z, MP_T * x, MP_T * y, int digs)
974 {
976   MP_T *t = nil_mp (p, digs);
977   if (MP_DIGIT (x, 1) == 0 && MP_DIGIT (y, 1) == 0) {
978     errno = EDOM;
979     A68_SP = pop_sp;
980     return NaN_MP;
981   } else {
982     BOOL_T flip = (BOOL_T) (MP_DIGIT (y, 1) < 0);
983     MP_DIGIT (y, 1) = ABS (MP_DIGIT (y, 1));
984     if (IS_ZERO_MP (x)) {
985       (void) mp_pi (p, z, MP_HALF_PI, digs);
986     } else {
987       BOOL_T flop = (BOOL_T) (MP_DIGIT (x, 1) <= 0);
988       MP_DIGIT (x, 1) = ABS (MP_DIGIT (x, 1));
989       (void) div_mp (p, z, y, x, digs);
990       (void) atan_mp (p, z, z, digs);
991       if (flop) {
992         (void) mp_pi (p, t, MP_PI, digs);
993         (void) sub_mp (p, z, t, z, digs);
994       }
995     }
996     if (flip) {
997       MP_DIGIT (z, 1) = -MP_DIGIT (z, 1);
998     }
999   }
1000   A68_SP = pop_sp;
1001   return z;
1002 }
1003
1004 //! @brief PROC (LONG REAL) LONG REAL csc
1005
1006 MP_T *csc_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1007 {
1008   sin_mp (p, z, x, digs);
1009   if (rec_mp (p, z, z, digs) == NaN_MP) {
1010     return NaN_MP;
1011   }
1012   return z;
1013 }
1014
1015 //! @brief PROC (LONG REAL) LONG REAL acsc
1016
1017 MP_T *acsc_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1018 {
1019   if (rec_mp (p, z, x, digs) == NaN_MP) {
1020     return NaN_MP;
1021   }
1022   return asin_mp (p, z, z, digs);
1023 }
1024
1025 //! @brief PROC (LONG REAL) LONG REAL sec
1026
1027 MP_T *sec_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1028 {
1029   cos_mp (p, z, x, digs);
1030   if (rec_mp (p, z, z, digs) == NaN_MP) {
1031     return NaN_MP;
1032   }
1033   return z;
1034 }
1035
1036 //! @brief PROC (LONG REAL) LONG REAL asec
1037
1038 MP_T *asec_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1039 {
1040   if (rec_mp (p, z, x, digs) == NaN_MP) {
1041     return NaN_MP;
1042   }
1043   return acos_mp (p, z, z, digs);
1044 }
1045
1046 //! @brief PROC (LONG REAL) LONG REAL cot
1047
1048 MP_T *cot_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1049 {
1050 // Use tan (x) = sin (x) / sqrt (1 - sin ^ 2 (x)).
1052   int gdigs = FUN_DIGITS (digs);
1053   BOOL_T negate;
1054   MP_T *pi = nil_mp (p, gdigs);
1055   MP_T *hpi = nil_mp (p, gdigs);
1056   (void) mp_pi (p, pi, MP_PI, gdigs);
1057   (void) mp_pi (p, hpi, MP_HALF_PI, gdigs);
1058   MP_T *x_g = len_mp (p, x, digs, gdigs);
1059   MP_T *y_g = nil_mp (p, gdigs);
1060   MP_T *sns = nil_mp (p, digs);
1061   MP_T *cns = nil_mp (p, digs);
1062 // Argument mod pi.
1063   (void) mod_mp (p, x_g, x_g, pi, gdigs);
1064   if (MP_DIGIT (x_g, 1) >= 0) {
1065     (void) sub_mp (p, y_g, x_g, hpi, gdigs);
1066     negate = (BOOL_T) (MP_DIGIT (y_g, 1) > 0);
1067   } else {
1068     (void) add_mp (p, y_g, x_g, hpi, gdigs);
1069     negate = (BOOL_T) (MP_DIGIT (y_g, 1) < 0);
1070   }
1071   (void) shorten_mp (p, x, digs, x_g, gdigs);
1072 // tan(x) = sin(x) / sqrt (1 - sin ** 2 (x)).
1073   (void) sin_mp (p, sns, x, digs);
1074   (void) mul_mp (p, cns, sns, sns, digs);
1075   (void) one_minus_mp (p, cns, cns, digs);
1076   (void) sqrt_mp (p, cns, cns, digs);
1077   if (div_mp (p, z, cns, sns, digs) == NaN_MP) {
1078     errno = EDOM;
1079     A68_SP = pop_sp;
1080     return NaN_MP;
1081   }
1082   A68_SP = pop_sp;
1083   if (negate) {
1084     MP_DIGIT (z, 1) = -MP_DIGIT (z, 1);
1085   }
1086   return z;
1087 }
1088
1089 //! @brief PROC (LONG REAL) LONG REAL arccot
1090
1091 MP_T *acot_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1092 {
1094   MP_T *f = nil_mp (p, digs);
1095   if (rec_mp (p, f, x, digs) == NaN_MP) {
1096     errno = EDOM;
1097     A68_SP = pop_sp;
1098     return NaN_MP;
1099   } else {
1100     (void) atan_mp (p, z, f, digs);
1101     A68_SP = pop_sp;
1102     return z;
1103   }
1104 }
1105
1106 //! @brief PROC (LONG REAL) LONG REAL sindg
1107
1108 MP_T *sindg_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1109 {
1111   MP_T *f = nil_mp (p, digs);
1112   MP_T *g = nil_mp (p, digs);
1113   (void) mp_pi (p, f, MP_PI_OVER_180, digs);
1114   (void) mul_mp (p, g, x, f, digs);
1115   (void) sin_mp (p, z, g, digs);
1116   A68_SP = pop_sp;
1117   return z;
1118 }
1119
1120 //! @brief PROC (LONG REAL) LONG REAL cosdg
1121
1122 MP_T *cosdg_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1123 {
1125   MP_T *f = nil_mp (p, digs);
1126   MP_T *g = nil_mp (p, digs);
1127   (void) mp_pi (p, f, MP_PI_OVER_180, digs);
1128   (void) mul_mp (p, g, x, f, digs);
1129   (void) cos_mp (p, z, g, digs);
1130   A68_SP = pop_sp;
1131   return z;
1132 }
1133
1134 //! @brief PROC (LONG REAL) LONG REAL tandg
1135
1136 MP_T *tandg_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1137 {
1139   MP_T *f = nil_mp (p, digs);
1140   MP_T *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 {
1158   MP_T *f = nil_mp (p, digs);
1159   MP_T *g = nil_mp (p, digs);
1160   (void) mp_pi (p, f, MP_PI_OVER_180, digs);
1161   (void) mul_mp (p, g, x, f, digs);
1162   if (cot_mp (p, z, g, digs) == NaN_MP) {
1163     errno = EDOM;
1164     A68_SP = pop_sp;
1165     return NaN_MP;
1166   } else {
1167     A68_SP = pop_sp;
1168     return z;
1169   }
1170 }
1171
1172 //! @brief PROC (LONG REAL) LONG REAL arcsindg
1173
1174 MP_T *asindg_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1175 {
1177   MP_T *f = nil_mp (p, digs);
1178   MP_T *g = nil_mp (p, digs);
1179   (void) asin_mp (p, f, x, digs);
1180   (void) mp_pi (p, g, MP_180_OVER_PI, digs);
1181   (void) mul_mp (p, z, f, g, digs);
1182   A68_SP = pop_sp;
1183   return z;
1184 }
1185
1186 //! @brief PROC (LONG REAL) LONG REAL arccosdg
1187
1188 MP_T *acosdg_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1189 {
1191   MP_T *f = nil_mp (p, digs);
1192   MP_T *g = nil_mp (p, digs);
1193   (void) acos_mp (p, f, x, digs);
1194   (void) mp_pi (p, g, MP_180_OVER_PI, digs);
1195   (void) mul_mp (p, z, f, g, digs);
1196   A68_SP = pop_sp;
1197   return z;
1198 }
1199
1200 //! @brief PROC (LONG REAL) LONG REAL arctandg
1201
1202 MP_T *atandg_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1203 {
1205   MP_T *f = nil_mp (p, digs);
1206   MP_T *g = nil_mp (p, digs);
1207   (void) atan_mp (p, f, x, digs);
1208   (void) mp_pi (p, g, MP_180_OVER_PI, digs);
1209   (void) mul_mp (p, z, f, g, digs);
1210   A68_SP = pop_sp;
1211   return z;
1212 }
1213
1214 //! @brief PROC (LONG REAL) LONG REAL arccotdg
1215
1216 MP_T *acotdg_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1217 {
1219   MP_T *f = nil_mp (p, digs);
1220   MP_T *g = nil_mp (p, digs);
1221   if (acot_mp (p, f, x, digs) == NaN_MP) {
1222     errno = EDOM;
1223     A68_SP = pop_sp;
1224     return NaN_MP;
1225   } else {
1226     (void) mp_pi (p, g, MP_180_OVER_PI, digs);
1227     (void) mul_mp (p, z, f, g, digs);
1228     A68_SP = pop_sp;
1229     return z;
1230   }
1231 }
1232
1233 //! @brief PROC (LONG REAL, LONG REAL) LONG REAL atan2dg
1234
1235 MP_T *atan2dg_mp (NODE_T * p, MP_T * z, MP_T * x, MP_T * y, int digs)
1236 {
1238   MP_T *f = nil_mp (p, digs);
1239   MP_T *g = nil_mp (p, digs);
1240   if (atan2_mp (p, f, x, y, digs) == NaN_MP) {
1241     errno = EDOM;
1242     A68_SP = pop_sp;
1243     return NaN_MP;
1244   } else {
1245     (void) mp_pi (p, g, MP_180_OVER_PI, digs);
1246     (void) mul_mp (p, z, f, g, digs);
1247     A68_SP = pop_sp;
1248     return z;
1249   }
1250 }
1251
1252 //! @brief PROC (LONG REAL) LONG REAL sinpi
1253
1254 MP_T *sinpi_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1255 {
1257   MP_T *f = nil_mp (p, digs);
1258   MP_T *g = lit_mp (p, 1, 0, digs);
1259   (void) mod_mp (p, f, x, g, digs);
1260   if (IS_ZERO_MP (f)) {
1261     SET_MP_ZERO (z, digs);
1262     A68_SP = pop_sp;
1263     return z;
1264   }
1265   (void) mp_pi (p, f, MP_PI, digs);
1266   (void) mul_mp (p, g, x, f, digs);
1267   (void) sin_mp (p, z, g, digs);
1268   A68_SP = pop_sp;
1269   return z;
1270 }
1271
1272 //! @brief PROC (LONG REAL) LONG REAL acospi
1273
1274 MP_T *cospi_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1275 {
1277   MP_T *f = nil_mp (p, digs);
1278   MP_T *g = lit_mp (p, 1, 0, digs);
1279   (void) mod_mp (p, f, x, g, digs);
1280   abs_mp (p, f, f, digs);
1281   SET_MP_HALF (g, digs);
1282   (void) sub_mp (p, g, f, g, digs);
1283   if (IS_ZERO_MP (g)) {
1284     SET_MP_ZERO (z, digs);
1285     A68_SP = pop_sp;
1286     return z;
1287   }
1288   (void) mp_pi (p, f, MP_PI, digs);
1289   (void) mul_mp (p, g, x, f, digs);
1290   (void) cos_mp (p, z, g, digs);
1291   A68_SP = pop_sp;
1292   return z;
1293 }
1294
1295 //! @brief PROC (LONG REAL) LONG REAL tanpi
1296
1297 MP_T *tanpi_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1298 {
1300   MP_T *f = nil_mp (p, digs);
1301   MP_T *g = lit_mp (p, 1, 0, digs);
1302   MP_T *h = nil_mp (p, digs);
1303   MP_T *half = nil_mp (p, digs);
1304   SET_MP_ONE (g, digs);
1305   (void) mod_mp (p, f, x, g, digs);
1306   if (IS_ZERO_MP (f)) {
1307     SET_MP_ZERO (z, digs);
1308     A68_SP = pop_sp;
1309     return z;
1310   }
1311   SET_MP_MINUS_HALF (half, digs);
1312   (void) sub_mp (p, h, f, half, digs);
1313   if (MP_DIGIT (h, 1) < 0) {
1314     (void) plus_one_mp (p, f, f, digs);
1315   } else {
1316     SET_MP_HALF (half, digs);
1317     (void) sub_mp (p, h, f, half, digs);
1318     if (MP_DIGIT (h, 1) < 0) {
1319       (void) minus_one_mp (p, f, f, digs);
1320     }
1321   }
1322   BOOL_T neg = MP_DIGIT (f, 1) < 0;
1323   (void) abs_mp (p, f, f, digs);
1324   SET_MP_QUART (g, digs);
1325   (void) sub_mp (p, g, f, g, digs);
1326   if (IS_ZERO_MP (g)) {
1327     if (neg) {
1328       SET_MP_MINUS_ONE (z, digs);
1329     } else {
1330       SET_MP_ONE (z, digs);
1331     }
1332     A68_SP = pop_sp;
1333     return z;
1334   }
1335   (void) mp_pi (p, f, MP_PI, digs);
1336   (void) mul_mp (p, g, x, f, digs);
1337   if (tan_mp (p, z, g, digs) == NaN_MP) {
1338     errno = EDOM;
1339     A68_SP = pop_sp;
1340     return NaN_MP;
1341   } else {
1342     A68_SP = pop_sp;
1343     return z;
1344   }
1345 }
1346
1347 //! @brief PROC (LONG REAL) LONG REAL cotpi
1348
1349 MP_T *cotpi_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1350 {
1352   MP_T *f = nil_mp (p, digs);
1353   MP_T *g = lit_mp (p, 1, 0, digs);
1354   MP_T *h = nil_mp (p, digs);
1355   MP_T *half = nil_mp (p, digs);
1356   (void) mod_mp (p, f, x, g, digs);
1357   if (IS_ZERO_MP (f)) {
1358     errno = EDOM;
1359     A68_SP = pop_sp;
1360     return NaN_MP;
1361   }
1362   SET_MP_MINUS_HALF (half, digs);
1363   (void) sub_mp (p, h, f, half, digs);
1364   if (MP_DIGIT (h, 1) < 0) {
1365     (void) plus_one_mp (p, f, f, digs);
1366   } else {
1367     SET_MP_HALF (half, digs);
1368     (void) sub_mp (p, h, f, half, digs);
1369     if (MP_DIGIT (h, 1) < 0) {
1370       (void) minus_one_mp (p, f, f, digs);
1371     }
1372   }
1373   BOOL_T neg = MP_DIGIT (f, 1) < 0;
1374   (void) abs_mp (p, f, f, digs);
1375   SET_MP_QUART (g, digs);
1376   (void) sub_mp (p, g, f, g, digs);
1377   if (IS_ZERO_MP (g)) {
1378     if (neg) {
1379       SET_MP_MINUS_ONE (z, digs);
1380     } else {
1381       SET_MP_ONE (z, digs);
1382     }
1383     A68_SP = pop_sp;
1384     return z;
1385   }
1386   (void) mp_pi (p, f, MP_PI, digs);
1387   (void) mul_mp (p, g, x, f, digs);
1388   if (cot_mp (p, z, g, digs) == NaN_MP) {
1389     errno = EDOM;
1390     A68_SP = pop_sp;
1391     return NaN_MP;
1392   } else {
1393     A68_SP = pop_sp;
1394     return z;
1395   }
1396 }
1397
1398 //! @brief LONG COMPLEX multiplication
1399
1400 MP_T *cmul_mp (NODE_T * p, MP_T * a, MP_T * b, MP_T * c, MP_T * d, int digs)
1401 {
1403   int gdigs = FUN_DIGITS (digs);
1404   MP_T *la = len_mp (p, a, digs, gdigs);
1405   MP_T *lb = len_mp (p, b, digs, gdigs);
1406   MP_T *lc = len_mp (p, c, digs, gdigs);
1407   MP_T *ld = len_mp (p, d, digs, gdigs);
1408   MP_T *ac = nil_mp (p, gdigs);
1409   MP_T *bd = nil_mp (p, gdigs);
1410   MP_T *ad = nil_mp (p, gdigs);
1411   MP_T *bc = nil_mp (p, gdigs);
1412   (void) mul_mp (p, ac, la, lc, gdigs);
1413   (void) mul_mp (p, bd, lb, ld, gdigs);
1414   (void) mul_mp (p, ad, la, ld, gdigs);
1415   (void) mul_mp (p, bc, lb, lc, gdigs);
1416   (void) sub_mp (p, la, ac, bd, gdigs);
1418   (void) shorten_mp (p, a, digs, la, gdigs);
1419   (void) shorten_mp (p, b, digs, lb, gdigs);
1420   A68_SP = pop_sp;
1421   return a;
1422 }
1423
1424 //! @brief LONG COMPLEX division
1425
1426 MP_T *cdiv_mp (NODE_T * p, MP_T * a, MP_T * b, MP_T * c, MP_T * d, int digs)
1427 {
1429   if (MP_DIGIT (c, 1) == (MP_T) 0 && MP_DIGIT (d, 1) == (MP_T) 0) {
1430     errno = ERANGE;
1431     return NaN_MP;
1432   }
1433   MP_T *q = nil_mp (p, digs);
1434   MP_T *r = nil_mp (p, digs);
1435   (void) move_mp (q, c, digs);
1436   (void) move_mp (r, d, digs);
1437   MP_DIGIT (q, 1) = ABS (MP_DIGIT (q, 1));
1438   MP_DIGIT (r, 1) = ABS (MP_DIGIT (r, 1));
1439   (void) sub_mp (p, q, q, r, digs);
1440   if (MP_DIGIT (q, 1) >= 0) {
1441     if (div_mp (p, q, d, c, digs) == NaN_MP) {
1442       errno = ERANGE;
1443       return NaN_MP;
1444     }
1445     (void) mul_mp (p, r, d, q, digs);
1446     (void) add_mp (p, r, r, c, digs);
1447     (void) mul_mp (p, c, b, q, digs);
1448     (void) add_mp (p, c, c, a, digs);
1449     (void) div_mp (p, c, c, r, digs);
1450     (void) mul_mp (p, d, a, q, digs);
1451     (void) sub_mp (p, d, b, d, digs);
1452     (void) div_mp (p, d, d, r, digs);
1453   } else {
1454     if (div_mp (p, q, c, d, digs) == NaN_MP) {
1455       errno = ERANGE;
1456       return NaN_MP;
1457     }
1458     (void) mul_mp (p, r, c, q, digs);
1459     (void) add_mp (p, r, r, d, digs);
1460     (void) mul_mp (p, c, a, q, digs);
1461     (void) add_mp (p, c, c, b, digs);
1462     (void) div_mp (p, c, c, r, digs);
1463     (void) mul_mp (p, d, b, q, digs);
1464     (void) sub_mp (p, d, d, a, digs);
1465     (void) div_mp (p, d, d, r, digs);
1466   }
1467   (void) move_mp (a, c, digs);
1468   (void) move_mp (b, d, digs);
1469   A68_SP = pop_sp;
1470   return a;
1471 }
1472
1473 //! @brief PROC (LONG COMPLEX) LONG COMPLEX csqrt
1474
1475 MP_T *csqrt_mp (NODE_T * p, MP_T * r, MP_T * i, int digs)
1476 {
1478   int gdigs = FUN_DIGITS (digs);
1479   MP_T *re = len_mp (p, r, digs, gdigs);
1480   MP_T *im = len_mp (p, i, digs, gdigs);
1481   if (IS_ZERO_MP (re) && IS_ZERO_MP (im)) {
1482     SET_MP_ZERO (re, gdigs);
1483     SET_MP_ZERO (im, gdigs);
1484   } else {
1485     MP_T *c1 = lit_mp (p, 1, 0, gdigs);
1486     MP_T *t = nil_mp (p, gdigs);
1487     MP_T *x = nil_mp (p, gdigs);
1488     MP_T *y = nil_mp (p, gdigs);
1489     MP_T *u = nil_mp (p, gdigs);
1490     MP_T *v = nil_mp (p, gdigs);
1491     MP_T *w = nil_mp (p, gdigs);
1492     SET_MP_ONE (c1, gdigs);
1493     (void) move_mp (x, re, gdigs);
1494     (void) move_mp (y, im, gdigs);
1495     MP_DIGIT (x, 1) = ABS (MP_DIGIT (x, 1));
1496     MP_DIGIT (y, 1) = ABS (MP_DIGIT (y, 1));
1497     (void) sub_mp (p, w, x, y, gdigs);
1498     if (MP_DIGIT (w, 1) >= 0) {
1499       (void) div_mp (p, t, y, x, gdigs);
1500       (void) mul_mp (p, v, t, t, gdigs);
1501       (void) add_mp (p, u, c1, v, gdigs);
1502       (void) sqrt_mp (p, v, u, gdigs);
1503       (void) add_mp (p, u, c1, v, gdigs);
1504       (void) half_mp (p, v, u, gdigs);
1505       (void) sqrt_mp (p, u, v, gdigs);
1506       (void) sqrt_mp (p, v, x, gdigs);
1507       (void) mul_mp (p, w, u, v, gdigs);
1508     } else {
1509       (void) div_mp (p, t, x, y, gdigs);
1510       (void) mul_mp (p, v, t, t, gdigs);
1511       (void) add_mp (p, u, c1, v, gdigs);
1512       (void) sqrt_mp (p, v, u, gdigs);
1513       (void) add_mp (p, u, t, v, gdigs);
1514       (void) half_mp (p, v, u, gdigs);
1515       (void) sqrt_mp (p, u, v, gdigs);
1516       (void) sqrt_mp (p, v, y, gdigs);
1517       (void) mul_mp (p, w, u, v, gdigs);
1518     }
1519     if (MP_DIGIT (re, 1) >= 0) {
1520       (void) move_mp (re, w, gdigs);
1521       (void) add_mp (p, u, w, w, gdigs);
1522       (void) div_mp (p, im, im, u, gdigs);
1523     } else {
1524       if (MP_DIGIT (im, 1) < 0) {
1525         MP_DIGIT (w, 1) = -MP_DIGIT (w, 1);
1526       }
1527       (void) add_mp (p, v, w, w, gdigs);
1528       (void) div_mp (p, re, im, v, gdigs);
1529       (void) move_mp (im, w, gdigs);
1530     }
1531   }
1532   (void) shorten_mp (p, r, digs, re, gdigs);
1533   (void) shorten_mp (p, i, digs, im, gdigs);
1534   A68_SP = pop_sp;
1535   return r;
1536 }
1537
1538 //! @brief PROC (LONG COMPLEX) LONG COMPLEX cexp
1539
1540 MP_T *cexp_mp (NODE_T * p, MP_T * r, MP_T * i, int digs)
1541 {
1543   int gdigs = FUN_DIGITS (digs);
1544   MP_T *re = len_mp (p, r, digs, gdigs);
1545   MP_T *im = len_mp (p, i, digs, gdigs);
1546   MP_T *u = nil_mp (p, gdigs);
1547   (void) exp_mp (p, u, re, gdigs);
1548   (void) cos_mp (p, re, im, gdigs);
1549   (void) sin_mp (p, im, im, gdigs);
1550   (void) mul_mp (p, re, re, u, gdigs);
1551   (void) mul_mp (p, im, im, u, gdigs);
1552   (void) shorten_mp (p, r, digs, re, gdigs);
1553   (void) shorten_mp (p, i, digs, im, gdigs);
1554   A68_SP = pop_sp;
1555   return r;
1556 }
1557
1558 //! @brief PROC (LONG COMPLEX) LONG COMPLEX cln
1559
1560 MP_T *cln_mp (NODE_T * p, MP_T * r, MP_T * i, int digs)
1561 {
1563   int gdigs = FUN_DIGITS (digs);
1564   MP_T *re = len_mp (p, r, digs, gdigs);
1565   MP_T *im = len_mp (p, i, digs, gdigs);
1566   MP_T *s = nil_mp (p, gdigs);
1567   MP_T *t = nil_mp (p, gdigs);
1568   MP_T *u = nil_mp (p, gdigs);
1569   MP_T *v = nil_mp (p, gdigs);
1570   (void) move_mp (u, re, gdigs);
1571   (void) move_mp (v, im, gdigs);
1572   (void) hypot_mp (p, s, u, v, gdigs);
1573   (void) move_mp (u, re, gdigs);
1574   (void) move_mp (v, im, gdigs);
1575   (void) atan2_mp (p, t, u, v, gdigs);
1576   (void) ln_mp (p, re, s, gdigs);
1577   (void) move_mp (im, t, gdigs);
1578   (void) shorten_mp (p, r, digs, re, gdigs);
1579   (void) shorten_mp (p, i, digs, im, gdigs);
1580   A68_SP = pop_sp;
1581   return r;
1582 }
1583
1584 //! @brief PROC (LONG COMPLEX) LONG COMPLEX csin
1585
1586 MP_T *csin_mp (NODE_T * p, MP_T * r, MP_T * i, int digs)
1587 {
1589   int gdigs = FUN_DIGITS (digs);
1590   MP_T *re = len_mp (p, r, digs, gdigs);
1591   MP_T *im = len_mp (p, i, digs, gdigs);
1592   MP_T *s = nil_mp (p, gdigs);
1593   MP_T *c = nil_mp (p, gdigs);
1594   MP_T *sh = nil_mp (p, gdigs);
1595   MP_T *ch = nil_mp (p, gdigs);
1596   if (IS_ZERO_MP (im)) {
1597     (void) sin_mp (p, re, re, gdigs);
1598     SET_MP_ZERO (im, gdigs);
1599   } else {
1600     (void) sin_mp (p, s, re, gdigs);
1601     (void) cos_mp (p, c, re, gdigs);
1602     (void) hyp_mp (p, sh, ch, im, gdigs);
1603     (void) mul_mp (p, re, s, ch, gdigs);
1604     (void) mul_mp (p, im, c, sh, gdigs);
1605   }
1606   (void) shorten_mp (p, r, digs, re, gdigs);
1607   (void) shorten_mp (p, i, digs, im, gdigs);
1608   A68_SP = pop_sp;
1609   return r;
1610 }
1611
1612 //! @brief PROC (LONG COMPLEX) LONG COMPLEX ccos
1613
1614 MP_T *ccos_mp (NODE_T * p, MP_T * r, MP_T * i, int digs)
1615 {
1617   int gdigs = FUN_DIGITS (digs);
1618   MP_T *re = len_mp (p, r, digs, gdigs);
1619   MP_T *im = len_mp (p, i, digs, gdigs);
1620   MP_T *s = nil_mp (p, gdigs);
1621   MP_T *c = nil_mp (p, gdigs);
1622   MP_T *sh = nil_mp (p, gdigs);
1623   MP_T *ch = nil_mp (p, gdigs);
1624   if (IS_ZERO_MP (im)) {
1625     (void) cos_mp (p, re, re, gdigs);
1626     SET_MP_ZERO (im, gdigs);
1627   } else {
1628     (void) sin_mp (p, s, re, gdigs);
1629     (void) cos_mp (p, c, re, gdigs);
1630     (void) hyp_mp (p, sh, ch, im, gdigs);
1631     MP_DIGIT (sh, 1) = -MP_DIGIT (sh, 1);
1632     (void) mul_mp (p, re, c, ch, gdigs);
1633     (void) mul_mp (p, im, s, sh, gdigs);
1634   }
1635   (void) shorten_mp (p, r, digs, re, gdigs);
1636   (void) shorten_mp (p, i, digs, im, gdigs);
1637   A68_SP = pop_sp;
1638   return r;
1639 }
1640
1641 //! @brief PROC (LONG COMPLEX) LONG COMPLEX ctan
1642
1643 MP_T *ctan_mp (NODE_T * p, MP_T * r, MP_T * i, int digs)
1644 {
1646   errno = 0;
1647   MP_T *s = nil_mp (p, digs);
1648   MP_T *t = nil_mp (p, digs);
1649   MP_T *u = nil_mp (p, digs);
1650   MP_T *v = nil_mp (p, digs);
1651   (void) move_mp (u, r, digs);
1652   (void) move_mp (v, i, digs);
1653   (void) csin_mp (p, u, v, digs);
1654   (void) move_mp (s, u, digs);
1655   (void) move_mp (t, v, digs);
1656   (void) move_mp (u, r, digs);
1657   (void) move_mp (v, i, digs);
1658   (void) ccos_mp (p, u, v, digs);
1659   (void) cdiv_mp (p, s, t, u, v, digs);
1660   (void) move_mp (r, s, digs);
1661   (void) move_mp (i, t, digs);
1662   A68_SP = pop_sp;
1663   return r;
1664 }
1665
1666 //! @brief PROC (LONG COMPLEX) LONG COMPLEX casin
1667
1668 MP_T *casin_mp (NODE_T * p, MP_T * r, MP_T * i, int digs)
1669 {
1671   int gdigs = FUN_DIGITS (digs);
1672   MP_T *re = len_mp (p, r, digs, gdigs);
1673   MP_T *im = len_mp (p, i, digs, gdigs);
1674   if (IS_ZERO_MP (im)) {
1675     BOOL_T neg = MP_DIGIT (re, 1) < 0;
1676     if (acos_mp (p, im, re, gdigs) == NaN_MP) {
1677       errno = 0;                // Ignore the acos error
1678       MP_DIGIT (re, 1) = ABS (MP_DIGIT (re, 1));
1679       (void) acosh_mp (p, im, re, gdigs);
1680     }
1681     (void) mp_pi (p, re, MP_HALF_PI, gdigs);
1682     if (neg) {
1683       MP_DIGIT (re, 1) = -MP_DIGIT (re, 1);
1684     }
1685   } else {
1686     MP_T *c1 = lit_mp (p, 1, 0, gdigs);
1687     MP_T *u = nil_mp (p, gdigs);
1688     MP_T *v = nil_mp (p, gdigs);
1689     MP_T *a = nil_mp (p, gdigs);
1690     MP_T *b = nil_mp (p, gdigs);
1691 // u=sqrt((r+1)^2+i^2), v=sqrt((r-1)^2+i^2).
1692     (void) add_mp (p, a, re, c1, gdigs);
1693     (void) sub_mp (p, b, re, c1, gdigs);
1694     (void) hypot_mp (p, u, a, im, gdigs);
1695     (void) hypot_mp (p, v, b, im, gdigs);
1696 // a=(u+v)/2, b=(u-v)/2.
1697     (void) add_mp (p, a, u, v, gdigs);
1698     (void) half_mp (p, a, a, gdigs);
1699     (void) sub_mp (p, b, u, v, gdigs);
1700     (void) half_mp (p, b, b, gdigs);
1701 // r=asin(b), i=ln(a+sqrt(a^2-1)).
1702     (void) mul_mp (p, u, a, a, gdigs);
1703     (void) sub_mp (p, u, u, c1, gdigs);
1704     (void) sqrt_mp (p, u, u, gdigs);
1705     (void) add_mp (p, u, a, u, gdigs);
1706     (void) ln_mp (p, im, u, gdigs);
1707     (void) asin_mp (p, re, b, gdigs);
1708   }
1709   (void) shorten_mp (p, r, digs, re, gdigs);
1710   (void) shorten_mp (p, i, digs, im, gdigs);
1711   A68_SP = pop_sp;
1712   return re;
1713 }
1714
1715 //! @brief PROC (LONG COMPLEX) LONG COMPLEX cacos
1716
1717 MP_T *cacos_mp (NODE_T * p, MP_T * r, MP_T * i, int digs)
1718 {
1720   int gdigs = FUN_DIGITS (digs);
1721   MP_T *re = len_mp (p, r, digs, gdigs);
1722   MP_T *im = len_mp (p, i, digs, gdigs);
1723   if (IS_ZERO_MP (im)) {
1724     BOOL_T neg = MP_DIGIT (re, 1) < 0;
1725     if (acos_mp (p, im, re, gdigs) == NaN_MP) {
1726       errno = 0;                // Ignore the acos error
1727       MP_DIGIT (re, 1) = ABS (MP_DIGIT (re, 1));
1728       (void) acosh_mp (p, im, re, gdigs);
1729       MP_DIGIT (im, 1) = -MP_DIGIT (im, 1);
1730     }
1731     if (neg) {
1732       (void) mp_pi (p, re, MP_PI, gdigs);
1733     } else {
1734       SET_MP_ZERO (re, gdigs);
1735     }
1736   } else {
1737     MP_T *c1 = lit_mp (p, 1, 0, gdigs);
1738     MP_T *u = nil_mp (p, gdigs);
1739     MP_T *v = nil_mp (p, gdigs);
1740     MP_T *a = nil_mp (p, gdigs);
1741     MP_T *b = nil_mp (p, gdigs);
1742 // u=sqrt((r+1)^2+i^2), v=sqrt((r-1)^2+i^2).
1743     (void) add_mp (p, a, re, c1, gdigs);
1744     (void) sub_mp (p, b, re, c1, gdigs);
1745     (void) hypot_mp (p, u, a, im, gdigs);
1746     (void) hypot_mp (p, v, b, im, gdigs);
1747 // a=(u+v)/2, b=(u-v)/2.
1748     (void) add_mp (p, a, u, v, gdigs);
1749     (void) half_mp (p, a, a, gdigs);
1750     (void) sub_mp (p, b, u, v, gdigs);
1751     (void) half_mp (p, b, b, gdigs);
1752 // r=acos(b), i=-ln(a+sqrt(a^2-1)).
1753     (void) mul_mp (p, u, a, a, gdigs);
1754     (void) sub_mp (p, u, u, c1, gdigs);
1755     (void) sqrt_mp (p, u, u, gdigs);
1756     (void) add_mp (p, u, a, u, gdigs);
1757     (void) ln_mp (p, im, u, gdigs);
1758     MP_DIGIT (im, 1) = -MP_DIGIT (im, 1);
1759     (void) acos_mp (p, re, b, gdigs);
1760   }
1761   (void) shorten_mp (p, r, digs, re, gdigs);
1762   (void) shorten_mp (p, i, digs, im, gdigs);
1763   A68_SP = pop_sp;
1764   return re;
1765 }
1766
1767 //! @brief PROC (LONG COMPLEX) LONG COMPLEX catan
1768
1769 MP_T *catan_mp (NODE_T * p, MP_T * r, MP_T * i, int digs)
1770 {
1772   int gdigs = FUN_DIGITS (digs);
1773   MP_T *re = len_mp (p, r, digs, gdigs);
1774   MP_T *im = len_mp (p, i, digs, gdigs);
1775   MP_T *u = nil_mp (p, gdigs);
1776   MP_T *v = nil_mp (p, gdigs);
1777   if (IS_ZERO_MP (im)) {
1778     (void) atan_mp (p, u, re, gdigs);
1779     SET_MP_ZERO (v, gdigs);
1780   } else {
1781     MP_T *c1 = lit_mp (p, 1, 0, gdigs);
1782     MP_T *a = nil_mp (p, gdigs);
1783     MP_T *b = nil_mp (p, gdigs);
1784 // a=sqrt(r^2+(i+1)^2), b=sqrt(r^2+(i-1)^2).
1785     (void) add_mp (p, a, im, c1, gdigs);
1786     (void) sub_mp (p, b, im, c1, gdigs);
1787     (void) hypot_mp (p, u, re, a, gdigs);
1788     (void) hypot_mp (p, v, re, b, gdigs);
1789 // im=ln(a/b)/4.
1790     (void) div_mp (p, u, u, v, gdigs);
1791     (void) ln_mp (p, u, u, gdigs);
1792     (void) half_mp (p, v, u, gdigs);
1793 // re=atan(2r/(1-r^2-i^2)).
1794     (void) mul_mp (p, a, re, re, gdigs);
1795     (void) mul_mp (p, b, im, im, gdigs);
1796     (void) sub_mp (p, u, c1, a, gdigs);
1797     (void) sub_mp (p, b, u, b, gdigs);
1798     (void) add_mp (p, a, re, re, gdigs);
1799     (void) div_mp (p, a, a, b, gdigs);
1800     (void) atan_mp (p, u, a, gdigs);
1801     (void) half_mp (p, u, u, gdigs);
1802   }
1803   (void) shorten_mp (p, r, digs, u, gdigs);
1804   (void) shorten_mp (p, i, digs, v, gdigs);
1805   A68_SP = pop_sp;
1806   return re;
1807 }
1808
1809 //! @brief PROC (LONG COMPLEX) LONG COMPLEX csinh
1810
1811 MP_T *csinh_mp (NODE_T * p, MP_T * r, MP_T * i, int digs)
1812 {
1814   int gdigs = FUN_DIGITS (digs);
1815   MP_T *u = len_mp (p, r, digs, gdigs);
1816   MP_T *v = len_mp (p, i, digs, gdigs);
1817   MP_T *s = nil_mp (p, gdigs);
1818   MP_T *t = nil_mp (p, gdigs);
1819 // sinh (z) =  -i sin (iz)
1820   SET_MP_ONE (t, gdigs);
1821   (void) cmul_mp (p, u, v, s, t, gdigs);
1822   (void) csin_mp (p, u, v, gdigs);
1823   SET_MP_ZERO (s, gdigs);
1824   SET_MP_MINUS_ONE (t, gdigs);
1825   (void) cmul_mp (p, u, v, s, t, gdigs);
1826 //
1827   (void) shorten_mp (p, r, digs, u, gdigs);
1828   (void) shorten_mp (p, i, digs, v, gdigs);
1829   A68_SP = pop_sp;
1830   return r;
1831 }
1832
1833 //! @brief PROC (LONG COMPLEX) LONG COMPLEX ccosh
1834
1835 MP_T *ccosh_mp (NODE_T * p, MP_T * r, MP_T * i, int digs)
1836 {
1838   int gdigs = FUN_DIGITS (digs);
1839   MP_T *u = len_mp (p, r, digs, gdigs);
1840   MP_T *v = len_mp (p, i, digs, gdigs);
1841   MP_T *s = nil_mp (p, gdigs);
1842   MP_T *t = nil_mp (p, gdigs);
1843 // cosh (z) =  cos (iz)
1844   SET_MP_ZERO (s, digs);
1845   SET_MP_ONE (t, gdigs);
1846   (void) cmul_mp (p, u, v, s, t, gdigs);
1847   (void) ccos_mp (p, u, v, gdigs);
1848 //
1849   (void) shorten_mp (p, r, digs, u, gdigs);
1850   (void) shorten_mp (p, i, digs, v, gdigs);
1851   A68_SP = pop_sp;
1852   return r;
1853 }
1854
1855 //! @brief PROC (LONG COMPLEX) LONG COMPLEX ctanh
1856
1857 MP_T *ctanh_mp (NODE_T * p, MP_T * r, MP_T * i, int digs)
1858 {
1860   int gdigs = FUN_DIGITS (digs);
1861   MP_T *u = len_mp (p, r, digs, gdigs);
1862   MP_T *v = len_mp (p, i, digs, gdigs);
1863   MP_T *s = nil_mp (p, gdigs);
1864   MP_T *t = nil_mp (p, gdigs);
1865 // tanh (z) =  -i tan (iz)
1866   SET_MP_ONE (t, gdigs);
1867   (void) cmul_mp (p, u, v, s, t, gdigs);
1868   (void) ctan_mp (p, u, v, gdigs);
1869   SET_MP_ZERO (u, gdigs);
1870   SET_MP_MINUS_ONE (v, gdigs);
1871   (void) cmul_mp (p, u, v, s, t, gdigs);
1872 //
1873   (void) shorten_mp (p, r, digs, u, gdigs);
1874   (void) shorten_mp (p, i, digs, v, gdigs);
1875   A68_SP = pop_sp;
1876   return r;
1877 }
1878
1879 //! @brief PROC (LONG COMPLEX) LONG COMPLEX casinh
1880
1881 MP_T *casinh_mp (NODE_T * p, MP_T * r, MP_T * i, int digs)
1882 {
1884   int gdigs = FUN_DIGITS (digs);
1885   MP_T *u = len_mp (p, r, digs, gdigs);
1886   MP_T *v = len_mp (p, i, digs, gdigs);
1887   MP_T *s = nil_mp (p, gdigs);
1888   MP_T *t = nil_mp (p, gdigs);
1889 // asinh (z) =  i asin (-iz)
1890   SET_MP_MINUS_ONE (t, gdigs);
1891   (void) cmul_mp (p, u, v, s, t, gdigs);
1892   (void) casin_mp (p, u, v, gdigs);
1893   SET_MP_ZERO (s, gdigs);
1894   SET_MP_ONE (t, gdigs);
1895   (void) cmul_mp (p, u, v, s, t, gdigs);
1896 //
1897   (void) shorten_mp (p, r, digs, u, gdigs);
1898   (void) shorten_mp (p, i, digs, v, gdigs);
1899   A68_SP = pop_sp;
1900   return r;
1901 }
1902
1903 //! @brief PROC (LONG COMPLEX) LONG COMPLEX cacosh
1904
1905 MP_T *cacosh_mp (NODE_T * p, MP_T * r, MP_T * i, int digs)
1906 {
1908   int gdigs = FUN_DIGITS (digs);
1909   MP_T *u = len_mp (p, r, digs, gdigs);
1910   MP_T *v = len_mp (p, i, digs, gdigs);
1911   MP_T *s = nil_mp (p, gdigs);
1912   MP_T *t = nil_mp (p, gdigs);
1913 // acosh (z) =  i * acos (z)
1914   (void) cacos_mp (p, u, v, gdigs);
1915   SET_MP_ONE (t, gdigs);
1916   (void) cmul_mp (p, u, v, s, t, gdigs);
1917 //
1918   (void) shorten_mp (p, r, digs, u, gdigs);
1919   (void) shorten_mp (p, i, digs, v, gdigs);
1920   A68_SP = pop_sp;
1921   return r;
1922 }
1923
1924 //! @brief PROC (LONG COMPLEX) LONG COMPLEX catanh
1925
1926 MP_T *catanh_mp (NODE_T * p, MP_T * r, MP_T * i, int digs)
1927 {
1929   int gdigs = FUN_DIGITS (digs);
1930   MP_T *u = len_mp (p, r, digs, gdigs);
1931   MP_T *v = len_mp (p, i, digs, gdigs);
1932   MP_T *s = nil_mp (p, gdigs);
1933   MP_T *t = nil_mp (p, gdigs);
1934 // atanh (z) =  i * atan (-iz)
1935   SET_MP_MINUS_ONE (t, gdigs);
1936   (void) cmul_mp (p, u, v, s, t, gdigs);
1937   (void) catan_mp (p, u, v, gdigs);
1938   SET_MP_ZERO (s, gdigs);
1939   SET_MP_ONE (t, gdigs);
1940   (void) cmul_mp (p, u, v, s, t, gdigs);
1941 //
1942   (void) shorten_mp (p, r, digs, u, gdigs);
1943   (void) shorten_mp (p, i, digs, v, gdigs);
1944   A68_SP = pop_sp;
1945   return r;
1946 }
```

 © 2023 J.M. van der Veer jmvdveer@xs4all.nl