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