double.c

     
   1  //! @file double.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 INT, LONG REAL and LONG BITS routines.
  25  
  26  #include "a68g.h"
  27  
  28  #if (A68_LEVEL >= 3)
  29  
  30  #include "a68g-genie.h"
  31  #include "a68g-prelude.h"
  32  #include "a68g-transput.h"
  33  #include "a68g-mp.h"
  34  #include "a68g-double.h"
  35  #include "a68g-lib.h"
  36  #include "a68g-numbers.h"
  37  #include "a68g-quad.h"
  38  
  39  // 128-bit REAL*16 stuff.
  40  
  41  #define RADIX (65536)
  42  #define RADIX_Q (65536.0q)
  43  #define CONST_2_UP_112_Q (5192296858534827628530496329220096.0q)
  44  #define DOUBLE_DIGITS MANT_DIGS (FLT128_MANT_DIG)
  45  
  46  #define IS_ZERO(u) (HW (u) == 0 && LW (u) == 0)
  47  #define EQ(u, v) (HW (u) == HW (v) && LW (u) == LW (v))
  48  #define GT(u, v) (HW (u) != HW (v) ? HW (u) > HW (v) : LW (u) > LW (v))
  49  #define GE(u, v) (HW (u) != HW (v) ? HW (u) >= HW (v) : LW (u) >= LW (v))
  50  
  51  DOUBLE_NUM_T double_ssub (NODE_T *, DOUBLE_NUM_T, DOUBLE_NUM_T);
  52  
  53  void m64to128 (DOUBLE_NUM_T * w, UNSIGNED_T u, UNSIGNED_T v)
  54  {
  55  // Knuth's 'M' algorithm.
  56  #define M (0xffffffff)
  57  #define N 32
  58    UNSIGNED_T k, t, w1, w2, w3;
  59    UNSIGNED_T hu = u >> N, lu = u & M, hv = v >> N, lv = v & M;
  60    t = lu * lv;
  61    w3 = t & M;
  62    k = t >> N;
  63    t = hu * lv + k;
  64    w2 = t & M;
  65    w1 = t >> N;
  66    t = lu * hv + w2;
  67    k = t >> N;
  68    HW (*w) = hu * hv + w1 + k;
  69    LW (*w) = (t << N) + w3;
  70  #undef M
  71  #undef N
  72  }
  73  
  74  void m128to128 (NODE_T * p, MOID_T * m, DOUBLE_NUM_T * w, DOUBLE_NUM_T u, DOUBLE_NUM_T v)
  75  {
  76  // Knuth's 'M' algorithm.
  77    DOUBLE_NUM_T w1, w2, w3;
  78    DOUBLE_NUM_T k, t, h;
  79    UNSIGNED_T hu = HW (u), lu = LW (u), hv = HW (v), lv = LW (v);
  80    if (lu == 0 || lv == 0) {
  81      set_lw (t, 0);
  82    } else {
  83      m64to128 (&t, lu, lv);
  84    }
  85    set_lw (w3, LW (t));
  86    set_lw (k, HW (t));
  87    if (hu == 0 || lv == 0) {
  88      set_lw (t, 0);
  89    } else {
  90      m64to128 (&t, hu, lv);
  91    }
  92    add_double (p, m, t, t, k);
  93    set_lw (w2, LW (t));
  94    set_lw (w1, HW (t));
  95    if (lu == 0 || hv == 0) {
  96      set_lw (t, 0);
  97    } else {
  98      m64to128 (&t, lu, hv);
  99    }
 100    add_double (p, m, t, t, w2);
 101    set_lw (k, HW (t));
 102    if (hu == 0 || hv == 0) {
 103      set_lw (h, 0);
 104    } else {
 105      m64to128 (&h, hu, hv);
 106    }
 107    add_double (p, m, h, h, w1);
 108    add_double (p, m, h, h, k);
 109    set_hw (*w, LW (t));
 110    add_double (p, m, *w, *w, w3);
 111    PRELUDE_ERROR (MODCHK (p, m, HW (h) != 0 || LW (h) != 0), p, ERROR_MATH, M_LONG_INT)
 112  }
 113  
 114  DOUBLE_NUM_T double_udiv (NODE_T * p, MOID_T * m, DOUBLE_NUM_T n, DOUBLE_NUM_T d, int mode)
 115  {
 116  // A bit naive long division.
 117    int k;
 118    UNSIGNED_T carry;
 119    DOUBLE_NUM_T q, r;
 120  // Special cases.
 121    PRELUDE_ERROR (IS_ZERO (d), p, ERROR_DIVISION_BY_ZERO, M_LONG_INT);
 122    if (IS_ZERO (n)) {
 123      if (mode == 0) {
 124        set_lw (q, 0);
 125        return q;
 126      } else {
 127        set_lw (r, 0);
 128        return r;
 129      }
 130    }
 131  // Would n and d be random, then ~50% of the divisions is trivial.
 132    if (EQ (n, d)) {
 133      if (mode == 0) {
 134        set_lw (q, 1);
 135        return q;
 136      } else {
 137        set_lw (r, 0);
 138        return r;
 139      }
 140    } else if (GT (d, n)) {
 141      if (mode == 0) {
 142        set_lw (q, 0);
 143        return q;
 144      } else {
 145        return n;
 146      }
 147    }
 148  // Halfword divide.
 149    if (HW (n) == 0 && HW (d) == 0) {
 150      if (mode == 0) {
 151        set_lw (q, LW (n) / LW (d));
 152        return q;
 153      } else {
 154        set_lw (r, LW (n) % LW (d));
 155        return r;
 156      }
 157    }
 158  // We now know that n and d both have > 64 bits.
 159  // Full divide.
 160    set_lw (q, 0);
 161    set_lw (r, 0);
 162    for (k = 128; k > 0; k--) {
 163      {
 164        carry = (LW (q) & D_SIGN) ? 0x1 : 0x0;
 165        LW (q) <<= 1;
 166        HW (q) = (HW (q) << 1) | carry;
 167      }
 168  // Left-shift r
 169      {
 170        carry = (LW (r) & D_SIGN) ? 0x1 : 0x0;
 171        LW (r) <<= 1;
 172        HW (r) = (HW (r) << 1) | carry;
 173      }
 174  // r[0] = n[k]
 175      {
 176        if (HW (n) & D_SIGN) {
 177          LW (r) |= 0x1;
 178        }
 179        carry = (LW (n) & D_SIGN) ? 0x1 : 0x0;
 180        LW (n) <<= 1;
 181        HW (n) = (HW (n) << 1) | carry;
 182      }
 183  // if r >= d
 184      if (GE (r, d)) {
 185  // r = r - d
 186        sub_double (p, m, r, r, d);
 187  // q[k] = 1
 188        LW (q) |= 0x1;
 189      }
 190    }
 191    if (mode == 0) {
 192      return q;
 193    } else {
 194      return r;
 195    }
 196  }
 197  
 198  DOUBLE_NUM_T double_uadd (NODE_T * p, MOID_T * m, DOUBLE_NUM_T u, DOUBLE_NUM_T v)
 199  {
 200    DOUBLE_NUM_T w;
 201    (void) p;
 202    add_double (p, m, w, u, v);
 203    return w;
 204  }
 205  
 206  DOUBLE_NUM_T double_usub (NODE_T * p, MOID_T * m, DOUBLE_NUM_T u, DOUBLE_NUM_T v)
 207  {
 208    DOUBLE_NUM_T w;
 209    (void) p;
 210    sub_double (p, m, w, u, v);
 211    return w;
 212  }
 213  
 214  DOUBLE_NUM_T double_umul (NODE_T * p, MOID_T * m, DOUBLE_NUM_T u, DOUBLE_NUM_T v)
 215  {
 216    DOUBLE_NUM_T w;
 217    m128to128 (p, m, &w, u, v);
 218    return w;
 219  }
 220  
 221  // Signed integer.
 222  
 223  DOUBLE_NUM_T double_sadd (NODE_T * p, DOUBLE_NUM_T u, DOUBLE_NUM_T v)
 224  {
 225    DOUBLE_NUM_T w;
 226    int neg_u = D_NEG (u), neg_v = D_NEG (v);
 227    set_lw (w, 0);
 228    if (neg_u) {
 229      u = neg_double_int (u);
 230    }
 231    if (neg_v) {
 232      v = neg_double_int (v);
 233    }
 234    if (!neg_u && !neg_v) {
 235      w = double_uadd (p, M_LONG_INT, u, v);
 236      PRELUDE_ERROR (D_NEG (w), p, ERROR_MATH, M_LONG_INT);
 237    } else if (neg_u && neg_v) {
 238      w = neg_double_int (double_sadd (p, u, v));
 239    } else if (neg_u) {
 240      w = double_ssub (p, v, u);
 241    } else if (neg_v) {
 242      w = double_ssub (p, u, v);
 243    }
 244    return w;
 245  }
 246  
 247  DOUBLE_NUM_T double_ssub (NODE_T * p, DOUBLE_NUM_T u, DOUBLE_NUM_T v)
 248  {
 249    DOUBLE_NUM_T w;
 250    int neg_u = D_NEG (u), neg_v = D_NEG (v);
 251    set_lw (w, 0);
 252    if (neg_u) {
 253      u = neg_double_int (u);
 254    }
 255    if (neg_v) {
 256      v = neg_double_int (v);
 257    }
 258    if (!neg_u && !neg_v) {
 259      if (D_LT (u, v)) {
 260        w = neg_double_int (double_usub (p, M_LONG_INT, v, u));
 261      } else {
 262        w = double_usub (p, M_LONG_INT, u, v);
 263      }
 264    } else if (neg_u && neg_v) {
 265      w = double_ssub (p, v, u);
 266    } else if (neg_u) {
 267      w = neg_double_int (double_sadd (p, u, v));
 268    } else if (neg_v) {
 269      w = double_sadd (p, u, v);
 270    }
 271    return w;
 272  }
 273  
 274  DOUBLE_NUM_T double_smul (NODE_T * p, DOUBLE_NUM_T u, DOUBLE_NUM_T v)
 275  {
 276    DOUBLE_NUM_T w;
 277    int neg_u = D_NEG (u), neg_v = D_NEG (v);
 278    if (neg_u) {
 279      u = neg_double_int (u);
 280    }
 281    if (neg_v) {
 282      v = neg_double_int (v);
 283    }
 284    w = double_umul (p, M_LONG_INT, u, v);
 285    if (neg_u != neg_v) {
 286      w = neg_double_int (w);
 287    }
 288    return w;
 289  }
 290  
 291  DOUBLE_NUM_T double_sdiv (NODE_T * p, DOUBLE_NUM_T u, DOUBLE_NUM_T v, int mode)
 292  {
 293    DOUBLE_NUM_T w;
 294    int neg_u = D_NEG (u), neg_v = D_NEG (v);
 295    if (neg_u) {
 296      u = neg_double_int (u);
 297    }
 298    if (neg_v) {
 299      v = neg_double_int (v);
 300    }
 301    w = double_udiv (p, M_LONG_INT, u, v, mode);
 302    if (mode == 0 && neg_u != neg_v) {
 303      w = neg_double_int (w);
 304    } else if (mode == 1 && D_NEG (w)) {
 305      w = double_sadd (p, w, v);
 306    }
 307    return w;
 308  }
 309  
 310  // Infinity.
 311  
 312  DOUBLE_T a68_divq (DOUBLE_T x, DOUBLE_T y)
 313  {
 314    return x / y;
 315  }
 316  
 317  DOUBLE_T a68_dposinf (void)
 318  {
 319    return a68_divq (+1.0, 0.0);
 320  }
 321  
 322  DOUBLE_T a68_dneginf (void)
 323  {
 324    return a68_divq (-1.0, 0.0);
 325  }
 326  
 327  //! @brief Sqrt (x^2 + y^2) that does not needlessly overflow.
 328  
 329  DOUBLE_T a68_double_hypot (DOUBLE_T x, DOUBLE_T y)
 330  {
 331    DOUBLE_T xabs = ABSQ (x), yabs = ABSQ (y), min, max;
 332    if (xabs < yabs) {
 333      min = xabs;
 334      max = yabs;
 335    } else {
 336      min = yabs;
 337      max = xabs;
 338    }
 339    if (min == 0.0q) {
 340      return max;
 341    } else {
 342      DOUBLE_T u = min / max;
 343      return max * sqrtq (1.0q + u * u);
 344    }
 345  }
 346  
 347  // Conversions.
 348  
 349  DOUBLE_NUM_T double_int_to_double_real (NODE_T * p, DOUBLE_NUM_T z)
 350  {
 351    DOUBLE_NUM_T w, radix;
 352    DOUBLE_T weight;
 353    int neg = D_NEG (z);
 354    if (neg) {
 355      z = abs_double_int (z);
 356    }
 357    w.f = 0.0q;
 358    set_lw (radix, RADIX);
 359    weight = 1.0q;
 360    while (!D_ZERO (z)) {
 361      DOUBLE_NUM_T digit;
 362      digit = double_udiv (p, M_LONG_INT, z, radix, 1);
 363      w.f = w.f + LW (digit) * weight;
 364      z = double_udiv (p, M_LONG_INT, z, radix, 0);
 365      weight = weight * RADIX_Q;
 366    }
 367    if (neg) {
 368      w.f = -w.f;
 369    }
 370    return w;
 371  }
 372  
 373  DOUBLE_NUM_T double_real_to_double_int (NODE_T * p, DOUBLE_NUM_T z)
 374  {
 375  // This routines looks a lot like "strtol". 
 376    DOUBLE_NUM_T sum, weight, radix;
 377    BOOL_T negative = (BOOL_T) (z.f < 0);
 378    z.f = fabsq (truncq (z.f));
 379    if (z.f > CONST_2_UP_112_Q) {
 380      errno = EDOM;
 381      MATH_RTE (p, errno != 0, M_LONG_REAL, NO_TEXT);
 382    }
 383    set_lw (sum, 0);
 384    set_lw (weight, 1);
 385    set_lw (radix, RADIX);
 386    while (z.f > 0) {
 387      DOUBLE_NUM_T term, digit, quot, rest;
 388      quot.f = truncq (z.f / RADIX_Q);
 389      rest.f = z.f - quot.f * RADIX_Q;
 390      z.f = quot.f;
 391      set_lw (digit, (INT_T) (rest.f));
 392      term = double_umul (p, M_LONG_INT, digit, weight);
 393      sum = double_uadd (p, M_LONG_INT, sum, term);
 394      if (z.f > 0.0q) {
 395        weight = double_umul (p, M_LONG_INT, weight, radix);
 396      }
 397    }
 398    if (negative) {
 399      return neg_double_int (sum);
 400    } else {
 401      return sum;
 402    }
 403  }
 404  
 405  //! @brief Value of LONG INT denotation
 406  
 407  int string_to_double_int (NODE_T * p, A68_LONG_INT * z, char *s)
 408  {
 409    int k, end, sign;
 410    DOUBLE_NUM_T weight, ten, sum;
 411    while (IS_SPACE (s[0])) {
 412      s++;
 413    }
 414  // Get the sign
 415    sign = (s[0] == '-' ? -1 : 1);
 416    if (s[0] == '+' || s[0] == '-') {
 417      s++;
 418    }
 419    end = 0;
 420    while (s[end] != '\0') {
 421      end++;
 422    }
 423    set_lw (sum, 0);
 424    set_lw (weight, 1);
 425    set_lw (ten, 10);
 426    for (k = end - 1; k >= 0; k--) {
 427      DOUBLE_NUM_T term;
 428      int digit = s[k] - '0';
 429      set_lw (term, digit);
 430      term = double_umul (p, M_LONG_INT, term, weight);
 431      sum = double_uadd (p, M_LONG_INT, sum, term);
 432      weight = double_umul (p, M_LONG_INT, weight, ten);
 433    }
 434    if (sign == -1) {
 435      HW (sum) = HW (sum) | D_SIGN;
 436    }
 437    VALUE (z) = sum;
 438    STATUS (z) = INIT_MASK;
 439    return A68_TRUE;
 440  }
 441  
 442  //! @brief LONG BITS value of LONG BITS denotation
 443  
 444  DOUBLE_NUM_T double_strtou (NODE_T * p, char *s)
 445  {
 446    int base = 0;
 447    DOUBLE_NUM_T z;
 448    char *radix = NO_TEXT;
 449    errno = 0;
 450    base = (int) a68_strtou (s, &radix, 10);
 451    if (base < 2 || base > 16) {
 452      diagnostic (A68_RUNTIME_ERROR, p, ERROR_INVALID_RADIX, base);
 453      exit_genie (p, A68_RUNTIME_ERROR);
 454    }
 455    set_lw (z, 0x0);
 456    if (radix != NO_TEXT && TO_UPPER (radix[0]) == TO_UPPER (RADIX_CHAR) && errno == 0) {
 457      DOUBLE_NUM_T w;
 458      char *q = radix;
 459      while (q[0] != NULL_CHAR) {
 460        q++;
 461      }
 462      set_lw (w, 1);
 463      while ((--q) != radix) {
 464        int digit = char_value (q[0]);
 465        if (digit < 0 && digit >= base) {
 466          diagnostic (A68_RUNTIME_ERROR, p, ERROR_IN_DENOTATION, M_LONG_BITS);
 467          exit_genie (p, A68_RUNTIME_ERROR);
 468        } else {
 469          DOUBLE_NUM_T v;
 470          set_lw (v, digit);
 471          v = double_umul (p, M_LONG_INT, v, w);
 472          z = double_uadd (p, M_LONG_INT, z, v);
 473          set_lw (v, base);
 474          w = double_umul (p, M_LONG_INT, w, v);
 475        }
 476      }
 477    } else {
 478      diagnostic (A68_RUNTIME_ERROR, p, ERROR_IN_DENOTATION, M_LONG_BITS);
 479      exit_genie (p, A68_RUNTIME_ERROR);
 480    }
 481    return (z);
 482  }
 483  
 484  //! @brief OP LENG = (BITS) LONG BITS
 485  
 486  void genie_lengthen_bits_to_double_bits (NODE_T * p)
 487  {
 488    A68_BITS k;
 489    DOUBLE_NUM_T d;
 490    POP_OBJECT (p, &k, A68_BITS);
 491    LW (d) = VALUE (&k);
 492    HW (d) = 0;
 493    PUSH_VALUE (p, d, A68_LONG_BITS);
 494  }
 495  
 496  //! @brief OP SHORTEN = (LONG BITS) BITS
 497  
 498  void genie_shorten_double_bits_to_bits (NODE_T * p)
 499  {
 500    A68_LONG_BITS k;
 501    DOUBLE_NUM_T j;
 502    POP_OBJECT (p, &k, A68_LONG_BITS);
 503    j = VALUE (&k);
 504    PRELUDE_ERROR (HW (j) != 0, p, ERROR_MATH, M_BITS);
 505    PUSH_VALUE (p, LW (j), A68_BITS);
 506  }
 507  
 508  //! @brief Convert to other radix, binary up to hexadecimal.
 509  
 510  BOOL_T convert_radix_double (NODE_T * p, DOUBLE_NUM_T z, int radix, int width)
 511  {
 512    DOUBLE_NUM_T w, rad;
 513    if (radix < 2 || radix > 16) {
 514      radix = 16;
 515    }
 516    set_lw (rad, radix);
 517    reset_transput_buffer (EDIT_BUFFER);
 518    if (width > 0) {
 519      while (width > 0) {
 520        w = double_udiv (p, M_LONG_INT, z, rad, 1);
 521        plusto_transput_buffer (p, digchar (LW (w)), EDIT_BUFFER);
 522        width--;
 523        z = double_udiv (p, M_LONG_INT, z, rad, 0);
 524      }
 525      return D_ZERO (z);
 526    } else if (width == 0) {
 527      do {
 528        w = double_udiv (p, M_LONG_INT, z, rad, 1);
 529        plusto_transput_buffer (p, digchar (LW (w)), EDIT_BUFFER);
 530        z = double_udiv (p, M_LONG_INT, z, rad, 0);
 531      }
 532      while (!D_ZERO (z));
 533      return A68_TRUE;
 534    } else {
 535      return A68_FALSE;
 536    }
 537  }
 538  
 539  //! @brief OP LENG = (LONG INT) LONG REAL
 540  
 541  void genie_widen_double_int_to_double_real (NODE_T * p)
 542  {
 543    A68_DOUBLE *z = (A68_DOUBLE *) STACK_TOP;
 544    EXECUTE_UNIT (SUB (p));
 545    VALUE (z) = double_int_to_double_real (p, VALUE (z));
 546  }
 547  
 548  //! @brief OP LENG = (REAL) LONG REAL
 549  
 550  DOUBLE_NUM_T dble_double_real (NODE_T * p, REAL_T z)
 551  {
 552  // Quick and dirty, only works with 64-bit INT_T.
 553    BOOL_T nega = (z < 0.0);
 554    REAL_T u = fabs (z);
 555    DOUBLE_NUM_T w;
 556    int expo = 0;
 557    standardise (&u, 1, REAL_DIG, &expo);
 558    u *= ten_up (REAL_DIG);
 559    expo -= REAL_DIG;
 560    set_lw (w, (INT_T) u);
 561    w = double_int_to_double_real (p, w);
 562    w.f *= ten_up_double (expo);
 563    if (nega) {
 564      w.f = -w.f;
 565    }
 566    return w;
 567  }
 568  
 569  //! @brief OP LENG = (REAL) LONG REAL
 570  
 571  void genie_lengthen_real_to_double_real (NODE_T * p)
 572  {
 573    A68_REAL z;
 574    POP_OBJECT (p, &z, A68_REAL);
 575    PUSH_VALUE (p, dble_double_real (p, VALUE (&z)), A68_LONG_REAL);
 576  }
 577  
 578  //! @brief OP SHORTEN = (LONG REAL) REAL
 579  
 580  void genie_shorten_double_real_to_real (NODE_T * p)
 581  {
 582    A68_LONG_REAL z;
 583    REAL_T w;
 584    POP_OBJECT (p, &z, A68_LONG_REAL);
 585    w = VALUE (&z).f;
 586    PUSH_VALUE (p, w, A68_REAL);
 587  }
 588  
 589  //! @brief Convert integer to multi-precison number.
 590  
 591  MP_T *double_int_to_mp (NODE_T * p, MP_T * z, DOUBLE_NUM_T k, int digs)
 592  {
 593    DOUBLE_NUM_T k2, radix;
 594    int n = 0, j, negative = D_NEG (k);
 595    if (negative) {
 596      k = neg_double_int (k);
 597    }
 598    set_lw (radix, MP_RADIX);
 599    k2 = k;
 600    do {
 601      k2 = double_udiv (p, M_LONG_INT, k2, radix, 0);
 602      if (!D_ZERO (k2)) {
 603        n++;
 604      }
 605    }
 606    while (!D_ZERO (k2));
 607    SET_MP_ZERO (z, digs);
 608    MP_EXPONENT (z) = (MP_T) n;
 609    for (j = 1 + n; j >= 1; j--) {
 610      DOUBLE_NUM_T term = double_udiv (p, M_LONG_INT, k, radix, 1);
 611      MP_DIGIT (z, j) = (MP_T) LW (term);
 612      k = double_udiv (p, M_LONG_INT, k, radix, 0);
 613    }
 614    MP_DIGIT (z, 1) = (negative ? -MP_DIGIT (z, 1) : MP_DIGIT (z, 1));
 615    check_mp_exp (p, z);
 616    return z;
 617  }
 618  
 619  //! @brief Convert multi-precision number to integer.
 620  
 621  DOUBLE_NUM_T mp_to_double_int (NODE_T * p, MP_T * z, int digs)
 622  {
 623  // This routines looks a lot like "strtol". 
 624    int j, expo = (int) MP_EXPONENT (z);
 625    DOUBLE_NUM_T sum, weight;
 626    set_lw (sum, 0);
 627    set_lw (weight, 1);
 628    BOOL_T negative;
 629    if (expo >= digs) {
 630      diagnostic (A68_RUNTIME_ERROR, p, ERROR_OUT_OF_BOUNDS, MOID (p));
 631      exit_genie (p, A68_RUNTIME_ERROR);
 632    }
 633    negative = (BOOL_T) (MP_DIGIT (z, 1) < 0);
 634    if (negative) {
 635      MP_DIGIT (z, 1) = -MP_DIGIT (z, 1);
 636    }
 637    for (j = 1 + expo; j >= 1; j--) {
 638      DOUBLE_NUM_T term, digit, radix;
 639      set_lw (digit, (MP_INT_T) MP_DIGIT (z, j));
 640      term = double_umul (p, M_LONG_INT, digit, weight);
 641      sum = double_uadd (p, M_LONG_INT, sum, term);
 642      set_lw (radix, MP_RADIX);
 643      weight = double_umul (p, M_LONG_INT, weight, radix);
 644    }
 645    if (negative) {
 646      return neg_double_int (sum);
 647    } else {
 648      return sum;
 649    }
 650  }
 651  
 652  //! @brief Convert real to multi-precison number.
 653  
 654  MP_T *double_real_to_mp (NODE_T * p, MP_T * z, DOUBLE_T x, int digs)
 655  {
 656    int j, k, sign_x, sum, weight;
 657    SET_MP_ZERO (z, digs);
 658    if (x == 0.0q) {
 659      return z;
 660    }
 661  // Small integers can be done better by int_to_mp.
 662    if (ABS (x) < MP_RADIX && truncq (x) == x) {
 663      return int_to_mp (p, z, (int) truncq (x), digs);
 664    }
 665    sign_x = SIGN (x);
 666  // Scale to [0, 0.1>.
 667    DOUBLE_T a = ABS (x);
 668    INT_T expo = (int) log10q (a);
 669    a /= ten_up_double (expo);
 670    expo--;
 671    if (a >= 1.0q) {
 672      a /= 10.0q;
 673      expo++;
 674    }
 675  // Transport digits of x to the mantissa of z.
 676    sum = 0;
 677    weight = (MP_RADIX / 10);
 678    for (k = 0, j = 1; a != 0.0q && j <= digs && k < DOUBLE_DIGITS; k++) {
 679      DOUBLE_T u = a * 10.0q;
 680      DOUBLE_T v = floorq (u);
 681      a = u - v;
 682      sum += weight * (int) v;
 683      weight /= 10;
 684      if (weight < 1) {
 685        MP_DIGIT (z, j++) = (MP_T) sum;
 686        sum = 0;
 687        weight = (MP_RADIX / 10);
 688      }
 689    }
 690  // Store the last digits.
 691    if (j <= digs) {
 692      MP_DIGIT (z, j) = (MP_T) sum;
 693    }
 694    (void) align_mp (z, &expo, digs);
 695    MP_EXPONENT (z) = (MP_T) expo;
 696    MP_DIGIT (z, 1) *= sign_x;
 697    check_mp_exp (p, z);
 698    return z;
 699  }
 700  
 701  //! @brief Convert multi-precision number to real.
 702  
 703  DOUBLE_T mp_to_double_real (NODE_T * p, MP_T * z, int digs)
 704  {
 705  // This routine looks a lot like "strtod".
 706    if (MP_EXPONENT (z) * (MP_T) LOG_MP_RADIX <= (MP_T) DOUBLE_MIN_10_EXP) {
 707      return 0;
 708    } else {
 709      DOUBLE_T terms[1 + MP_MAX_DIGITS];
 710      DOUBLE_T weight = ten_up_double ((int) (MP_EXPONENT (z) * LOG_MP_RADIX));
 711      int lim = MIN (digs, MP_MAX_DIGITS);
 712      for (unt k = 1; k <= lim; k++) {
 713        terms[k] = ABS (MP_DIGIT (z, k)) * weight;
 714        weight /= MP_RADIX;
 715      }
 716  // Sum terms from small to large.
 717      DOUBLE_T sum = 0;
 718      for (unt k = lim; k >= 1; k--) {
 719        sum += terms[k];
 720      }
 721      CHECK_DOUBLE_REAL (p, sum);
 722      return MP_DIGIT (z, 1) >= 0 ? sum : -sum;
 723    }
 724  }
 725  
 726  DOUBLE_T inverf_double_real (DOUBLE_T z)
 727  {
 728    if (fabsq (z) >= 1.0q) {
 729      errno = EDOM;
 730      return z;
 731    } else {
 732  // Newton-Raphson.
 733      DOUBLE_T f = sqrtq (M_PIq) / 2, g, x = z;
 734      int its = 10;
 735      x = dble (a68_inverf ((REAL_T) x)).f;
 736      do {
 737        g = x;
 738        x -= f * (erfq (x) - z) / expq (-(x * x));
 739      } while (its-- > 0 && errno == 0 && fabsq (x - g) > (3 * FLT128_EPSILON));
 740      return x;
 741    }
 742  }
 743  
 744  //! @brief OP LENG = (LONG REAL) LONG LONG REAL
 745  
 746  void genie_lengthen_double_real_to_mp (NODE_T * p)
 747  {
 748    int digs = DIGITS (M_LONG_LONG_REAL);
 749    A68_LONG_REAL x;
 750    POP_OBJECT (p, &x, A68_LONG_REAL);
 751    MP_T *z = nil_mp (p, digs);
 752    (void) double_real_to_mp (p, z, VALUE (&x).f, digs);
 753    MP_STATUS (z) = (MP_T) INIT_MASK;
 754  }
 755  
 756  //! @brief OP SHORTEN = (LONG LONG REAL) LONG REAL
 757  
 758  void genie_shorten_mp_to_double_real (NODE_T * p)
 759  {
 760    MOID_T *mode = LHS_MODE (p);
 761    int digs = DIGITS (mode), size = SIZE (mode);
 762    MP_T *z;
 763    DOUBLE_NUM_T d;
 764    DECREMENT_STACK_POINTER (p, size);
 765    z = (MP_T *) STACK_TOP;
 766    MP_STATUS (z) = (MP_T) INIT_MASK;
 767    d.f = mp_to_double_real (p, z, digs);
 768    PUSH_VALUE (p, d, A68_LONG_REAL);
 769  }
 770  
 771  //! @brief OP SHORTEN = (LONG LONG COMPLEX) LONG COMPLEX
 772  
 773  void genie_shorten_long_mp_complex_to_double_compl (NODE_T * p)
 774  {
 775    int digs = DIGITS (M_LONG_LONG_REAL), size = SIZE (M_LONG_LONG_REAL);
 776    MP_T *b = (MP_T *) STACK_OFFSET (-size);
 777    MP_T *a = (MP_T *) STACK_OFFSET (-2 * size);
 778    DECREMENT_STACK_POINTER (p, 2 * size);
 779    DOUBLE_NUM_T u, v;
 780    u.f = mp_to_double_real (p, a, digs);
 781    v.f = mp_to_double_real (p, b, digs);
 782    PUSH_VALUE (p, u, A68_LONG_REAL);
 783    PUSH_VALUE (p, v, A68_LONG_REAL);
 784  }
 785  
 786  //! @brief OP LENG = (LONG INT) LONG LONG INT
 787  
 788  void genie_lengthen_double_int_to_mp (NODE_T * p)
 789  {
 790    int digs = DIGITS (M_LONG_LONG_INT);
 791    A68_LONG_INT k;
 792    POP_OBJECT (p, &k, A68_LONG_INT);
 793    MP_T *z = nil_mp (p, digs);
 794    (void) double_int_to_mp (p, z, VALUE (&k), digs);
 795    MP_STATUS (z) = (MP_T) INIT_MASK;
 796  }
 797  
 798  //! @brief OP SHORTEN = (LONG LONG INT) LONG INT
 799  
 800  void genie_shorten_mp_to_double_int (NODE_T * p)
 801  {
 802    MOID_T *mode = LHS_MODE (p);
 803    int digs = DIGITS (mode), size = SIZE (mode);
 804    MP_T *z;
 805    DECREMENT_STACK_POINTER (p, size);
 806    z = (MP_T *) STACK_TOP;
 807    MP_STATUS (z) = (MP_T) INIT_MASK;
 808    PUSH_VALUE (p, mp_to_double_int (p, z, digs), A68_LONG_INT);
 809  }
 810  
 811  //! @brief OP LENG = (INT) LONG INT
 812  
 813  void genie_lengthen_int_to_double_int (NODE_T * p)
 814  {
 815    A68_INT k;
 816    INT_T v;
 817    DOUBLE_NUM_T d;
 818    POP_OBJECT (p, &k, A68_INT);
 819    v = VALUE (&k);
 820    if (v >= 0) {
 821      LW (d) = v;
 822      HW (d) = 0;
 823    } else {
 824      LW (d) = -v;
 825      HW (d) = D_SIGN;
 826    }
 827    PUSH_VALUE (p, d, A68_LONG_INT);
 828  }
 829  
 830  //! @brief OP SHORTEN = (LONG INT) INT
 831  
 832  void genie_shorten_long_int_to_int (NODE_T * p)
 833  {
 834    A68_LONG_INT k;
 835    DOUBLE_NUM_T j;
 836    POP_OBJECT (p, &k, A68_LONG_INT);
 837    j = VALUE (&k);
 838    PRELUDE_ERROR (HW (j) != 0 && HW (j) != D_SIGN, p, ERROR_MATH, M_INT);
 839    PRELUDE_ERROR (LW (j) & D_SIGN, p, ERROR_MATH, M_INT);
 840    if (D_NEG (j)) {
 841      PUSH_VALUE (p, -LW (j), A68_INT);
 842    } else {
 843      PUSH_VALUE (p, LW (j), A68_INT);
 844    }
 845  }
 846  
 847  // Constants.
 848  
 849  //! @brief PROC long max int = LONG INT
 850  
 851  void genie_double_max_int (NODE_T * p)
 852  {
 853    DOUBLE_NUM_T d;
 854    HW (d) = 0x7fffffffffffffffLL;
 855    LW (d) = 0xffffffffffffffffLL;
 856    PUSH_VALUE (p, d, A68_LONG_INT);
 857  }
 858  
 859  //! @brief PROC long max bits = LONG BITS
 860  
 861  void genie_double_max_bits (NODE_T * p)
 862  {
 863    DOUBLE_NUM_T d;
 864    HW (d) = 0xffffffffffffffffLL;
 865    LW (d) = 0xffffffffffffffffLL;
 866    PUSH_VALUE (p, d, A68_LONG_INT);
 867  }
 868  
 869  //! @brief LONG REAL max long real
 870  
 871  void genie_double_max_real (NODE_T * p)
 872  {
 873    DOUBLE_NUM_T d;
 874    d.f = FLT128_MAX;
 875    PUSH_VALUE (p, d, A68_LONG_REAL);
 876  }
 877  
 878  //! @brief LONG REAL min long real
 879  
 880  void genie_double_min_real (NODE_T * p)
 881  {
 882    DOUBLE_NUM_T d;
 883    d.f = FLT128_MIN;
 884    PUSH_VALUE (p, d, A68_LONG_REAL);
 885  }
 886  
 887  //! @brief LONG REAL small long real
 888  
 889  void genie_double_small_real (NODE_T * p)
 890  {
 891    DOUBLE_NUM_T d;
 892    d.f = FLT128_EPSILON;
 893    PUSH_VALUE (p, d, A68_LONG_REAL);
 894  }
 895  
 896  //! @brief PROC long pi = LON REAL
 897  
 898  void genie_pi_double (NODE_T * p)
 899  {
 900    DOUBLE_NUM_T w;
 901    w.f = M_PIq;
 902    PUSH_VALUE (p, w, A68_LONG_INT);
 903  }
 904  
 905  // MONADs and DYADs
 906  
 907  //! @brief OP SIGN = (LONG INT) INT
 908  
 909  void genie_sign_double_int (NODE_T * p)
 910  {
 911    A68_LONG_INT k;
 912    POP_OBJECT (p, &k, A68_LONG_INT);
 913    PUSH_VALUE (p, sign_double_int (VALUE (&k)), A68_INT);
 914  }
 915  
 916  //! @brief OP ABS = (LONG INT) LONG INT
 917  
 918  void genie_abs_double_int (NODE_T * p)
 919  {
 920    A68_LONG_INT *k;
 921    POP_OPERAND_ADDRESS (p, k, A68_LONG_INT);
 922    VALUE (k) = abs_double_int (VALUE (k));
 923  }
 924  
 925  //! @brief OP ODD = (LONG INT) BOOL
 926  
 927  void genie_odd_double_int (NODE_T * p)
 928  {
 929    A68_LONG_INT j;
 930    DOUBLE_NUM_T w;
 931    POP_OBJECT (p, &j, A68_LONG_INT);
 932    w = abs_double_int (VALUE (&j));
 933    if (LW (w) & 0x1) {
 934      PUSH_VALUE (p, A68_TRUE, A68_BOOL);
 935    } else {
 936      PUSH_VALUE (p, A68_FALSE, A68_BOOL);
 937    }
 938  }
 939  
 940  //! @brief OP - = (LONG INT) LONG INT
 941  
 942  void genie_minus_double_int (NODE_T * p)
 943  {
 944    A68_LONG_INT *k;
 945    POP_OPERAND_ADDRESS (p, k, A68_LONG_INT);
 946    VALUE (k) = neg_double_int (VALUE (k));
 947  }
 948  
 949  //! @brief OP ROUND = (LONG REAL) LONG INT
 950  
 951  void genie_round_double_real (NODE_T * p)
 952  {
 953    A68_LONG_REAL x;
 954    DOUBLE_NUM_T u;
 955    POP_OBJECT (p, &x, A68_LONG_REAL);
 956    u = VALUE (&x);
 957    if (u.f < 0.0q) {
 958      u.f = u.f - 0.5q;
 959    } else {
 960      u.f = u.f + 0.5q;
 961    }
 962    PUSH_VALUE (p, double_real_to_double_int (p, u), A68_LONG_INT);
 963  }
 964  
 965  //! @brief OP ENTIER = (LONG REAL) LONG INT
 966  
 967  void genie_entier_double_real (NODE_T * p)
 968  {
 969    A68_LONG_REAL x;
 970    DOUBLE_NUM_T u;
 971    POP_OBJECT (p, &x, A68_LONG_REAL);
 972    u = VALUE (&x);
 973    u.f = floorq (u.f);
 974    PUSH_VALUE (p, double_real_to_double_int (p, u), A68_LONG_INT);
 975  }
 976  
 977  //! @brief OP + = (LONG INT, LONG INT) LONG INT
 978  
 979  void genie_add_double_int (NODE_T * p)
 980  {
 981    A68_LONG_INT i, j;
 982    POP_OBJECT (p, &j, A68_LONG_INT);
 983    POP_OBJECT (p, &i, A68_LONG_INT);
 984    PUSH_VALUE (p, double_sadd (p, VALUE (&i), VALUE (&j)), A68_LONG_INT);
 985  }
 986  
 987  //! @brief OP - = (LONG INT, LONG INT) LONG INT
 988  
 989  void genie_sub_double_int (NODE_T * p)
 990  {
 991    A68_LONG_INT i, j;
 992    POP_OBJECT (p, &j, A68_LONG_INT);
 993    POP_OBJECT (p, &i, A68_LONG_INT);
 994    PUSH_VALUE (p, double_ssub (p, VALUE (&i), VALUE (&j)), A68_LONG_INT);
 995  }
 996  
 997  //! @brief OP * = (LONG INT, LONG INT) LONG INT
 998  
 999  void genie_mul_double_int (NODE_T * p)
1000  {
1001    A68_LONG_INT i, j;
1002    POP_OBJECT (p, &j, A68_LONG_INT);
1003    POP_OBJECT (p, &i, A68_LONG_INT);
1004    PUSH_VALUE (p, double_smul (p, VALUE (&i), VALUE (&j)), A68_LONG_INT);
1005  }
1006  
1007  //! @brief OP / = (LONG INT, LONG INT) LONG INT
1008  
1009  void genie_over_double_int (NODE_T * p)
1010  {
1011    A68_LONG_INT i, j;
1012    POP_OBJECT (p, &j, A68_LONG_INT);
1013    POP_OBJECT (p, &i, A68_LONG_INT);
1014    PRELUDE_ERROR (D_ZERO (VALUE (&j)), p, ERROR_DIVISION_BY_ZERO, M_LONG_INT);
1015    PUSH_VALUE (p, double_sdiv (p, VALUE (&i), VALUE (&j), 0), A68_LONG_INT);
1016  }
1017  
1018  //! @brief OP MOD = (LONG INT, LONG INT) LONG INT
1019  
1020  void genie_mod_double_int (NODE_T * p)
1021  {
1022    A68_LONG_INT i, j;
1023    POP_OBJECT (p, &j, A68_LONG_INT);
1024    POP_OBJECT (p, &i, A68_LONG_INT);
1025    PRELUDE_ERROR (D_ZERO (VALUE (&j)), p, ERROR_DIVISION_BY_ZERO, M_LONG_INT);
1026    PUSH_VALUE (p, double_sdiv (p, VALUE (&i), VALUE (&j), 1), A68_LONG_INT);
1027  }
1028  
1029  //! @brief OP / = (LONG INT, LONG INT) LONG REAL
1030  
1031  void genie_div_double_int (NODE_T * p)
1032  {
1033    A68_LONG_INT i, j;
1034    DOUBLE_NUM_T u, v, w;
1035    POP_OBJECT (p, &j, A68_LONG_INT);
1036    POP_OBJECT (p, &i, A68_LONG_INT);
1037    PRELUDE_ERROR (D_ZERO (VALUE (&j)), p, ERROR_DIVISION_BY_ZERO, M_LONG_INT);
1038    v = double_int_to_double_real (p, VALUE (&j));
1039    u = double_int_to_double_real (p, VALUE (&i));
1040    w.f = u.f / v.f;
1041    PUSH_VALUE (p, w, A68_LONG_REAL);
1042  }
1043  
1044  //! @brief OP ** = (LONG INT, INT) INT
1045  
1046  void genie_pow_double_int_int (NODE_T * p)
1047  {
1048    A68_LONG_INT i;
1049    A68_INT j;
1050    UNSIGNED_T expo, top;
1051    DOUBLE_NUM_T mult, prod;
1052    POP_OBJECT (p, &j, A68_INT);
1053    PRELUDE_ERROR (VALUE (&j) < 0, p, ERROR_EXPONENT_INVALID, M_INT);
1054    top = (UNSIGNED_T) VALUE (&j);
1055    POP_OBJECT (p, &i, A68_LONG_INT);
1056    set_lw (prod, 1);
1057    mult = VALUE (&i);
1058    expo = 1;
1059    while (expo <= top) {
1060      if (expo & top) {
1061        prod = double_smul (p, prod, mult);
1062      }
1063      expo <<= 1;
1064      if (expo <= top) {
1065        mult = double_smul (p, mult, mult);
1066      }
1067    }
1068    PUSH_VALUE (p, prod, A68_LONG_INT);
1069  }
1070  
1071  //! @brief OP - = (LONG REAL) LONG REAL
1072  
1073  void genie_minus_double_real (NODE_T * p)
1074  {
1075    A68_LONG_REAL *u;
1076    POP_OPERAND_ADDRESS (p, u, A68_LONG_REAL);
1077    VALUE (u).f = -(VALUE (u).f);
1078  }
1079  
1080  //! @brief OP ABS = (LONG REAL) LONG REAL
1081  
1082  void genie_abs_double_real (NODE_T * p)
1083  {
1084    A68_LONG_REAL *u;
1085    POP_OPERAND_ADDRESS (p, u, A68_LONG_REAL);
1086    VALUE (u).f = fabsq (VALUE (u).f);
1087  }
1088  
1089  //! @brief OP SIGN = (LONG REAL) INT
1090  
1091  void genie_sign_double_real (NODE_T * p)
1092  {
1093    A68_LONG_REAL u;
1094    POP_OBJECT (p, &u, A68_LONG_REAL);
1095    PUSH_VALUE (p, sign_double_real (VALUE (&u)), A68_INT);
1096  }
1097  
1098  //! @brief OP ** = (LONG REAL, INT) INT
1099  
1100  void genie_pow_double_real_int (NODE_T * p)
1101  {
1102    A68_LONG_REAL z;
1103    A68_INT j;
1104    INT_T top;
1105    UNSIGNED_T expo;
1106    DOUBLE_NUM_T mult, prod;
1107    int negative;
1108    POP_OBJECT (p, &j, A68_INT);
1109    top = (UNSIGNED_T) VALUE (&j);
1110    POP_OBJECT (p, &z, A68_LONG_INT);
1111    prod.f = 1.0q;
1112    mult.f = VALUE (&z).f;
1113    if (top < 0) {
1114      top = -top;
1115      negative = A68_TRUE;
1116    } else {
1117      negative = A68_FALSE;
1118    }
1119    expo = 1;
1120    while (expo <= top) {
1121      if (expo & top) {
1122        prod.f = prod.f * mult.f;
1123        CHECK_DOUBLE_REAL (p, prod.f);
1124      }
1125      expo <<= 1;
1126      if (expo <= top) {
1127        mult.f = mult.f * mult.f;
1128        CHECK_DOUBLE_REAL (p, mult.f);
1129      }
1130    }
1131    if (negative) {
1132      prod.f = 1.0q / prod.f;
1133    }
1134    PUSH_VALUE (p, prod, A68_LONG_REAL);
1135  }
1136  
1137  //! @brief OP ** = (LONG REAL, LONG REAL) LONG REAL
1138  
1139  void genie_pow_double_real (NODE_T * p)
1140  {
1141    A68_LONG_REAL x, y;
1142    DOUBLE_T z = 0.0q;
1143    POP_OBJECT (p, &y, A68_LONG_REAL);
1144    POP_OBJECT (p, &x, A68_LONG_REAL);
1145    errno = 0;
1146    PRELUDE_ERROR (VALUE (&x).f < 0.0q, p, ERROR_INVALID_ARGUMENT, M_LONG_REAL);
1147    if (VALUE (&x).f == 0.0q) {
1148      if (VALUE (&y).f < 0.0q) {
1149        errno = ERANGE;
1150        MATH_RTE (p, errno != 0, M_LONG_REAL, NO_TEXT);
1151      } else {
1152        z = (VALUE (&y).f == 0.0q ? 1.0q : 0.0q);
1153      }
1154    } else {
1155      z = expq (VALUE (&y).f * logq (VALUE (&x).f));
1156      MATH_RTE (p, errno != 0, M_LONG_REAL, NO_TEXT);
1157    }
1158    PUSH_VALUE (p, dble (z), A68_LONG_REAL);
1159  }
1160  
1161  //! @brief OP + = (LONG REAL, LONG REAL) LONG REAL
1162  
1163  void genie_add_double_real (NODE_T * p)
1164  {
1165    A68_LONG_REAL u, v;
1166    DOUBLE_NUM_T w;
1167    POP_OBJECT (p, &v, A68_LONG_REAL);
1168    POP_OBJECT (p, &u, A68_LONG_REAL);
1169    w.f = VALUE (&u).f + VALUE (&v).f;
1170    CHECK_DOUBLE_REAL (p, w.f);
1171    PUSH_VALUE (p, w, A68_LONG_REAL);
1172  }
1173  
1174  //! @brief OP - = (LONG REAL, LONG REAL) LONG REAL
1175  
1176  void genie_sub_double_real (NODE_T * p)
1177  {
1178    A68_LONG_REAL u, v;
1179    DOUBLE_NUM_T w;
1180    POP_OBJECT (p, &v, A68_LONG_REAL);
1181    POP_OBJECT (p, &u, A68_LONG_REAL);
1182    w.f = VALUE (&u).f - VALUE (&v).f;
1183    CHECK_DOUBLE_REAL (p, w.f);
1184    PUSH_VALUE (p, w, A68_LONG_REAL);
1185  }
1186  
1187  //! @brief OP * = (LONG REAL, LONG REAL) LONG REAL
1188  
1189  void genie_mul_double_real (NODE_T * p)
1190  {
1191    A68_LONG_REAL u, v;
1192    DOUBLE_NUM_T w;
1193    POP_OBJECT (p, &v, A68_LONG_REAL);
1194    POP_OBJECT (p, &u, A68_LONG_REAL);
1195    w.f = VALUE (&u).f * VALUE (&v).f;
1196    CHECK_DOUBLE_REAL (p, w.f);
1197    PUSH_VALUE (p, w, A68_LONG_REAL);
1198  }
1199  
1200  //! @brief OP / = (LONG REAL, LONG REAL) LONG REAL
1201  
1202  void genie_over_double_real (NODE_T * p)
1203  {
1204    A68_LONG_REAL u, v;
1205    DOUBLE_NUM_T w;
1206    POP_OBJECT (p, &v, A68_LONG_REAL);
1207    POP_OBJECT (p, &u, A68_LONG_REAL);
1208    PRELUDE_ERROR (VALUE (&v).f == 0.0q, p, ERROR_DIVISION_BY_ZERO, M_LONG_REAL);
1209    w.f = VALUE (&u).f / VALUE (&v).f;
1210    PUSH_VALUE (p, w, A68_LONG_REAL);
1211  }
1212  
1213  //! @brief OP +:= = (REF LONG INT, LONG INT) REF LONG INT
1214  
1215  void genie_plusab_double_int (NODE_T * p)
1216  {
1217    genie_f_and_becomes (p, M_REF_LONG_INT, genie_add_double_int);
1218  }
1219  
1220  //! @brief OP -:= = (REF LONG INT, LONG INT) REF LONG INT
1221  
1222  void genie_minusab_double_int (NODE_T * p)
1223  {
1224    genie_f_and_becomes (p, M_REF_LONG_INT, genie_sub_double_int);
1225  }
1226  
1227  //! @brief OP *:= = (REF LONG INT, LONG INT) REF LONG INT
1228  
1229  void genie_timesab_double_int (NODE_T * p)
1230  {
1231    genie_f_and_becomes (p, M_REF_LONG_INT, genie_mul_double_int);
1232  }
1233  
1234  //! @brief OP %:= = (REF LONG INT, LONG INT) REF LONG INT
1235  
1236  void genie_overab_double_int (NODE_T * p)
1237  {
1238    genie_f_and_becomes (p, M_REF_LONG_INT, genie_over_double_int);
1239  }
1240  
1241  //! @brief OP %*:= = (REF LONG INT, LONG INT) REF LONG INT
1242  
1243  void genie_modab_double_int (NODE_T * p)
1244  {
1245    genie_f_and_becomes (p, M_REF_LONG_INT, genie_mod_double_int);
1246  }
1247  
1248  //! @brief OP +:= = (REF LONG REAL, LONG REAL) REF LONG REAL
1249  
1250  void genie_plusab_double_real (NODE_T * p)
1251  {
1252    genie_f_and_becomes (p, M_REF_LONG_REAL, genie_add_double_real);
1253  }
1254  
1255  //! @brief OP -:= = (REF LONG REAL, LONG REAL) REF LONG REAL
1256  
1257  void genie_minusab_double_real (NODE_T * p)
1258  {
1259    genie_f_and_becomes (p, M_REF_LONG_REAL, genie_sub_double_real);
1260  }
1261  
1262  //! @brief OP *:= = (REF LONG REAL, LONG REAL) REF LONG REAL
1263  
1264  void genie_timesab_double_real (NODE_T * p)
1265  {
1266    genie_f_and_becomes (p, M_REF_LONG_REAL, genie_mul_double_real);
1267  }
1268  
1269  //! @brief OP /:= = (REF LONG REAL, LONG REAL) REF LONG REAL
1270  
1271  void genie_divab_double_real (NODE_T * p)
1272  {
1273    genie_f_and_becomes (p, M_REF_LONG_REAL, genie_over_double_real);
1274  }
1275  
1276  // OP (LONG INT, LONG INT) BOOL.
1277  
1278  #define A68_CMP_INT(n, OP)\
1279  void n (NODE_T * p) {\
1280    A68_LONG_INT i, j;\
1281    int k;\
1282    POP_OBJECT (p, &j, A68_LONG_INT);\
1283    POP_OBJECT (p, &i, A68_LONG_INT);\
1284    k = sign_double_int (double_ssub (p, VALUE (&i), VALUE (&j)));\
1285    PUSH_VALUE (p, (BOOL_T) (k OP 0), A68_BOOL);\
1286    }
1287  A68_CMP_INT (genie_eq_double_int, ==)
1288    A68_CMP_INT (genie_ne_double_int, !=)
1289    A68_CMP_INT (genie_lt_double_int, <)
1290    A68_CMP_INT (genie_gt_double_int, >)
1291    A68_CMP_INT (genie_le_double_int, <=)
1292    A68_CMP_INT (genie_ge_double_int, >=)
1293  // OP (LONG REAL, LONG REAL) BOOL.
1294  #define A68_CMP_REAL(n, OP)\
1295  void n (NODE_T * p) {\
1296    A68_LONG_REAL i, j;\
1297    POP_OBJECT (p, &j, A68_LONG_REAL);\
1298    POP_OBJECT (p, &i, A68_LONG_REAL);\
1299    PUSH_VALUE (p, (BOOL_T) (VALUE (&i).f OP VALUE (&j).f), A68_BOOL);\
1300    }
1301    A68_CMP_REAL (genie_eq_double_real, ==)
1302    A68_CMP_REAL (genie_ne_double_real, !=)
1303    A68_CMP_REAL (genie_lt_double_real, <)
1304    A68_CMP_REAL (genie_gt_double_real, >)
1305    A68_CMP_REAL (genie_le_double_real, <=)
1306    A68_CMP_REAL (genie_ge_double_real, >=)
1307  //! @brief OP NOT = (LONG BITS) LONG BITS
1308       void genie_not_double_bits (NODE_T * p)
1309  {
1310    A68_LONG_BITS i;
1311    DOUBLE_NUM_T w;
1312    POP_OBJECT (p, &i, A68_LONG_BITS);
1313    HW (w) = ~HW (VALUE (&i));
1314    LW (w) = ~LW (VALUE (&i));
1315    PUSH_VALUE (p, w, A68_LONG_BITS);
1316  }
1317  
1318  //! @brief OP = = (LONG BITS, LONG BITS) BOOL.
1319  
1320  void genie_eq_double_bits (NODE_T * p)
1321  {
1322    A68_LONG_BITS i, j;
1323    BOOL_T u, v;
1324    POP_OBJECT (p, &j, A68_LONG_BITS);
1325    POP_OBJECT (p, &i, A68_LONG_BITS);
1326    u = HW (VALUE (&i)) == HW (VALUE (&j));
1327    v = LW (VALUE (&i)) == LW (VALUE (&j));
1328    PUSH_VALUE (p, (BOOL_T) (u & v ? A68_TRUE : A68_FALSE), A68_BOOL);
1329  }
1330  
1331  //! @brief OP ~= = (LONG BITS, LONG BITS) BOOL.
1332  
1333  void genie_ne_double_bits (NODE_T * p)
1334  {
1335    A68_LONG_BITS i, j;
1336    BOOL_T u, v;
1337    POP_OBJECT (p, &j, A68_LONG_BITS);    // (i ~= j) == ~ (i = j) 
1338    POP_OBJECT (p, &i, A68_LONG_BITS);
1339    u = HW (VALUE (&i)) == HW (VALUE (&j));
1340    v = LW (VALUE (&i)) == LW (VALUE (&j));
1341    PUSH_VALUE (p, (BOOL_T) (u & v ? A68_FALSE : A68_TRUE), A68_BOOL);
1342  }
1343  
1344  //! @brief OP <= = (LONG BITS, LONG BITS) BOOL
1345  
1346  void genie_le_double_bits (NODE_T * p)
1347  {
1348    A68_LONG_BITS i, j;
1349    BOOL_T u, v;
1350    POP_OBJECT (p, &j, A68_LONG_BITS);
1351    POP_OBJECT (p, &i, A68_LONG_BITS);
1352    u = (HW (VALUE (&i)) | HW (VALUE (&j))) == HW (VALUE (&j));
1353    v = (LW (VALUE (&i)) | LW (VALUE (&j))) == LW (VALUE (&j));
1354    PUSH_VALUE (p, (BOOL_T) (u & v ? A68_TRUE : A68_FALSE), A68_BOOL);
1355  }
1356  
1357  //! @brief OP > = (LONG BITS, LONG BITS) BOOL
1358  
1359  void genie_gt_double_bits (NODE_T * p)
1360  {
1361    A68_LONG_BITS i, j;
1362    BOOL_T u, v;
1363    POP_OBJECT (p, &j, A68_LONG_BITS);    // (i > j) == ! (i <= j)
1364    POP_OBJECT (p, &i, A68_LONG_BITS);
1365    u = (HW (VALUE (&i)) | HW (VALUE (&j))) == HW (VALUE (&j));
1366    v = (LW (VALUE (&i)) | LW (VALUE (&j))) == LW (VALUE (&j));
1367    PUSH_VALUE (p, (BOOL_T) (u & v ? A68_FALSE : A68_TRUE), A68_BOOL);
1368  }
1369  
1370  //! @brief OP >= = (LONG BITS, LONG BITS) BOOL
1371  
1372  void genie_ge_double_bits (NODE_T * p)
1373  {
1374    A68_LONG_BITS i, j;
1375    BOOL_T u, v;
1376    POP_OBJECT (p, &j, A68_LONG_BITS);    // (i >= j) == (j <= i)
1377    POP_OBJECT (p, &i, A68_LONG_BITS);
1378    u = (HW (VALUE (&i)) | HW (VALUE (&j))) == HW (VALUE (&i));
1379    v = (LW (VALUE (&i)) | LW (VALUE (&j))) == LW (VALUE (&i));
1380    PUSH_VALUE (p, (BOOL_T) (u & v ? A68_TRUE : A68_FALSE), A68_BOOL);
1381  }
1382  
1383  //! @brief OP < = (LONG BITS, LONG BITS) BOOL
1384  
1385  void genie_lt_double_bits (NODE_T * p)
1386  {
1387    A68_LONG_BITS i, j;
1388    BOOL_T u, v;
1389    POP_OBJECT (p, &j, A68_LONG_BITS);    // (i < j) == ! (i >= j)
1390    POP_OBJECT (p, &i, A68_LONG_BITS);
1391    u = (HW (VALUE (&i)) | HW (VALUE (&j))) == HW (VALUE (&i));
1392    v = (LW (VALUE (&i)) | LW (VALUE (&j))) == LW (VALUE (&i));
1393    PUSH_VALUE (p, (BOOL_T) (u & v ? A68_FALSE : A68_TRUE), A68_BOOL);
1394  }
1395  
1396  //! @brief PROC long bits pack = ([] BOOL) BITS
1397  
1398  void genie_double_bits_pack (NODE_T * p)
1399  {
1400    A68_REF z;
1401    DOUBLE_NUM_T w;
1402    A68_ARRAY *arr;
1403    A68_TUPLE *tup;
1404    int size;
1405    POP_REF (p, &z);
1406    CHECK_REF (p, z, M_ROW_BOOL);
1407    GET_DESCRIPTOR (arr, tup, &z);
1408    size = ROW_SIZE (tup);
1409    PRELUDE_ERROR (size < 0 || size > BITS_WIDTH, p, ERROR_OUT_OF_BOUNDS, M_ROW_BOOL);
1410    set_lw (w, 0x0);
1411    if (ROW_SIZE (tup) > 0) {
1412      UNSIGNED_T bit = 0x0;
1413      int k, n = 0;
1414      BYTE_T *base = DEREF (BYTE_T, &ARRAY (arr));
1415      for (k = UPB (tup); k >= LWB (tup); k--) {
1416        A68_BOOL *boo = (A68_BOOL *) & (base[INDEX_1_DIM (arr, tup, k)]);
1417        CHECK_INIT (p, INITIALISED (boo), M_BOOL);
1418        if (n == 0 || n == BITS_WIDTH) {
1419          bit = 0x1;
1420        }
1421        if (VALUE (boo)) {
1422          if (n > BITS_WIDTH) {
1423            LW (w) |= bit;
1424          } else {
1425            HW (w) |= bit;
1426          };
1427        }
1428        n++;
1429        bit <<= 1;
1430      }
1431    }
1432    PUSH_VALUE (p, w, A68_LONG_BITS);
1433  }
1434  
1435  //! @brief OP AND = (LONG BITS, LONG BITS) LONG BITS
1436  
1437  void genie_and_double_bits (NODE_T * p)
1438  {
1439    A68_LONG_BITS i, j;
1440    DOUBLE_NUM_T w;
1441    POP_OBJECT (p, &j, A68_LONG_BITS);
1442    POP_OBJECT (p, &i, A68_LONG_BITS);
1443    HW (w) = HW (VALUE (&i)) & HW (VALUE (&j));
1444    LW (w) = LW (VALUE (&i)) & LW (VALUE (&j));
1445    PUSH_VALUE (p, w, A68_LONG_BITS);
1446  }
1447  
1448  //! @brief OP OR = (LONG BITS, LONG BITS) LONG BITS
1449  
1450  void genie_or_double_bits (NODE_T * p)
1451  {
1452    A68_LONG_BITS i, j;
1453    DOUBLE_NUM_T w;
1454    POP_OBJECT (p, &j, A68_LONG_BITS);
1455    POP_OBJECT (p, &i, A68_LONG_BITS);
1456    HW (w) = HW (VALUE (&i)) | HW (VALUE (&j));
1457    LW (w) = LW (VALUE (&i)) | LW (VALUE (&j));
1458    PUSH_VALUE (p, w, A68_LONG_BITS);
1459  }
1460  
1461  //! @brief OP XOR = (LONG BITS, LONG BITS) LONG BITS
1462  
1463  void genie_xor_double_bits (NODE_T * p)
1464  {
1465    A68_LONG_BITS i, j;
1466    DOUBLE_NUM_T w;
1467    POP_OBJECT (p, &j, A68_LONG_BITS);
1468    POP_OBJECT (p, &i, A68_LONG_BITS);
1469    HW (w) = HW (VALUE (&i)) ^ HW (VALUE (&j));
1470    LW (w) = LW (VALUE (&i)) ^ LW (VALUE (&j));
1471    PUSH_VALUE (p, w, A68_LONG_BITS);
1472  }
1473  
1474  //! @brief OP + = (LONG BITS, LONG BITS) LONG BITS
1475  
1476  void genie_add_double_bits (NODE_T * p)
1477  {
1478    A68_LONG_BITS i, j;
1479    DOUBLE_NUM_T w;
1480    POP_OBJECT (p, &j, A68_LONG_BITS);
1481    POP_OBJECT (p, &i, A68_LONG_BITS);
1482    add_double (p, M_LONG_BITS, w, VALUE (&i), VALUE (&j));
1483    PUSH_VALUE (p, w, A68_LONG_BITS);
1484  }
1485  
1486  //! @brief OP - = (LONG BITS, LONG BITS) LONG BITS
1487  
1488  void genie_sub_double_bits (NODE_T * p)
1489  {
1490    A68_LONG_BITS i, j;
1491    DOUBLE_NUM_T w;
1492    POP_OBJECT (p, &j, A68_LONG_BITS);
1493    POP_OBJECT (p, &i, A68_LONG_BITS);
1494    sub_double (p, M_LONG_BITS, w, VALUE (&i), VALUE (&j));
1495    PUSH_VALUE (p, w, A68_LONG_BITS);
1496  }
1497  
1498  //! @brief OP * = (LONG BITS, LONG BITS) LONG BITS
1499  
1500  void genie_times_double_bits (NODE_T * p)
1501  {
1502    A68_LONG_BITS i, j;
1503    DOUBLE_NUM_T w;
1504    POP_OBJECT (p, &j, A68_LONG_BITS);
1505    POP_OBJECT (p, &i, A68_LONG_BITS);
1506    w = double_umul (p, M_LONG_BITS, VALUE (&i), VALUE (&j));
1507    PUSH_VALUE (p, w, A68_LONG_BITS);
1508  }
1509  
1510  //! @brief OP OVER = (LONG BITS, LONG BITS) LONG BITS
1511  
1512  void genie_over_double_bits (NODE_T * p)
1513  {
1514    A68_LONG_BITS i, j;
1515    DOUBLE_NUM_T w;
1516    POP_OBJECT (p, &j, A68_LONG_BITS);
1517    POP_OBJECT (p, &i, A68_LONG_BITS);
1518    w = double_udiv (p, M_LONG_BITS, VALUE (&i), VALUE (&j), 0);
1519    PUSH_VALUE (p, w, A68_LONG_BITS);
1520  }
1521  
1522  //! @brief OP MOD = (LONG BITS, LONG BITS) LONG BITS
1523  
1524  void genie_mod_double_bits (NODE_T * p)
1525  {
1526    A68_LONG_BITS i, j;
1527    DOUBLE_NUM_T w;
1528    POP_OBJECT (p, &j, A68_LONG_BITS);
1529    POP_OBJECT (p, &i, A68_LONG_BITS);
1530    w = double_udiv (p, M_LONG_BITS, VALUE (&i), VALUE (&j), 1);
1531    PUSH_VALUE (p, w, A68_LONG_BITS);
1532  }
1533  
1534  //! @brief OP ELEM = (INT, LONG BITS) BOOL
1535  
1536  void genie_elem_double_bits (NODE_T * p)
1537  {
1538    A68_LONG_BITS j;
1539    A68_INT i;
1540    int k, n;
1541    UNSIGNED_T mask = 0x1, *w;
1542    POP_OBJECT (p, &j, A68_LONG_BITS);
1543    POP_OBJECT (p, &i, A68_INT);
1544    k = VALUE (&i);
1545    PRELUDE_ERROR (k < 1 || k > LONG_BITS_WIDTH, p, ERROR_OUT_OF_BOUNDS, M_INT);
1546    if (k <= BITS_WIDTH) {
1547      w = &(LW (VALUE (&j)));
1548    } else {
1549      w = &(HW (VALUE (&j)));
1550      k -= BITS_WIDTH;
1551    }
1552    for (n = 0; n < (BITS_WIDTH - VALUE (&i)); n++) {
1553      mask = mask << 1;
1554    }
1555    PUSH_VALUE (p, (BOOL_T) ((*w & mask) ? A68_TRUE : A68_FALSE), A68_BOOL);
1556  }
1557  
1558  //! @brief OP SET = (INT, LONG BITS) LONG BITS
1559  
1560  void genie_set_double_bits (NODE_T * p)
1561  {
1562    A68_LONG_BITS j;
1563    A68_INT i;
1564    int k, n;
1565    UNSIGNED_T mask = 0x1, *w;
1566    POP_OBJECT (p, &j, A68_LONG_BITS);
1567    POP_OBJECT (p, &i, A68_INT);
1568    k = VALUE (&i);
1569    PRELUDE_ERROR (k < 1 || k > LONG_BITS_WIDTH, p, ERROR_OUT_OF_BOUNDS, M_INT);
1570    if (k <= BITS_WIDTH) {
1571      w = &(LW (VALUE (&j)));
1572    } else {
1573      w = &(HW (VALUE (&j)));
1574      k -= BITS_WIDTH;
1575    }
1576    for (n = 0; n < (BITS_WIDTH - VALUE (&i)); n++) {
1577      mask = mask << 1;
1578    }
1579    (*w) |= mask;
1580    PUSH_OBJECT (p, j, A68_LONG_BITS);
1581  }
1582  
1583  //! @brief OP CLEAR = (INT, LONG BITS) LONG BITS
1584  
1585  void genie_clear_double_bits (NODE_T * p)
1586  {
1587    A68_LONG_BITS j;
1588    A68_INT i;
1589    int k, n;
1590    UNSIGNED_T mask = 0x1, *w;
1591    POP_OBJECT (p, &j, A68_LONG_BITS);
1592    POP_OBJECT (p, &i, A68_INT);
1593    k = VALUE (&i);
1594    PRELUDE_ERROR (k < 1 || k > LONG_BITS_WIDTH, p, ERROR_OUT_OF_BOUNDS, M_INT);
1595    if (k <= BITS_WIDTH) {
1596      w = &(LW (VALUE (&j)));
1597    } else {
1598      w = &(HW (VALUE (&j)));
1599      k -= BITS_WIDTH;
1600    }
1601    for (n = 0; n < (BITS_WIDTH - VALUE (&i)); n++) {
1602      mask = mask << 1;
1603    }
1604    (*w) &= ~mask;
1605    PUSH_OBJECT (p, j, A68_LONG_BITS);
1606  }
1607  
1608  //! @brief OP SHL = (LONG BITS, INT) LONG BITS
1609  
1610  void genie_shl_double_bits (NODE_T * p)
1611  {
1612    A68_LONG_BITS i;
1613    A68_INT j;
1614    DOUBLE_NUM_T *w;
1615    int k, n;
1616    POP_OBJECT (p, &j, A68_INT);
1617    POP_OBJECT (p, &i, A68_LONG_BITS);
1618    w = &VALUE (&i);
1619    k = VALUE (&j);
1620    if (VALUE (&j) >= 0) {
1621      for (n = 0; n < k; n++) {
1622        UNSIGNED_T carry = ((LW (*w) & D_SIGN) ? 0x1 : 0x0);
1623        PRELUDE_ERROR (MODCHK (p, M_LONG_BITS, HW (*w) | D_SIGN), p, ERROR_MATH, M_LONG_BITS);
1624        HW (*w) = (HW (*w) << 1) | carry;
1625        LW (*w) = (LW (*w) << 1);
1626      }
1627    } else {
1628      k = -k;
1629      for (n = 0; n < k; n++) {
1630        UNSIGNED_T carry = ((HW (*w) & 0x1) ? D_SIGN : 0x0);
1631        HW (*w) = (HW (*w) >> 1);
1632        LW (*w) = (LW (*w) >> 1) | carry;
1633      }
1634    }
1635    PUSH_OBJECT (p, i, A68_LONG_BITS);
1636  }
1637  
1638  //! @brief OP SHR = (LONG BITS, INT) LONG BITS
1639  
1640  void genie_shr_double_bits (NODE_T * p)
1641  {
1642    A68_INT *j;
1643    POP_OPERAND_ADDRESS (p, j, A68_INT);
1644    VALUE (j) = -VALUE (j);
1645    genie_shl_double_bits (p);    // Conform RR
1646  }
1647  
1648  //! @brief OP ROL = (LONG BITS, INT) LONG BITS
1649  
1650  void genie_rol_double_bits (NODE_T * p)
1651  {
1652    A68_LONG_BITS i;
1653    A68_INT j;
1654    DOUBLE_NUM_T *w = &VALUE (&i);
1655    int k, n;
1656    POP_OBJECT (p, &j, A68_INT);
1657    POP_OBJECT (p, &i, A68_LONG_BITS);
1658    k = VALUE (&j);
1659    if (k >= 0) {
1660      for (n = 0; n < k; n++) {
1661        UNSIGNED_T carry = ((HW (*w) & D_SIGN) ? 0x1 : 0x0);
1662        UNSIGNED_T carry_between = ((LW (*w) & D_SIGN) ? 0x1 : 0x0);
1663        HW (*w) = (HW (*w) << 1) | carry_between;
1664        LW (*w) = (LW (*w) << 1) | carry;
1665      }
1666    } else {
1667      k = -k;
1668      for (n = 0; n < k; n++) {
1669        UNSIGNED_T carry = ((LW (*w) & 0x1) ? D_SIGN : 0x0);
1670        UNSIGNED_T carry_between = ((HW (*w) & 0x1) ? D_SIGN : 0x0);
1671        HW (*w) = (HW (*w) >> 1) | carry;
1672        LW (*w) = (LW (*w) >> 1) | carry_between;
1673      }
1674    }
1675    PUSH_OBJECT (p, i, A68_LONG_BITS);
1676  }
1677  
1678  //! @brief OP ROR = (LONG BITS, INT) LONG BITS
1679  
1680  void genie_ror_double_bits (NODE_T * p)
1681  {
1682    A68_INT *j;
1683    POP_OPERAND_ADDRESS (p, j, A68_INT);
1684    VALUE (j) = -VALUE (j);
1685    genie_rol_double_bits (p);    // Conform RR
1686  }
1687  
1688  //! @brief OP BIN = (LONG INT) LONG BITS
1689  
1690  void genie_bin_double_int (NODE_T * p)
1691  {
1692    A68_LONG_INT i;
1693    POP_OBJECT (p, &i, A68_LONG_INT);
1694  // RR does not convert negative numbers
1695    if (D_NEG (VALUE (&i))) {
1696      errno = EDOM;
1697      diagnostic (A68_RUNTIME_ERROR, p, ERROR_OUT_OF_BOUNDS, M_BITS);
1698      exit_genie (p, A68_RUNTIME_ERROR);
1699    }
1700    PUSH_OBJECT (p, i, A68_LONG_BITS);
1701  }
1702  
1703  //! @brief OP +* = (LONG REAL, LONG REAL) LONG COMPLEX
1704  
1705  void genie_i_double_compl (NODE_T * p)
1706  {
1707    (void) p;
1708  }
1709  
1710  //! @brief OP SHORTEN = (LONG COMPLEX) COMPLEX
1711  
1712  void genie_shorten_double_compl_to_complex (NODE_T * p)
1713  {
1714    A68_LONG_REAL re, im;
1715    REAL_T w;
1716    POP_OBJECT (p, &im, A68_LONG_REAL);
1717    POP_OBJECT (p, &re, A68_LONG_REAL);
1718    w = VALUE (&re).f;
1719    PUSH_VALUE (p, w, A68_REAL);
1720    w = VALUE (&im).f;
1721    PUSH_VALUE (p, w, A68_REAL);
1722  }
1723  
1724  //! @brief OP LENG = (LONG COMPLEX) LONG LONG COMPLEX
1725  
1726  void genie_lengthen_double_compl_to_long_mp_complex (NODE_T * p)
1727  {
1728    int digs = DIGITS (M_LONG_LONG_REAL);
1729    A68_LONG_REAL re, im;
1730    POP_OBJECT (p, &im, A68_LONG_REAL);
1731    POP_OBJECT (p, &re, A68_LONG_REAL);
1732    MP_T *z = nil_mp (p, digs);
1733    (void) double_real_to_mp (p, z, VALUE (&re).f, digs);
1734    MP_STATUS (z) = (MP_T) INIT_MASK;
1735    z = nil_mp (p, digs);
1736    (void) double_real_to_mp (p, z, VALUE (&im).f, digs);
1737    MP_STATUS (z) = (MP_T) INIT_MASK;
1738  }
1739  
1740  //! @brief OP +* = (LONG INT, LONG INT) LONG COMPLEX
1741  
1742  void genie_i_int_double_compl (NODE_T * p)
1743  {
1744    A68_LONG_INT re, im;
1745    POP_OBJECT (p, &im, A68_LONG_INT);
1746    POP_OBJECT (p, &re, A68_LONG_INT);
1747    PUSH_VALUE (p, double_int_to_double_real (p, VALUE (&re)), A68_LONG_REAL);
1748    PUSH_VALUE (p, double_int_to_double_real (p, VALUE (&im)), A68_LONG_REAL);
1749  }
1750  
1751  //! @brief OP RE = (LONG COMPLEX) LONG REAL
1752  
1753  void genie_re_double_compl (NODE_T * p)
1754  {
1755    DECREMENT_STACK_POINTER (p, SIZE (M_LONG_REAL));
1756  }
1757  
1758  //! @brief OP IM = (LONG COMPLEX) LONG REAL
1759  
1760  void genie_im_double_compl (NODE_T * p)
1761  {
1762    A68_LONG_REAL re, im;
1763    POP_OBJECT (p, &im, A68_LONG_REAL);
1764    POP_OBJECT (p, &re, A68_LONG_REAL);
1765    PUSH_OBJECT (p, im, A68_LONG_REAL);
1766  }
1767  
1768  //! @brief OP - = (LONG COMPLEX) LONG COMPLEX
1769  
1770  void genie_minus_double_compl (NODE_T * p)
1771  {
1772    A68_LONG_REAL re, im;
1773    POP_OBJECT (p, &im, A68_LONG_REAL);
1774    POP_OBJECT (p, &re, A68_LONG_REAL);
1775    VALUE (&re).f = -VALUE (&re).f;
1776    VALUE (&im).f = -VALUE (&im).f;
1777    PUSH_OBJECT (p, im, A68_LONG_REAL);
1778    PUSH_OBJECT (p, re, A68_LONG_REAL);
1779  }
1780  
1781  //! @brief OP ABS = (LONG COMPLEX) LONG REAL
1782  
1783  void genie_abs_double_compl (NODE_T * p)
1784  {
1785    A68_LONG_REAL re, im;
1786    POP_LONG_COMPLEX (p, &re, &im);
1787    PUSH_VALUE (p, dble (a68_double_hypot (VALUE (&re).f, VALUE (&im).f)), A68_LONG_REAL);
1788  }
1789  
1790  //! @brief OP ARG = (LONG COMPLEX) LONG REAL
1791  
1792  void genie_arg_double_compl (NODE_T * p)
1793  {
1794    A68_LONG_REAL re, im;
1795    POP_LONG_COMPLEX (p, &re, &im);
1796    PRELUDE_ERROR (VALUE (&re).f == 0.0q && VALUE (&im).f == 0.0q, p, ERROR_INVALID_ARGUMENT, M_LONG_COMPLEX);
1797    PUSH_VALUE (p, dble (atan2q (VALUE (&im).f, VALUE (&re).f)), A68_LONG_REAL);
1798  }
1799  
1800  //! @brief OP CONJ = (LONG COMPLEX) LONG COMPLEX
1801  
1802  void genie_conj_double_compl (NODE_T * p)
1803  {
1804    A68_LONG_REAL im;
1805    POP_OBJECT (p, &im, A68_LONG_REAL);
1806    VALUE (&im).f = -VALUE (&im).f;
1807    PUSH_OBJECT (p, im, A68_LONG_REAL);
1808  }
1809  
1810  //! @brief OP + = (COMPLEX, COMPLEX) COMPLEX
1811  
1812  void genie_add_double_compl (NODE_T * p)
1813  {
1814    A68_LONG_REAL re_x, im_x, re_y, im_y;
1815    POP_LONG_COMPLEX (p, &re_y, &im_y);
1816    POP_LONG_COMPLEX (p, &re_x, &im_x);
1817    VALUE (&re_x).f += VALUE (&re_y).f;
1818    VALUE (&im_x).f += VALUE (&im_y).f;
1819    CHECK_DOUBLE_COMPLEX (p, VALUE (&im_x).f, VALUE (&im_y).f);
1820    PUSH_OBJECT (p, re_x, A68_LONG_REAL);
1821    PUSH_OBJECT (p, im_x, A68_LONG_REAL);
1822  }
1823  
1824  //! @brief OP - = (COMPLEX, COMPLEX) COMPLEX
1825  
1826  void genie_sub_double_compl (NODE_T * p)
1827  {
1828    A68_LONG_REAL re_x, im_x, re_y, im_y;
1829    POP_LONG_COMPLEX (p, &re_y, &im_y);
1830    POP_LONG_COMPLEX (p, &re_x, &im_x);
1831    VALUE (&re_x).f -= VALUE (&re_y).f;
1832    VALUE (&im_x).f -= VALUE (&im_y).f;
1833    CHECK_DOUBLE_COMPLEX (p, VALUE (&im_x).f, VALUE (&im_y).f);
1834    PUSH_OBJECT (p, re_x, A68_LONG_REAL);
1835    PUSH_OBJECT (p, im_x, A68_LONG_REAL);
1836  }
1837  
1838  //! @brief OP * = (COMPLEX, COMPLEX) COMPLEX
1839  
1840  void genie_mul_double_compl (NODE_T * p)
1841  {
1842    A68_LONG_REAL re_x, im_x, re_y, im_y;
1843    DOUBLE_T re, im;
1844    POP_LONG_COMPLEX (p, &re_y, &im_y);
1845    POP_LONG_COMPLEX (p, &re_x, &im_x);
1846    re = VALUE (&re_x).f * VALUE (&re_y).f - VALUE (&im_x).f * VALUE (&im_y).f;
1847    im = VALUE (&im_x).f * VALUE (&re_y).f + VALUE (&re_x).f * VALUE (&im_y).f;
1848    CHECK_DOUBLE_COMPLEX (p, VALUE (&im_x).f, VALUE (&im_y).f);
1849    PUSH_VALUE (p, dble (re), A68_LONG_REAL);
1850    PUSH_VALUE (p, dble (im), A68_LONG_REAL);
1851  }
1852  
1853  //! @brief OP / = (COMPLEX, COMPLEX) COMPLEX
1854  
1855  void genie_div_double_compl (NODE_T * p)
1856  {
1857    A68_LONG_REAL re_x, im_x, re_y, im_y;
1858    DOUBLE_T re = 0.0, im = 0.0;
1859    POP_LONG_COMPLEX (p, &re_y, &im_y);
1860    POP_LONG_COMPLEX (p, &re_x, &im_x);
1861    PRELUDE_ERROR (VALUE (&re_y).f == 0.0q && VALUE (&im_y).f == 0.0q, p, ERROR_DIVISION_BY_ZERO, M_LONG_COMPLEX);
1862    if (ABSQ (VALUE (&re_y).f) >= ABSQ (VALUE (&im_y).f)) {
1863      DOUBLE_T r = VALUE (&im_y).f / VALUE (&re_y).f, den = VALUE (&re_y).f + r * VALUE (&im_y).f;
1864      re = (VALUE (&re_x).f + r * VALUE (&im_x).f) / den;
1865      im = (VALUE (&im_x).f - r * VALUE (&re_x).f) / den;
1866    } else {
1867      DOUBLE_T r = VALUE (&re_y).f / VALUE (&im_y).f, den = VALUE (&im_y).f + r * VALUE (&re_y).f;
1868      re = (VALUE (&re_x).f * r + VALUE (&im_x).f) / den;
1869      im = (VALUE (&im_x).f * r - VALUE (&re_x).f) / den;
1870    }
1871    PUSH_VALUE (p, dble (re), A68_LONG_REAL);
1872    PUSH_VALUE (p, dble (im), A68_LONG_REAL);
1873  }
1874  
1875  //! @brief OP ** = (LONG COMPLEX, INT) LONG COMPLEX
1876  
1877  void genie_pow_double_compl_int (NODE_T * p)
1878  {
1879    A68_LONG_REAL re_x, im_x;
1880    DOUBLE_T re_y, im_y, re_z, im_z;
1881    A68_INT j;
1882    INT_T expo;
1883    BOOL_T negative;
1884    POP_OBJECT (p, &j, A68_INT);
1885    POP_LONG_COMPLEX (p, &re_x, &im_x);
1886    re_z = 1.0q;
1887    im_z = 0.0q;
1888    re_y = VALUE (&re_x).f;
1889    im_y = VALUE (&im_x).f;
1890    expo = 1;
1891    negative = (BOOL_T) (VALUE (&j) < 0);
1892    if (negative) {
1893      VALUE (&j) = -VALUE (&j);
1894    }
1895    while ((UNSIGNED_T) expo <= (UNSIGNED_T) (VALUE (&j))) {
1896      DOUBLE_T z;
1897      if (expo & VALUE (&j)) {
1898        z = re_z * re_y - im_z * im_y;
1899        im_z = re_z * im_y + im_z * re_y;
1900        re_z = z;
1901      }
1902      z = re_y * re_y - im_y * im_y;
1903      im_y = im_y * re_y + re_y * im_y;
1904      re_y = z;
1905      CHECK_DOUBLE_COMPLEX (p, re_y, im_y);
1906      CHECK_DOUBLE_COMPLEX (p, re_z, im_z);
1907      expo <<= 1;
1908    }
1909    if (negative) {
1910      PUSH_VALUE (p, dble (1.0q), A68_LONG_REAL);
1911      PUSH_VALUE (p, dble (0.0q), A68_LONG_REAL);
1912      PUSH_VALUE (p, dble (re_z), A68_LONG_REAL);
1913      PUSH_VALUE (p, dble (im_z), A68_LONG_REAL);
1914      genie_div_double_compl (p);
1915    } else {
1916      PUSH_VALUE (p, dble (re_z), A68_LONG_REAL);
1917      PUSH_VALUE (p, dble (im_z), A68_LONG_REAL);
1918    }
1919  }
1920  
1921  //! @brief OP = = (COMPLEX, COMPLEX) BOOL
1922  
1923  void genie_eq_double_compl (NODE_T * p)
1924  {
1925    A68_LONG_REAL re_x, im_x, re_y, im_y;
1926    POP_LONG_COMPLEX (p, &re_y, &im_y);
1927    POP_LONG_COMPLEX (p, &re_x, &im_x);
1928    PUSH_VALUE (p, (BOOL_T) ((VALUE (&re_x).f == VALUE (&re_y).f) && (VALUE (&im_x).f == VALUE (&im_y).f)), A68_BOOL);
1929  }
1930  
1931  //! @brief OP /= = (COMPLEX, COMPLEX) BOOL
1932  
1933  void genie_ne_double_compl (NODE_T * p)
1934  {
1935    A68_LONG_REAL re_x, im_x, re_y, im_y;
1936    POP_LONG_COMPLEX (p, &re_y, &im_y);
1937    POP_LONG_COMPLEX (p, &re_x, &im_x);
1938    PUSH_VALUE (p, (BOOL_T) ! ((VALUE (&re_x).f == VALUE (&re_y).f) && (VALUE (&im_x).f == VALUE (&im_y).f)), A68_BOOL);
1939  }
1940  
1941  //! @brief OP +:= = (REF COMPLEX, COMPLEX) REF COMPLEX
1942  
1943  void genie_plusab_double_compl (NODE_T * p)
1944  {
1945    genie_f_and_becomes (p, M_REF_LONG_COMPLEX, genie_add_double_compl);
1946  }
1947  
1948  //! @brief OP -:= = (REF COMPLEX, COMPLEX) REF COMPLEX
1949  
1950  void genie_minusab_double_compl (NODE_T * p)
1951  {
1952    genie_f_and_becomes (p, M_REF_LONG_COMPLEX, genie_sub_double_compl);
1953  }
1954  
1955  //! @brief OP *:= = (REF COMPLEX, COMPLEX) REF COMPLEX
1956  
1957  void genie_timesab_double_compl (NODE_T * p)
1958  {
1959    genie_f_and_becomes (p, M_REF_LONG_COMPLEX, genie_mul_double_compl);
1960  }
1961  
1962  //! @brief OP /:= = (REF COMPLEX, COMPLEX) REF COMPLEX
1963  
1964  void genie_divab_double_compl (NODE_T * p)
1965  {
1966    genie_f_and_becomes (p, M_REF_LONG_COMPLEX, genie_div_double_compl);
1967  }
1968  
1969  //! @brief OP LENG = (COMPLEX) LONG COMPLEX 
1970  
1971  void genie_lengthen_complex_to_double_compl (NODE_T * p)
1972  {
1973    A68_REAL i;
1974    POP_OBJECT (p, &i, A68_REAL);
1975    genie_lengthen_real_to_double_real (p);
1976    PUSH_OBJECT (p, i, A68_REAL);
1977    genie_lengthen_real_to_double_real (p);
1978  }
1979  
1980  // Functions
1981  
1982  #define CD_FUNCTION(name, fun)\
1983  void name (NODE_T * p) {\
1984    A68_LONG_REAL *x;\
1985    POP_OPERAND_ADDRESS (p, x, A68_LONG_REAL);\
1986    errno = 0;\
1987    VALUE (x).f = fun (VALUE (x).f);\
1988    MATH_RTE (p, errno != 0, M_LONG_REAL, NO_TEXT);\
1989  }
1990  
1991  CD_FUNCTION (genie_acos_double_real, acosq);
1992  CD_FUNCTION (genie_acosh_double_real, acoshq);
1993  CD_FUNCTION (genie_asinh_double_real, asinhq);
1994  CD_FUNCTION (genie_atanh_double_real, atanhq);
1995  CD_FUNCTION (genie_asin_double_real, asinq);
1996  CD_FUNCTION (genie_atan_double_real, atanq);
1997  CD_FUNCTION (genie_cosh_double_real, coshq);
1998  CD_FUNCTION (genie_cos_double_real, cosq);
1999  CD_FUNCTION (genie_curt_double_real, cbrtq);
2000  CD_FUNCTION (genie_exp_double_real, expq);
2001  CD_FUNCTION (genie_ln_double_real, logq);
2002  CD_FUNCTION (genie_log_double_real, log10q);
2003  CD_FUNCTION (genie_sinh_double_real, sinhq);
2004  CD_FUNCTION (genie_sin_double_real, sinq);
2005  CD_FUNCTION (genie_sqrt_double_real, sqrtq);
2006  CD_FUNCTION (genie_tanh_double_real, tanhq);
2007  CD_FUNCTION (genie_tan_double_real, tanq);
2008  CD_FUNCTION (genie_erf_double_real, erfq);
2009  CD_FUNCTION (genie_erfc_double_real, erfcq);
2010  CD_FUNCTION (genie_lngamma_double_real, lgammaq);
2011  CD_FUNCTION (genie_gamma_double_real, tgammaq);
2012  CD_FUNCTION (genie_csc_double_real, csc_double_real);
2013  CD_FUNCTION (genie_acsc_double_real, acsc_double_real);
2014  CD_FUNCTION (genie_sec_double_real, sec_double_real);
2015  CD_FUNCTION (genie_asec_double_real, asec_double_real);
2016  CD_FUNCTION (genie_cot_double_real, cot_double_real);
2017  CD_FUNCTION (genie_acot_double_real, acot_double_real);
2018  CD_FUNCTION (genie_sindg_double_real, sindg_double_real);
2019  CD_FUNCTION (genie_cosdg_double_real, cosdg_double_real);
2020  CD_FUNCTION (genie_tandg_double_real, tandg_double_real);
2021  CD_FUNCTION (genie_asindg_double_real, asindg_double_real);
2022  CD_FUNCTION (genie_acosdg_double_real, acosdg_double_real);
2023  CD_FUNCTION (genie_atandg_double_real, atandg_double_real);
2024  CD_FUNCTION (genie_cotdg_double_real, cotdg_double_real);
2025  CD_FUNCTION (genie_acotdg_double_real, acotdg_double_real);
2026  CD_FUNCTION (genie_sinpi_double_real, sinpi_double_real);
2027  CD_FUNCTION (genie_cospi_double_real, cospi_double_real);
2028  CD_FUNCTION (genie_tanpi_double_real, tanpi_double_real);
2029  CD_FUNCTION (genie_cotpi_double_real, cotpi_double_real);
2030  
2031  //! @brief PROC long arctan2 = (LONG REAL) LONG REAL
2032  
2033  void genie_atan2_double_real (NODE_T * p)
2034  {
2035    A68_LONG_REAL x, y;
2036    POP_OBJECT (p, &y, A68_LONG_REAL);
2037    POP_OBJECT (p, &x, A68_LONG_REAL);
2038    errno = 0;
2039    PRELUDE_ERROR (VALUE (&x).f == 0.0q && VALUE (&y).f == 0.0q, p, ERROR_INVALID_ARGUMENT, M_LONG_REAL);
2040    VALUE (&x).f = a68_atan2 (VALUE (&y).f, VALUE (&x).f);
2041    PRELUDE_ERROR (errno != 0, p, ERROR_MATH_EXCEPTION, NO_TEXT);
2042    PUSH_OBJECT (p, x, A68_LONG_REAL);
2043  }
2044  
2045  //! @brief PROC long arctan2dg = (LONG REAL) LONG REAL
2046  
2047  void genie_atan2dg_double_real (NODE_T * p)
2048  {
2049    A68_LONG_REAL x, y;
2050    POP_OBJECT (p, &y, A68_LONG_REAL);
2051    POP_OBJECT (p, &x, A68_LONG_REAL);
2052    errno = 0;
2053    PRELUDE_ERROR (VALUE (&x).f == 0.0q && VALUE (&y).f == 0.0q, p, ERROR_INVALID_ARGUMENT, M_LONG_REAL);
2054    VALUE (&x).f = CONST_180_OVER_PI_Q * a68_atan2 (VALUE (&y).f, VALUE (&x).f);
2055    PRELUDE_ERROR (errno != 0, p, ERROR_MATH_EXCEPTION, NO_TEXT);
2056    PUSH_OBJECT (p, x, A68_LONG_REAL);
2057  }
2058  
2059  //! @brief PROC (LONG REAL) LONG REAL inverf
2060  
2061  void genie_inverf_double_real (NODE_T * _p_)
2062  {
2063    A68_LONG_REAL x;
2064    DOUBLE_T y, z;
2065    A68 (f_entry) = _p_;
2066    POP_OBJECT (_p_, &x, A68_LONG_REAL);
2067    errno = 0;
2068    y = VALUE (&x).f;
2069    z = inverf_double_real (y);
2070    MATH_RTE (_p_, errno != 0, M_LONG_REAL, NO_TEXT);
2071    CHECK_DOUBLE_REAL (_p_, z);
2072    PUSH_VALUE (_p_, dble (z), A68_LONG_REAL);
2073  }
2074  
2075  //! @brief PROC (LONG REAL) LONG REAL inverfc
2076  
2077  void genie_inverfc_double_real (NODE_T * p)
2078  {
2079    A68_LONG_REAL *u;
2080    POP_OPERAND_ADDRESS (p, u, A68_LONG_REAL);
2081    VALUE (u).f = 1.0q - (VALUE (u).f);
2082    genie_inverf_double_real (p);
2083  }
2084  
2085  #define _re_ (VALUE (&re).f)
2086  #define _im_ (VALUE (&im).f)
2087  
2088  #define CD_C_FUNCTION(p, g)\
2089    A68_LONG_REAL re, im;\
2090    DOUBLE_COMPLEX_T z;\
2091    POP_OBJECT (p, &im, A68_LONG_REAL);\
2092    POP_OBJECT (p, &re, A68_LONG_REAL);\
2093    errno = 0;\
2094    z = VALUE (&re).f + VALUE (&im).f * _Complex_I;\
2095    z = g (z);\
2096    PUSH_VALUE (p, dble ((DOUBLE_T) crealq (z)), A68_LONG_REAL);\
2097    PUSH_VALUE (p, dble ((DOUBLE_T) cimagq (z)), A68_LONG_REAL);\
2098    MATH_RTE (p, errno != 0, M_COMPLEX, NO_TEXT);
2099  
2100  //! @brief PROC long csqrt = (LONG COMPLEX) LONG COMPLEX
2101  
2102  void genie_sqrt_double_compl (NODE_T * p)
2103  {
2104    CD_C_FUNCTION (p, csqrtq);
2105  }
2106  
2107  //! @brief PROC long csin = (LONG COMPLEX) LONG COMPLEX
2108  
2109  void genie_sin_double_compl (NODE_T * p)
2110  {
2111    CD_C_FUNCTION (p, csinq);
2112  }
2113  
2114  //! @brief PROC long ccos = (LONG COMPLEX) LONG COMPLEX
2115  
2116  void genie_cos_double_compl (NODE_T * p)
2117  {
2118    CD_C_FUNCTION (p, ccosq);
2119  }
2120  
2121  //! @brief PROC long ctan = (LONG COMPLEX) LONG COMPLEX
2122  
2123  void genie_tan_double_compl (NODE_T * p)
2124  {
2125    CD_C_FUNCTION (p, ctanq);
2126  }
2127  
2128  //! @brief PROC long casin = (LONG COMPLEX) LONG COMPLEX
2129  
2130  void genie_asin_double_compl (NODE_T * p)
2131  {
2132    CD_C_FUNCTION (p, casinq);
2133  }
2134  
2135  //! @brief PROC long cacos = (LONG COMPLEX) LONG COMPLEX
2136  
2137  void genie_acos_double_compl (NODE_T * p)
2138  {
2139    CD_C_FUNCTION (p, cacosq);
2140  }
2141  
2142  //! @brief PROC long catan = (LONG COMPLEX) LONG COMPLEX
2143  
2144  void genie_atan_double_compl (NODE_T * p)
2145  {
2146    CD_C_FUNCTION (p, catanq);
2147  }
2148  
2149  //! @brief PROC long cexp = (LONG COMPLEX) LONG COMPLEX
2150  
2151  void genie_exp_double_compl (NODE_T * p)
2152  {
2153    CD_C_FUNCTION (p, cexpq);
2154  }
2155  
2156  //! @brief PROC long cln = (LONG COMPLEX) LONG COMPLEX
2157  
2158  void genie_ln_double_compl (NODE_T * p)
2159  {
2160    CD_C_FUNCTION (p, clogq);
2161  }
2162  
2163  //! @brief PROC long csinh = (LONG COMPLEX) LONG COMPLEX
2164  
2165  void genie_sinh_double_compl (NODE_T * p)
2166  {
2167    CD_C_FUNCTION (p, csinhq);
2168  }
2169  
2170  //! @brief PROC long ccosh = (LONG COMPLEX) LONG COMPLEX
2171  
2172  void genie_cosh_double_compl (NODE_T * p)
2173  {
2174    CD_C_FUNCTION (p, ccoshq);
2175  }
2176  
2177  //! @brief PROC long ctanh = (LONG COMPLEX) LONG COMPLEX
2178  
2179  void genie_tanh_double_compl (NODE_T * p)
2180  {
2181    CD_C_FUNCTION (p, ctanhq);
2182  }
2183  
2184  //! @brief PROC long casinh = (LONG COMPLEX) LONG COMPLEX
2185  
2186  void genie_asinh_double_compl (NODE_T * p)
2187  {
2188    CD_C_FUNCTION (p, casinhq);
2189  }
2190  
2191  //! @brief PROC long cacosh = (LONG COMPLEX) LONG COMPLEX
2192  
2193  void genie_acosh_double_compl (NODE_T * p)
2194  {
2195    CD_C_FUNCTION (p, cacoshq);
2196  }
2197  
2198  //! @brief PROC long catanh = (LONG COMPLEX) LONG COMPLEX
2199  
2200  void genie_atanh_double_compl (NODE_T * p)
2201  {
2202    CD_C_FUNCTION (p, catanhq);
2203  }
2204  
2205  #undef _re_
2206  #undef _im_
2207  
2208  //! @brief PROC next long random = LONG REAL
2209  
2210  void genie_next_random_double_real (NODE_T * p)
2211  {
2212  // This is 'real width' digits only.
2213    genie_next_random (p);
2214    genie_lengthen_real_to_double_real (p);
2215  }
2216  
2217  #define CALL(g, x, y) {\
2218    ADDR_T pop_sp = A68_SP;\
2219    A68_LONG_REAL *z = (A68_LONG_REAL *) STACK_TOP;\
2220    DOUBLE_NUM_T _w_;\
2221    _w_.f = (x);\
2222    PUSH_VALUE (_p_, _w_, A68_LONG_REAL);\
2223    genie_call_procedure (_p_, M_PROC_LONG_REAL_LONG_REAL, M_PROC_LONG_REAL_LONG_REAL, M_PROC_LONG_REAL_LONG_REAL, &(g), pop_sp, pop_fp);\
2224    (y) = VALUE (z).f;\
2225    A68_SP = pop_sp;\
2226  }
2227  
2228  //! @brief Transform string into real-16.
2229  
2230  DOUBLE_T string_to_double_real (char *s, char **end)
2231  {
2232    int i, dot = -1, pos = 0, pow = 0, expo;
2233    DOUBLE_T sum, W, y[FLT128_DIG];
2234    errno = 0;
2235    for (i = 0; i < FLT128_DIG; i++) {
2236      y[i] = 0.0q;
2237    }
2238    while (IS_SPACE (s[0])) {
2239      s++;
2240    }
2241  // Scan mantissa digits and put them into "y".
2242    if (s[0] == '-') {
2243      W = -1.0q;
2244    } else {
2245      W = 1.0q;
2246    }
2247    if (s[0] == '+' || s[0] == '-') {
2248      s++;
2249    }
2250    while (s[0] == '0') {
2251      s++;
2252    }
2253    while (pow < FLT128_DIG && s[pos] != NULL_CHAR && (IS_DIGIT (s[pos]) || s[pos] == POINT_CHAR)) {
2254      if (s[pos] == POINT_CHAR) {
2255        dot = pos;
2256      } else {
2257        int val = (int) s[pos] - (int) '0';
2258        y[pow] = W * val;
2259        W /= 10.0q;
2260        pow++;
2261      }
2262      pos++;
2263    }
2264    (*end) = &(s[pos]);
2265  // Sum from low to high to preserve precision.
2266    sum = 0.0q;
2267    for (i = FLT128_DIG - 1; i >= 0; i--) {
2268      sum = sum + y[i];
2269    }
2270  // See if there is an exponent.
2271    if (s[pos] != NULL_CHAR && TO_UPPER (s[pos]) == TO_UPPER (EXPONENT_CHAR)) {
2272      expo = (int) strtol (&(s[++pos]), end, 10);
2273    } else {
2274      expo = 0;
2275    }
2276  // Standardise.
2277    if (dot >= 0) {
2278      expo += dot - 1;
2279    } else {
2280      expo += pow - 1;
2281    }
2282    while (sum != 0.0q && fabsq (sum) < 1.0q) {
2283      sum *= 10.0q;
2284      expo -= 1;
2285    }
2286  //
2287    if (errno == 0) {
2288      return sum * ten_up_double (expo);
2289    } else {
2290      return 0.0q;
2291    }
2292  }
2293  
2294  void genie_beta_inc_cf_double_real (NODE_T * p)
2295  {
2296    A68_LONG_REAL x, s, t;
2297    POP_OBJECT (p, &x, A68_LONG_REAL);
2298    POP_OBJECT (p, &t, A68_LONG_REAL);
2299    POP_OBJECT (p, &s, A68_LONG_REAL);
2300    errno = 0;
2301    PUSH_VALUE (p, dble (a68_beta_inc_double_real (VALUE (&s).f, VALUE (&t).f, VALUE (&x).f)), A68_LONG_REAL);
2302    MATH_RTE (p, errno != 0, M_LONG_REAL, NO_TEXT);
2303  }
2304  
2305  void genie_beta_double_real (NODE_T * p)
2306  {
2307    A68_LONG_REAL a, b;
2308    POP_OBJECT (p, &b, A68_LONG_REAL);
2309    POP_OBJECT (p, &a, A68_LONG_REAL);
2310    errno = 0;
2311    PUSH_VALUE (p, dble (expq (lgammaq (VALUE (&a).f) + lgammaq (VALUE (&b).f) - lgammaq (VALUE (&a).f + VALUE (&b).f))), A68_LONG_REAL);
2312    MATH_RTE (p, errno != 0, M_LONG_REAL, NO_TEXT);
2313  }
2314  
2315  void genie_ln_beta_double_real (NODE_T * p)
2316  {
2317    A68_LONG_REAL a, b;
2318    POP_OBJECT (p, &b, A68_LONG_REAL);
2319    POP_OBJECT (p, &a, A68_LONG_REAL);
2320    errno = 0;
2321    PUSH_VALUE (p, dble (lgammaq (VALUE (&a).f) + lgammaq (VALUE (&b).f) - lgammaq (VALUE (&a).f + VALUE (&b).f)), A68_LONG_REAL);
2322    MATH_RTE (p, errno != 0, M_LONG_REAL, NO_TEXT);
2323  }
2324  
2325  // LONG REAL infinity
2326  
2327  void genie_infinity_double_real (NODE_T * p)
2328  {
2329    PUSH_VALUE (p, dble (a68_posinf ()), A68_LONG_REAL);
2330  }
2331  
2332  // LONG REAL minus infinity
2333  
2334  void genie_minus_infinity_double_real (NODE_T * p)
2335  {
2336    PUSH_VALUE (p, dble (a68_dneginf ()), A68_LONG_REAL);
2337  }
2338  
2339  #endif