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  }