mp.c

     
   1  //! @file mp.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 INT, REAL and COMPLEX routines.
  25  
  26  // Multiprecision calculations are useful in these cases:
  27  //
  28  // Ill-conditioned linear systems
  29  // Summation of large series
  30  // Long-time or large-scale simulations
  31  // Small-scale phenomena
  32  // 'Experimental mathematics'
  33  //
  34  // The routines in this library follow algorithms as described in the
  35  // literature, notably
  36  // 
  37  //   D.M. Smith, "Efficient Multiple-Precision Evaluation of Elementary Functions"
  38  //   Mathematics of Computation 52 (1989) 131-134
  39  // 
  40  //   D.M. Smith, "A Multiple-Precision Division Algorithm"
  41  //   Mathematics of Computation 66 (1996) 157-163
  42  //
  43  //   The GNU MPFR library documentation
  44  // 
  45  // Multiprecision libraries are (freely) available, but this one is particularly 
  46  // designed to work with Algol68G. It implements following modes:
  47  // 
  48  //    LONG INT, LONG REAL, LONG COMPLEX, LONG BITS
  49  //    LONG LONG INT, LONG LONG REAL, LONG LONG COMPLEX, LONG LONG BITS
  50  //
  51  // Note that recent implementations of GCC make available 64-bit LONG INT and 
  52  // 128-bit LONG REAL. This suits many multiprecision needs already. 
  53  // On such platforms, below code is used for LONG LONG modes only.
  54  // Now that 64-bit integers are commonplace, this library has been adapted to
  55  // exploit them. Having some more digits per word gives performance gain.
  56  // This is implemented through macros to keep the library compatible with
  57  // old platforms with 32-bit integers and 64-bit doubles.
  58  // 
  59  // Currently, LONG modes have a fixed precision, and LONG LONG modes have
  60  // user-definable precision. Precisions span about 30 decimal digits for
  61  // LONG modes up to (default) about 60 decimal digits for LONG LONG modes, a
  62  // range that is thought to be adequate for most multiprecision applications.
  63  //
  64  // Although the maximum length of a number is in principle unbound, this 
  65  // implementation is not designed for more than a few hundred decimal places. 
  66  // At higher precisions, expect a performance penalty with respect to
  67  // state of the art implementations that may for instance use convolution for
  68  // multiplication. 
  69  // 
  70  // This library takes a sloppy approach towards LONG INT and LONG BITS which are
  71  // implemented as LONG REAL and truncated where appropriate. This keeps the code
  72  // short at the penalty of some performance loss.
  73  // 
  74  // As is common practice, mp numbers are represented by a row of digits
  75  // in a large base. Layout of a mp number "z" is:
  76  // 
  77  //    MP_T *z;
  78  //
  79  //    MP_STATUS (z)        Status word
  80  //    MP_EXPONENT (z)      Exponent with base MP_RADIX
  81  //    MP_DIGIT (z, 1 .. N) Digits 1 .. N
  82  // 
  83  // Note that this library assumes IEEE 754 compatible implementation of
  84  // type "double". It also assumes a 32- (or 64-) bit type "int".
  85  // 
  86  // Most legacy multiple precision libraries stored numbers as [] int*4.
  87  // However, since division and multiplication are O(N ** 2) operations, it is
  88  // advantageous to keep the base as high as possible. Modern computers handle
  89  // doubles at similar or better speed as integers, therefore this library
  90  // opts for storing numbers as [] words were a word is real*8 (legacy) or 
  91  // int*8 (on f.i. ix86 processors that have real*10), trading space for speed.
  92  // 
  93  // Set a base such that "base^2" can be exactly represented by a word.
  94  // To facilitate transput, we require a base that is a power of 10.
  95  // 
  96  // If we choose the base right then in multiplication and division we do not need
  97  // to normalise intermediate results at each step since a number of additions
  98  // can be made before overflow occurs. That is why we specify "MAX_REPR_INT".
  99  // 
 100  // Mind that the precision of a mp number is at worst just
 101  // (LONG_MP_DIGITS - 1) * LOG_MP_RADIX + 1, since the most significant mp digit
 102  // is also in range [0 .. MP_RADIX>. Do not specify less than 2 digits.
 103  // 
 104  // Since this software is distributed without any warranty, it is your
 105  // responsibility to validate the behaviour of the routines and their accuracy
 106  // using the source code provided. See the GNU General Public License for details.
 107  
 108  #include "a68g.h"
 109  #include "a68g-genie.h"
 110  #include "a68g-prelude.h"
 111  #include "a68g-double.h"
 112  #include "a68g-mp.h"
 113  
 114  // Internal mp constants.
 115  
 116  //! @brief Set number of digits for long long numbers.
 117  
 118  void set_long_mp_digits (int n)
 119  {
 120    A68_MP (varying_mp_digits) = n;
 121  }
 122  
 123  //! @brief Convert precision to digits for long long number.
 124  
 125  int width_to_mp_digits (int n)
 126  {
 127    return 1 + n / LOG_MP_RADIX;
 128  }
 129  
 130  //! @brief Unformatted write of z to stdout; debugging routine.
 131  
 132  #if !defined (BUILD_WIN32)
 133  
 134  void raw_write_mp (char *str, MP_T * z, int digs)
 135  {
 136    fprintf (stdout, "\n(%d digits)%s", digs, str);
 137    for (unt i = 1; i <= digs; i++) {
 138  #if (A68_LEVEL >= 3)
 139      fprintf (stdout, " %09lld", (MP_INT_T) MP_DIGIT (z, i));
 140  #else
 141      fprintf (stdout, " %07d", (MP_INT_T) MP_DIGIT (z, i));
 142  #endif
 143    }
 144    fprintf (stdout, " E" A68_LD, (MP_INT_T) MP_EXPONENT (z));
 145    fprintf (stdout, " S" A68_LD, (MP_INT_T) MP_STATUS (z));
 146    fprintf (stdout, "\n");
 147    ASSERT (fflush (stdout) == 0);
 148  }
 149  
 150  #endif
 151  
 152  //! @brief Whether z is a valid representation for its mode.
 153  
 154  BOOL_T check_mp_int (MP_T * z, MOID_T * m)
 155  {
 156    if (m == M_LONG_INT || m == M_LONG_BITS) {
 157      return (BOOL_T) ((MP_EXPONENT (z) >= (MP_T) 0) && (MP_EXPONENT (z) < (MP_T) LONG_MP_DIGITS));
 158    } else if (m == M_LONG_LONG_INT || m == M_LONG_LONG_BITS) {
 159      return (BOOL_T) ((MP_EXPONENT (z) >= (MP_T) 0) && (MP_EXPONENT (z) < (MP_T) A68_MP (varying_mp_digits)));
 160    }
 161    return A68_FALSE;
 162  }
 163  
 164  //! @brief |x|
 165  
 166  MP_T *abs_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
 167  {
 168    (void) p;
 169    if (x != z) {
 170      (void) move_mp (z, x, digs);
 171    }
 172    MP_DIGIT (z, 1) = ABS (MP_DIGIT (z, 1));
 173    MP_STATUS (z) = (MP_T) INIT_MASK;
 174    return z;
 175  }
 176  
 177  //! @brief -x
 178  
 179  MP_T *minus_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
 180  {
 181    (void) p;
 182    if (x != z) {
 183      (void) move_mp (z, x, digs);
 184    }
 185    MP_DIGIT (z, 1) = -(MP_DIGIT (z, 1));
 186    MP_STATUS (z) = (MP_T) INIT_MASK;
 187    return z;
 188  }
 189  
 190  //! @brief 1 - x
 191  
 192  MP_T *one_minus_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
 193  {
 194    ADDR_T pop_sp = A68_SP;
 195    (void) sub_mp (p, z, mp_one (digs), x, digs);
 196    A68_SP = pop_sp;
 197    return z;
 198  }
 199  
 200  //! @brief x - 1
 201  
 202  MP_T *minus_one_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
 203  {
 204    ADDR_T pop_sp = A68_SP;
 205    (void) sub_mp (p, z, x, mp_one (digs), digs);
 206    A68_SP = pop_sp;
 207    return z;
 208  }
 209  
 210  //! @brief x + 1
 211  
 212  MP_T *plus_one_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
 213  {
 214    ADDR_T pop_sp = A68_SP;
 215    (void) add_mp (p, z, x, mp_one (digs), digs);
 216    A68_SP = pop_sp;
 217    return z;
 218  }
 219  
 220  //! @brief Test whether x = y.
 221  
 222  BOOL_T same_mp (NODE_T * p, MP_T * x, MP_T * y, int digs)
 223  {
 224    (void) p;
 225    if ((MP_STATUS (x) == MP_STATUS (y)) && (MP_EXPONENT (x) == MP_EXPONENT (y))) {
 226      for (unt k = digs; k >= 1; k--) {
 227        if (MP_DIGIT (x, k) != MP_DIGIT (y, k)) {
 228          return A68_FALSE;
 229        }
 230      }
 231      return A68_TRUE;
 232    } else {
 233      return A68_FALSE;
 234    }
 235  }
 236  
 237  //! @brief Align 10-base z in a MP_RADIX mantissa.
 238  
 239  MP_T *align_mp (MP_T * z, INT_T * expo, int digs)
 240  {
 241    INT_T shift;
 242    if (*expo >= 0) {
 243      shift = LOG_MP_RADIX - (*expo) % LOG_MP_RADIX - 1;
 244      (*expo) /= LOG_MP_RADIX;
 245    } else {
 246      shift = (-(*expo) - 1) % LOG_MP_RADIX;
 247      (*expo) = ((*expo) + 1) / LOG_MP_RADIX;
 248      (*expo)--;
 249    }
 250  // Optimising below code does not make the library noticeably faster.
 251    for (unt i = 1; i <= shift; i++) {
 252      INT_T carry = 0;
 253      for (INT_T j = 1; j <= digs; j++) {
 254        MP_INT_T k = ((MP_INT_T) MP_DIGIT (z, j)) % 10;
 255        MP_DIGIT (z, j) = (MP_T) ((MP_INT_T) (MP_DIGIT (z, j) / 10) + carry * (MP_RADIX / 10));
 256        carry = k;
 257      }
 258    }
 259    return z;
 260  }
 261  
 262  //! @brief Transform string into multi-precision number.
 263  
 264  MP_T *strtomp (NODE_T * p, MP_T * z, char *s, int digs)
 265  {
 266    BOOL_T ok = A68_TRUE;
 267    errno = 0;
 268    SET_MP_ZERO (z, digs);
 269    while (IS_SPACE (s[0])) {
 270      s++;
 271    }
 272  // Get the sign.
 273    int sign = (s[0] == '-' ? -1 : 1);
 274    if (s[0] == '+' || s[0] == '-') {
 275      s++;
 276    }
 277  // Scan mantissa digs and put them into "z".
 278    while (s[0] == '0') {
 279      s++;
 280    }
 281    int i = 0, dig = 1;
 282    INT_T sum = 0, dot = -1, one = -1, pow = 0, W = MP_RADIX / 10;
 283    while (s[i] != NULL_CHAR && dig <= digs && (IS_DIGIT (s[i]) || s[i] == POINT_CHAR)) {
 284      if (s[i] == POINT_CHAR) {
 285        dot = i;
 286      } else {
 287        int value = (int) s[i] - (int) '0';
 288        if (one < 0 && value > 0) {
 289          one = pow;
 290        }
 291        sum += W * value;
 292        if (one >= 0) {
 293          W /= 10;
 294        }
 295        pow++;
 296        if (W < 1) {
 297          MP_DIGIT (z, dig++) = (MP_T) sum;
 298          sum = 0;
 299          W = MP_RADIX / 10;
 300        }
 301      }
 302      i++;
 303    }
 304  // Store the last digs.
 305    if (dig <= digs) {
 306      MP_DIGIT (z, dig++) = (MP_T) sum;
 307    }
 308  // See if there is an exponent.
 309    INT_T expo;
 310    if (s[i] != NULL_CHAR && TO_UPPER (s[i]) == TO_UPPER (EXPONENT_CHAR)) {
 311      char *end;
 312      expo = (int) strtol (&(s[++i]), &end, 10);
 313      ok = (BOOL_T) (end[0] == NULL_CHAR);
 314    } else {
 315      expo = 0;
 316      ok = (BOOL_T) (s[i] == NULL_CHAR);
 317    }
 318  // Calculate effective exponent.
 319    if (dot >= 0) {
 320      if (one > dot) {
 321        expo -= one - dot + 1;
 322      } else {
 323        expo += dot - 1;
 324      }
 325    } else {
 326      expo += pow - 1;
 327    }
 328    (void) align_mp (z, &expo, digs);
 329    MP_EXPONENT (z) = (MP_DIGIT (z, 1) == 0 ? 0 : (MP_T) expo);
 330    MP_DIGIT (z, 1) *= sign;
 331    check_mp_exp (p, z);
 332    if (errno == 0 && ok) {
 333      return z;
 334    } else {
 335      return NaN_MP;
 336    }
 337  }
 338  
 339  //! @brief Convert integer to multi-precison number.
 340  
 341  MP_T *int_to_mp (NODE_T * p, MP_T * z, INT_T k, int digs)
 342  {
 343    int sign_k = 1;
 344    if (k < 0) {
 345      k = -k;
 346      sign_k = -1;
 347    }
 348    int m = k, n = 0;
 349    while ((m /= MP_RADIX) != 0) {
 350      n++;
 351    }
 352    set_mp (z, 0, n, digs);
 353    for (unt j = 1 + n; j >= 1; j--) {
 354      MP_DIGIT (z, j) = (MP_T) (k % MP_RADIX);
 355      k /= MP_RADIX;
 356    }
 357    MP_DIGIT (z, 1) = sign_k * MP_DIGIT (z, 1);
 358    check_mp_exp (p, z);
 359    return z;
 360  }
 361  
 362  //! @brief Convert unt to multi-precison number.
 363  
 364  MP_T *unt_to_mp (NODE_T * p, MP_T * z, UNSIGNED_T k, int digs)
 365  {
 366    int m = k, n = 0;
 367    while ((m /= MP_RADIX) != 0) {
 368      n++;
 369    }
 370    set_mp (z, 0, n, digs);
 371    for (unt j = 1 + n; j >= 1; j--) {
 372      MP_DIGIT (z, j) = (MP_T) (k % MP_RADIX);
 373      k /= MP_RADIX;
 374    }
 375    check_mp_exp (p, z);
 376    return z;
 377  }
 378  
 379  //! @brief Convert multi-precision number to integer.
 380  
 381  INT_T mp_to_int (NODE_T * p, MP_T * z, int digs)
 382  {
 383  // This routines looks a lot like "strtol". 
 384    INT_T expo = (int) MP_EXPONENT (z), sum = 0, weight = 1;
 385    if (expo >= digs) {
 386      diagnostic (A68_RUNTIME_ERROR, p, ERROR_OUT_OF_BOUNDS, MOID (p));
 387      exit_genie (p, A68_RUNTIME_ERROR);
 388    }
 389    BOOL_T negative = (BOOL_T) (MP_DIGIT (z, 1) < 0);
 390    if (negative) {
 391      MP_DIGIT (z, 1) = -MP_DIGIT (z, 1);
 392    }
 393    for (int j = 1 + expo; j >= 1; j--) {
 394      if ((MP_INT_T) MP_DIGIT (z, j) > A68_MAX_INT / weight) {
 395        diagnostic (A68_RUNTIME_ERROR, p, ERROR_OUT_OF_BOUNDS, M_INT);
 396        exit_genie (p, A68_RUNTIME_ERROR);
 397      }
 398      INT_T term = (MP_INT_T) MP_DIGIT (z, j) * weight;
 399      if (sum > A68_MAX_INT - term) {
 400        diagnostic (A68_RUNTIME_ERROR, p, ERROR_OUT_OF_BOUNDS, M_INT);
 401        exit_genie (p, A68_RUNTIME_ERROR);
 402      }
 403      sum += term;
 404      weight *= MP_RADIX;
 405    }
 406    return negative ? -sum : sum;
 407  }
 408  
 409  //! @brief Convert REAL_T to multi-precison number.
 410  
 411  MP_T *real_to_mp (NODE_T * p, MP_T * z, REAL_T x, int digs)
 412  {
 413    SET_MP_ZERO (z, digs);
 414    if (x == 0.0) {
 415      return z;
 416    }
 417  // Small integers can be done better by int_to_mp.
 418    if (ABS (x) < MP_RADIX && trunc (x) == x) {
 419      return int_to_mp (p, z, (INT_T) trunc (x), digs);
 420    }
 421    int sign_x = SIGN (x);
 422  // Scale to [0, 0.1>.
 423    REAL_T a = ABS (x);
 424    INT_T expo = (int) log10 (a);
 425    a /= ten_up (expo);
 426    expo--;
 427    if (a >= 1) {
 428      a /= 10;
 429      expo++;
 430    }
 431  // Transport digs of x to the mantissa of z.
 432    INT_T sum = 0, weight = (MP_RADIX / 10);
 433    int j = 1;
 434    for (unt k = 0; a != 0.0 && j <= digs && k < REAL_DIGITS; k++) {
 435      REAL_T u = a * 10;
 436      REAL_T v = floor (u);
 437      a = u - v;
 438      sum += weight * (INT_T) v;
 439      weight /= 10;
 440      if (weight < 1) {
 441        MP_DIGIT (z, j++) = (MP_T) sum;
 442        sum = 0;
 443        weight = (MP_RADIX / 10);
 444      }
 445    }
 446  // Store the last digs.
 447    if (j <= digs) {
 448      MP_DIGIT (z, j) = (MP_T) sum;
 449    }
 450    (void) align_mp (z, &expo, digs);
 451    MP_EXPONENT (z) = (MP_T) expo;
 452    MP_DIGIT (z, 1) *= sign_x;
 453    check_mp_exp (p, z);
 454    return z;
 455  }
 456  
 457  //! @brief Convert multi-precision number to real.
 458  
 459  REAL_T mp_to_real (NODE_T * p, MP_T * z, int digs)
 460  {
 461  // This routine looks a lot like "strtod".
 462    (void) p;
 463    if (MP_EXPONENT (z) * (MP_T) LOG_MP_RADIX <= (MP_T) REAL_MIN_10_EXP) {
 464      return 0;
 465    } else {
 466      REAL_T terms[1 + MP_MAX_DIGITS];
 467      REAL_T weight = ten_up ((int) (MP_EXPONENT (z) * LOG_MP_RADIX));
 468      unt lim = MIN (digs, MP_MAX_DIGITS);
 469      for (unt k = 1; k <= lim; k++) {
 470        terms[k] = ABS (MP_DIGIT (z, k)) * weight;
 471        weight /= MP_RADIX;
 472      }
 473  // Sum terms from small to large.
 474      REAL_T sum = 0;
 475      for (unt k = lim; k >= 1; k--) {
 476        sum += terms[k];
 477      }
 478      CHECK_REAL (p, sum);
 479      return MP_DIGIT (z, 1) >= 0 ? sum : -sum;
 480    }
 481  }
 482  
 483  //! @brief Normalise positive intermediate, fast.
 484  
 485  static inline void norm_mp_light (MP_T * w, int k, int digs)
 486  {
 487  // Bring every digit back to [0 .. MP_RADIX>.
 488    MP_T *z = &MP_DIGIT (w, digs);
 489    for (unt j = digs; j >= k; j--, z--) {
 490      if (z[0] >= MP_RADIX) {
 491        z[0] -= (MP_T) MP_RADIX;
 492        z[-1] += 1;
 493      } else if (z[0] < 0) {
 494        z[0] += (MP_T) MP_RADIX;
 495        z[-1] -= 1;
 496      }
 497    }
 498  }
 499  
 500  //! @brief Normalise positive intermediate.
 501  
 502  static inline void norm_mp (MP_T * w, int k, int digs)
 503  {
 504  // Bring every digit back to [0 .. MP_RADIX>.
 505    unt j; MP_T *z;
 506    for (j = digs, z = &MP_DIGIT (w, digs); j >= k; j--, z--) {
 507      if (z[0] >= (MP_T) MP_RADIX) {
 508        MP_T carry = (MP_T) ((MP_INT_T) (z[0] / (MP_T) MP_RADIX));
 509        z[0] -= carry * (MP_T) MP_RADIX;
 510        z[-1] += carry;
 511      } else if (z[0] < (MP_T) 0) {
 512        MP_T carry = (MP_T) 1 + (MP_T) ((MP_INT_T) ((-z[0] - 1) / (MP_T) MP_RADIX));
 513        z[0] += carry * (MP_T) MP_RADIX;
 514        z[-1] -= carry;
 515      }
 516    }
 517  }
 518  
 519  //! @brief Round multi-precision number.
 520  
 521  static inline void round_internal_mp (MP_T * z, MP_T * w, int digs)
 522  {
 523  // Assume that w has precision of at least 2 + digs.
 524    int last = (MP_DIGIT (w, 1) == 0 ? 2 + digs : 1 + digs);
 525    if (MP_DIGIT (w, last) >= MP_RADIX / 2) {
 526      MP_DIGIT (w, last - 1) += 1;
 527    }
 528    if (MP_DIGIT (w, last - 1) >= MP_RADIX) {
 529      norm_mp (w, 2, last); // Hardly ever happens - no need to optimise
 530    }
 531    if (MP_DIGIT (w, 1) == 0) {
 532      (void) move_mp_part (&MP_DIGIT (z, 1), &MP_DIGIT (w, 2), digs);
 533      MP_EXPONENT (z) = MP_EXPONENT (w) - 1;
 534    } else if (z != w) {
 535      (void) move_mp_part (&MP_EXPONENT (z), &MP_EXPONENT (w), (1 + digs));
 536    }
 537  // Zero is zero is zero.
 538    if (MP_DIGIT (z, 1) == 0) {
 539      MP_EXPONENT (z) = (MP_T) 0;
 540    }
 541  }
 542  
 543  //! @brief Truncate at decimal point.
 544  
 545  MP_T *trunc_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
 546  {
 547    if (MP_EXPONENT (x) < 0) {
 548      SET_MP_ZERO (z, digs);
 549    } else if (MP_EXPONENT (x) >= (MP_T) digs) {
 550      errno = EDOM;
 551      diagnostic (A68_RUNTIME_ERROR, p, ERROR_OUT_OF_BOUNDS, (IS (MOID (p), PROC_SYMBOL) ? SUB_MOID (p) : MOID (p)));
 552      exit_genie (p, A68_RUNTIME_ERROR);
 553    } else {
 554      (void) move_mp (z, x, digs);
 555      for (int k = (int) (MP_EXPONENT (x) + 2); k <= digs; k++) {
 556        MP_DIGIT (z, k) = (MP_T) 0;
 557      }
 558    }
 559    return z;
 560  }
 561  
 562  //! @brief Floor - largest integer smaller than x.
 563  
 564  MP_T *floor_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
 565  {
 566    if (MP_EXPONENT (x) < 0) {
 567      SET_MP_ZERO (z, digs);
 568    } else if (MP_EXPONENT (x) >= (MP_T) digs) {
 569      errno = EDOM;
 570      diagnostic (A68_RUNTIME_ERROR, p, ERROR_OUT_OF_BOUNDS, (IS (MOID (p), PROC_SYMBOL) ? SUB_MOID (p) : MOID (p)));
 571      exit_genie (p, A68_RUNTIME_ERROR);
 572    } else {
 573      (void) move_mp (z, x, digs);
 574      for (int k = (int) (MP_EXPONENT (x) + 2); k <= digs; k++) {
 575        MP_DIGIT (z, k) = (MP_T) 0;
 576      }
 577    }
 578    if (MP_DIGIT (x, 1) < 0 && ! same_mp (p, z, x, digs)) {
 579      (void) minus_one_mp (p, z, z, digs);
 580    }
 581    return z;
 582  }
 583  
 584  BOOL_T is_int_mp (NODE_T *p, MP_T *z, int digs)
 585  {
 586    ADDR_T pop_sp = A68_SP;
 587    MP_T *y = nil_mp (p, digs);
 588    trunc_mp (p, y, z, digs);
 589    BOOL_T tst = same_mp (p, y, z, digs);
 590    A68_SP = pop_sp;
 591    return tst;
 592  }
 593  
 594  //! @brief Shorten and round.
 595  
 596  MP_T *shorten_mp (NODE_T * p, MP_T * z, int digs, MP_T * x, int digs_x)
 597  {
 598    if (digs > digs_x) {
 599      return lengthen_mp (p, z, digs, x, digs_x);
 600    } else if (digs == digs_x) {
 601      return move_mp (z, x, digs);
 602    } else {
 603  // Reserve extra digs for proper rounding.
 604      ADDR_T pop_sp = A68_SP;
 605      int digs_h = digs + 2;
 606      BOOL_T negative = (BOOL_T) (MP_DIGIT (x, 1) < 0);
 607      MP_T *w = nil_mp (p, digs_h);
 608      if (negative) {
 609        MP_DIGIT (x, 1) = -MP_DIGIT (x, 1);
 610      }
 611      MP_STATUS (w) = (MP_T) 0;
 612      MP_EXPONENT (w) = MP_EXPONENT (x) + 1;
 613      MP_DIGIT (w, 1) = (MP_T) 0;
 614      (void) move_mp_part (&MP_DIGIT (w, 2), &MP_DIGIT (x, 1), digs + 1);
 615      round_internal_mp (z, w, digs);
 616      if (negative) {
 617        MP_DIGIT (z, 1) = -MP_DIGIT (z, 1);
 618      }
 619      A68_SP = pop_sp;
 620      return z;
 621    }
 622  }
 623  
 624  //! @brief Lengthen x and assign to z.
 625  
 626  MP_T *lengthen_mp (NODE_T * p, MP_T * z, int digs_z, MP_T * x, int digs_x)
 627  {
 628    if (digs_z < digs_x) {
 629      return shorten_mp (p, z, digs_z, x, digs_x);
 630    } else if (digs_z == digs_x) {
 631      return move_mp (z, x, digs_z);
 632    } else {
 633      if (z != x) {
 634        (void) move_mp_part (&MP_DIGIT (z, 1), &MP_DIGIT (x, 1), digs_x);
 635        MP_EXPONENT (z) = MP_EXPONENT (x);
 636        MP_STATUS (z) = MP_STATUS (x);
 637      }
 638      for (unt j = 1 + digs_x; j <= digs_z; j++) {
 639        MP_DIGIT (z, j) = (MP_T) 0;
 640      }
 641    }
 642    return z;
 643  }
 644  
 645  //! @brief Set "z" to the sum of positive "x" and positive "y".
 646  
 647  MP_T *add_mp (NODE_T * p, MP_T * z, MP_T * x, MP_T * y, int digs)
 648  {
 649    MP_STATUS (z) = (MP_T) INIT_MASK;
 650  // Trivial cases.
 651    if (MP_DIGIT (x, 1) == (MP_T) 0) {
 652      (void) move_mp (z, y, digs);
 653      return z;
 654    } else if (MP_DIGIT (y, 1) == 0) {
 655      (void) move_mp (z, x, digs);
 656      return z;
 657    }
 658  // We want positive arguments.
 659    ADDR_T pop_sp = A68_SP;
 660    MP_T x_1 = MP_DIGIT (x, 1), y_1 = MP_DIGIT (y, 1);
 661    MP_DIGIT (x, 1) = ABS (x_1);
 662    MP_DIGIT (y, 1) = ABS (y_1);
 663    if (x_1 >= 0 && y_1 < 0) {
 664      (void) sub_mp (p, z, x, y, digs);
 665    } else if (x_1 < 0 && y_1 >= 0) {
 666      (void) sub_mp (p, z, y, x, digs);
 667    } else if (x_1 < 0 && y_1 < 0) {
 668      (void) add_mp (p, z, x, y, digs);
 669      MP_DIGIT (z, 1) = -MP_DIGIT (z, 1);
 670    } else {
 671  // Add.
 672      int digs_h = 2 + digs;
 673      MP_T *w = nil_mp (p, digs_h);
 674      if (MP_EXPONENT (x) == MP_EXPONENT (y)) {
 675        MP_EXPONENT (w) = (MP_T) 1 + MP_EXPONENT (x);
 676        for (unt j = 1; j <= digs; j++) {
 677          MP_DIGIT (w, j + 1) = MP_DIGIT (x, j) + MP_DIGIT (y, j);
 678        }
 679        MP_DIGIT (w, digs_h) = (MP_T) 0;
 680      } else if (MP_EXPONENT (x) > MP_EXPONENT (y)) {
 681        int shl_y = (int) MP_EXPONENT (x) - (int) MP_EXPONENT (y);
 682        MP_EXPONENT (w) = (MP_T) 1 + MP_EXPONENT (x);
 683        for (unt j = 1; j < digs_h; j++) {
 684          int i_y = j - shl_y;
 685          MP_T x_j = (j > digs ? 0 : MP_DIGIT (x, j));
 686          MP_T y_j = (i_y <= 0 || i_y > digs ? 0 : MP_DIGIT (y, i_y));
 687          MP_DIGIT (w, j + 1) = x_j + y_j;
 688        }
 689      } else {
 690        int shl_x = (int) MP_EXPONENT (y) - (int) MP_EXPONENT (x);
 691        MP_EXPONENT (w) = (MP_T) 1 + MP_EXPONENT (y);
 692        for (unt j = 1; j < digs_h; j++) {
 693          int i_x = j - shl_x;
 694          MP_T x_j = (i_x <= 0 || i_x > digs ? 0 : MP_DIGIT (x, i_x));
 695          MP_T y_j = (j > digs ? 0 : MP_DIGIT (y, j));
 696          MP_DIGIT (w, j + 1) = x_j + y_j;
 697        }
 698      }
 699      norm_mp_light (w, 2, digs_h);
 700      round_internal_mp (z, w, digs);
 701      check_mp_exp (p, z);
 702    }
 703  // Restore and exit.
 704    A68_SP = pop_sp;
 705    MP_T z_1 = MP_DIGIT (z, 1);
 706    MP_DIGIT (x, 1) = x_1;
 707    MP_DIGIT (y, 1) = y_1;
 708    MP_DIGIT (z, 1) = z_1;        // In case z IS x OR z IS y
 709    return z;
 710  }
 711  
 712  //! @brief Set "z" to the difference of positive "x" and positive "y".
 713  
 714  MP_T *sub_mp (NODE_T * p, MP_T * z, MP_T * x, MP_T * y, int digs)
 715  {
 716    MP_STATUS (z) = (MP_T) INIT_MASK;
 717  // Trivial cases.
 718    if (MP_DIGIT (x, 1) == (MP_T) 0) {
 719      (void) move_mp (z, y, digs);
 720      MP_DIGIT (z, 1) = -MP_DIGIT (z, 1);
 721      return z;
 722    } else if (MP_DIGIT (y, 1) == (MP_T) 0) {
 723      (void) move_mp (z, x, digs);
 724      return z;
 725    }
 726  // We want positive arguments.
 727    ADDR_T pop_sp = A68_SP;
 728    MP_T x_1 = MP_DIGIT (x, 1), y_1 = MP_DIGIT (y, 1);
 729    MP_DIGIT (x, 1) = ABS (x_1);
 730    MP_DIGIT (y, 1) = ABS (y_1);
 731    if (x_1 >= 0 && y_1 < 0) {
 732      (void) add_mp (p, z, x, y, digs);
 733    } else if (x_1 < 0 && y_1 >= 0) {
 734      (void) add_mp (p, z, y, x, digs);
 735      MP_DIGIT (z, 1) = -MP_DIGIT (z, 1);
 736    } else if (x_1 < 0 && y_1 < 0) {
 737      (void) sub_mp (p, z, y, x, digs);
 738    } else {
 739  // Subtract.
 740      BOOL_T negative = A68_FALSE;
 741      int fnz, digs_h = 2 + digs;
 742      MP_T *w = nil_mp (p, digs_h);
 743      if (MP_EXPONENT (x) == MP_EXPONENT (y)) {
 744        MP_EXPONENT (w) = (MP_T) 1 + MP_EXPONENT (x);
 745        for (unt j = 1; j <= digs; j++) {
 746          MP_DIGIT (w, j + 1) = MP_DIGIT (x, j) - MP_DIGIT (y, j);
 747        }
 748        MP_DIGIT (w, digs_h) = (MP_T) 0;
 749      } else if (MP_EXPONENT (x) > MP_EXPONENT (y)) {
 750        int shl_y = (int) MP_EXPONENT (x) - (int) MP_EXPONENT (y);
 751        MP_EXPONENT (w) = (MP_T) 1 + MP_EXPONENT (x);
 752        for (unt j = 1; j < digs_h; j++) {
 753          int i_y = j - shl_y;
 754          MP_T x_j = (j > digs ? 0 : MP_DIGIT (x, j));
 755          MP_T y_j = (i_y <= 0 || i_y > digs ? 0 : MP_DIGIT (y, i_y));
 756          MP_DIGIT (w, j + 1) = x_j - y_j;
 757        }
 758      } else {
 759        int shl_x = (int) MP_EXPONENT (y) - (int) MP_EXPONENT (x);
 760        MP_EXPONENT (w) = (MP_T) 1 + MP_EXPONENT (y);
 761        for (unt j = 1; j < digs_h; j++) {
 762          int i_x = j - shl_x;
 763          MP_T x_j = (i_x <= 0 || i_x > digs ? 0 : MP_DIGIT (x, i_x));
 764          MP_T y_j = (j > digs ? 0 : MP_DIGIT (y, j));
 765          MP_DIGIT (w, j + 1) = x_j - y_j;
 766        }
 767      }
 768  // Correct if we subtract large from small.
 769      if (MP_DIGIT (w, 2) <= 0) {
 770        fnz = -1;
 771        for (unt j = 2; j <= digs_h && fnz < 0; j++) {
 772          if (MP_DIGIT (w, j) != 0) {
 773            fnz = j;
 774          }
 775        }
 776        negative = (BOOL_T) (MP_DIGIT (w, fnz) < 0);
 777        if (negative) {
 778          for (unt j = fnz; j <= digs_h; j++) {
 779            MP_DIGIT (w, j) = -MP_DIGIT (w, j);
 780          }
 781        }
 782      }
 783  // Normalise.
 784      norm_mp_light (w, 2, digs_h);
 785      fnz = -1;
 786      for (unt j = 1; j <= digs_h && fnz < 0; j++) {
 787        if (MP_DIGIT (w, j) != 0) {
 788          fnz = j;
 789        }
 790      }
 791      if (fnz > 1) {
 792        int j2 = fnz - 1;
 793        for (int j = 1; j <= digs_h - j2; j++) {
 794          MP_DIGIT (w, j) = MP_DIGIT (w, j + j2);
 795          MP_DIGIT (w, j + j2) = (MP_T) 0;
 796        }
 797        MP_EXPONENT (w) -= j2;
 798      }
 799  // Round.
 800      round_internal_mp (z, w, digs);
 801      if (negative) {
 802        MP_DIGIT (z, 1) = -MP_DIGIT (z, 1);
 803      }
 804      check_mp_exp (p, z);
 805    }
 806  // Restore and exit.
 807    A68_SP = pop_sp;
 808    MP_T z_1 = MP_DIGIT (z, 1);
 809    MP_DIGIT (x, 1) = x_1;
 810    MP_DIGIT (y, 1) = y_1;
 811    MP_DIGIT (z, 1) = z_1;        // In case z IS x OR z IS y
 812    return z;
 813  }
 814  
 815  //! @brief Set "z" to the product of "x" and "y".
 816  
 817  MP_T *mul_mp (NODE_T * p, MP_T * z, MP_T * x, MP_T * y, int digs)
 818  {
 819    if (IS_ZERO_MP (x) || IS_ZERO_MP (y)) {
 820      SET_MP_ZERO (z, digs);
 821      return z;
 822    }
 823  // Grammar school algorithm with intermittent normalisation.
 824    ADDR_T pop_sp = A68_SP;
 825    int digs_h = 2 + digs;
 826    MP_T x_1 = MP_DIGIT (x, 1), y_1 = MP_DIGIT (y, 1);
 827    MP_DIGIT (x, 1) = ABS (x_1);
 828    MP_DIGIT (y, 1) = ABS (y_1);
 829    MP_STATUS (z) = (MP_T) INIT_MASK;
 830    MP_T *w = lit_mp (p, 0, MP_EXPONENT (x) + MP_EXPONENT (y) + 1, digs_h);
 831    int oflow = (int) FLOOR_MP ((MP_REAL_T) MAX_REPR_INT / (2 * MP_REAL_RADIX * MP_REAL_RADIX)) - 1;
 832    for (unt i = digs; i >= 1; i--) {
 833      MP_T yi = MP_DIGIT (y, i);
 834      if (yi != 0) {
 835        int k = digs_h - i;
 836        int j = (k > digs ? digs : k);
 837        MP_T *u = &MP_DIGIT (w, i + j), *v = &MP_DIGIT (x, j);
 838        if ((digs - i + 1) % oflow == 0) {
 839          norm_mp (w, 2, digs_h);
 840        }
 841        while (j-- >= 1) {
 842          (u--)[0] += yi * (v--)[0];
 843        }
 844      }
 845    }
 846    norm_mp (w, 2, digs_h);
 847    round_internal_mp (z, w, digs);
 848  // Restore and exit.
 849    A68_SP = pop_sp;
 850    MP_T z_1 = MP_DIGIT (z, 1);
 851    MP_DIGIT (x, 1) = x_1;
 852    MP_DIGIT (y, 1) = y_1;
 853    MP_DIGIT (z, 1) = ((x_1 * y_1) >= 0 ? z_1 : -z_1);
 854    check_mp_exp (p, z);
 855    return z;
 856  }
 857  
 858  //! @brief Set "z" to the quotient of "x" and "y".
 859  
 860  MP_T *div_mp (NODE_T * p, MP_T * z, MP_T * x, MP_T * y, int digs)
 861  {
 862  // This routine is based on
 863  // 
 864  //    D. M. Smith, "A Multiple-Precision Division Algorithm"
 865  //    Mathematics of Computation 66 (1996) 157-163.
 866  // 
 867  // This is O(N^2) but runs faster than straightforward methods by skipping
 868  // most of the intermediate normalisation and recovering from wrong
 869  // guesses without separate correction steps.
 870  //
 871  // Depending on application, div_mp cost is circa 3 times that of mul_mp.
 872  // Therefore Newton-Raphson division makes no sense here.
 873  //
 874    if (IS_ZERO_MP (y)) {
 875      errno = ERANGE;
 876      return NaN_MP;
 877    }
 878  // Determine normalisation interval assuming that q < 2b in each step.
 879  #if (A68_LEVEL <= 2)
 880    int oflow = (int) FLOOR_MP ((MP_REAL_T) MAX_REPR_INT / (3 * MP_REAL_RADIX * MP_REAL_RADIX)) - 1;
 881  #else
 882    int oflow = (int) FLOOR_MP ((MP_REAL_T) MAX_REPR_INT / (2 * MP_REAL_RADIX * MP_REAL_RADIX)) - 1;
 883  #endif
 884  //
 885    MP_T x_1 = MP_DIGIT (x, 1), y_1 = MP_DIGIT (y, 1);
 886    MP_DIGIT (x, 1) = ABS (x_1);
 887    MP_DIGIT (y, 1) = ABS (y_1);
 888    MP_STATUS (z) = (MP_T) INIT_MASK;
 889  // Slight optimisation when the denominator has few digits.
 890    int nzdigs = digs;
 891    while (MP_DIGIT (y, nzdigs) == 0 && nzdigs > 1) {
 892      nzdigs--;
 893    }
 894    if (nzdigs == 1 && MP_EXPONENT (y) == 0) {
 895      (void) div_mp_digit (p, z, x, MP_DIGIT (y, 1), digs);
 896      MP_T z_1 = MP_DIGIT (z, 1);
 897      MP_DIGIT (x, 1) = x_1;
 898      MP_DIGIT (y, 1) = y_1;
 899      MP_DIGIT (z, 1) = ((x_1 * y_1) >= 0 ? z_1 : -z_1);
 900      check_mp_exp (p, z);
 901      return z;
 902    }
 903  // Working nominator in which the quotient develops.
 904    ADDR_T pop_sp = A68_SP;
 905    int wdigs = 4 + digs;
 906    MP_T *w = lit_mp (p, 0, MP_EXPONENT (x) - MP_EXPONENT (y), wdigs);
 907    (void) move_mp_part (&MP_DIGIT (w, 2), &MP_DIGIT (x, 1), digs);
 908  // Estimate the denominator. For small MP_RADIX add: MP_DIGIT (y, 4) / MP_REAL_RADIX.
 909    MP_REAL_T den = (MP_DIGIT (y, 1) * MP_REAL_RADIX + MP_DIGIT (y, 2)) * MP_REAL_RADIX + MP_DIGIT (y, 3);
 910    MP_T *t = &MP_DIGIT (w, 2);
 911    for (unt k = 1, len = digs + 2, first = 3; k <= digs + 2; k++, len++, first++, t++) {
 912  // Estimate quotient digit.
 913      MP_REAL_T q, nom = ((t[-1] * MP_REAL_RADIX + t[0]) * MP_REAL_RADIX + t[1]) * MP_REAL_RADIX + (wdigs >= (first + 2) ? t[2] : 0);
 914      if (nom == 0) {
 915        q = 0;
 916      } else {
 917  // Correct the nominator.
 918        q = (MP_T) (MP_INT_T) (nom / den);
 919        int lim = MINIMUM (len, wdigs);
 920        if (nzdigs <= lim - first + 1) {
 921          lim = first + nzdigs - 1;
 922        }
 923        MP_T *u = t, *v = &MP_DIGIT (y, 1);
 924        for (unt j = first; j <= lim; j++) {
 925          (u++)[0] -= q * (v++)[0];
 926        }
 927      }
 928      t[0] += t[-1] * MP_RADIX;
 929      t[-1] = q;
 930      if (k % oflow == 0 || k == digs + 2) {
 931        norm_mp (w, first, wdigs);
 932      }
 933    }
 934    norm_mp (w, 2, digs);
 935    round_internal_mp (z, w, digs);
 936  // Restore and exit.
 937    A68_SP = pop_sp;
 938    MP_T z_1 = MP_DIGIT (z, 1);
 939    MP_DIGIT (x, 1) = x_1;
 940    MP_DIGIT (y, 1) = y_1;
 941    MP_DIGIT (z, 1) = ((x_1 * y_1) >= 0 ? z_1 : -z_1);
 942    check_mp_exp (p, z);
 943    return z;
 944  }
 945  
 946  //! @brief Set "z" to the integer quotient of "x" and "y".
 947  
 948  MP_T *over_mp (NODE_T * p, MP_T * z, MP_T * x, MP_T * y, int digs)
 949  {
 950    if (MP_DIGIT (y, 1) == 0) {
 951      errno = ERANGE;
 952      return NaN_MP;
 953    }
 954    int digs_g = FUN_DIGITS (digs);
 955    ADDR_T pop_sp = A68_SP;
 956    MP_T *x_g = len_mp (p, x, digs, digs_g);
 957    MP_T *y_g = len_mp (p, y, digs, digs_g);
 958    MP_T *z_g = nil_mp (p, digs_g);
 959    (void) div_mp (p, z_g, x_g, y_g, digs_g);
 960    trunc_mp (p, z_g, z_g, digs_g);
 961    (void) shorten_mp (p, z, digs, z_g, digs_g);
 962    MP_STATUS (z) = (MP_T) INIT_MASK;
 963  // Restore and exit.
 964    A68_SP = pop_sp;
 965    return z;
 966  }
 967  
 968  //! @brief Set "z" to x mod y.
 969  
 970  MP_T *mod_mp (NODE_T * p, MP_T * z, MP_T * x, MP_T * y, int digs)
 971  {
 972    if (MP_DIGIT (y, 1) == 0) {
 973      errno = EDOM;
 974      return NaN_MP;
 975    }
 976    int digs_g = FUN_DIGITS (digs);
 977    ADDR_T pop_sp = A68_SP;
 978    MP_T *x_g = len_mp (p, x, digs, digs_g);
 979    MP_T *y_g = len_mp (p, y, digs, digs_g);
 980    MP_T *z_g = nil_mp (p, digs_g);
 981  // x mod y = x - y * trunc (x / y).
 982    (void) over_mp (p, z_g, x_g, y_g, digs_g);
 983    (void) mul_mp (p, z_g, y_g, z_g, digs_g);
 984    (void) sub_mp (p, z_g, x_g, z_g, digs_g);
 985    (void) shorten_mp (p, z, digs, z_g, digs_g);
 986  // Restore and exit.
 987    A68_SP = pop_sp;
 988    return z;
 989  }
 990  
 991  //! @brief Set "z" to the product of x and digit y.
 992  
 993  MP_T *mul_mp_digit (NODE_T * p, MP_T * z, MP_T * x, MP_T y, int digs)
 994  {
 995  // This is an O(N) routine for multiplication by a short value.
 996    MP_T x_1 = MP_DIGIT (x, 1);
 997    int digs_h = 2 + digs;
 998    ADDR_T pop_sp = A68_SP;
 999    MP_DIGIT (x, 1) = ABS (x_1);
1000    MP_STATUS (z) = (MP_T) INIT_MASK;
1001    MP_T y_1 = y;
1002    y = ABS (y_1);
1003    if (y == 2) {
1004      (void) add_mp (p, z, x, x, digs);
1005    } else {
1006      MP_T *w = lit_mp (p, 0, MP_EXPONENT (x) + 1, digs_h);
1007      MP_T *u = &MP_DIGIT (w, 1 + digs), *v = &MP_DIGIT (x, digs);
1008      int j = digs;
1009      while (j-- >= 1) {
1010        (u--)[0] += y * (v--)[0];
1011      }
1012      norm_mp (w, 2, digs_h);
1013      round_internal_mp (z, w, digs);
1014    }
1015  // Restore and exit.
1016    A68_SP = pop_sp;
1017    MP_T z_1 = MP_DIGIT (z, 1);
1018    MP_DIGIT (x, 1) = x_1;
1019    MP_DIGIT (z, 1) = ((x_1 * y_1) >= 0 ? z_1 : -z_1);
1020    check_mp_exp (p, z);
1021    return z;
1022  }
1023  
1024  //! @brief Set "z" to x/2.
1025  
1026  MP_T *half_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1027  {
1028    MP_T *w, z_1, x_1 = MP_DIGIT (x, 1), *u, *v;
1029    int j, digs_h = 2 + digs;
1030    ADDR_T pop_sp = A68_SP;
1031    MP_DIGIT (x, 1) = ABS (x_1);
1032    MP_STATUS (z) = (MP_T) INIT_MASK;
1033  // Calculate x * 0.5.
1034    w = lit_mp (p, 0, MP_EXPONENT (x), digs_h);
1035    j = digs;
1036    u = &MP_DIGIT (w, 1 + digs);
1037    v = &MP_DIGIT (x, digs);
1038    while (j-- >= 1) {
1039      (u--)[0] += (MP_RADIX / 2) * (v--)[0];
1040    }
1041    norm_mp (w, 2, digs_h);
1042    round_internal_mp (z, w, digs);
1043  // Restore and exit.
1044    A68_SP = pop_sp;
1045    z_1 = MP_DIGIT (z, 1);
1046    MP_DIGIT (x, 1) = x_1;
1047    MP_DIGIT (z, 1) = (x_1 >= 0 ? z_1 : -z_1);
1048    check_mp_exp (p, z);
1049    return z;
1050  }
1051  
1052  //! @brief Set "z" to x/10.
1053  
1054  MP_T *tenth_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1055  {
1056    MP_T *w, z_1, x_1 = MP_DIGIT (x, 1), *u, *v;
1057    int j, digs_h = 2 + digs;
1058    ADDR_T pop_sp = A68_SP;
1059    MP_DIGIT (x, 1) = ABS (x_1);
1060    MP_STATUS (z) = (MP_T) INIT_MASK;
1061  // Calculate x * 0.1.
1062    w = lit_mp (p, 0, MP_EXPONENT (x), digs_h);
1063    j = digs;
1064    u = &MP_DIGIT (w, 1 + digs);
1065    v = &MP_DIGIT (x, digs);
1066    while (j-- >= 1) {
1067      (u--)[0] += (MP_RADIX / 10) * (v--)[0];
1068    }
1069    norm_mp (w, 2, digs_h);
1070    round_internal_mp (z, w, digs);
1071  // Restore and exit.
1072    A68_SP = pop_sp;
1073    z_1 = MP_DIGIT (z, 1);
1074    MP_DIGIT (x, 1) = x_1;
1075    MP_DIGIT (z, 1) = (x_1 >= 0 ? z_1 : -z_1);
1076    check_mp_exp (p, z);
1077    return z;
1078  }
1079  
1080  //! @brief Set "z" to the quotient of x and digit y.
1081  
1082  MP_T *div_mp_digit (NODE_T * p, MP_T * z, MP_T * x, MP_T y, int digs)
1083  {
1084    if (y == 0) {
1085      errno = ERANGE;
1086      return NaN_MP;
1087    }
1088  // Determine normalisation interval assuming that q < 2b in each step.
1089  #if (A68_LEVEL <= 2)
1090    int oflow = (int) FLOOR_MP ((MP_REAL_T) MAX_REPR_INT / (3 * MP_REAL_RADIX * MP_REAL_RADIX)) - 1;
1091  #else
1092    int oflow = (int) FLOOR_MP ((MP_REAL_T) MAX_REPR_INT / (2 * MP_REAL_RADIX * MP_REAL_RADIX)) - 1;
1093  #endif
1094  // Work with positive operands.
1095    ADDR_T pop_sp = A68_SP;
1096    MP_T x_1 = MP_DIGIT (x, 1), y_1 = y;
1097    MP_DIGIT (x, 1) = ABS (x_1);
1098    MP_STATUS (z) = (MP_T) INIT_MASK;
1099    y = ABS (y_1);
1100  //
1101    if (y == 2) {
1102      (void) half_mp (p, z, x, digs);
1103    } else if (y == 10) {
1104      (void) tenth_mp (p, z, x, digs);
1105    } else {
1106      int wdigs = 4 + digs;
1107      MP_T *w = lit_mp (p, 0, MP_EXPONENT (x), wdigs);
1108      (void) move_mp_part (&MP_DIGIT (w, 2), &MP_DIGIT (x, 1), digs);
1109  // Estimate the denominator.
1110      MP_REAL_T den = (MP_REAL_T) y * MP_REAL_RADIX * MP_REAL_RADIX;
1111      MP_T *t = &MP_DIGIT (w, 2);
1112      for (unt k = 1, first = 3; k <= digs + 2; k++, first++, t++) {
1113  // Estimate quotient digit and correct.
1114        MP_REAL_T nom = ((t[-1] * MP_REAL_RADIX + t[0]) * MP_REAL_RADIX + t[1]) * MP_REAL_RADIX + (wdigs >= (first + 2) ? t[2] : 0);
1115        MP_REAL_T q = (MP_T) (MP_INT_T) (nom / den);
1116        t[0] += t[-1] * MP_RADIX - q * y;
1117        t[-1] = q;
1118        if (k % oflow == 0 || k == digs + 2) {
1119          norm_mp (w, first, wdigs);
1120        }
1121      }
1122      norm_mp (w, 2, digs);
1123      round_internal_mp (z, w, digs);
1124    }
1125  // Restore and exit.
1126    A68_SP = pop_sp;
1127    MP_T z_1 = MP_DIGIT (z, 1);
1128    MP_DIGIT (x, 1) = x_1;
1129    MP_DIGIT (z, 1) = ((x_1 * y_1) >= 0 ? z_1 : -z_1);
1130    check_mp_exp (p, z);
1131    return z;
1132  }
1133  
1134  //! @brief Set "z" to the integer quotient of "x" and "y".
1135  
1136  MP_T *over_mp_digit (NODE_T * p, MP_T * z, MP_T * x, MP_T y, int digs)
1137  {
1138    if (y == 0) {
1139      errno = ERANGE;
1140      return NaN_MP;
1141    }
1142    int digs_g = FUN_DIGITS (digs);
1143    ADDR_T pop_sp = A68_SP;
1144    MP_T *x_g = len_mp (p, x, digs, digs_g);
1145    MP_T *z_g = nil_mp (p, digs_g);
1146    (void) div_mp_digit (p, z_g, x_g, y, digs_g);
1147    trunc_mp (p, z_g, z_g, digs_g);
1148    (void) shorten_mp (p, z, digs, z_g, digs_g);
1149  // Restore and exit.
1150    A68_SP = pop_sp;
1151    return z;
1152  }
1153  
1154  //! @brief Set "z" to the reciprocal of "x".
1155  
1156  MP_T *rec_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1157  {
1158    if (IS_ZERO_MP (x)) {
1159      errno = ERANGE;
1160      return NaN_MP;
1161    }
1162    ADDR_T pop_sp = A68_SP;
1163    (void) div_mp (p, z, mp_one (digs), x, digs);
1164    A68_SP = pop_sp;
1165    return z;
1166  }
1167  
1168  //! @brief LONG REAL long pi
1169  
1170  void genie_pi_mp (NODE_T * p)
1171  {
1172    int digs = DIGITS (MOID (p));
1173    MP_T *z = nil_mp (p, digs);
1174    (void) mp_pi (p, z, MP_PI, digs);
1175    MP_STATUS (z) = (MP_T) INIT_MASK;
1176  }
1177  
1178  //! @brief Set "z" to "x" ** "n".
1179  
1180  MP_T *pow_mp_int (NODE_T * p, MP_T * z, MP_T * x, INT_T n, int digs)
1181  {
1182    ADDR_T pop_sp = A68_SP;
1183    int bit, digs_g = FUN_DIGITS (digs);
1184    BOOL_T negative;
1185    MP_T *x_g = len_mp (p, x, digs, digs_g);
1186    MP_T *z_g = lit_mp (p, 1, 0, digs_g);
1187    negative = (BOOL_T) (n < 0);
1188    if (negative) {
1189      n = -n;
1190    }
1191    bit = 1;
1192    while ((int) bit <= (int) n) {
1193      if (n & bit) {
1194        (void) mul_mp (p, z_g, z_g, x_g, digs_g);
1195      }
1196      (void) mul_mp (p, x_g, x_g, x_g, digs_g);
1197      bit <<= 1;
1198    }
1199    (void) shorten_mp (p, z, digs, z_g, digs_g);
1200    A68_SP = pop_sp;
1201    if (negative) {
1202      (void) rec_mp (p, z, z, digs);
1203    }
1204    check_mp_exp (p, z);
1205    return z;
1206  }
1207  
1208  //! @brief Set "z" to "x" ** "y".
1209  
1210  MP_T *pow_mp (NODE_T * p, MP_T * z, MP_T * x, MP_T * y, int digs)
1211  {
1212    PRELUDE_ERROR (ln_mp (p, z, x, digs) == NaN_MP, p, ERROR_INVALID_ARGUMENT, MOID (p));
1213    (void) mul_mp (p, z, y, z, digs);
1214    (void) exp_mp (p, z, z, digs);
1215    return z;
1216  }
1217  
1218  //! @brief Set "z" to 10 ** "n".
1219  
1220  MP_T *ten_up_mp (NODE_T * p, MP_T * z, int n, int digs)
1221  {
1222  #if (A68_LEVEL >= 3)
1223    static MP_T y[LOG_MP_RADIX] = { 1, 10, 100, 1000, 10000, 100000, 1000000, 10000000, 100000000 };
1224  #else
1225    static MP_T y[LOG_MP_RADIX] = { 1, 10, 100, 1000, 10000, 100000, 1000000 };
1226  #endif
1227    if (n >= 0) {
1228      set_mp (z, y[n % LOG_MP_RADIX], n / LOG_MP_RADIX, digs);
1229    } else {
1230      set_mp (z, y[(LOG_MP_RADIX + n % LOG_MP_RADIX) % LOG_MP_RADIX], (n + 1) / LOG_MP_RADIX - 1, digs);
1231    }
1232    check_mp_exp (p, z);
1233    return z;
1234  }
1235  
1236  //! @brief Comparison of "x" and "y".
1237  
1238  void eq_mp (NODE_T * p, A68_BOOL * z, MP_T * x, MP_T * y, int digs)
1239  {
1240    ADDR_T pop_sp = A68_SP;
1241    MP_T *v = nil_mp (p, digs);
1242    (void) sub_mp (p, v, x, y, digs);
1243    STATUS (z) = INIT_MASK;
1244    VALUE (z) = (MP_DIGIT (v, 1) == 0 ? A68_TRUE : A68_FALSE);
1245    A68_SP = pop_sp;
1246  }
1247  
1248  //! @brief Comparison of "x" and "y".
1249  
1250  void ne_mp (NODE_T * p, A68_BOOL * z, MP_T * x, MP_T * y, int digs)
1251  {
1252    ADDR_T pop_sp = A68_SP;
1253    MP_T *v = nil_mp (p, digs);
1254    (void) sub_mp (p, v, x, y, digs);
1255    STATUS (z) = INIT_MASK;
1256    VALUE (z) = (MP_DIGIT (v, 1) != 0 ? A68_TRUE : A68_FALSE);
1257    A68_SP = pop_sp;
1258  }
1259  
1260  //! @brief Comparison of "x" and "y".
1261  
1262  void lt_mp (NODE_T * p, A68_BOOL * z, MP_T * x, MP_T * y, int digs)
1263  {
1264    ADDR_T pop_sp = A68_SP;
1265    MP_T *v = nil_mp (p, digs);
1266    (void) sub_mp (p, v, x, y, digs);
1267    STATUS (z) = INIT_MASK;
1268    VALUE (z) = (MP_DIGIT (v, 1) < 0 ? A68_TRUE : A68_FALSE);
1269    A68_SP = pop_sp;
1270  }
1271  
1272  //! @brief Comparison of "x" and "y".
1273  
1274  void le_mp (NODE_T * p, A68_BOOL * z, MP_T * x, MP_T * y, int digs)
1275  {
1276    ADDR_T pop_sp = A68_SP;
1277    MP_T *v = nil_mp (p, digs);
1278    (void) sub_mp (p, v, x, y, digs);
1279    STATUS (z) = INIT_MASK;
1280    VALUE (z) = (MP_DIGIT (v, 1) <= 0 ? A68_TRUE : A68_FALSE);
1281    A68_SP = pop_sp;
1282  }
1283  
1284  //! @brief Comparison of "x" and "y".
1285  
1286  void gt_mp (NODE_T * p, A68_BOOL * z, MP_T * x, MP_T * y, int digs)
1287  {
1288    ADDR_T pop_sp = A68_SP;
1289    MP_T *v = nil_mp (p, digs);
1290    (void) sub_mp (p, v, x, y, digs);
1291    STATUS (z) = INIT_MASK;
1292    VALUE (z) = (MP_DIGIT (v, 1) > 0 ? A68_TRUE : A68_FALSE);
1293    A68_SP = pop_sp;
1294  }
1295  
1296  //! @brief Comparison of "x" and "y".
1297  
1298  void ge_mp (NODE_T * p, A68_BOOL * z, MP_T * x, MP_T * y, int digs)
1299  {
1300    ADDR_T pop_sp = A68_SP;
1301    MP_T *v = nil_mp (p, digs);
1302    (void) sub_mp (p, v, x, y, digs);
1303    STATUS (z) = INIT_MASK;
1304    VALUE (z) = (MP_DIGIT (v, 1) >= 0 ? A68_TRUE : A68_FALSE);
1305    A68_SP = pop_sp;
1306  }
1307  
1308  //! @brief round (x).
1309  
1310  MP_T *round_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1311  {
1312    MP_T *y = nil_mp (p, digs);
1313    SET_MP_HALF (y, digs);
1314    if (MP_DIGIT (x, 1) >= 0) {
1315      (void) add_mp (p, z, x, y, digs);
1316      (void) trunc_mp (p, z, z, digs);
1317    } else {
1318      (void) sub_mp (p, z, x, y, digs);
1319      (void) trunc_mp (p, z, z, digs);
1320    }
1321    MP_STATUS (z) = (MP_T) INIT_MASK;
1322    return z;
1323  }
1324  
1325  //! @brief Entier (x).
1326  
1327  MP_T *entier_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
1328  {
1329    if (MP_DIGIT (x, 1) >= 0) {
1330      (void) trunc_mp (p, z, x, digs);
1331    } else {
1332      MP_T *y = nil_mp (p, digs);
1333      (void) move_mp (y, z, digs);
1334      (void) trunc_mp (p, z, x, digs);
1335      (void) sub_mp (p, y, y, z, digs);
1336      if (MP_DIGIT (y, 1) != 0) {
1337        SET_MP_ONE (y, digs);
1338        (void) sub_mp (p, z, z, y, digs);
1339      }
1340    }
1341    MP_STATUS (z) = (MP_T) INIT_MASK;
1342    return z;
1343  }