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