mp-complex.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 COMPLEX math functions.
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 LONG COMPLEX multiplication
33
34  MP_T *cmul_mp (NODE_T * p, MP_T * a, MP_T * b, MP_T * c, MP_T * d, int digs)
35  {
36    ADDR_T pop_sp = A68_SP;
37    int gdigs = FUN_DIGITS (digs);
38    MP_T *la = len_mp (p, a, digs, gdigs);
39    MP_T *lb = len_mp (p, b, digs, gdigs);
40    MP_T *lc = len_mp (p, c, digs, gdigs);
41    MP_T *ld = len_mp (p, d, digs, gdigs);
42    MP_T *ac = nil_mp (p, gdigs);
43    MP_T *bd = nil_mp (p, gdigs);
44    MP_T *ad = nil_mp (p, gdigs);
45    MP_T *bc = nil_mp (p, gdigs);
46    (void) mul_mp (p, ac, la, lc, gdigs);
47    (void) mul_mp (p, bd, lb, ld, gdigs);
48    (void) mul_mp (p, ad, la, ld, gdigs);
49    (void) mul_mp (p, bc, lb, lc, gdigs);
50    (void) sub_mp (p, la, ac, bd, gdigs);
51    (void) add_mp (p, lb, ad, bc, gdigs);
52    (void) shorten_mp (p, a, digs, la, gdigs);
53    (void) shorten_mp (p, b, digs, lb, gdigs);
54    A68_SP = pop_sp;
55    return a;
56  }
57
58  //! @brief LONG COMPLEX division
59
60  MP_T *cdiv_mp (NODE_T * p, MP_T * a, MP_T * b, MP_T * c, MP_T * d, int digs)
61  {
62    ADDR_T pop_sp = A68_SP;
63    if (MP_DIGIT (c, 1) == (MP_T) 0 && MP_DIGIT (d, 1) == (MP_T) 0) {
64      errno = ERANGE;
65      return NaN_MP;
66    }
67    MP_T *q = nil_mp (p, digs);
68    MP_T *r = nil_mp (p, digs);
69    (void) move_mp (q, c, digs);
70    (void) move_mp (r, d, digs);
71    MP_DIGIT (q, 1) = ABS (MP_DIGIT (q, 1));
72    MP_DIGIT (r, 1) = ABS (MP_DIGIT (r, 1));
73    (void) sub_mp (p, q, q, r, digs);
74    if (MP_DIGIT (q, 1) >= 0) {
75      if (div_mp (p, q, d, c, digs) == NaN_MP) {
76        errno = ERANGE;
77        return NaN_MP;
78      }
79      (void) mul_mp (p, r, d, q, digs);
80      (void) add_mp (p, r, r, c, digs);
81      (void) mul_mp (p, c, b, q, digs);
82      (void) add_mp (p, c, c, a, digs);
83      (void) div_mp (p, c, c, r, digs);
84      (void) mul_mp (p, d, a, q, digs);
85      (void) sub_mp (p, d, b, d, digs);
86      (void) div_mp (p, d, d, r, digs);
87    } else {
88      if (div_mp (p, q, c, d, digs) == NaN_MP) {
89        errno = ERANGE;
90        return NaN_MP;
91      }
92      (void) mul_mp (p, r, c, q, digs);
93      (void) add_mp (p, r, r, d, digs);
94      (void) mul_mp (p, c, a, q, digs);
95      (void) add_mp (p, c, c, b, digs);
96      (void) div_mp (p, c, c, r, digs);
97      (void) mul_mp (p, d, b, q, digs);
98      (void) sub_mp (p, d, d, a, digs);
99      (void) div_mp (p, d, d, r, digs);
100    }
101    (void) move_mp (a, c, digs);
102    (void) move_mp (b, d, digs);
103    A68_SP = pop_sp;
104    return a;
105  }
106
107  //! @brief PROC (LONG COMPLEX) LONG COMPLEX csqrt
108
109  MP_T *csqrt_mp (NODE_T * p, MP_T * r, MP_T * i, int digs)
110  {
111    ADDR_T pop_sp = A68_SP;
112    int gdigs = FUN_DIGITS (digs);
113    MP_T *re = len_mp (p, r, digs, gdigs);
114    MP_T *im = len_mp (p, i, digs, gdigs);
115    if (IS_ZERO_MP (re) && IS_ZERO_MP (im)) {
116      SET_MP_ZERO (re, gdigs);
117      SET_MP_ZERO (im, gdigs);
118    } else {
119      MP_T *c1 = lit_mp (p, 1, 0, gdigs);
120      MP_T *t = nil_mp (p, gdigs);
121      MP_T *x = nil_mp (p, gdigs);
122      MP_T *y = nil_mp (p, gdigs);
123      MP_T *u = nil_mp (p, gdigs);
124      MP_T *v = nil_mp (p, gdigs);
125      MP_T *w = nil_mp (p, gdigs);
126      SET_MP_ONE (c1, gdigs);
127      (void) move_mp (x, re, gdigs);
128      (void) move_mp (y, im, gdigs);
129      MP_DIGIT (x, 1) = ABS (MP_DIGIT (x, 1));
130      MP_DIGIT (y, 1) = ABS (MP_DIGIT (y, 1));
131      (void) sub_mp (p, w, x, y, gdigs);
132      if (MP_DIGIT (w, 1) >= 0) {
133        (void) div_mp (p, t, y, x, gdigs);
134        (void) mul_mp (p, v, t, t, gdigs);
135        (void) add_mp (p, u, c1, v, gdigs);
136        (void) sqrt_mp (p, v, u, gdigs);
137        (void) add_mp (p, u, c1, v, gdigs);
138        (void) half_mp (p, v, u, gdigs);
139        (void) sqrt_mp (p, u, v, gdigs);
140        (void) sqrt_mp (p, v, x, gdigs);
141        (void) mul_mp (p, w, u, v, gdigs);
142      } else {
143        (void) div_mp (p, t, x, y, gdigs);
144        (void) mul_mp (p, v, t, t, gdigs);
145        (void) add_mp (p, u, c1, v, gdigs);
146        (void) sqrt_mp (p, v, u, gdigs);
147        (void) add_mp (p, u, t, v, gdigs);
148        (void) half_mp (p, v, u, gdigs);
149        (void) sqrt_mp (p, u, v, gdigs);
150        (void) sqrt_mp (p, v, y, gdigs);
151        (void) mul_mp (p, w, u, v, gdigs);
152      }
153      if (MP_DIGIT (re, 1) >= 0) {
154        (void) move_mp (re, w, gdigs);
155        (void) add_mp (p, u, w, w, gdigs);
156        (void) div_mp (p, im, im, u, gdigs);
157      } else {
158        if (MP_DIGIT (im, 1) < 0) {
159          MP_DIGIT (w, 1) = -MP_DIGIT (w, 1);
160        }
161        (void) add_mp (p, v, w, w, gdigs);
162        (void) div_mp (p, re, im, v, gdigs);
163        (void) move_mp (im, w, gdigs);
164      }
165    }
166    (void) shorten_mp (p, r, digs, re, gdigs);
167    (void) shorten_mp (p, i, digs, im, gdigs);
168    A68_SP = pop_sp;
169    return r;
170  }
171
172  //! @brief PROC (LONG COMPLEX) LONG COMPLEX cexp
173
174  MP_T *cexp_mp (NODE_T * p, MP_T * r, MP_T * i, int digs)
175  {
176    ADDR_T pop_sp = A68_SP;
177    int gdigs = FUN_DIGITS (digs);
178    MP_T *re = len_mp (p, r, digs, gdigs);
179    MP_T *im = len_mp (p, i, digs, gdigs);
180    MP_T *u = nil_mp (p, gdigs);
181    (void) exp_mp (p, u, re, gdigs);
182    (void) cos_mp (p, re, im, gdigs);
183    (void) sin_mp (p, im, im, gdigs);
184    (void) mul_mp (p, re, re, u, gdigs);
185    (void) mul_mp (p, im, im, u, gdigs);
186    (void) shorten_mp (p, r, digs, re, gdigs);
187    (void) shorten_mp (p, i, digs, im, gdigs);
188    A68_SP = pop_sp;
189    return r;
190  }
191
192  //! @brief PROC (LONG COMPLEX) LONG COMPLEX cln
193
194  MP_T *cln_mp (NODE_T * p, MP_T * r, MP_T * i, int digs)
195  {
196    ADDR_T pop_sp = A68_SP;
197    int gdigs = FUN_DIGITS (digs);
198    MP_T *re = len_mp (p, r, digs, gdigs);
199    MP_T *im = len_mp (p, i, digs, gdigs);
200    MP_T *s = nil_mp (p, gdigs);
201    MP_T *t = nil_mp (p, gdigs);
202    MP_T *u = nil_mp (p, gdigs);
203    MP_T *v = nil_mp (p, gdigs);
204    (void) move_mp (u, re, gdigs);
205    (void) move_mp (v, im, gdigs);
206    (void) hypot_mp (p, s, u, v, gdigs);
207    (void) move_mp (u, re, gdigs);
208    (void) move_mp (v, im, gdigs);
209    (void) atan2_mp (p, t, u, v, gdigs);
210    (void) ln_mp (p, re, s, gdigs);
211    (void) move_mp (im, t, gdigs);
212    (void) shorten_mp (p, r, digs, re, gdigs);
213    (void) shorten_mp (p, i, digs, im, gdigs);
214    A68_SP = pop_sp;
215    return r;
216  }
217
218  //! @brief PROC (LONG COMPLEX) LONG COMPLEX csin
219
220  MP_T *csin_mp (NODE_T * p, MP_T * r, MP_T * i, int digs)
221  {
222    ADDR_T pop_sp = A68_SP;
223    int gdigs = FUN_DIGITS (digs);
224    MP_T *re = len_mp (p, r, digs, gdigs);
225    MP_T *im = len_mp (p, i, digs, gdigs);
226    MP_T *s = nil_mp (p, gdigs);
227    MP_T *c = nil_mp (p, gdigs);
228    MP_T *sh = nil_mp (p, gdigs);
229    MP_T *ch = nil_mp (p, gdigs);
230    if (IS_ZERO_MP (im)) {
231      (void) sin_mp (p, re, re, gdigs);
232      SET_MP_ZERO (im, gdigs);
233    } else {
234      (void) sin_mp (p, s, re, gdigs);
235      (void) cos_mp (p, c, re, gdigs);
236      (void) hyp_mp (p, sh, ch, im, gdigs);
237      (void) mul_mp (p, re, s, ch, gdigs);
238      (void) mul_mp (p, im, c, sh, gdigs);
239    }
240    (void) shorten_mp (p, r, digs, re, gdigs);
241    (void) shorten_mp (p, i, digs, im, gdigs);
242    A68_SP = pop_sp;
243    return r;
244  }
245
246  //! @brief PROC (LONG COMPLEX) LONG COMPLEX ccos
247
248  MP_T *ccos_mp (NODE_T * p, MP_T * r, MP_T * i, int digs)
249  {
250    ADDR_T pop_sp = A68_SP;
251    int gdigs = FUN_DIGITS (digs);
252    MP_T *re = len_mp (p, r, digs, gdigs);
253    MP_T *im = len_mp (p, i, digs, gdigs);
254    MP_T *s = nil_mp (p, gdigs);
255    MP_T *c = nil_mp (p, gdigs);
256    MP_T *sh = nil_mp (p, gdigs);
257    MP_T *ch = nil_mp (p, gdigs);
258    if (IS_ZERO_MP (im)) {
259      (void) cos_mp (p, re, re, gdigs);
260      SET_MP_ZERO (im, gdigs);
261    } else {
262      (void) sin_mp (p, s, re, gdigs);
263      (void) cos_mp (p, c, re, gdigs);
264      (void) hyp_mp (p, sh, ch, im, gdigs);
265      MP_DIGIT (sh, 1) = -MP_DIGIT (sh, 1);
266      (void) mul_mp (p, re, c, ch, gdigs);
267      (void) mul_mp (p, im, s, sh, gdigs);
268    }
269    (void) shorten_mp (p, r, digs, re, gdigs);
270    (void) shorten_mp (p, i, digs, im, gdigs);
271    A68_SP = pop_sp;
272    return r;
273  }
274
275  //! @brief PROC (LONG COMPLEX) LONG COMPLEX ctan
276
277  MP_T *ctan_mp (NODE_T * p, MP_T * r, MP_T * i, int digs)
278  {
279    ADDR_T pop_sp = A68_SP;
280    errno = 0;
281    MP_T *s = nil_mp (p, digs);
282    MP_T *t = nil_mp (p, digs);
283    MP_T *u = nil_mp (p, digs);
284    MP_T *v = nil_mp (p, digs);
285    (void) move_mp (u, r, digs);
286    (void) move_mp (v, i, digs);
287    (void) csin_mp (p, u, v, digs);
288    (void) move_mp (s, u, digs);
289    (void) move_mp (t, v, digs);
290    (void) move_mp (u, r, digs);
291    (void) move_mp (v, i, digs);
292    (void) ccos_mp (p, u, v, digs);
293    (void) cdiv_mp (p, s, t, u, v, digs);
294    (void) move_mp (r, s, digs);
295    (void) move_mp (i, t, digs);
296    A68_SP = pop_sp;
297    return r;
298  }
299
300  //! @brief PROC (LONG COMPLEX) LONG COMPLEX casin
301
302  MP_T *casin_mp (NODE_T * p, MP_T * r, MP_T * i, int digs)
303  {
304    ADDR_T pop_sp = A68_SP;
305    int gdigs = FUN_DIGITS (digs);
306    MP_T *re = len_mp (p, r, digs, gdigs);
307    MP_T *im = len_mp (p, i, digs, gdigs);
308    BOOL_T negim = MP_DIGIT (im, 1) < 0;
309    if (IS_ZERO_MP (im)) {
310      BOOL_T neg = MP_DIGIT (re, 1) < 0;
311      if (acos_mp (p, im, re, gdigs) == NaN_MP) {
312        errno = 0;                // Ignore the acos error
313        MP_DIGIT (re, 1) = ABS (MP_DIGIT (re, 1));
314        (void) acosh_mp (p, im, re, gdigs);
315      }
316      (void) mp_pi (p, re, MP_HALF_PI, gdigs);
317      if (neg) {
318        MP_DIGIT (re, 1) = -MP_DIGIT (re, 1);
319      }
320    } else {
321      MP_T *c1 = lit_mp (p, 1, 0, gdigs);
322      MP_T *u = nil_mp (p, gdigs);
323      MP_T *v = nil_mp (p, gdigs);
324      MP_T *a = nil_mp (p, gdigs);
325      MP_T *b = nil_mp (p, gdigs);
326  // u=sqrt((r+1)^2+i^2), v=sqrt((r-1)^2+i^2).
327      (void) add_mp (p, a, re, c1, gdigs);
328      (void) sub_mp (p, b, re, c1, gdigs);
329      (void) hypot_mp (p, u, a, im, gdigs);
330      (void) hypot_mp (p, v, b, im, gdigs);
331  // a=(u+v)/2, b=(u-v)/2.
332      (void) add_mp (p, a, u, v, gdigs);
333      (void) half_mp (p, a, a, gdigs);
334      (void) sub_mp (p, b, u, v, gdigs);
335      (void) half_mp (p, b, b, gdigs);
336  // r=asin(b), i=ln(a+sqrt(a^2-1)).
337      (void) mul_mp (p, u, a, a, gdigs);
338      (void) sub_mp (p, u, u, c1, gdigs);
339      (void) sqrt_mp (p, u, u, gdigs);
340      (void) add_mp (p, u, a, u, gdigs);
341      (void) ln_mp (p, im, u, gdigs);
342      (void) asin_mp (p, re, b, gdigs);
343    }
344    if (negim) {
345      MP_DIGIT (im, 1) = -MP_DIGIT (im, 1);
346    }
347    (void) shorten_mp (p, r, digs, re, gdigs);
348    (void) shorten_mp (p, i, digs, im, gdigs);
349    A68_SP = pop_sp;
350    return re;
351  }
352
353  //! @brief PROC (LONG COMPLEX) LONG COMPLEX cacos
354
355  MP_T *cacos_mp (NODE_T * p, MP_T * r, MP_T * i, int digs)
356  {
357    ADDR_T pop_sp = A68_SP;
358    int gdigs = FUN_DIGITS (digs);
359    MP_T *re = len_mp (p, r, digs, gdigs);
360    MP_T *im = len_mp (p, i, digs, gdigs);
361    BOOL_T negim = MP_DIGIT (im, 1) < 0;
362    if (IS_ZERO_MP (im)) {
363      BOOL_T neg = MP_DIGIT (re, 1) < 0;
364      if (acos_mp (p, im, re, gdigs) == NaN_MP) {
365        errno = 0;                // Ignore the acos error
366        MP_DIGIT (re, 1) = ABS (MP_DIGIT (re, 1));
367        (void) acosh_mp (p, im, re, gdigs);
368        MP_DIGIT (im, 1) = -MP_DIGIT (im, 1);
369      }
370      if (neg) {
371        (void) mp_pi (p, re, MP_PI, gdigs);
372      } else {
373        SET_MP_ZERO (re, gdigs);
374      }
375    } else {
376      MP_T *c1 = lit_mp (p, 1, 0, gdigs);
377      MP_T *u = nil_mp (p, gdigs);
378      MP_T *v = nil_mp (p, gdigs);
379      MP_T *a = nil_mp (p, gdigs);
380      MP_T *b = nil_mp (p, gdigs);
381  // u=sqrt((r+1)^2+i^2), v=sqrt((r-1)^2+i^2).
382      (void) add_mp (p, a, re, c1, gdigs);
383      (void) sub_mp (p, b, re, c1, gdigs);
384      (void) hypot_mp (p, u, a, im, gdigs);
385      (void) hypot_mp (p, v, b, im, gdigs);
386  // a=(u+v)/2, b=(u-v)/2.
387      (void) add_mp (p, a, u, v, gdigs);
388      (void) half_mp (p, a, a, gdigs);
389      (void) sub_mp (p, b, u, v, gdigs);
390      (void) half_mp (p, b, b, gdigs);
391  // r=acos(b), i=-ln(a+sqrt(a^2-1)).
392      (void) mul_mp (p, u, a, a, gdigs);
393      (void) sub_mp (p, u, u, c1, gdigs);
394      (void) sqrt_mp (p, u, u, gdigs);
395      (void) add_mp (p, u, a, u, gdigs);
396      (void) ln_mp (p, im, u, gdigs);
397      (void) acos_mp (p, re, b, gdigs);
398    }
399    if (!negim) {
400      MP_DIGIT (im, 1) = -MP_DIGIT (im, 1);
401    }
402    (void) shorten_mp (p, r, digs, re, gdigs);
403    (void) shorten_mp (p, i, digs, im, gdigs);
404    A68_SP = pop_sp;
405    return re;
406  }
407
408  //! @brief PROC (LONG COMPLEX) LONG COMPLEX catan
409
410  MP_T *catan_mp (NODE_T * p, MP_T * r, MP_T * i, int digs)
411  {
412    ADDR_T pop_sp = A68_SP;
413    int gdigs = FUN_DIGITS (digs);
414    MP_T *re = len_mp (p, r, digs, gdigs);
415    MP_T *im = len_mp (p, i, digs, gdigs);
416    MP_T *u = nil_mp (p, gdigs);
417    MP_T *v = nil_mp (p, gdigs);
418    if (IS_ZERO_MP (im)) {
419      (void) atan_mp (p, u, re, gdigs);
420      SET_MP_ZERO (v, gdigs);
421    } else {
422      MP_T *c1 = lit_mp (p, 1, 0, gdigs);
423      MP_T *a = nil_mp (p, gdigs);
424      MP_T *b = nil_mp (p, gdigs);
425  // a=sqrt(r^2+(i+1)^2), b=sqrt(r^2+(i-1)^2).
426      (void) add_mp (p, a, im, c1, gdigs);
427      (void) sub_mp (p, b, im, c1, gdigs);
428      (void) hypot_mp (p, u, re, a, gdigs);
429      (void) hypot_mp (p, v, re, b, gdigs);
430  // im=ln(a/b).
431      (void) div_mp (p, u, u, v, gdigs);
432      (void) ln_mp (p, v, u, gdigs);
433  // re=atan(2r/(1-r^2-i^2)).
434      (void) mul_mp (p, a, re, re, gdigs);
435      (void) mul_mp (p, b, im, im, gdigs);
436      (void) add_mp (p, a, a, b, gdigs);
437      (void) sub_mp (p, u, c1, a, gdigs);
438      if (IS_ZERO_MP (u)) {
439        (void) mp_pi (p, u, MP_HALF_PI, gdigs);
440      } else {
441        int neg = MP_DIGIT (u, 1) < 0;
442        (void) add_mp (p, a, re, re, gdigs);
443        (void) div_mp (p, a, a, u, gdigs);
444        (void) atan_mp (p, u, a, gdigs);
445        if (neg) {
446          (void) mp_pi (p, a, MP_PI, gdigs);
447          if (MP_DIGIT (re, 1) < 0) {
448            (void) sub_mp (p, u, u, a, gdigs);
449          } else {
450            (void) add_mp (p, u, u, a, gdigs);
451          }
452        }
453      }
454    }
455    (void) shorten_mp (p, r, digs, u, gdigs);
456    (void) shorten_mp (p, i, digs, v, gdigs);
457    (void) half_mp (p, r, r, digs);
458    (void) half_mp (p, i, i, digs);
459    A68_SP = pop_sp;
460    return re;
461  }
462
463  //! @brief PROC (LONG COMPLEX) LONG COMPLEX csinh
464
465  MP_T *csinh_mp (NODE_T * p, MP_T * r, MP_T * i, int digs)
466  {
467    ADDR_T pop_sp = A68_SP;
468    int gdigs = FUN_DIGITS (digs);
469    MP_T *u = len_mp (p, r, digs, gdigs);
470    MP_T *v = len_mp (p, i, digs, gdigs);
471    MP_T *s = nil_mp (p, gdigs);
472    MP_T *t = nil_mp (p, gdigs);
473  // sinh (z) =  -i sin (iz)
474    SET_MP_ONE (t, gdigs);
475    (void) cmul_mp (p, u, v, s, t, gdigs);
476    (void) csin_mp (p, u, v, gdigs);
477    SET_MP_ZERO (s, gdigs);
478    SET_MP_MINUS_ONE (t, gdigs);
479    (void) cmul_mp (p, u, v, s, t, gdigs);
480  //
481    (void) shorten_mp (p, r, digs, u, gdigs);
482    (void) shorten_mp (p, i, digs, v, gdigs);
483    A68_SP = pop_sp;
484    return r;
485  }
486
487  //! @brief PROC (LONG COMPLEX) LONG COMPLEX ccosh
488
489  MP_T *ccosh_mp (NODE_T * p, MP_T * r, MP_T * i, int digs)
490  {
491    ADDR_T pop_sp = A68_SP;
492    int gdigs = FUN_DIGITS (digs);
493    MP_T *u = len_mp (p, r, digs, gdigs);
494    MP_T *v = len_mp (p, i, digs, gdigs);
495    MP_T *s = nil_mp (p, gdigs);
496    MP_T *t = nil_mp (p, gdigs);
497  // cosh (z) =  cos (iz)
498    SET_MP_ZERO (s, digs);
499    SET_MP_ONE (t, gdigs);
500    (void) cmul_mp (p, u, v, s, t, gdigs);
501    (void) ccos_mp (p, u, v, gdigs);
502  //
503    (void) shorten_mp (p, r, digs, u, gdigs);
504    (void) shorten_mp (p, i, digs, v, gdigs);
505    A68_SP = pop_sp;
506    return r;
507  }
508
509  //! @brief PROC (LONG COMPLEX) LONG COMPLEX ctanh
510
511  MP_T *ctanh_mp (NODE_T * p, MP_T * r, MP_T * i, int digs)
512  {
513    ADDR_T pop_sp = A68_SP;
514    int gdigs = FUN_DIGITS (digs);
515    MP_T *u = len_mp (p, r, digs, gdigs);
516    MP_T *v = len_mp (p, i, digs, gdigs);
517    MP_T *s = nil_mp (p, gdigs);
518    MP_T *t = nil_mp (p, gdigs);
519  // tanh (z) =  -i tan (iz)
520    SET_MP_ONE (t, gdigs);
521    (void) cmul_mp (p, u, v, s, t, gdigs);
522    (void) ctan_mp (p, u, v, gdigs);
523    SET_MP_ZERO (u, gdigs);
524    SET_MP_MINUS_ONE (v, gdigs);
525    (void) cmul_mp (p, u, v, s, t, gdigs);
526  //
527    (void) shorten_mp (p, r, digs, u, gdigs);
528    (void) shorten_mp (p, i, digs, v, gdigs);
529    A68_SP = pop_sp;
530    return r;
531  }
532
533  //! @brief PROC (LONG COMPLEX) LONG COMPLEX casinh
534
535  MP_T *casinh_mp (NODE_T * p, MP_T * r, MP_T * i, int digs)
536  {
537    ADDR_T pop_sp = A68_SP;
538    int gdigs = FUN_DIGITS (digs);
539    MP_T *u = len_mp (p, r, digs, gdigs);
540    MP_T *v = len_mp (p, i, digs, gdigs);
541    MP_T *s = nil_mp (p, gdigs);
542    MP_T *t = nil_mp (p, gdigs);
543  // asinh (z) =  i asin (-iz)
544    SET_MP_MINUS_ONE (t, gdigs);
545    (void) cmul_mp (p, u, v, s, t, gdigs);
546    (void) casin_mp (p, u, v, gdigs);
547    SET_MP_ZERO (s, gdigs);
548    SET_MP_ONE (t, gdigs);
549    (void) cmul_mp (p, u, v, s, t, gdigs);
550  //
551    (void) shorten_mp (p, r, digs, u, gdigs);
552    (void) shorten_mp (p, i, digs, v, gdigs);
553    A68_SP = pop_sp;
554    return r;
555  }
556
557  //! @brief PROC (LONG COMPLEX) LONG COMPLEX cacosh
558
559  MP_T *cacosh_mp (NODE_T * p, MP_T * r, MP_T * i, int digs)
560  {
561    ADDR_T pop_sp = A68_SP;
562    int gdigs = FUN_DIGITS (digs);
563    MP_T *u = len_mp (p, r, digs, gdigs);
564    MP_T *v = len_mp (p, i, digs, gdigs);
565    MP_T *s = nil_mp (p, gdigs);
566    MP_T *t = nil_mp (p, gdigs);
567  // acosh (z) =  i * acos (z)
568    (void) cacos_mp (p, u, v, gdigs);
569    SET_MP_ZERO (s, gdigs);
570    SET_MP_ONE (t, gdigs);
571    (void) cmul_mp (p, u, v, s, t, gdigs);
572  //
573    (void) shorten_mp (p, r, digs, u, gdigs);
574    (void) shorten_mp (p, i, digs, v, gdigs);
575    A68_SP = pop_sp;
576    return r;
577  }
578
579  //! @brief PROC (LONG COMPLEX) LONG COMPLEX catanh
580
581  MP_T *catanh_mp (NODE_T * p, MP_T * r, MP_T * i, int digs)
582  {
583    ADDR_T pop_sp = A68_SP;
584    int gdigs = FUN_DIGITS (digs);
585    MP_T *u = len_mp (p, r, digs, gdigs);
586    MP_T *v = len_mp (p, i, digs, gdigs);
587    MP_T *s = nil_mp (p, gdigs);
588    MP_T *t = nil_mp (p, gdigs);
589  // atanh (z) =  i * atan (-iz)
590    SET_MP_MINUS_ONE (t, gdigs);
591    (void) cmul_mp (p, u, v, s, t, gdigs);
592    (void) catan_mp (p, u, v, gdigs);
593    SET_MP_ZERO (s, gdigs);
594    SET_MP_ONE (t, gdigs);
595    (void) cmul_mp (p, u, v, s, t, gdigs);
596  //
597    (void) shorten_mp (p, r, digs, u, gdigs);
598    (void) shorten_mp (p, i, digs, v, gdigs);
599    A68_SP = pop_sp;
600    return r;
601  }

```