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