mp-math.c

     
   1  //! @file mp-math.c
   2  //! @author J. Marcel van der Veer
   3  //!
   4  //! @section Copyright
   5  //!
   6  //! This file is part of Algol68G - an Algol 68 compiler-interpreter.
   7  //! Copyright 2001-2023 J. Marcel van der Veer [algol68g@xs4all.nl].
   8  //!
   9  //! @section License
  10  //!
  11  //! This program is free software; you can redistribute it and/or modify it 
  12  //! under the terms of the GNU General Public License as published by the 
  13  //! Free Software Foundation; either version 3 of the License, or 
  14  //! (at your option) any later version.
  15  //!
  16  //! This program is distributed in the hope that it will be useful, but 
  17  //! WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 
  18  //! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for 
  19  //! more details. You should have received a copy of the GNU General Public 
  20  //! License along with this program. If not, see [http://www.gnu.org/licenses/].
  21  
  22  //! @section Synopsis
  23  //!
  24  //! [LONG] LONG REAL math functions.
  25  
  26  #include "a68g.h"
  27  #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  {
  52    ADDR_T pop_sp = A68_SP;
  53    if (MP_DIGIT (x, 1) == 0) {
  54      A68_SP = pop_sp;
  55      SET_MP_ZERO (z, digs);
  56      return z;
  57    }
  58    if (MP_DIGIT (x, 1) < 0) {
  59      A68_SP = pop_sp;
  60      errno = EDOM;
  61      return NaN_MP;
  62    }
  63    int gdigs = FUN_DIGITS (digs), hdigs;
  64    BOOL_T reciprocal = A68_FALSE;
  65    MP_T *z_g = nil_mp (p, gdigs);
  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_double_real (p, x_g, gdigs);
  82      (void) double_real_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<n+1> = (x<n> + a / x<n>) / 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  {
 110    ADDR_T pop_sp = A68_SP;
 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_double_real (p, x_g, gdigs);
 141      (void) double_real_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<n+1> = (2 x<n> + a / x<n> ^ 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).
 175    ADDR_T pop_sp = A68_SP;
 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);
 214    ADDR_T pop_sp = A68_SP;
 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);
 289    ADDR_T pop_sp = A68_SP;
 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  {
 351    ADDR_T pop_sp = A68_SP;
 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  {
 374    ADDR_T pop_sp = A68_SP;
 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.
 398    ADDR_T pop_sp = A68_SP;
 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<n+1> = x<n> - 1 + a / exp (x<n>).
 445      MP_T *tmp = nil_mp (p, gdigs);
 446  // Construct an estimate.
 447  #if (A68_LEVEL >= 3)
 448      (void) double_real_to_mp (p, z_g, logq (mp_to_double_real (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  {
 484    ADDR_T pop_sp = A68_SP;
 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  {
 500    ADDR_T pop_sp = A68_SP;
 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  {
 525    ADDR_T pop_sp = A68_SP;
 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 {
 544      ADDR_T pop_sp = A68_SP;
 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  {
 575    ADDR_T pop_sp = A68_SP;
 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  {
 590    ADDR_T pop_sp = A68_SP;
 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  {
 616    ADDR_T pop_sp = A68_SP;
 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  {
 632    ADDR_T pop_sp = A68_SP;
 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.
 653    ADDR_T pop_sp = A68_SP;
 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.
 752    ADDR_T pop_sp = A68_SP;
 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)).
 773    ADDR_T pop_sp = A68_SP;
 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  {
 815    ADDR_T pop_sp = A68_SP;
 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  {
 848    ADDR_T pop_sp = A68_SP;
 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.
 886    ADDR_T pop_sp = A68_SP;
 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<n+1> = x<n> - cos (x<n>) * (sin (x<n>) - a cos (x<n>)).
 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) double_real_to_mp (p, z_g, atanq (mp_to_double_real (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  {
 975    ADDR_T pop_sp = A68_SP;
 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)).
1051    ADDR_T pop_sp = A68_SP;
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  {
1093    ADDR_T pop_sp = A68_SP;
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  {
1110    ADDR_T pop_sp = A68_SP;
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  {
1124    ADDR_T pop_sp = A68_SP;
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  {
1138    ADDR_T pop_sp = A68_SP;
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  {
1157    ADDR_T pop_sp = A68_SP;
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  {
1176    ADDR_T pop_sp = A68_SP;
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  {
1190    ADDR_T pop_sp = A68_SP;
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  {
1204    ADDR_T pop_sp = A68_SP;
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  {
1218    ADDR_T pop_sp = A68_SP;
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  {
1237    ADDR_T pop_sp = A68_SP;
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  {
1256    ADDR_T pop_sp = A68_SP;
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  {
1276    ADDR_T pop_sp = A68_SP;
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  {
1299    ADDR_T pop_sp = A68_SP;
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  {
1351    ADDR_T pop_sp = A68_SP;
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  }