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  
  38  // 128-bit REAL support.
  39  
  40  // Conversions.
  41  
  42  DOUBLE_NUM_T double_int_to_double (NODE_T * p, DOUBLE_NUM_T z)
  43  {
  44    int neg = D_NEG (z);
  45    if (neg) {
  46      z = abs_double_int (z);
  47    }
  48    DOUBLE_NUM_T w, radix;
  49    w.f = 0.0q;
  50    set_lw (radix, RADIX);
  51    DOUBLE_T weight = 1.0q;
  52    while (!D_ZERO (z)) {
  53      DOUBLE_NUM_T digit;
  54      digit = double_udiv (p, M_LONG_INT, z, radix, 1);
  55      w.f = w.f + LW (digit) * weight;
  56      z = double_udiv (p, M_LONG_INT, z, radix, 0);
  57      weight = weight * RADIX_Q;
  58    }
  59    if (neg) {
  60      w.f = -w.f;
  61    }
  62    return w;
  63  }
  64  
  65  DOUBLE_NUM_T double_to_double_int (NODE_T * p, DOUBLE_NUM_T z)
  66  {
  67  // This routines looks a lot like "strtol". 
  68    BOOL_T negative = (BOOL_T) (z.f < 0);
  69    z.f = fabs_double (trunc_double (z.f));
  70    if (z.f > CONST_2_UP_112_Q) {
  71      errno = EDOM;
  72      MATH_RTE (p, errno != 0, M_LONG_REAL, NO_TEXT);
  73    }
  74    DOUBLE_NUM_T sum, weight, radix;
  75    set_lw (sum, 0);
  76    set_lw (weight, 1);
  77    set_lw (radix, RADIX);
  78    while (z.f > 0) {
  79      DOUBLE_NUM_T term, digit, quot, rest;
  80      quot.f = trunc_double (z.f / RADIX_Q);
  81      rest.f = z.f - quot.f * RADIX_Q;
  82      z.f = quot.f;
  83      set_lw (digit, (INT_T) (rest.f));
  84      term = double_umul (p, M_LONG_INT, digit, weight);
  85      sum = double_uadd (p, M_LONG_INT, sum, term);
  86      if (z.f > 0.0q) {
  87        weight = double_umul (p, M_LONG_INT, weight, radix);
  88      }
  89    }
  90    if (negative) {
  91      return neg_double_int (sum);
  92    } else {
  93      return sum;
  94    }
  95  }
  96  
  97  //! @brief Value of LONG INT denotation
  98  
  99  int string_to_double_int (NODE_T * p, A68_LONG_INT * z, char *s)
 100  {
 101    while (IS_SPACE (s[0])) {
 102      s++;
 103    }
 104  // Get the sign
 105    int sign = (s[0] == '-' ? -1 : 1);
 106    if (s[0] == '+' || s[0] == '-') {
 107      s++;
 108    }
 109    int end = 0;
 110    while (s[end] != '\0') {
 111      end++;
 112    }
 113    DOUBLE_NUM_T weight, ten, sum;
 114    set_lw (sum, 0);
 115    set_lw (weight, 1);
 116    set_lw (ten, 10);
 117    for (int k = end - 1; k >= 0; k--) {
 118      DOUBLE_NUM_T term;
 119      int digit = s[k] - '0';
 120      set_lw (term, digit);
 121      term = double_umul (p, M_LONG_INT, term, weight);
 122      sum = double_uadd (p, M_LONG_INT, sum, term);
 123      weight = double_umul (p, M_LONG_INT, weight, ten);
 124    }
 125    if (sign == -1) {
 126      HW (sum) = HW (sum) | D_SIGN;
 127    }
 128    VALUE (z) = sum;
 129    STATUS (z) = INIT_MASK;
 130    return A68_TRUE;
 131  }
 132  
 133  //! @brief LONG BITS value of LONG BITS denotation
 134  
 135  DOUBLE_NUM_T double_strtou (NODE_T * p, char *s)
 136  {
 137    errno = 0;
 138    char *radix = NO_TEXT;
 139    int base = (int) a68_strtou (s, &radix, 10);
 140    if (base < 2 || base > 16) {
 141      diagnostic (A68_RUNTIME_ERROR, p, ERROR_INVALID_RADIX, base);
 142      exit_genie (p, A68_RUNTIME_ERROR);
 143    }
 144    DOUBLE_NUM_T z;
 145    set_lw (z, 0x0);
 146    if (radix != NO_TEXT && TO_UPPER (radix[0]) == TO_UPPER (RADIX_CHAR) && errno == 0) {
 147      DOUBLE_NUM_T w;
 148      char *q = radix;
 149      while (q[0] != NULL_CHAR) {
 150        q++;
 151      }
 152      set_lw (w, 1);
 153      while ((--q) != radix) {
 154        int digit = char_value (q[0]);
 155        if (digit < 0 && digit >= base) {
 156          diagnostic (A68_RUNTIME_ERROR, p, ERROR_IN_DENOTATION, M_LONG_BITS);
 157          exit_genie (p, A68_RUNTIME_ERROR);
 158        } else {
 159          DOUBLE_NUM_T v;
 160          set_lw (v, digit);
 161          v = double_umul (p, M_LONG_INT, v, w);
 162          z = double_uadd (p, M_LONG_INT, z, v);
 163          set_lw (v, base);
 164          w = double_umul (p, M_LONG_INT, w, v);
 165        }
 166      }
 167    } else {
 168      diagnostic (A68_RUNTIME_ERROR, p, ERROR_IN_DENOTATION, M_LONG_BITS);
 169      exit_genie (p, A68_RUNTIME_ERROR);
 170    }
 171    return (z);
 172  }
 173  
 174  //! @brief OP LENG = (BITS) LONG BITS
 175  
 176  void genie_lengthen_bits_to_double_bits (NODE_T * p)
 177  {
 178    A68_BITS k;
 179    POP_OBJECT (p, &k, A68_BITS);
 180    DOUBLE_NUM_T d;
 181    LW (d) = VALUE (&k);
 182    HW (d) = 0;
 183    PUSH_VALUE (p, d, A68_LONG_BITS);
 184  }
 185  
 186  //! @brief OP SHORTEN = (LONG BITS) BITS
 187  
 188  void genie_shorten_double_bits_to_bits (NODE_T * p)
 189  {
 190    A68_LONG_BITS k;
 191    POP_OBJECT (p, &k, A68_LONG_BITS);
 192    DOUBLE_NUM_T j = VALUE (&k);
 193    PRELUDE_ERROR (HW (j) != 0, p, ERROR_MATH, M_BITS);
 194    PUSH_VALUE (p, LW (j), A68_BITS);
 195  }
 196  
 197  //! @brief Convert to other radix, binary up to hexadecimal.
 198  
 199  BOOL_T convert_radix_double (NODE_T * p, DOUBLE_NUM_T z, int radix, int width)
 200  {
 201    if (radix < 2 || radix > 16) {
 202      radix = 16;
 203    }
 204    DOUBLE_NUM_T w, rad;
 205    set_lw (rad, radix);
 206    reset_transput_buffer (EDIT_BUFFER);
 207    if (width > 0) {
 208      while (width > 0) {
 209        w = double_udiv (p, M_LONG_INT, z, rad, 1);
 210        plusto_transput_buffer (p, digchar (LW (w)), EDIT_BUFFER);
 211        width--;
 212        z = double_udiv (p, M_LONG_INT, z, rad, 0);
 213      }
 214      return D_ZERO (z);
 215    } else if (width == 0) {
 216      do {
 217        w = double_udiv (p, M_LONG_INT, z, rad, 1);
 218        plusto_transput_buffer (p, digchar (LW (w)), EDIT_BUFFER);
 219        z = double_udiv (p, M_LONG_INT, z, rad, 0);
 220      }
 221      while (!D_ZERO (z));
 222      return A68_TRUE;
 223    } else {
 224      return A68_FALSE;
 225    }
 226  }
 227  
 228  //! @brief OP LENG = (LONG INT) LONG REAL
 229  
 230  void genie_widen_double_int_to_double (NODE_T * p)
 231  {
 232    A68_DOUBLE *z = (A68_DOUBLE *) STACK_TOP;
 233    GENIE_UNIT (SUB (p));
 234    VALUE (z) = double_int_to_double (p, VALUE (z));
 235  }
 236  
 237  //! @brief OP LENG = (REAL) LONG REAL
 238  
 239  DOUBLE_NUM_T dble_double (NODE_T * p, REAL_T z)
 240  {
 241  // Quick and dirty, only works with 64-bit INT_T.
 242    BOOL_T nega = (z < 0.0);
 243    REAL_T u = fabs (z);
 244    int expo = 0;
 245    standardise (&u, 1, A68_REAL_DIG, &expo);
 246    u *= ten_up (A68_REAL_DIG);
 247    expo -= A68_REAL_DIG;
 248    DOUBLE_NUM_T w;
 249    set_lw (w, (INT_T) u);
 250    w = double_int_to_double (p, w);
 251    w.f *= ten_up_double (expo);
 252    if (nega) {
 253      w.f = -w.f;
 254    }
 255    return w;
 256  }
 257  
 258  //! @brief OP LENG = (REAL) LONG REAL
 259  
 260  void genie_lengthen_real_to_double (NODE_T * p)
 261  {
 262    A68_REAL z;
 263    POP_OBJECT (p, &z, A68_REAL);
 264    PUSH_VALUE (p, dble_double (p, VALUE (&z)), A68_LONG_REAL);
 265  }
 266  
 267  //! @brief OP SHORTEN = (LONG REAL) REAL
 268  
 269  void genie_shorten_double_to_real (NODE_T * p)
 270  {
 271    A68_LONG_REAL z;
 272    POP_OBJECT (p, &z, A68_LONG_REAL);
 273    REAL_T w = VALUE (&z).f;
 274    PUSH_VALUE (p, w, A68_REAL);
 275  }
 276  
 277  //! @brief Convert integer to multi-precison number.
 278  
 279  MP_T *double_int_to_mp (NODE_T * p, MP_T * z, DOUBLE_NUM_T k, int digs)
 280  {
 281    int negative = D_NEG (k);
 282    if (negative) {
 283      k = neg_double_int (k);
 284    }
 285    DOUBLE_NUM_T radix;
 286    set_lw (radix, MP_RADIX);
 287    DOUBLE_NUM_T k2 = k;
 288    int n = 0;
 289    do {
 290      k2 = double_udiv (p, M_LONG_INT, k2, radix, 0);
 291      if (!D_ZERO (k2)) {
 292        n++;
 293      }
 294    }
 295    while (!D_ZERO (k2));
 296    SET_MP_ZERO (z, digs);
 297    MP_EXPONENT (z) = (MP_T) n;
 298    for (int j = 1 + n; j >= 1; j--) {
 299      DOUBLE_NUM_T term = double_udiv (p, M_LONG_INT, k, radix, 1);
 300      MP_DIGIT (z, j) = (MP_T) LW (term);
 301      k = double_udiv (p, M_LONG_INT, k, radix, 0);
 302    }
 303    MP_DIGIT (z, 1) = (negative ? -MP_DIGIT (z, 1) : MP_DIGIT (z, 1));
 304    check_mp_exp (p, z);
 305    return z;
 306  }
 307  
 308  //! @brief Convert multi-precision number to integer.
 309  
 310  DOUBLE_NUM_T mp_to_double_int (NODE_T * p, MP_T * z, int digs)
 311  {
 312  // This routines looks a lot like "strtol". 
 313    int expo = (int) MP_EXPONENT (z);
 314    DOUBLE_NUM_T sum, weight;
 315    set_lw (sum, 0);
 316    set_lw (weight, 1);
 317    BOOL_T negative;
 318    if (expo >= digs) {
 319      diagnostic (A68_RUNTIME_ERROR, p, ERROR_OUT_OF_BOUNDS, MOID (p));
 320      exit_genie (p, A68_RUNTIME_ERROR);
 321    }
 322    negative = (BOOL_T) (MP_DIGIT (z, 1) < 0);
 323    if (negative) {
 324      MP_DIGIT (z, 1) = -MP_DIGIT (z, 1);
 325    }
 326    for (int j = 1 + expo; j >= 1; j--) {
 327      DOUBLE_NUM_T term, digit, radix;
 328      set_lw (digit, (MP_INT_T) MP_DIGIT (z, j));
 329      term = double_umul (p, M_LONG_INT, digit, weight);
 330      sum = double_uadd (p, M_LONG_INT, sum, term);
 331      set_lw (radix, MP_RADIX);
 332      weight = double_umul (p, M_LONG_INT, weight, radix);
 333    }
 334    if (negative) {
 335      return neg_double_int (sum);
 336    } else {
 337      return sum;
 338    }
 339  }
 340  
 341  //! @brief Convert real to multi-precison number.
 342  
 343  MP_T *double_to_mp (NODE_T * p, MP_T * z, DOUBLE_T x, int digs)
 344  {
 345    SET_MP_ZERO (z, digs);
 346    if (x == 0.0q) {
 347      return z;
 348    }
 349  // Small integers can be done better by int_to_mp.
 350    if (ABS (x) < MP_RADIX && trunc_double (x) == x) {
 351      return int_to_mp (p, z, (int) trunc_double (x), digs);
 352    }
 353    int sign_x = SIGN (x);
 354  // Scale to [0, 0.1>.
 355    DOUBLE_T a = ABS (x);
 356    INT_T expo = (int) log10_double (a);
 357    a /= ten_up_double (expo);
 358    expo--;
 359    if (a >= 1.0q) {
 360      a /= 10.0q;
 361      expo++;
 362    }
 363  // Transport digits of x to the mantissa of z.
 364    int j = 1, sum = 0, weight = (MP_RADIX / 10);
 365    for (int k = 0; a != 0.0q && j <= digs && k < A68_DOUBLE_MAN; k++) {
 366      DOUBLE_T u = a * 10.0q;
 367      DOUBLE_T v = floor_double (u);
 368      a = u - v;
 369      sum += weight * (int) v;
 370      weight /= 10;
 371      if (weight < 1) {
 372        MP_DIGIT (z, j++) = (MP_T) sum;
 373        sum = 0;
 374        weight = (MP_RADIX / 10);
 375      }
 376    }
 377  // Store the last digits.
 378    if (j <= digs) {
 379      MP_DIGIT (z, j) = (MP_T) sum;
 380    }
 381    (void) align_mp (z, &expo, digs);
 382    MP_EXPONENT (z) = (MP_T) expo;
 383    MP_DIGIT (z, 1) *= sign_x;
 384    check_mp_exp (p, z);
 385    return z;
 386  }
 387  
 388  //! @brief Convert multi-precision number to real.
 389  
 390  DOUBLE_T mp_to_double (NODE_T * p, MP_T * z, int digs)
 391  {
 392  // This routine looks a lot like "strtod".
 393    if (MP_EXPONENT (z) * (MP_T) LOG_MP_RADIX <= (MP_T) A68_DOUBLE_MIN_EXP) {
 394      return 0;
 395    } else {
 396      DOUBLE_T weight = ten_up_double ((int) (MP_EXPONENT (z) * LOG_MP_RADIX));
 397      int lim = MIN (digs, MP_MAX_DIGITS);
 398      DOUBLE_T terms[1 + MP_MAX_DIGITS];
 399      for (int k = 1; k <= lim; k++) {
 400        terms[k] = ABS (MP_DIGIT (z, k)) * weight;
 401        weight /= MP_RADIX;
 402      }
 403  // Sum terms from small to large.
 404      DOUBLE_T sum = 0;
 405      for (int k = lim; k >= 1; k--) {
 406        sum += terms[k];
 407      }
 408      CHECK_DOUBLE_REAL (p, sum);
 409      return MP_DIGIT (z, 1) >= 0 ? sum : -sum;
 410    }
 411  }
 412  
 413  DOUBLE_T inverf_double (DOUBLE_T z)
 414  {
 415    if (fabs_double (z) >= 1.0q) {
 416      errno = EDOM;
 417      return z;
 418    } else {
 419  // Newton-Raphson.
 420      DOUBLE_T f = sqrt_double (M_PIq) / 2, g, x = z;
 421      int its = 10;
 422      x = dble (a68_inverf_real ((REAL_T) x)).f;
 423      do {
 424        g = x;
 425        x -= f * (erf_double (x) - z) / exp_double (-(x * x));
 426      } while (its-- > 0 && errno == 0 && fabs_double (x - g) > (3 * A68_DOUBLE_EPS));
 427      return x;
 428    }
 429  }
 430  
 431  //! @brief OP LENG = (LONG REAL) LONG LONG REAL
 432  
 433  void genie_lengthen_double_to_mp (NODE_T * p)
 434  {
 435    int digs = DIGITS (M_LONG_LONG_REAL);
 436    A68_LONG_REAL x;
 437    POP_OBJECT (p, &x, A68_LONG_REAL);
 438    MP_T *z = nil_mp (p, digs);
 439    (void) double_to_mp (p, z, VALUE (&x).f, digs);
 440    MP_STATUS (z) = (MP_T) INIT_MASK;
 441  }
 442  
 443  //! @brief OP SHORTEN = (LONG LONG REAL) LONG REAL
 444  
 445  void genie_shorten_mp_to_double (NODE_T * p)
 446  {
 447    MOID_T *mode = LHS_MODE (p);
 448    int digs = DIGITS (mode), size = SIZE (mode);
 449    DOUBLE_NUM_T d;
 450    DECREMENT_STACK_POINTER (p, size);
 451    MP_T *z = (MP_T *) STACK_TOP;
 452    MP_STATUS (z) = (MP_T) INIT_MASK;
 453    d.f = mp_to_double (p, z, digs);
 454    PUSH_VALUE (p, d, A68_LONG_REAL);
 455  }
 456  
 457  //! @brief OP SHORTEN = (LONG LONG COMPLEX) LONG COMPLEX
 458  
 459  void genie_shorten_long_mp_complex_to_double_compl (NODE_T * p)
 460  {
 461    int digs = DIGITS (M_LONG_LONG_REAL), size = SIZE (M_LONG_LONG_REAL);
 462    MP_T *b = (MP_T *) STACK_OFFSET (-size);
 463    MP_T *a = (MP_T *) STACK_OFFSET (-2 * size);
 464    DECREMENT_STACK_POINTER (p, 2 * size);
 465    DOUBLE_NUM_T u, v;
 466    u.f = mp_to_double (p, a, digs);
 467    v.f = mp_to_double (p, b, digs);
 468    PUSH_VALUE (p, u, A68_LONG_REAL);
 469    PUSH_VALUE (p, v, A68_LONG_REAL);
 470  }
 471  
 472  //! @brief OP LENG = (LONG INT) LONG LONG INT
 473  
 474  void genie_lengthen_double_int_to_mp (NODE_T * p)
 475  {
 476    int digs = DIGITS (M_LONG_LONG_INT);
 477    A68_LONG_INT k;
 478    POP_OBJECT (p, &k, A68_LONG_INT);
 479    MP_T *z = nil_mp (p, digs);
 480    (void) double_int_to_mp (p, z, VALUE (&k), digs);
 481    MP_STATUS (z) = (MP_T) INIT_MASK;
 482  }
 483  
 484  //! @brief OP SHORTEN = (LONG LONG INT) LONG INT
 485  
 486  void genie_shorten_mp_to_double_int (NODE_T * p)
 487  {
 488    MOID_T *mode = LHS_MODE (p);
 489    int digs = DIGITS (mode), size = SIZE (mode);
 490    DECREMENT_STACK_POINTER (p, size);
 491    MP_T *z = (MP_T *) STACK_TOP;
 492    MP_STATUS (z) = (MP_T) INIT_MASK;
 493    PUSH_VALUE (p, mp_to_double_int (p, z, digs), A68_LONG_INT);
 494  }
 495  
 496  //! @brief OP LENG = (INT) LONG INT
 497  
 498  void genie_lengthen_int_to_double_int (NODE_T * p)
 499  {
 500    A68_INT k;
 501    POP_OBJECT (p, &k, A68_INT);
 502    INT_T v = VALUE (&k);
 503    DOUBLE_NUM_T d;
 504    if (v >= 0) {
 505      LW (d) = v;
 506      HW (d) = 0;
 507    } else {
 508      LW (d) = -v;
 509      HW (d) = D_SIGN;
 510    }
 511    PUSH_VALUE (p, d, A68_LONG_INT);
 512  }
 513  
 514  //! @brief OP SHORTEN = (LONG INT) INT
 515  
 516  void genie_shorten_long_int_to_int (NODE_T * p)
 517  {
 518    A68_LONG_INT k;
 519    POP_OBJECT (p, &k, A68_LONG_INT);
 520    DOUBLE_NUM_T j = VALUE (&k);
 521    PRELUDE_ERROR (HW (j) != 0 && HW (j) != D_SIGN, p, ERROR_MATH, M_INT);
 522    PRELUDE_ERROR (LW (j) & D_SIGN, p, ERROR_MATH, M_INT);
 523    if (D_NEG (j)) {
 524      PUSH_VALUE (p, -LW (j), A68_INT);
 525    } else {
 526      PUSH_VALUE (p, LW (j), A68_INT);
 527    }
 528  }
 529  
 530  // Constants.
 531  
 532  //! @brief PROC long max int = LONG INT
 533  
 534  void genie_double_max_int (NODE_T * p)
 535  {
 536    DOUBLE_NUM_T d;
 537    HW (d) = 0x7fffffffffffffffLL;
 538    LW (d) = 0xffffffffffffffffLL;
 539    PUSH_VALUE (p, d, A68_LONG_INT);
 540  }
 541  
 542  //! @brief PROC long max bits = LONG BITS
 543  
 544  void genie_double_max_bits (NODE_T * p)
 545  {
 546    DOUBLE_NUM_T d;
 547    HW (d) = 0xffffffffffffffffLL;
 548    LW (d) = 0xffffffffffffffffLL;
 549    PUSH_VALUE (p, d, A68_LONG_INT);
 550  }
 551  
 552  //! @brief LONG REAL max long real
 553  
 554  void genie_double_max_real (NODE_T * p)
 555  {
 556    DOUBLE_NUM_T d;
 557    d.f = A68_DOUBLE_MAX;
 558    PUSH_VALUE (p, d, A68_LONG_REAL);
 559  }
 560  
 561  //! @brief LONG REAL min long real
 562  
 563  void genie_double_min_real (NODE_T * p)
 564  {
 565    DOUBLE_NUM_T d;
 566    d.f = A68_DOUBLE_MIN;
 567    PUSH_VALUE (p, d, A68_LONG_REAL);
 568  }
 569  
 570  //! @brief LONG REAL small long real
 571  
 572  void genie_double_small_real (NODE_T * p)
 573  {
 574    DOUBLE_NUM_T d;
 575    d.f = A68_DOUBLE_EPS;
 576    PUSH_VALUE (p, d, A68_LONG_REAL);
 577  }
 578  
 579  //! @brief PROC long pi = LON REAL
 580  
 581  void genie_pi_double (NODE_T * p)
 582  {
 583    DOUBLE_NUM_T w;
 584    w.f = M_PIq;
 585    PUSH_VALUE (p, w, A68_LONG_INT);
 586  }
 587  
 588  // MONADs and DYADs
 589  
 590  //! @brief OP SIGN = (LONG INT) INT
 591  
 592  void genie_sign_double_int (NODE_T * p)
 593  {
 594    A68_LONG_INT k;
 595    POP_OBJECT (p, &k, A68_LONG_INT);
 596    PUSH_VALUE (p, sign_double_int (VALUE (&k)), A68_INT);
 597  }
 598  
 599  //! @brief OP ABS = (LONG INT) LONG INT
 600  
 601  void genie_abs_double_int (NODE_T * p)
 602  {
 603    A68_LONG_INT *k;
 604    POP_OPERAND_ADDRESS (p, k, A68_LONG_INT);
 605    VALUE (k) = abs_double_int (VALUE (k));
 606  }
 607  
 608  //! @brief OP ODD = (LONG INT) BOOL
 609  
 610  void genie_odd_double_int (NODE_T * p)
 611  {
 612    A68_LONG_INT j;
 613    POP_OBJECT (p, &j, A68_LONG_INT);
 614    DOUBLE_NUM_T w = abs_double_int (VALUE (&j));
 615    if (LW (w) & 0x1) {
 616      PUSH_VALUE (p, A68_TRUE, A68_BOOL);
 617    } else {
 618      PUSH_VALUE (p, A68_FALSE, A68_BOOL);
 619    }
 620  }
 621  
 622  //! @brief OP - = (LONG INT) LONG INT
 623  
 624  void genie_minus_double_int (NODE_T * p)
 625  {
 626    A68_LONG_INT *k;
 627    POP_OPERAND_ADDRESS (p, k, A68_LONG_INT);
 628    VALUE (k) = neg_double_int (VALUE (k));
 629  }
 630  
 631  //! @brief OP ROUND = (LONG REAL) LONG INT
 632  
 633  void genie_round_double (NODE_T * p)
 634  {
 635    A68_LONG_REAL x;
 636    POP_OBJECT (p, &x, A68_LONG_REAL);
 637    DOUBLE_NUM_T u = VALUE (&x);
 638    if (u.f < 0.0q) {
 639      u.f = u.f - 0.5q;
 640    } else {
 641      u.f = u.f + 0.5q;
 642    }
 643    PUSH_VALUE (p, double_to_double_int (p, u), A68_LONG_INT);
 644  }
 645  
 646  //! @brief OP ENTIER = (LONG REAL) LONG INT
 647  
 648  void genie_entier_double (NODE_T * p)
 649  {
 650    A68_LONG_REAL x;
 651    POP_OBJECT (p, &x, A68_LONG_REAL);
 652    DOUBLE_NUM_T u = VALUE (&x);
 653    u.f = floor_double (u.f);
 654    PUSH_VALUE (p, double_to_double_int (p, u), A68_LONG_INT);
 655  }
 656  
 657  //! @brief OP + = (LONG INT, LONG INT) LONG INT
 658  
 659  void genie_add_double_int (NODE_T * p)
 660  {
 661    A68_LONG_INT i, j;
 662    POP_OBJECT (p, &j, A68_LONG_INT);
 663    POP_OBJECT (p, &i, A68_LONG_INT);
 664    PUSH_VALUE (p, double_sadd (p, VALUE (&i), VALUE (&j)), A68_LONG_INT);
 665  }
 666  
 667  //! @brief OP - = (LONG INT, LONG INT) LONG INT
 668  
 669  void genie_sub_double_int (NODE_T * p)
 670  {
 671    A68_LONG_INT i, j;
 672    POP_OBJECT (p, &j, A68_LONG_INT);
 673    POP_OBJECT (p, &i, A68_LONG_INT);
 674    PUSH_VALUE (p, double_ssub (p, VALUE (&i), VALUE (&j)), A68_LONG_INT);
 675  }
 676  
 677  //! @brief OP * = (LONG INT, LONG INT) LONG INT
 678  
 679  void genie_mul_double_int (NODE_T * p)
 680  {
 681    A68_LONG_INT i, j;
 682    POP_OBJECT (p, &j, A68_LONG_INT);
 683    POP_OBJECT (p, &i, A68_LONG_INT);
 684    PUSH_VALUE (p, double_smul (p, VALUE (&i), VALUE (&j)), A68_LONG_INT);
 685  }
 686  
 687  //! @brief OP / = (LONG INT, LONG INT) LONG INT
 688  
 689  void genie_over_double_int (NODE_T * p)
 690  {
 691    A68_LONG_INT i, j;
 692    POP_OBJECT (p, &j, A68_LONG_INT);
 693    POP_OBJECT (p, &i, A68_LONG_INT);
 694    PRELUDE_ERROR (D_ZERO (VALUE (&j)), p, ERROR_DIVISION_BY_ZERO, M_LONG_INT);
 695    PUSH_VALUE (p, double_sdiv (p, VALUE (&i), VALUE (&j), 0), A68_LONG_INT);
 696  }
 697  
 698  //! @brief OP MOD = (LONG INT, LONG INT) LONG INT
 699  
 700  void genie_mod_double_int (NODE_T * p)
 701  {
 702    A68_LONG_INT i, j;
 703    POP_OBJECT (p, &j, A68_LONG_INT);
 704    POP_OBJECT (p, &i, A68_LONG_INT);
 705    PRELUDE_ERROR (D_ZERO (VALUE (&j)), p, ERROR_DIVISION_BY_ZERO, M_LONG_INT);
 706    PUSH_VALUE (p, double_sdiv (p, VALUE (&i), VALUE (&j), 1), A68_LONG_INT);
 707  }
 708  
 709  //! @brief OP / = (LONG INT, LONG INT) LONG REAL
 710  
 711  void genie_div_double_int (NODE_T * p)
 712  {
 713    A68_LONG_INT i, j;
 714    POP_OBJECT (p, &j, A68_LONG_INT);
 715    POP_OBJECT (p, &i, A68_LONG_INT);
 716    PRELUDE_ERROR (D_ZERO (VALUE (&j)), p, ERROR_DIVISION_BY_ZERO, M_LONG_INT);
 717    DOUBLE_NUM_T u, v, w;
 718    v = double_int_to_double (p, VALUE (&j));
 719    u = double_int_to_double (p, VALUE (&i));
 720    w.f = u.f / v.f;
 721    PUSH_VALUE (p, w, A68_LONG_REAL);
 722  }
 723  
 724  //! @brief OP ** = (LONG INT, INT) INT
 725  
 726  void genie_pow_double_int_int (NODE_T * p)
 727  {
 728    A68_LONG_INT i; A68_INT j;
 729    POP_OBJECT (p, &j, A68_INT);
 730    PRELUDE_ERROR (VALUE (&j) < 0, p, ERROR_EXPONENT_INVALID, M_INT);
 731    POP_OBJECT (p, &i, A68_LONG_INT);
 732    DOUBLE_NUM_T mult = VALUE (&i), prod;
 733    set_lw (prod, 1);
 734    UNSIGNED_T top = (UNSIGNED_T) VALUE (&j), expo = 1;
 735    while (expo <= top) {
 736      if (expo & top) {
 737        prod = double_smul (p, prod, mult);
 738      }
 739      expo <<= 1;
 740      if (expo <= top) {
 741        mult = double_smul (p, mult, mult);
 742      }
 743    }
 744    PUSH_VALUE (p, prod, A68_LONG_INT);
 745  }
 746  
 747  //! @brief OP - = (LONG REAL) LONG REAL
 748  
 749  void genie_minus_double (NODE_T * p)
 750  {
 751    A68_LONG_REAL *u;
 752    POP_OPERAND_ADDRESS (p, u, A68_LONG_REAL);
 753    VALUE (u).f = -(VALUE (u).f);
 754  }
 755  
 756  //! @brief OP ABS = (LONG REAL) LONG REAL
 757  
 758  void genie_abs_double (NODE_T * p)
 759  {
 760    A68_LONG_REAL *u;
 761    POP_OPERAND_ADDRESS (p, u, A68_LONG_REAL);
 762    VALUE (u).f = fabs_double (VALUE (u).f);
 763  }
 764  
 765  //! @brief OP SIGN = (LONG REAL) INT
 766  
 767  void genie_sign_double (NODE_T * p)
 768  {
 769    A68_LONG_REAL u;
 770    POP_OBJECT (p, &u, A68_LONG_REAL);
 771    PUSH_VALUE (p, sign_double (VALUE (&u)), A68_INT);
 772  }
 773  
 774  //! @brief OP ** = (LONG REAL, INT) INT
 775  
 776  void genie_pow_double_int (NODE_T * p)
 777  {
 778    A68_INT j;
 779    POP_OBJECT (p, &j, A68_INT);
 780    INT_T top = (INT_T) VALUE (&j);
 781    A68_LONG_REAL z;
 782    POP_OBJECT (p, &z, A68_LONG_INT);
 783    DOUBLE_NUM_T mult, prod;
 784    prod.f = 1.0q;
 785    mult.f = VALUE (&z).f;
 786    int negative;
 787    if (top < 0) {
 788      top = -top;
 789      negative = A68_TRUE;
 790    } else {
 791      negative = A68_FALSE;
 792    }
 793    UNSIGNED_T expo = 1;
 794    while (expo <= top) {
 795      if (expo & top) {
 796        prod.f = prod.f * mult.f;
 797        CHECK_DOUBLE_REAL (p, prod.f);
 798      }
 799      expo <<= 1;
 800      if (expo <= top) {
 801        mult.f = mult.f * mult.f;
 802        CHECK_DOUBLE_REAL (p, mult.f);
 803      }
 804    }
 805    if (negative) {
 806      prod.f = 1.0q / prod.f;
 807    }
 808    PUSH_VALUE (p, prod, A68_LONG_REAL);
 809  }
 810  
 811  //! @brief OP ** = (LONG REAL, LONG REAL) LONG REAL
 812  
 813  void genie_pow_double (NODE_T * p)
 814  {
 815    A68_LONG_REAL x, y;
 816    POP_OBJECT (p, &y, A68_LONG_REAL);
 817    POP_OBJECT (p, &x, A68_LONG_REAL);
 818    errno = 0;
 819    PRELUDE_ERROR (VALUE (&x).f < 0.0q, p, ERROR_INVALID_ARGUMENT, M_LONG_REAL);
 820    DOUBLE_T z = 0.0q;
 821    if (VALUE (&x).f == 0.0q) {
 822      if (VALUE (&y).f < 0.0q) {
 823        errno = ERANGE;
 824        MATH_RTE (p, errno != 0, M_LONG_REAL, NO_TEXT);
 825      } else {
 826        z = (VALUE (&y).f == 0.0q ? 1.0q : 0.0q);
 827      }
 828    } else {
 829      z = exp_double (VALUE (&y).f * log_double (VALUE (&x).f));
 830      MATH_RTE (p, errno != 0, M_LONG_REAL, NO_TEXT);
 831    }
 832    PUSH_VALUE (p, dble (z), A68_LONG_REAL);
 833  }
 834  
 835  //! @brief OP + = (LONG REAL, LONG REAL) LONG REAL
 836  
 837  void genie_add_double (NODE_T * p)
 838  {
 839    A68_LONG_REAL u, v;
 840    POP_OBJECT (p, &v, A68_LONG_REAL);
 841    POP_OBJECT (p, &u, A68_LONG_REAL);
 842    DOUBLE_NUM_T w;
 843    w.f = VALUE (&u).f + VALUE (&v).f;
 844    CHECK_DOUBLE_REAL (p, w.f);
 845    PUSH_VALUE (p, w, A68_LONG_REAL);
 846  }
 847  
 848  //! @brief OP - = (LONG REAL, LONG REAL) LONG REAL
 849  
 850  void genie_sub_double (NODE_T * p)
 851  {
 852    A68_LONG_REAL u, v;
 853    POP_OBJECT (p, &v, A68_LONG_REAL);
 854    POP_OBJECT (p, &u, A68_LONG_REAL);
 855    DOUBLE_NUM_T w;
 856    w.f = VALUE (&u).f - VALUE (&v).f;
 857    CHECK_DOUBLE_REAL (p, w.f);
 858    PUSH_VALUE (p, w, A68_LONG_REAL);
 859  }
 860  
 861  //! @brief OP * = (LONG REAL, LONG REAL) LONG REAL
 862  
 863  void genie_mul_double (NODE_T * p)
 864  {
 865    A68_LONG_REAL u, v;
 866    POP_OBJECT (p, &v, A68_LONG_REAL);
 867    POP_OBJECT (p, &u, A68_LONG_REAL);
 868    DOUBLE_NUM_T w;
 869    w.f = VALUE (&u).f * VALUE (&v).f;
 870    CHECK_DOUBLE_REAL (p, w.f);
 871    PUSH_VALUE (p, w, A68_LONG_REAL);
 872  }
 873  
 874  //! @brief OP / = (LONG REAL, LONG REAL) LONG REAL
 875  
 876  void genie_over_double (NODE_T * p)
 877  {
 878    A68_LONG_REAL u, v;
 879    POP_OBJECT (p, &v, A68_LONG_REAL);
 880    POP_OBJECT (p, &u, A68_LONG_REAL);
 881    PRELUDE_ERROR (VALUE (&v).f == 0.0q, p, ERROR_DIVISION_BY_ZERO, M_LONG_REAL);
 882    DOUBLE_NUM_T w;
 883    w.f = VALUE (&u).f / VALUE (&v).f;
 884    PUSH_VALUE (p, w, A68_LONG_REAL);
 885  }
 886  
 887  //! @brief OP +:= = (REF LONG INT, LONG INT) REF LONG INT
 888  
 889  void genie_plusab_double_int (NODE_T * p)
 890  {
 891    genie_f_and_becomes (p, M_REF_LONG_INT, genie_add_double_int);
 892  }
 893  
 894  //! @brief OP -:= = (REF LONG INT, LONG INT) REF LONG INT
 895  
 896  void genie_minusab_double_int (NODE_T * p)
 897  {
 898    genie_f_and_becomes (p, M_REF_LONG_INT, genie_sub_double_int);
 899  }
 900  
 901  //! @brief OP *:= = (REF LONG INT, LONG INT) REF LONG INT
 902  
 903  void genie_timesab_double_int (NODE_T * p)
 904  {
 905    genie_f_and_becomes (p, M_REF_LONG_INT, genie_mul_double_int);
 906  }
 907  
 908  //! @brief OP %:= = (REF LONG INT, LONG INT) REF LONG INT
 909  
 910  void genie_overab_double_int (NODE_T * p)
 911  {
 912    genie_f_and_becomes (p, M_REF_LONG_INT, genie_over_double_int);
 913  }
 914  
 915  //! @brief OP %*:= = (REF LONG INT, LONG INT) REF LONG INT
 916  
 917  void genie_modab_double_int (NODE_T * p)
 918  {
 919    genie_f_and_becomes (p, M_REF_LONG_INT, genie_mod_double_int);
 920  }
 921  
 922  //! @brief OP +:= = (REF LONG REAL, LONG REAL) REF LONG REAL
 923  
 924  void genie_plusab_double (NODE_T * p)
 925  {
 926    genie_f_and_becomes (p, M_REF_LONG_REAL, genie_add_double);
 927  }
 928  
 929  //! @brief OP -:= = (REF LONG REAL, LONG REAL) REF LONG REAL
 930  
 931  void genie_minusab_double (NODE_T * p)
 932  {
 933    genie_f_and_becomes (p, M_REF_LONG_REAL, genie_sub_double);
 934  }
 935  
 936  //! @brief OP *:= = (REF LONG REAL, LONG REAL) REF LONG REAL
 937  
 938  void genie_timesab_double (NODE_T * p)
 939  {
 940    genie_f_and_becomes (p, M_REF_LONG_REAL, genie_mul_double);
 941  }
 942  
 943  //! @brief OP /:= = (REF LONG REAL, LONG REAL) REF LONG REAL
 944  
 945  void genie_divab_double (NODE_T * p)
 946  {
 947    genie_f_and_becomes (p, M_REF_LONG_REAL, genie_over_double);
 948  }
 949  
 950  // OP (LONG INT, LONG INT) BOOL.
 951  
 952  #define A68_CMP_INT(n, OP)\
 953  void n (NODE_T * p) {\
 954    A68_LONG_INT i, j;\
 955    POP_OBJECT (p, &j, A68_LONG_INT);\
 956    POP_OBJECT (p, &i, A68_LONG_INT);\
 957    int k = sign_double_int (double_ssub (p, VALUE (&i), VALUE (&j)));\
 958    PUSH_VALUE (p, (BOOL_T) (k OP 0), A68_BOOL);\
 959    }
 960  
 961  A68_CMP_INT (genie_eq_double_int, ==);
 962  A68_CMP_INT (genie_ne_double_int, !=);
 963  A68_CMP_INT (genie_lt_double_int, <);
 964  A68_CMP_INT (genie_gt_double_int, >);
 965  A68_CMP_INT (genie_le_double_int, <=);
 966  A68_CMP_INT (genie_ge_double_int, >=);
 967  
 968  // OP (LONG REAL, LONG REAL) BOOL.
 969  #define A68_CMP_REAL(n, OP)\
 970  void n (NODE_T * p) {\
 971    A68_LONG_REAL i, j;\
 972    POP_OBJECT (p, &j, A68_LONG_REAL);\
 973    POP_OBJECT (p, &i, A68_LONG_REAL);\
 974    PUSH_VALUE (p, (BOOL_T) (VALUE (&i).f OP VALUE (&j).f), A68_BOOL);\
 975    }
 976  
 977  A68_CMP_REAL (genie_eq_double, ==);
 978  A68_CMP_REAL (genie_ne_double, !=);
 979  A68_CMP_REAL (genie_lt_double, <);
 980  A68_CMP_REAL (genie_gt_double, >);
 981  A68_CMP_REAL (genie_le_double, <=);
 982  A68_CMP_REAL (genie_ge_double, >=);
 983  
 984  //! @brief OP NOT = (LONG BITS) LONG BITS
 985  
 986  void genie_not_double_bits (NODE_T * p)
 987  {
 988    A68_LONG_BITS i;
 989    POP_OBJECT (p, &i, A68_LONG_BITS);
 990    DOUBLE_NUM_T w;
 991    HW (w) = ~HW (VALUE (&i));
 992    LW (w) = ~LW (VALUE (&i));
 993    PUSH_VALUE (p, w, A68_LONG_BITS);
 994  }
 995  
 996  //! @brief OP = = (LONG BITS, LONG BITS) BOOL.
 997  
 998  void genie_eq_double_bits (NODE_T * p)
 999  {
1000    A68_LONG_BITS i, j;
1001    POP_OBJECT (p, &j, A68_LONG_BITS);
1002    POP_OBJECT (p, &i, A68_LONG_BITS);
1003    BOOL_T u = HW (VALUE (&i)) == HW (VALUE (&j));
1004    BOOL_T v = LW (VALUE (&i)) == LW (VALUE (&j));
1005    PUSH_VALUE (p, (BOOL_T) (u & v ? A68_TRUE : A68_FALSE), A68_BOOL);
1006  }
1007  
1008  //! @brief OP ~= = (LONG BITS, LONG BITS) BOOL.
1009  
1010  void genie_ne_double_bits (NODE_T * p)
1011  {
1012    A68_LONG_BITS i, j;
1013    POP_OBJECT (p, &j, A68_LONG_BITS);    // (i ~= j) == ~ (i = j) 
1014    POP_OBJECT (p, &i, A68_LONG_BITS);
1015    BOOL_T u = HW (VALUE (&i)) == HW (VALUE (&j));
1016    BOOL_T v = LW (VALUE (&i)) == LW (VALUE (&j));
1017    PUSH_VALUE (p, (BOOL_T) (u & v ? A68_FALSE : A68_TRUE), A68_BOOL);
1018  }
1019  
1020  //! @brief OP <= = (LONG BITS, LONG BITS) BOOL
1021  
1022  void genie_le_double_bits (NODE_T * p)
1023  {
1024    A68_LONG_BITS i, j;
1025    POP_OBJECT (p, &j, A68_LONG_BITS);
1026    POP_OBJECT (p, &i, A68_LONG_BITS);
1027    BOOL_T u = (HW (VALUE (&i)) | HW (VALUE (&j))) == HW (VALUE (&j));
1028    BOOL_T v = (LW (VALUE (&i)) | LW (VALUE (&j))) == LW (VALUE (&j));
1029    PUSH_VALUE (p, (BOOL_T) (u & v ? A68_TRUE : A68_FALSE), A68_BOOL);
1030  }
1031  
1032  //! @brief OP > = (LONG BITS, LONG BITS) BOOL
1033  
1034  void genie_gt_double_bits (NODE_T * p)
1035  {
1036    A68_LONG_BITS i, j;
1037    POP_OBJECT (p, &j, A68_LONG_BITS);    // (i > j) == ! (i <= j)
1038    POP_OBJECT (p, &i, A68_LONG_BITS);
1039    BOOL_T u = (HW (VALUE (&i)) | HW (VALUE (&j))) == HW (VALUE (&j));
1040    BOOL_T v = (LW (VALUE (&i)) | LW (VALUE (&j))) == LW (VALUE (&j));
1041    PUSH_VALUE (p, (BOOL_T) (u & v ? A68_FALSE : A68_TRUE), A68_BOOL);
1042  }
1043  
1044  //! @brief OP >= = (LONG BITS, LONG BITS) BOOL
1045  
1046  void genie_ge_double_bits (NODE_T * p)
1047  {
1048    A68_LONG_BITS i, j;
1049    POP_OBJECT (p, &j, A68_LONG_BITS);    // (i >= j) == (j <= i)
1050    POP_OBJECT (p, &i, A68_LONG_BITS);
1051    BOOL_T u = (HW (VALUE (&i)) | HW (VALUE (&j))) == HW (VALUE (&i));
1052    BOOL_T v = (LW (VALUE (&i)) | LW (VALUE (&j))) == LW (VALUE (&i));
1053    PUSH_VALUE (p, (BOOL_T) (u & v ? A68_TRUE : A68_FALSE), A68_BOOL);
1054  }
1055  
1056  //! @brief OP < = (LONG BITS, LONG BITS) BOOL
1057  
1058  void genie_lt_double_bits (NODE_T * p)
1059  {
1060    A68_LONG_BITS i, j;
1061    POP_OBJECT (p, &j, A68_LONG_BITS);    // (i < j) == ! (i >= j)
1062    POP_OBJECT (p, &i, A68_LONG_BITS);
1063    BOOL_T u = (HW (VALUE (&i)) | HW (VALUE (&j))) == HW (VALUE (&i));
1064    BOOL_T v = (LW (VALUE (&i)) | LW (VALUE (&j))) == LW (VALUE (&i));
1065    PUSH_VALUE (p, (BOOL_T) (u & v ? A68_FALSE : A68_TRUE), A68_BOOL);
1066  }
1067  
1068  //! @brief PROC long bits pack = ([] BOOL) BITS
1069  
1070  void genie_double_bits_pack (NODE_T * p)
1071  {
1072    A68_REF z;
1073    POP_REF (p, &z);
1074    CHECK_REF (p, z, M_ROW_BOOL);
1075    A68_ARRAY *arr; A68_TUPLE *tup;
1076    GET_DESCRIPTOR (arr, tup, &z);
1077    int size = ROW_SIZE (tup);
1078    PRELUDE_ERROR (size < 0 || size > A68_BITS_WIDTH, p, ERROR_OUT_OF_BOUNDS, M_ROW_BOOL);
1079    DOUBLE_NUM_T w;
1080    set_lw (w, 0x0);
1081    if (ROW_SIZE (tup) > 0) {
1082      UNSIGNED_T bit = 0x0;
1083      BYTE_T *base = DEREF (BYTE_T, &ARRAY (arr));
1084      int n = 0;
1085      for (int k = UPB (tup); k >= LWB (tup); k--) {
1086        A68_BOOL *boo = (A68_BOOL *) & (base[INDEX_1_DIM (arr, tup, k)]);
1087        CHECK_INIT (p, INITIALISED (boo), M_BOOL);
1088        if (n == 0 || n == A68_BITS_WIDTH) {
1089          bit = 0x1;
1090        }
1091        if (VALUE (boo)) {
1092          if (n > A68_BITS_WIDTH) {
1093            LW (w) |= bit;
1094          } else {
1095            HW (w) |= bit;
1096          };
1097        }
1098        n++;
1099        bit <<= 1;
1100      }
1101    }
1102    PUSH_VALUE (p, w, A68_LONG_BITS);
1103  }
1104  
1105  //! @brief OP AND = (LONG BITS, LONG BITS) LONG BITS
1106  
1107  void genie_and_double_bits (NODE_T * p)
1108  {
1109    A68_LONG_BITS i, j;
1110    POP_OBJECT (p, &j, A68_LONG_BITS);
1111    POP_OBJECT (p, &i, A68_LONG_BITS);
1112    DOUBLE_NUM_T w;
1113    HW (w) = HW (VALUE (&i)) & HW (VALUE (&j));
1114    LW (w) = LW (VALUE (&i)) & LW (VALUE (&j));
1115    PUSH_VALUE (p, w, A68_LONG_BITS);
1116  }
1117  
1118  //! @brief OP OR = (LONG BITS, LONG BITS) LONG BITS
1119  
1120  void genie_or_double_bits (NODE_T * p)
1121  {
1122    A68_LONG_BITS i, j;
1123    POP_OBJECT (p, &j, A68_LONG_BITS);
1124    POP_OBJECT (p, &i, A68_LONG_BITS);
1125    DOUBLE_NUM_T w;
1126    HW (w) = HW (VALUE (&i)) | HW (VALUE (&j));
1127    LW (w) = LW (VALUE (&i)) | LW (VALUE (&j));
1128    PUSH_VALUE (p, w, A68_LONG_BITS);
1129  }
1130  
1131  //! @brief OP XOR = (LONG BITS, LONG BITS) LONG BITS
1132  
1133  void genie_xor_double_bits (NODE_T * p)
1134  {
1135    A68_LONG_BITS i, j;
1136    POP_OBJECT (p, &j, A68_LONG_BITS);
1137    POP_OBJECT (p, &i, A68_LONG_BITS);
1138    DOUBLE_NUM_T w;
1139    HW (w) = HW (VALUE (&i)) ^ HW (VALUE (&j));
1140    LW (w) = LW (VALUE (&i)) ^ LW (VALUE (&j));
1141    PUSH_VALUE (p, w, A68_LONG_BITS);
1142  }
1143  
1144  //! @brief OP + = (LONG BITS, LONG BITS) LONG BITS
1145  
1146  void genie_add_double_bits (NODE_T * p)
1147  {
1148    A68_LONG_BITS i, j;
1149    POP_OBJECT (p, &j, A68_LONG_BITS);
1150    POP_OBJECT (p, &i, A68_LONG_BITS);
1151    DOUBLE_NUM_T w;
1152    add_double (p, M_LONG_BITS, w, VALUE (&i), VALUE (&j));
1153    PUSH_VALUE (p, w, A68_LONG_BITS);
1154  }
1155  
1156  //! @brief OP - = (LONG BITS, LONG BITS) LONG BITS
1157  
1158  void genie_sub_double_bits (NODE_T * p)
1159  {
1160    A68_LONG_BITS i, j;
1161    POP_OBJECT (p, &j, A68_LONG_BITS);
1162    POP_OBJECT (p, &i, A68_LONG_BITS);
1163    DOUBLE_NUM_T w;
1164    sub_double (p, M_LONG_BITS, w, VALUE (&i), VALUE (&j));
1165    PUSH_VALUE (p, w, A68_LONG_BITS);
1166  }
1167  
1168  //! @brief OP * = (LONG BITS, LONG BITS) LONG BITS
1169  
1170  void genie_times_double_bits (NODE_T * p)
1171  {
1172    A68_LONG_BITS i, j;
1173    POP_OBJECT (p, &j, A68_LONG_BITS);
1174    POP_OBJECT (p, &i, A68_LONG_BITS);
1175    DOUBLE_NUM_T w = double_umul (p, M_LONG_BITS, VALUE (&i), VALUE (&j));
1176    PUSH_VALUE (p, w, A68_LONG_BITS);
1177  }
1178  
1179  //! @brief OP OVER = (LONG BITS, LONG BITS) LONG BITS
1180  
1181  void genie_over_double_bits (NODE_T * p)
1182  {
1183    A68_LONG_BITS i, j;
1184    POP_OBJECT (p, &j, A68_LONG_BITS);
1185    POP_OBJECT (p, &i, A68_LONG_BITS);
1186    DOUBLE_NUM_T w = double_udiv (p, M_LONG_BITS, VALUE (&i), VALUE (&j), 0);
1187    PUSH_VALUE (p, w, A68_LONG_BITS);
1188  }
1189  
1190  //! @brief OP MOD = (LONG BITS, LONG BITS) LONG BITS
1191  
1192  void genie_mod_double_bits (NODE_T * p)
1193  {
1194    A68_LONG_BITS i, j;
1195    DOUBLE_NUM_T w;
1196    POP_OBJECT (p, &j, A68_LONG_BITS);
1197    POP_OBJECT (p, &i, A68_LONG_BITS);
1198    w = double_udiv (p, M_LONG_BITS, VALUE (&i), VALUE (&j), 1);
1199    PUSH_VALUE (p, w, A68_LONG_BITS);
1200  }
1201  
1202  //! @brief OP ELEM = (INT, LONG BITS) BOOL
1203  
1204  void genie_elem_double_bits (NODE_T * p)
1205  {
1206    A68_LONG_BITS j; A68_INT i;
1207    POP_OBJECT (p, &j, A68_LONG_BITS);
1208    POP_OBJECT (p, &i, A68_INT);
1209    int k = VALUE (&i);
1210    PRELUDE_ERROR (k < 1 || k > A68_LONG_BITS_WIDTH, p, ERROR_OUT_OF_BOUNDS, M_INT);
1211    UNSIGNED_T mask = 0x1, *w;
1212    if (k <= A68_BITS_WIDTH) {
1213      w = &(LW (VALUE (&j)));
1214    } else {
1215      w = &(HW (VALUE (&j)));
1216      k -= A68_BITS_WIDTH;
1217    }
1218    for (int n = 0; n < (A68_BITS_WIDTH - VALUE (&i)); n++) {
1219      mask = mask << 1;
1220    }
1221    PUSH_VALUE (p, (BOOL_T) ((*w & mask) ? A68_TRUE : A68_FALSE), A68_BOOL);
1222  }
1223  
1224  //! @brief OP SET = (INT, LONG BITS) LONG BITS
1225  
1226  void genie_set_double_bits (NODE_T * p)
1227  {
1228    A68_LONG_BITS j; A68_INT i;
1229    POP_OBJECT (p, &j, A68_LONG_BITS);
1230    POP_OBJECT (p, &i, A68_INT);
1231    int k = VALUE (&i);
1232    PRELUDE_ERROR (k < 1 || k > A68_LONG_BITS_WIDTH, p, ERROR_OUT_OF_BOUNDS, M_INT);
1233    UNSIGNED_T mask = 0x1, *w;
1234    if (k <= A68_BITS_WIDTH) {
1235      w = &(LW (VALUE (&j)));
1236    } else {
1237      w = &(HW (VALUE (&j)));
1238      k -= A68_BITS_WIDTH;
1239    }
1240    for (int n = 0; n < (A68_BITS_WIDTH - VALUE (&i)); n++) {
1241      mask = mask << 1;
1242    }
1243    (*w) |= mask;
1244    PUSH_OBJECT (p, j, A68_LONG_BITS);
1245  }
1246  
1247  //! @brief OP CLEAR = (INT, LONG BITS) LONG BITS
1248  
1249  void genie_clear_double_bits (NODE_T * p)
1250  {
1251    A68_LONG_BITS j; A68_INT i;
1252    POP_OBJECT (p, &j, A68_LONG_BITS);
1253    POP_OBJECT (p, &i, A68_INT);
1254    int k = VALUE (&i);
1255    PRELUDE_ERROR (k < 1 || k > A68_LONG_BITS_WIDTH, p, ERROR_OUT_OF_BOUNDS, M_INT);
1256    UNSIGNED_T mask = 0x1, *w;
1257    if (k <= A68_BITS_WIDTH) {
1258      w = &(LW (VALUE (&j)));
1259    } else {
1260      w = &(HW (VALUE (&j)));
1261      k -= A68_BITS_WIDTH;
1262    }
1263    for (int n = 0; n < (A68_BITS_WIDTH - VALUE (&i)); n++) {
1264      mask = mask << 1;
1265    }
1266    (*w) &= ~mask;
1267    PUSH_OBJECT (p, j, A68_LONG_BITS);
1268  }
1269  
1270  //! @brief OP SHL = (LONG BITS, INT) LONG BITS
1271  
1272  void genie_shl_double_bits (NODE_T * p)
1273  {
1274    A68_LONG_BITS i; A68_INT j;
1275    POP_OBJECT (p, &j, A68_INT);
1276    POP_OBJECT (p, &i, A68_LONG_BITS);
1277    DOUBLE_NUM_T *w = &VALUE (&i);
1278    int k = VALUE (&j);
1279    if (VALUE (&j) >= 0) {
1280      for (int n = 0; n < k; n++) {
1281        UNSIGNED_T carry = ((LW (*w) & D_SIGN) ? 0x1 : 0x0);
1282        PRELUDE_ERROR (MODCHK (p, M_LONG_BITS, HW (*w) | D_SIGN), p, ERROR_MATH, M_LONG_BITS);
1283        HW (*w) = (HW (*w) << 1) | carry;
1284        LW (*w) = (LW (*w) << 1);
1285      }
1286    } else {
1287      k = -k;
1288      for (int n = 0; n < k; n++) {
1289        UNSIGNED_T carry = ((HW (*w) & 0x1) ? D_SIGN : 0x0);
1290        HW (*w) = (HW (*w) >> 1);
1291        LW (*w) = (LW (*w) >> 1) | carry;
1292      }
1293    }
1294    PUSH_OBJECT (p, i, A68_LONG_BITS);
1295  }
1296  
1297  //! @brief OP SHR = (LONG BITS, INT) LONG BITS
1298  
1299  void genie_shr_double_bits (NODE_T * p)
1300  {
1301    A68_INT *j;
1302    POP_OPERAND_ADDRESS (p, j, A68_INT);
1303    VALUE (j) = -VALUE (j);
1304    genie_shl_double_bits (p);    // Conform RR
1305  }
1306  
1307  //! @brief OP ROL = (LONG BITS, INT) LONG BITS
1308  
1309  void genie_rol_double_bits (NODE_T * p)
1310  {
1311    A68_LONG_BITS i; A68_INT j;
1312    DOUBLE_NUM_T *w = &VALUE (&i);
1313    POP_OBJECT (p, &j, A68_INT);
1314    POP_OBJECT (p, &i, A68_LONG_BITS);
1315    int k = VALUE (&j);
1316    if (k >= 0) {
1317      for (int n = 0; n < k; n++) {
1318        UNSIGNED_T carry = ((HW (*w) & D_SIGN) ? 0x1 : 0x0);
1319        UNSIGNED_T carry_between = ((LW (*w) & D_SIGN) ? 0x1 : 0x0);
1320        HW (*w) = (HW (*w) << 1) | carry_between;
1321        LW (*w) = (LW (*w) << 1) | carry;
1322      }
1323    } else {
1324      k = -k;
1325      for (int n = 0; n < k; n++) {
1326        UNSIGNED_T carry = ((LW (*w) & 0x1) ? D_SIGN : 0x0);
1327        UNSIGNED_T carry_between = ((HW (*w) & 0x1) ? D_SIGN : 0x0);
1328        HW (*w) = (HW (*w) >> 1) | carry;
1329        LW (*w) = (LW (*w) >> 1) | carry_between;
1330      }
1331    }
1332    PUSH_OBJECT (p, i, A68_LONG_BITS);
1333  }
1334  
1335  //! @brief OP ROR = (LONG BITS, INT) LONG BITS
1336  
1337  void genie_ror_double_bits (NODE_T * p)
1338  {
1339    A68_INT *j;
1340    POP_OPERAND_ADDRESS (p, j, A68_INT);
1341    VALUE (j) = -VALUE (j);
1342    genie_rol_double_bits (p);    // Conform RR
1343  }
1344  
1345  //! @brief OP BIN = (LONG INT) LONG BITS
1346  
1347  void genie_bin_double_int (NODE_T * p)
1348  {
1349    A68_LONG_INT i;
1350    POP_OBJECT (p, &i, A68_LONG_INT);
1351  // RR does not convert negative numbers
1352    if (D_NEG (VALUE (&i))) {
1353      errno = EDOM;
1354      diagnostic (A68_RUNTIME_ERROR, p, ERROR_OUT_OF_BOUNDS, M_BITS);
1355      exit_genie (p, A68_RUNTIME_ERROR);
1356    }
1357    PUSH_OBJECT (p, i, A68_LONG_BITS);
1358  }
1359  
1360  //! @brief OP +* = (LONG REAL, LONG REAL) LONG COMPLEX
1361  
1362  void genie_i_double_compl (NODE_T * p)
1363  {
1364    (void) p;
1365  }
1366  
1367  //! @brief OP SHORTEN = (LONG COMPLEX) COMPLEX
1368  
1369  void genie_shorten_double_compl_to_complex (NODE_T * p)
1370  {
1371    A68_LONG_REAL re, im;
1372    POP_OBJECT (p, &im, A68_LONG_REAL);
1373    POP_OBJECT (p, &re, A68_LONG_REAL);
1374    REAL_T w = VALUE (&re).f;
1375    PUSH_VALUE (p, w, A68_REAL);
1376    w = VALUE (&im).f;
1377    PUSH_VALUE (p, w, A68_REAL);
1378  }
1379  
1380  //! @brief OP LENG = (LONG COMPLEX) LONG LONG COMPLEX
1381  
1382  void genie_lengthen_double_compl_to_long_mp_complex (NODE_T * p)
1383  {
1384    int digs = DIGITS (M_LONG_LONG_REAL);
1385    A68_LONG_REAL re, im;
1386    POP_OBJECT (p, &im, A68_LONG_REAL);
1387    POP_OBJECT (p, &re, A68_LONG_REAL);
1388    MP_T *z = nil_mp (p, digs);
1389    (void) double_to_mp (p, z, VALUE (&re).f, digs);
1390    MP_STATUS (z) = (MP_T) INIT_MASK;
1391    z = nil_mp (p, digs);
1392    (void) double_to_mp (p, z, VALUE (&im).f, digs);
1393    MP_STATUS (z) = (MP_T) INIT_MASK;
1394  }
1395  
1396  //! @brief OP +* = (LONG INT, LONG INT) LONG COMPLEX
1397  
1398  void genie_i_int_double_compl (NODE_T * p)
1399  {
1400    A68_LONG_INT re, im;
1401    POP_OBJECT (p, &im, A68_LONG_INT);
1402    POP_OBJECT (p, &re, A68_LONG_INT);
1403    PUSH_VALUE (p, double_int_to_double (p, VALUE (&re)), A68_LONG_REAL);
1404    PUSH_VALUE (p, double_int_to_double (p, VALUE (&im)), A68_LONG_REAL);
1405  }
1406  
1407  //! @brief OP RE = (LONG COMPLEX) LONG REAL
1408  
1409  void genie_re_double_compl (NODE_T * p)
1410  {
1411    DECREMENT_STACK_POINTER (p, SIZE (M_LONG_REAL));
1412  }
1413  
1414  //! @brief OP IM = (LONG COMPLEX) LONG REAL
1415  
1416  void genie_im_double_compl (NODE_T * p)
1417  {
1418    A68_LONG_REAL re, im;
1419    POP_OBJECT (p, &im, A68_LONG_REAL);
1420    POP_OBJECT (p, &re, A68_LONG_REAL);
1421    PUSH_OBJECT (p, im, A68_LONG_REAL);
1422  }
1423  
1424  //! @brief OP - = (LONG COMPLEX) LONG COMPLEX
1425  
1426  void genie_minus_double_compl (NODE_T * p)
1427  {
1428    A68_LONG_REAL re, im;
1429    POP_OBJECT (p, &im, A68_LONG_REAL);
1430    POP_OBJECT (p, &re, A68_LONG_REAL);
1431    VALUE (&re).f = -VALUE (&re).f;
1432    VALUE (&im).f = -VALUE (&im).f;
1433    PUSH_OBJECT (p, im, A68_LONG_REAL);
1434    PUSH_OBJECT (p, re, A68_LONG_REAL);
1435  }
1436  
1437  //! @brief OP ABS = (LONG COMPLEX) LONG REAL
1438  
1439  void genie_abs_double_compl (NODE_T * p)
1440  {
1441    A68_LONG_REAL re, im;
1442    POP_LONG_COMPLEX (p, &re, &im);
1443    PUSH_VALUE (p, dble (a68_hypot_double (VALUE (&re).f, VALUE (&im).f)), A68_LONG_REAL);
1444  }
1445  
1446  //! @brief OP ARG = (LONG COMPLEX) LONG REAL
1447  
1448  void genie_arg_double_compl (NODE_T * p)
1449  {
1450    A68_LONG_REAL re, im;
1451    POP_LONG_COMPLEX (p, &re, &im);
1452    PRELUDE_ERROR (VALUE (&re).f == 0.0q && VALUE (&im).f == 0.0q, p, ERROR_INVALID_ARGUMENT, M_LONG_COMPLEX);
1453    PUSH_VALUE (p, dble (atan2_double (VALUE (&im).f, VALUE (&re).f)), A68_LONG_REAL);
1454  }
1455  
1456  //! @brief OP CONJ = (LONG COMPLEX) LONG COMPLEX
1457  
1458  void genie_conj_double_compl (NODE_T * p)
1459  {
1460    A68_LONG_REAL im;
1461    POP_OBJECT (p, &im, A68_LONG_REAL);
1462    VALUE (&im).f = -VALUE (&im).f;
1463    PUSH_OBJECT (p, im, A68_LONG_REAL);
1464  }
1465  
1466  //! @brief OP + = (COMPLEX, COMPLEX) COMPLEX
1467  
1468  void genie_add_double_compl (NODE_T * p)
1469  {
1470    A68_LONG_REAL re_x, im_x, re_y, im_y;
1471    POP_LONG_COMPLEX (p, &re_y, &im_y);
1472    POP_LONG_COMPLEX (p, &re_x, &im_x);
1473    VALUE (&re_x).f += VALUE (&re_y).f;
1474    VALUE (&im_x).f += VALUE (&im_y).f;
1475    CHECK_DOUBLE_COMPLEX (p, VALUE (&im_x).f, VALUE (&im_y).f);
1476    PUSH_OBJECT (p, re_x, A68_LONG_REAL);
1477    PUSH_OBJECT (p, im_x, A68_LONG_REAL);
1478  }
1479  
1480  //! @brief OP - = (COMPLEX, COMPLEX) COMPLEX
1481  
1482  void genie_sub_double_compl (NODE_T * p)
1483  {
1484    A68_LONG_REAL re_x, im_x, re_y, im_y;
1485    POP_LONG_COMPLEX (p, &re_y, &im_y);
1486    POP_LONG_COMPLEX (p, &re_x, &im_x);
1487    VALUE (&re_x).f -= VALUE (&re_y).f;
1488    VALUE (&im_x).f -= VALUE (&im_y).f;
1489    CHECK_DOUBLE_COMPLEX (p, VALUE (&im_x).f, VALUE (&im_y).f);
1490    PUSH_OBJECT (p, re_x, A68_LONG_REAL);
1491    PUSH_OBJECT (p, im_x, A68_LONG_REAL);
1492  }
1493  
1494  //! @brief OP * = (COMPLEX, COMPLEX) COMPLEX
1495  
1496  void genie_mul_double_compl (NODE_T * p)
1497  {
1498    A68_LONG_REAL re_x, im_x, re_y, im_y;
1499    POP_LONG_COMPLEX (p, &re_y, &im_y);
1500    POP_LONG_COMPLEX (p, &re_x, &im_x);
1501    DOUBLE_T re = VALUE (&re_x).f * VALUE (&re_y).f - VALUE (&im_x).f * VALUE (&im_y).f;
1502    DOUBLE_T im = VALUE (&im_x).f * VALUE (&re_y).f + VALUE (&re_x).f * VALUE (&im_y).f;
1503    CHECK_DOUBLE_COMPLEX (p, VALUE (&im_x).f, VALUE (&im_y).f);
1504    PUSH_VALUE (p, dble (re), A68_LONG_REAL);
1505    PUSH_VALUE (p, dble (im), A68_LONG_REAL);
1506  }
1507  
1508  //! @brief OP / = (COMPLEX, COMPLEX) COMPLEX
1509  
1510  void genie_div_double_compl (NODE_T * p)
1511  {
1512    A68_LONG_REAL re_x, im_x, re_y, im_y;
1513    DOUBLE_T re = 0.0, im = 0.0;
1514    POP_LONG_COMPLEX (p, &re_y, &im_y);
1515    POP_LONG_COMPLEX (p, &re_x, &im_x);
1516    PRELUDE_ERROR (VALUE (&re_y).f == 0.0q && VALUE (&im_y).f == 0.0q, p, ERROR_DIVISION_BY_ZERO, M_LONG_COMPLEX);
1517    if (ABSQ (VALUE (&re_y).f) >= ABSQ (VALUE (&im_y).f)) {
1518      DOUBLE_T r = VALUE (&im_y).f / VALUE (&re_y).f, den = VALUE (&re_y).f + r * VALUE (&im_y).f;
1519      re = (VALUE (&re_x).f + r * VALUE (&im_x).f) / den;
1520      im = (VALUE (&im_x).f - r * VALUE (&re_x).f) / den;
1521    } else {
1522      DOUBLE_T r = VALUE (&re_y).f / VALUE (&im_y).f, den = VALUE (&im_y).f + r * VALUE (&re_y).f;
1523      re = (VALUE (&re_x).f * r + VALUE (&im_x).f) / den;
1524      im = (VALUE (&im_x).f * r - VALUE (&re_x).f) / den;
1525    }
1526    PUSH_VALUE (p, dble (re), A68_LONG_REAL);
1527    PUSH_VALUE (p, dble (im), A68_LONG_REAL);
1528  }
1529  
1530  //! @brief OP ** = (LONG COMPLEX, INT) LONG COMPLEX
1531  
1532  void genie_pow_double_compl_int (NODE_T * p)
1533  {
1534    A68_LONG_REAL re_x, im_x;
1535    A68_INT j;
1536    BOOL_T negative;
1537    POP_OBJECT (p, &j, A68_INT);
1538    POP_LONG_COMPLEX (p, &re_x, &im_x);
1539    DOUBLE_T re_z = 1.0q, im_z = 0.0q;
1540    DOUBLE_T re_y = VALUE (&re_x).f, im_y = VALUE (&im_x).f;
1541    negative = (BOOL_T) (VALUE (&j) < 0);
1542    if (negative) {
1543      VALUE (&j) = -VALUE (&j);
1544    }
1545    INT_T expo = 1;
1546    while ((UNSIGNED_T) expo <= (UNSIGNED_T) (VALUE (&j))) {
1547      DOUBLE_T z;
1548      if (expo & VALUE (&j)) {
1549        z = re_z * re_y - im_z * im_y;
1550        im_z = re_z * im_y + im_z * re_y;
1551        re_z = z;
1552      }
1553      z = re_y * re_y - im_y * im_y;
1554      im_y = im_y * re_y + re_y * im_y;
1555      re_y = z;
1556      CHECK_DOUBLE_COMPLEX (p, re_y, im_y);
1557      CHECK_DOUBLE_COMPLEX (p, re_z, im_z);
1558      expo <<= 1;
1559    }
1560    if (negative) {
1561      PUSH_VALUE (p, dble (1.0q), A68_LONG_REAL);
1562      PUSH_VALUE (p, dble (0.0q), A68_LONG_REAL);
1563      PUSH_VALUE (p, dble (re_z), A68_LONG_REAL);
1564      PUSH_VALUE (p, dble (im_z), A68_LONG_REAL);
1565      genie_div_double_compl (p);
1566    } else {
1567      PUSH_VALUE (p, dble (re_z), A68_LONG_REAL);
1568      PUSH_VALUE (p, dble (im_z), A68_LONG_REAL);
1569    }
1570  }
1571  
1572  //! @brief OP = = (COMPLEX, COMPLEX) BOOL
1573  
1574  void genie_eq_double_compl (NODE_T * p)
1575  {
1576    A68_LONG_REAL re_x, im_x, re_y, im_y;
1577    POP_LONG_COMPLEX (p, &re_y, &im_y);
1578    POP_LONG_COMPLEX (p, &re_x, &im_x);
1579    PUSH_VALUE (p, (BOOL_T) ((VALUE (&re_x).f == VALUE (&re_y).f) && (VALUE (&im_x).f == VALUE (&im_y).f)), A68_BOOL);
1580  }
1581  
1582  //! @brief OP /= = (COMPLEX, COMPLEX) BOOL
1583  
1584  void genie_ne_double_compl (NODE_T * p)
1585  {
1586    A68_LONG_REAL re_x, im_x, re_y, im_y;
1587    POP_LONG_COMPLEX (p, &re_y, &im_y);
1588    POP_LONG_COMPLEX (p, &re_x, &im_x);
1589    PUSH_VALUE (p, (BOOL_T) ! ((VALUE (&re_x).f == VALUE (&re_y).f) && (VALUE (&im_x).f == VALUE (&im_y).f)), A68_BOOL);
1590  }
1591  
1592  //! @brief OP +:= = (REF COMPLEX, COMPLEX) REF COMPLEX
1593  
1594  void genie_plusab_double_compl (NODE_T * p)
1595  {
1596    genie_f_and_becomes (p, M_REF_LONG_COMPLEX, genie_add_double_compl);
1597  }
1598  
1599  //! @brief OP -:= = (REF COMPLEX, COMPLEX) REF COMPLEX
1600  
1601  void genie_minusab_double_compl (NODE_T * p)
1602  {
1603    genie_f_and_becomes (p, M_REF_LONG_COMPLEX, genie_sub_double_compl);
1604  }
1605  
1606  //! @brief OP *:= = (REF COMPLEX, COMPLEX) REF COMPLEX
1607  
1608  void genie_timesab_double_compl (NODE_T * p)
1609  {
1610    genie_f_and_becomes (p, M_REF_LONG_COMPLEX, genie_mul_double_compl);
1611  }
1612  
1613  //! @brief OP /:= = (REF COMPLEX, COMPLEX) REF COMPLEX
1614  
1615  void genie_divab_double_compl (NODE_T * p)
1616  {
1617    genie_f_and_becomes (p, M_REF_LONG_COMPLEX, genie_div_double_compl);
1618  }
1619  
1620  //! @brief OP LENG = (COMPLEX) LONG COMPLEX 
1621  
1622  void genie_lengthen_complex_to_double_compl (NODE_T * p)
1623  {
1624    A68_REAL i;
1625    POP_OBJECT (p, &i, A68_REAL);
1626    genie_lengthen_real_to_double (p);
1627    PUSH_OBJECT (p, i, A68_REAL);
1628    genie_lengthen_real_to_double (p);
1629  }
1630  
1631  // Functions
1632  
1633  #define CD_FUNCTION(name, fun)\
1634  void name (NODE_T * p) {\
1635    A68_LONG_REAL *x;\
1636    POP_OPERAND_ADDRESS (p, x, A68_LONG_REAL);\
1637    errno = 0;\
1638    VALUE (x).f = fun (VALUE (x).f);\
1639    MATH_RTE (p, errno != 0, M_LONG_REAL, NO_TEXT);\
1640  }
1641  
1642  CD_FUNCTION (genie_acos_double, acos_double);
1643  CD_FUNCTION (genie_acosh_double, acosh_double);
1644  CD_FUNCTION (genie_asinh_double, asinh_double);
1645  CD_FUNCTION (genie_atanh_double, atanh_double);
1646  CD_FUNCTION (genie_asin_double, asin_double);
1647  CD_FUNCTION (genie_atan_double, atan_double);
1648  CD_FUNCTION (genie_cosh_double, cosh_double);
1649  CD_FUNCTION (genie_cos_double, cos_double);
1650  CD_FUNCTION (genie_curt_double, cbrt_double);
1651  CD_FUNCTION (genie_exp_double, exp_double);
1652  CD_FUNCTION (genie_ln_double, log_double);
1653  CD_FUNCTION (genie_log_double, log10_double);
1654  CD_FUNCTION (genie_sinh_double, sinh_double);
1655  CD_FUNCTION (genie_sin_double, sin_double);
1656  CD_FUNCTION (genie_sqrt_double, sqrt_double);
1657  CD_FUNCTION (genie_tanh_double, tanh_double);
1658  CD_FUNCTION (genie_tan_double, tan_double);
1659  CD_FUNCTION (genie_erf_double, erf_double);
1660  CD_FUNCTION (genie_erfc_double, erfc_double);
1661  CD_FUNCTION (genie_lngamma_double, lgamma_double);
1662  CD_FUNCTION (genie_gamma_double, tgamma_double);
1663  CD_FUNCTION (genie_csc_double, a68_csc_double);
1664  CD_FUNCTION (genie_cscdg_double, a68_cscdg_double);
1665  CD_FUNCTION (genie_acsc_double, a68_acsc_double);
1666  CD_FUNCTION (genie_acscdg_double, a68_acscdg_double);
1667  CD_FUNCTION (genie_sec_double, a68_sec_double);
1668  CD_FUNCTION (genie_secdg_double, a68_secdg_double);
1669  CD_FUNCTION (genie_asec_double, a68_asec_double);
1670  CD_FUNCTION (genie_asecdg_double, a68_asecdg_double);
1671  CD_FUNCTION (genie_cot_double, a68_cot_double);
1672  CD_FUNCTION (genie_acot_double, a68_acot_double);
1673  CD_FUNCTION (genie_sindg_double, a68_sindg_double);
1674  CD_FUNCTION (genie_cas_double, a68_cas_double);
1675  CD_FUNCTION (genie_cosdg_double, a68_cosdg_double);
1676  CD_FUNCTION (genie_tandg_double, a68_tandg_double);
1677  CD_FUNCTION (genie_asindg_double, a68_asindg_double);
1678  CD_FUNCTION (genie_acosdg_double, a68_acosdg_double);
1679  CD_FUNCTION (genie_atandg_double, a68_atandg_double);
1680  CD_FUNCTION (genie_cotdg_double, a68_cotdg_double);
1681  CD_FUNCTION (genie_acotdg_double, a68_acotdg_double);
1682  CD_FUNCTION (genie_sinpi_double, a68_sinpi_double);
1683  CD_FUNCTION (genie_cospi_double, a68_cospi_double);
1684  CD_FUNCTION (genie_tanpi_double, a68_tanpi_double);
1685  CD_FUNCTION (genie_cotpi_double, a68_cotpi_double);
1686  
1687  //! @brief PROC long arctan2 = (LONG REAL) LONG REAL
1688  
1689  void genie_atan2_double (NODE_T * p)
1690  {
1691    A68_LONG_REAL x, y;
1692    POP_OBJECT (p, &y, A68_LONG_REAL);
1693    POP_OBJECT (p, &x, A68_LONG_REAL);
1694    errno = 0;
1695    PRELUDE_ERROR (VALUE (&x).f == 0.0q && VALUE (&y).f == 0.0q, p, ERROR_INVALID_ARGUMENT, M_LONG_REAL);
1696    VALUE (&x).f = a68_atan2_real (VALUE (&y).f, VALUE (&x).f);
1697    PRELUDE_ERROR (errno != 0, p, ERROR_MATH_EXCEPTION, NO_TEXT);
1698    PUSH_OBJECT (p, x, A68_LONG_REAL);
1699  }
1700  
1701  //! @brief PROC long arctan2dg = (LONG REAL) LONG REAL
1702  
1703  void genie_atan2dg_double (NODE_T * p)
1704  {
1705    A68_LONG_REAL x, y;
1706    POP_OBJECT (p, &y, A68_LONG_REAL);
1707    POP_OBJECT (p, &x, A68_LONG_REAL);
1708    errno = 0;
1709    PRELUDE_ERROR (VALUE (&x).f == 0.0q && VALUE (&y).f == 0.0q, p, ERROR_INVALID_ARGUMENT, M_LONG_REAL);
1710    VALUE (&x).f = CONST_180_OVER_PI_Q * a68_atan2_real (VALUE (&y).f, VALUE (&x).f);
1711    PRELUDE_ERROR (errno != 0, p, ERROR_MATH_EXCEPTION, NO_TEXT);
1712    PUSH_OBJECT (p, x, A68_LONG_REAL);
1713  }
1714  
1715  //! @brief PROC (LONG REAL) LONG REAL inverf
1716  
1717  void genie_inverf_double (NODE_T * _p_)
1718  {
1719    A68_LONG_REAL x;
1720    DOUBLE_T y, z;
1721    POP_OBJECT (_p_, &x, A68_LONG_REAL);
1722    errno = 0;
1723    y = VALUE (&x).f;
1724    z = inverf_double (y);
1725    MATH_RTE (_p_, errno != 0, M_LONG_REAL, NO_TEXT);
1726    CHECK_DOUBLE_REAL (_p_, z);
1727    PUSH_VALUE (_p_, dble (z), A68_LONG_REAL);
1728  }
1729  
1730  //! @brief PROC (LONG REAL) LONG REAL inverfc
1731  
1732  void genie_inverfc_double (NODE_T * p)
1733  {
1734    A68_LONG_REAL *u;
1735    POP_OPERAND_ADDRESS (p, u, A68_LONG_REAL);
1736    VALUE (u).f = 1.0q - (VALUE (u).f);
1737    genie_inverf_double (p);
1738  }
1739  
1740  #define CD_C_FUNCTION(p, g)\
1741    A68_LONG_REAL re, im;\
1742    DOUBLE_COMPLEX_T z;\
1743    POP_OBJECT (p, &im, A68_LONG_REAL);\
1744    POP_OBJECT (p, &re, A68_LONG_REAL);\
1745    errno = 0;\
1746    z = VALUE (&re).f + VALUE (&im).f * _Complex_I;\
1747    z = g (z);\
1748    PUSH_VALUE (p, dble ((DOUBLE_T) creal_double (z)), A68_LONG_REAL);\
1749    PUSH_VALUE (p, dble ((DOUBLE_T) cimag_double (z)), A68_LONG_REAL);\
1750    MATH_RTE (p, errno != 0, M_COMPLEX, NO_TEXT);
1751  
1752  //! @brief PROC long csqrt = (LONG COMPLEX) LONG COMPLEX
1753  
1754  void genie_sqrt_double_compl (NODE_T * p)
1755  {
1756    CD_C_FUNCTION (p, csqrt_double);
1757  }
1758  
1759  //! @brief PROC long csin = (LONG COMPLEX) LONG COMPLEX
1760  
1761  void genie_sin_double_compl (NODE_T * p)
1762  {
1763    CD_C_FUNCTION (p, csin_double);
1764  }
1765  
1766  //! @brief PROC long ccos = (LONG COMPLEX) LONG COMPLEX
1767  
1768  void genie_cos_double_compl (NODE_T * p)
1769  {
1770    CD_C_FUNCTION (p, ccos_double);
1771  }
1772  
1773  //! @brief PROC long ctan = (LONG COMPLEX) LONG COMPLEX
1774  
1775  void genie_tan_double_compl (NODE_T * p)
1776  {
1777    CD_C_FUNCTION (p, ctan_double);
1778  }
1779  
1780  //! @brief PROC long casin = (LONG COMPLEX) LONG COMPLEX
1781  
1782  void genie_asin_double_compl (NODE_T * p)
1783  {
1784    CD_C_FUNCTION (p, casin_double);
1785  }
1786  
1787  //! @brief PROC long cacos = (LONG COMPLEX) LONG COMPLEX
1788  
1789  void genie_acos_double_compl (NODE_T * p)
1790  {
1791    CD_C_FUNCTION (p, cacos_double);
1792  }
1793  
1794  //! @brief PROC long catan = (LONG COMPLEX) LONG COMPLEX
1795  
1796  void genie_atan_double_compl (NODE_T * p)
1797  {
1798    CD_C_FUNCTION (p, catan_double);
1799  }
1800  
1801  //! @brief PROC long cexp = (LONG COMPLEX) LONG COMPLEX
1802  
1803  void genie_exp_double_compl (NODE_T * p)
1804  {
1805    CD_C_FUNCTION (p, cexp_double);
1806  }
1807  
1808  //! @brief PROC long cln = (LONG COMPLEX) LONG COMPLEX
1809  
1810  void genie_ln_double_compl (NODE_T * p)
1811  {
1812    CD_C_FUNCTION (p, clog_double);
1813  }
1814  
1815  //! @brief PROC long csinh = (LONG COMPLEX) LONG COMPLEX
1816  
1817  void genie_sinh_double_compl (NODE_T * p)
1818  {
1819    CD_C_FUNCTION (p, csinh_double);
1820  }
1821  
1822  //! @brief PROC long ccosh = (LONG COMPLEX) LONG COMPLEX
1823  
1824  void genie_cosh_double_compl (NODE_T * p)
1825  {
1826    CD_C_FUNCTION (p, ccosh_double);
1827  }
1828  
1829  //! @brief PROC long ctanh = (LONG COMPLEX) LONG COMPLEX
1830  
1831  void genie_tanh_double_compl (NODE_T * p)
1832  {
1833    CD_C_FUNCTION (p, ctanh_double);
1834  }
1835  
1836  //! @brief PROC long casinh = (LONG COMPLEX) LONG COMPLEX
1837  
1838  void genie_asinh_double_compl (NODE_T * p)
1839  {
1840    CD_C_FUNCTION (p, casinh_double);
1841  }
1842  
1843  //! @brief PROC long cacosh = (LONG COMPLEX) LONG COMPLEX
1844  
1845  void genie_acosh_double_compl (NODE_T * p)
1846  {
1847    CD_C_FUNCTION (p, cacosh_double);
1848  }
1849  
1850  //! @brief PROC long catanh = (LONG COMPLEX) LONG COMPLEX
1851  
1852  void genie_atanh_double_compl (NODE_T * p)
1853  {
1854    CD_C_FUNCTION (p, catanh_double);
1855  }
1856  
1857  //! @brief PROC next long random = LONG REAL
1858  
1859  void genie_next_random_double (NODE_T * p)
1860  {
1861  // This is 'real width' digits only.
1862    genie_next_random (p);
1863    genie_lengthen_real_to_double (p);
1864  }
1865  
1866  #define CALL(g, x, y) {\
1867    ADDR_T pop_sp = A68_SP;\
1868    A68_LONG_REAL *z = (A68_LONG_REAL *) STACK_TOP;\
1869    DOUBLE_NUM_T _w_;\
1870    _w_.f = (x);\
1871    PUSH_VALUE (_p_, _w_, A68_LONG_REAL);\
1872    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);\
1873    (y) = VALUE (z).f;\
1874    A68_SP = pop_sp;\
1875  }
1876  
1877  //! @brief Transform string into real-16.
1878  
1879  DOUBLE_T string_to_double (char *s, char **end)
1880  {
1881    errno = 0;
1882    DOUBLE_T y[A68_DOUBLE_DIG];
1883    for (int i = 0; i < A68_DOUBLE_DIG; i++) {
1884      y[i] = 0.0q;
1885    }
1886    while (IS_SPACE (s[0])) {
1887      s++;
1888    }
1889  // Scan mantissa digits and put them into "y".
1890    DOUBLE_T W;
1891    if (s[0] == '-') {
1892      W = -1.0q;
1893    } else {
1894      W = 1.0q;
1895    }
1896    if (s[0] == '+' || s[0] == '-') {
1897      s++;
1898    }
1899    while (s[0] == '0') {
1900      s++;
1901    }
1902    int dot = -1, pos = 0, pow = 0;
1903    while (pow < A68_DOUBLE_DIG && s[pos] != NULL_CHAR && (IS_DIGIT (s[pos]) || s[pos] == POINT_CHAR)) {
1904      if (s[pos] == POINT_CHAR) {
1905        dot = pos;
1906      } else {
1907        int val = (int) s[pos] - (int) '0';
1908        y[pow] = W * val;
1909        W /= 10.0q;
1910        pow++;
1911      }
1912      pos++;
1913    }
1914    (*end) = &(s[pos]);
1915  // Sum from low to high to preserve precision.
1916    DOUBLE_T sum = 0.0q;
1917    for (int i = A68_DOUBLE_DIG - 1; i >= 0; i--) {
1918      sum = sum + y[i];
1919    }
1920  // See if there is an exponent.
1921    int expo;
1922    if (s[pos] != NULL_CHAR && TO_UPPER (s[pos]) == TO_UPPER (EXPONENT_CHAR)) {
1923      expo = (int) strtol (&(s[++pos]), end, 10);
1924    } else {
1925      expo = 0;
1926    }
1927  // Standardise.
1928    if (dot >= 0) {
1929      expo += dot - 1;
1930    } else {
1931      expo += pow - 1;
1932    }
1933    while (sum != 0.0q && fabs_double (sum) < 1.0q) {
1934      sum *= 10.0q;
1935      expo -= 1;
1936    }
1937    if (errno == 0) {
1938      return sum * ten_up_double (expo);
1939    } else {
1940      return 0.0q;
1941    }
1942  }
1943  
1944  void genie_beta_inc_cf_double (NODE_T * p)
1945  {
1946    A68_LONG_REAL x, s, t;
1947    POP_OBJECT (p, &x, A68_LONG_REAL);
1948    POP_OBJECT (p, &t, A68_LONG_REAL);
1949    POP_OBJECT (p, &s, A68_LONG_REAL);
1950    errno = 0;
1951    PUSH_VALUE (p, dble (a68_beta_inc_double (VALUE (&s).f, VALUE (&t).f, VALUE (&x).f)), A68_LONG_REAL);
1952    MATH_RTE (p, errno != 0, M_LONG_REAL, NO_TEXT);
1953  }
1954  
1955  void genie_beta_double (NODE_T * p)
1956  {
1957    A68_LONG_REAL a, b;
1958    POP_OBJECT (p, &b, A68_LONG_REAL);
1959    POP_OBJECT (p, &a, A68_LONG_REAL);
1960    errno = 0;
1961    PUSH_VALUE (p, dble (exp_double (lgamma_double (VALUE (&a).f) + lgamma_double (VALUE (&b).f) - lgamma_double (VALUE (&a).f + VALUE (&b).f))), A68_LONG_REAL);
1962    MATH_RTE (p, errno != 0, M_LONG_REAL, NO_TEXT);
1963  }
1964  
1965  void genie_ln_beta_double (NODE_T * p)
1966  {
1967    A68_LONG_REAL a, b;
1968    POP_OBJECT (p, &b, A68_LONG_REAL);
1969    POP_OBJECT (p, &a, A68_LONG_REAL);
1970    errno = 0;
1971    PUSH_VALUE (p, dble (lgamma_double (VALUE (&a).f) + lgamma_double (VALUE (&b).f) - lgamma_double (VALUE (&a).f + VALUE (&b).f)), A68_LONG_REAL);
1972    MATH_RTE (p, errno != 0, M_LONG_REAL, NO_TEXT);
1973  }
1974  
1975  // LONG REAL infinity
1976  
1977  void genie_infinity_double (NODE_T * p)
1978  {
1979    PUSH_VALUE (p, dble (a68_posinf_real ()), A68_LONG_REAL);
1980  }
1981  
1982  // LONG REAL minus infinity
1983  
1984  void genie_minus_infinity_double (NODE_T * p)
1985  {
1986    PUSH_VALUE (p, dble (a68_neginf_double ()), A68_LONG_REAL);
1987  }
1988  
1989  #endif