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 MP_STATUS (z) = (MP_T) INIT_MASK;
495 }
496
497 //! @brief OP SHORTEN = (LONG LONG REAL) LONG REAL
498
499 void genie_shorten_mp_to_double (NODE_T * p)
500 {
501 MOID_T *mode = LHS_MODE (p);
502 int digs = DIGITS (mode), size = SIZE (mode);
503 DOUBLE_NUM_T d;
504 DECREMENT_STACK_POINTER (p, size);
505 MP_T *z = (MP_T *) STACK_TOP;
506 MP_STATUS (z) = (MP_T) INIT_MASK;
507 d.f = mp_to_double (p, z, digs);
508 PUSH_VALUE (p, d, A68G_LONG_REAL);
509 }
510
511 //! @brief OP SHORTEN = (LONG LONG COMPLEX) LONG COMPLEX
512
513 void genie_shorten_long_mp_complex_to_double_compl (NODE_T * p)
514 {
515 int digs = DIGITS (M_LONG_LONG_REAL), size = SIZE (M_LONG_LONG_REAL);
516 MP_T *b = (MP_T *) STACK_OFFSET (-size);
517 MP_T *a = (MP_T *) STACK_OFFSET (-2 * size);
518 DECREMENT_STACK_POINTER (p, 2 * size);
519 DOUBLE_NUM_T u, v;
520 u.f = mp_to_double (p, a, digs);
521 v.f = mp_to_double (p, b, digs);
522 PUSH_VALUE (p, u, A68G_LONG_REAL);
523 PUSH_VALUE (p, v, A68G_LONG_REAL);
524 }
525
526 //! @brief OP LENG = (LONG INT) LONG LONG INT
527
528 void genie_lengthen_double_int_to_mp (NODE_T * p)
529 {
530 int digs = DIGITS (M_LONG_LONG_INT);
531 A68G_LONG_INT k;
532 POP_OBJECT (p, &k, A68G_LONG_INT);
533 MP_T *z = nil_mp (p, digs);
534 (void) double_int_to_mp (p, z, VALUE (&k), digs);
535 MP_STATUS (z) = (MP_T) INIT_MASK;
536 }
537
538 //! @brief OP SHORTEN = (LONG LONG INT) LONG INT
539
540 void genie_shorten_mp_to_double_int (NODE_T * p)
541 {
542 MOID_T *mode = LHS_MODE (p);
543 int digs = DIGITS (mode), size = SIZE (mode);
544 DECREMENT_STACK_POINTER (p, size);
545 MP_T *z = (MP_T *) STACK_TOP;
546 MP_STATUS (z) = (MP_T) INIT_MASK;
547 PUSH_VALUE (p, mp_to_double_int (p, z, digs), A68G_LONG_INT);
548 }
549
550 //! @brief OP LENG = (INT) LONG INT
551
552 void genie_lengthen_int_to_double_int (NODE_T * p)
553 {
554 A68G_INT k;
555 POP_OBJECT (p, &k, A68G_INT);
556 INT_T v = VALUE (&k);
557 DOUBLE_NUM_T d;
558 if (v >= 0) {
559 LW (d) = v;
560 HW (d) = 0;
561 } else {
562 LW (d) = -v;
563 HW (d) = D_SIGN;
564 }
565 PUSH_VALUE (p, d, A68G_LONG_INT);
566 }
567
568 //! @brief OP SHORTEN = (LONG INT) INT
569
570 void genie_shorten_long_int_to_int (NODE_T * p)
571 {
572 A68G_LONG_INT k;
573 POP_OBJECT (p, &k, A68G_LONG_INT);
574 DOUBLE_NUM_T j = VALUE (&k);
575 PRELUDE_ERROR (HW (j) != 0 && HW (j) != D_SIGN, p, ERROR_MATH, M_INT);
576 PRELUDE_ERROR (LW (j) & D_SIGN, p, ERROR_MATH, M_INT);
577 if (D_NEG (j)) {
578 PUSH_VALUE (p, -LW (j), A68G_INT);
579 } else {
580 PUSH_VALUE (p, LW (j), A68G_INT);
581 }
582 }
583
584 // Constants.
585
586 //! @brief PROC long max int = LONG INT
587
588 void genie_double_max_int (NODE_T * p)
589 {
590 DOUBLE_NUM_T d;
591 HW (d) = 0x7fffffffffffffffLL;
592 LW (d) = 0xffffffffffffffffLL;
593 PUSH_VALUE (p, d, A68G_LONG_INT);
594 }
595
596 //! @brief PROC long max bits = LONG BITS
597
598 void genie_double_max_bits (NODE_T * p)
599 {
600 DOUBLE_NUM_T d;
601 HW (d) = 0xffffffffffffffffLL;
602 LW (d) = 0xffffffffffffffffLL;
603 PUSH_VALUE (p, d, A68G_LONG_INT);
604 }
605
606 //! @brief LONG REAL max long real
607
608 void genie_double_max_real (NODE_T * p)
609 {
610 DOUBLE_NUM_T d;
611 d.f = A68G_DOUBLE_MAX;
612 PUSH_VALUE (p, d, A68G_LONG_REAL);
613 }
614
615 //! @brief LONG REAL min long real
616
617 void genie_double_min_real (NODE_T * p)
618 {
619 DOUBLE_NUM_T d;
620 d.f = A68G_DOUBLE_MIN;
621 PUSH_VALUE (p, d, A68G_LONG_REAL);
622 }
623
624 //! @brief LONG REAL small long real
625
626 void genie_double_small_real (NODE_T * p)
627 {
628 DOUBLE_NUM_T d;
629 d.f = A68G_DOUBLE_EPS;
630 PUSH_VALUE (p, d, A68G_LONG_REAL);
631 }
632
633 //! @brief PROC long pi = LON REAL
634
635 void genie_pi_double (NODE_T * p)
636 {
637 DOUBLE_NUM_T w;
638 w.f = M_PIq;
639 PUSH_VALUE (p, w, A68G_LONG_INT);
640 }
641
642 // MONADs and DYADs
643
644 //! @brief OP SIGN = (LONG INT) INT
645
646 void genie_sign_double_int (NODE_T * p)
647 {
648 A68G_LONG_INT k;
649 POP_OBJECT (p, &k, A68G_LONG_INT);
650 PUSH_VALUE (p, sign_double_int (VALUE (&k)), A68G_INT);
651 }
652
653 //! @brief OP ABS = (LONG INT) LONG INT
654
655 void genie_abs_double_int (NODE_T * p)
656 {
657 A68G_LONG_INT *k;
658 POP_OPERAND_ADDRESS (p, k, A68G_LONG_INT);
659 VALUE (k) = abs_double_int (VALUE (k));
660 }
661
662 //! @brief OP ODD = (LONG INT) BOOL
663
664 void genie_odd_double_int (NODE_T * p)
665 {
666 A68G_LONG_INT j;
667 POP_OBJECT (p, &j, A68G_LONG_INT);
668 DOUBLE_NUM_T w = abs_double_int (VALUE (&j));
669 if (LW (w) & 0x1) {
670 PUSH_VALUE (p, A68G_TRUE, A68G_BOOL);
671 } else {
672 PUSH_VALUE (p, A68G_FALSE, A68G_BOOL);
673 }
674 }
675
676 //! @brief OP - = (LONG INT) LONG INT
677
678 void genie_minus_double_int (NODE_T * p)
679 {
680 A68G_LONG_INT *k;
681 POP_OPERAND_ADDRESS (p, k, A68G_LONG_INT);
682 VALUE (k) = neg_double_int (VALUE (k));
683 }
684
685 //! @brief OP + = (LONG INT, LONG INT) LONG INT
686
687 void genie_add_double_int (NODE_T * p)
688 {
689 A68G_LONG_INT i, j;
690 POP_OBJECT (p, &j, A68G_LONG_INT);
691 POP_OBJECT (p, &i, A68G_LONG_INT);
692 PUSH_VALUE (p, double_sadd (p, VALUE (&i), VALUE (&j)), A68G_LONG_INT);
693 }
694
695 //! @brief OP - = (LONG INT, LONG INT) LONG INT
696
697 void genie_sub_double_int (NODE_T * p)
698 {
699 A68G_LONG_INT i, j;
700 POP_OBJECT (p, &j, A68G_LONG_INT);
701 POP_OBJECT (p, &i, A68G_LONG_INT);
702 PUSH_VALUE (p, double_ssub (p, VALUE (&i), VALUE (&j)), A68G_LONG_INT);
703 }
704
705 //! @brief OP * = (LONG INT, LONG INT) LONG INT
706
707 void genie_mul_double_int (NODE_T * p)
708 {
709 A68G_LONG_INT i, j;
710 POP_OBJECT (p, &j, A68G_LONG_INT);
711 POP_OBJECT (p, &i, A68G_LONG_INT);
712 PUSH_VALUE (p, double_smul (p, VALUE (&i), VALUE (&j)), A68G_LONG_INT);
713 }
714
715 //! @brief OP / = (LONG INT, LONG INT) LONG INT
716
717 void genie_over_double_int (NODE_T * p)
718 {
719 A68G_LONG_INT i, j;
720 POP_OBJECT (p, &j, A68G_LONG_INT);
721 POP_OBJECT (p, &i, A68G_LONG_INT);
722 PRELUDE_ERROR (D_ZERO (VALUE (&j)), p, ERROR_DIVISION_BY_ZERO, M_LONG_INT);
723 PUSH_VALUE (p, double_sdiv (p, VALUE (&i), VALUE (&j), 0), A68G_LONG_INT);
724 }
725
726 //! @brief OP MOD = (LONG INT, LONG INT) LONG INT
727
728 void genie_mod_double_int (NODE_T * p)
729 {
730 A68G_LONG_INT i, j;
731 POP_OBJECT (p, &j, A68G_LONG_INT);
732 POP_OBJECT (p, &i, A68G_LONG_INT);
733 PRELUDE_ERROR (D_ZERO (VALUE (&j)), p, ERROR_DIVISION_BY_ZERO, M_LONG_INT);
734 PUSH_VALUE (p, double_sdiv (p, VALUE (&i), VALUE (&j), 1), A68G_LONG_INT);
735 }
736
737 //! @brief OP / = (LONG INT, LONG INT) LONG REAL
738
739 void genie_div_double_int (NODE_T * p)
740 {
741 A68G_LONG_INT i, j;
742 POP_OBJECT (p, &j, A68G_LONG_INT);
743 POP_OBJECT (p, &i, A68G_LONG_INT);
744 PRELUDE_ERROR (D_ZERO (VALUE (&j)), p, ERROR_DIVISION_BY_ZERO, M_LONG_INT);
745 DOUBLE_NUM_T u, v, w;
746 v = double_int_to_double (p, VALUE (&j));
747 u = double_int_to_double (p, VALUE (&i));
748 w.f = u.f / v.f;
749 PUSH_VALUE (p, w, A68G_LONG_REAL);
750 }
751
752 //! @brief OP ** = (LONG INT, INT) INT
753
754 void genie_pow_double_int_int (NODE_T * p)
755 {
756 A68G_LONG_INT i; A68G_INT j;
757 POP_OBJECT (p, &j, A68G_INT);
758 PRELUDE_ERROR (VALUE (&j) < 0, p, ERROR_EXPONENT_INVALID, M_INT);
759 POP_OBJECT (p, &i, A68G_LONG_INT);
760 DOUBLE_NUM_T mult = VALUE (&i), prod;
761 set_lw (prod, 1);
762 UNSIGNED_T top = (UNSIGNED_T) VALUE (&j), expo = 1;
763 while (expo <= top) {
764 if (expo & top) {
765 prod = double_smul (p, prod, mult);
766 }
767 expo <<= 1;
768 if (expo <= top) {
769 mult = double_smul (p, mult, mult);
770 }
771 }
772 PUSH_VALUE (p, prod, A68G_LONG_INT);
773 }
774
775 //! @brief OP - = (LONG REAL) LONG REAL
776
777 void genie_minus_double (NODE_T * p)
778 {
779 A68G_LONG_REAL *u;
780 POP_OPERAND_ADDRESS (p, u, A68G_LONG_REAL);
781 VALUE (u).f = -(VALUE (u).f);
782 }
783
784 //! @brief OP ABS = (LONG REAL) LONG REAL
785
786 void genie_abs_double (NODE_T * p)
787 {
788 A68G_LONG_REAL *u;
789 POP_OPERAND_ADDRESS (p, u, A68G_LONG_REAL);
790 VALUE (u).f = fabs_double (VALUE (u).f);
791 }
792
793 //! @brief OP SIGN = (LONG REAL) INT
794
795 void genie_sign_double (NODE_T * p)
796 {
797 A68G_LONG_REAL u;
798 POP_OBJECT (p, &u, A68G_LONG_REAL);
799 PUSH_VALUE (p, sign_double (VALUE (&u)), A68G_INT);
800 }
801
802 //! @brief OP ** = (LONG REAL, INT) INT
803
804 void genie_pow_double_int (NODE_T * p)
805 {
806 A68G_INT j;
807 POP_OBJECT (p, &j, A68G_INT);
808 INT_T top = (INT_T) VALUE (&j);
809 A68G_LONG_REAL z;
810 POP_OBJECT (p, &z, A68G_LONG_INT);
811 DOUBLE_NUM_T mult, prod;
812 prod.f = 1.0q;
813 mult.f = VALUE (&z).f;
814 int negative;
815 if (top < 0) {
816 top = -top;
817 negative = A68G_TRUE;
818 } else {
819 negative = A68G_FALSE;
820 }
821 UNSIGNED_T expo = 1;
822 while (expo <= top) {
823 if (expo & top) {
824 prod.f = prod.f * mult.f;
825 CHECK_DOUBLE_REAL (p, prod.f);
826 }
827 expo <<= 1;
828 if (expo <= top) {
829 mult.f = mult.f * mult.f;
830 CHECK_DOUBLE_REAL (p, mult.f);
831 }
832 }
833 if (negative) {
834 prod.f = 1.0q / prod.f;
835 }
836 PUSH_VALUE (p, prod, A68G_LONG_REAL);
837 }
838
839 //! @brief OP ** = (LONG REAL, LONG REAL) LONG REAL
840
841 void genie_pow_double (NODE_T * p)
842 {
843 A68G_LONG_REAL x, y;
844 POP_OBJECT (p, &y, A68G_LONG_REAL);
845 POP_OBJECT (p, &x, A68G_LONG_REAL);
846 errno = 0;
847 PRELUDE_ERROR (VALUE (&x).f < 0.0q, p, ERROR_INVALID_ARGUMENT, M_LONG_REAL);
848 DOUBLE_T z = 0.0q;
849 if (VALUE (&x).f == 0.0q) {
850 if (VALUE (&y).f < 0.0q) {
851 errno = ERANGE;
852 MATH_RTE (p, errno != 0, M_LONG_REAL, NO_TEXT);
853 } else {
854 z = (VALUE (&y).f == 0.0q ? 1.0q : 0.0q);
855 }
856 } else {
857 z = exp_double (VALUE (&y).f * log_double (VALUE (&x).f));
858 MATH_RTE (p, errno != 0, M_LONG_REAL, NO_TEXT);
859 }
860 PUSH_VALUE (p, dble (z), A68G_LONG_REAL);
861 }
862
863 //! @brief OP + = (LONG REAL, LONG REAL) LONG REAL
864
865 void genie_add_double (NODE_T * p)
866 {
867 A68G_LONG_REAL u, v;
868 POP_OBJECT (p, &v, A68G_LONG_REAL);
869 POP_OBJECT (p, &u, A68G_LONG_REAL);
870 DOUBLE_NUM_T w;
871 w.f = VALUE (&u).f + VALUE (&v).f;
872 CHECK_DOUBLE_REAL (p, w.f);
873 PUSH_VALUE (p, w, A68G_LONG_REAL);
874 }
875
876 //! @brief OP - = (LONG REAL, LONG REAL) LONG REAL
877
878 void genie_sub_double (NODE_T * p)
879 {
880 A68G_LONG_REAL u, v;
881 POP_OBJECT (p, &v, A68G_LONG_REAL);
882 POP_OBJECT (p, &u, A68G_LONG_REAL);
883 DOUBLE_NUM_T w;
884 w.f = VALUE (&u).f - VALUE (&v).f;
885 CHECK_DOUBLE_REAL (p, w.f);
886 PUSH_VALUE (p, w, A68G_LONG_REAL);
887 }
888
889 //! @brief OP * = (LONG REAL, LONG REAL) LONG REAL
890
891 void genie_mul_double (NODE_T * p)
892 {
893 A68G_LONG_REAL u, v;
894 POP_OBJECT (p, &v, A68G_LONG_REAL);
895 POP_OBJECT (p, &u, A68G_LONG_REAL);
896 DOUBLE_NUM_T w;
897 w.f = VALUE (&u).f * VALUE (&v).f;
898 CHECK_DOUBLE_REAL (p, w.f);
899 PUSH_VALUE (p, w, A68G_LONG_REAL);
900 }
901
902 //! @brief OP / = (LONG REAL, LONG REAL) LONG REAL
903
904 void genie_over_double (NODE_T * p)
905 {
906 A68G_LONG_REAL u, v;
907 POP_OBJECT (p, &v, A68G_LONG_REAL);
908 POP_OBJECT (p, &u, A68G_LONG_REAL);
909 PRELUDE_ERROR (VALUE (&v).f == 0.0q, p, ERROR_DIVISION_BY_ZERO, M_LONG_REAL);
910 DOUBLE_NUM_T w;
911 w.f = VALUE (&u).f / VALUE (&v).f;
912 PUSH_VALUE (p, w, A68G_LONG_REAL);
913 }
914
915 //! @brief OP +:= = (REF LONG INT, LONG INT) REF LONG INT
916
917 void genie_plusab_double_int (NODE_T * p)
918 {
919 genie_f_and_becomes (p, M_REF_LONG_INT, genie_add_double_int);
920 }
921
922 //! @brief OP -:= = (REF LONG INT, LONG INT) REF LONG INT
923
924 void genie_minusab_double_int (NODE_T * p)
925 {
926 genie_f_and_becomes (p, M_REF_LONG_INT, genie_sub_double_int);
927 }
928
929 //! @brief OP *:= = (REF LONG INT, LONG INT) REF LONG INT
930
931 void genie_timesab_double_int (NODE_T * p)
932 {
933 genie_f_and_becomes (p, M_REF_LONG_INT, genie_mul_double_int);
934 }
935
936 //! @brief OP %:= = (REF LONG INT, LONG INT) REF LONG INT
937
938 void genie_overab_double_int (NODE_T * p)
939 {
940 genie_f_and_becomes (p, M_REF_LONG_INT, genie_over_double_int);
941 }
942
943 //! @brief OP %*:= = (REF LONG INT, LONG INT) REF LONG INT
944
945 void genie_modab_double_int (NODE_T * p)
946 {
947 genie_f_and_becomes (p, M_REF_LONG_INT, genie_mod_double_int);
948 }
949
950 //! @brief OP +:= = (REF LONG REAL, LONG REAL) REF LONG REAL
951
952 void genie_plusab_double (NODE_T * p)
953 {
954 genie_f_and_becomes (p, M_REF_LONG_REAL, genie_add_double);
955 }
956
957 //! @brief OP -:= = (REF LONG REAL, LONG REAL) REF LONG REAL
958
959 void genie_minusab_double (NODE_T * p)
960 {
961 genie_f_and_becomes (p, M_REF_LONG_REAL, genie_sub_double);
962 }
963
964 //! @brief OP *:= = (REF LONG REAL, LONG REAL) REF LONG REAL
965
966 void genie_timesab_double (NODE_T * p)
967 {
968 genie_f_and_becomes (p, M_REF_LONG_REAL, genie_mul_double);
969 }
970
971 //! @brief OP /:= = (REF LONG REAL, LONG REAL) REF LONG REAL
972
973 void genie_divab_double (NODE_T * p)
974 {
975 genie_f_and_becomes (p, M_REF_LONG_REAL, genie_over_double);
976 }
977
978 // OP (LONG INT, LONG INT) BOOL.
979
980 #define A68G_CMP_INT(n, OP)\
981 void n (NODE_T * p) {\
982 A68G_LONG_INT i, j;\
983 POP_OBJECT (p, &j, A68G_LONG_INT);\
984 POP_OBJECT (p, &i, A68G_LONG_INT);\
985 DOUBLE_NUM_T w = double_ssub (p, VALUE (&i), VALUE (&j));\
986 int k = sign_double_int (w);\
987 PUSH_VALUE (p, (BOOL_T) (k OP 0), A68G_BOOL);\
988 }
989
990 A68G_CMP_INT (genie_eq_double_int, ==);
991 A68G_CMP_INT (genie_ne_double_int, !=);
992 A68G_CMP_INT (genie_lt_double_int, <);
993 A68G_CMP_INT (genie_gt_double_int, >);
994 A68G_CMP_INT (genie_le_double_int, <=);
995 A68G_CMP_INT (genie_ge_double_int, >=);
996
997 // OP (LONG REAL, LONG REAL) BOOL.
998 #define A68G_CMP_REAL(n, OP)\
999 void n (NODE_T * p) {\
1000 A68G_LONG_REAL i, j;\
1001 POP_OBJECT (p, &j, A68G_LONG_REAL);\
1002 POP_OBJECT (p, &i, A68G_LONG_REAL);\
1003 PUSH_VALUE (p, (BOOL_T) (VALUE (&i).f OP VALUE (&j).f), A68G_BOOL);\
1004 }
1005
1006 A68G_CMP_REAL (genie_eq_double, ==);
1007 A68G_CMP_REAL (genie_ne_double, !=);
1008 A68G_CMP_REAL (genie_lt_double, <);
1009 A68G_CMP_REAL (genie_gt_double, >);
1010 A68G_CMP_REAL (genie_le_double, <=);
1011 A68G_CMP_REAL (genie_ge_double, >=);
1012
1013 //! @brief OP NOT = (LONG BITS) LONG BITS
1014
1015 void genie_not_double_bits (NODE_T * p)
1016 {
1017 A68G_LONG_BITS i;
1018 POP_OBJECT (p, &i, A68G_LONG_BITS);
1019 DOUBLE_NUM_T w;
1020 HW (w) = ~HW (VALUE (&i));
1021 LW (w) = ~LW (VALUE (&i));
1022 PUSH_VALUE (p, w, A68G_LONG_BITS);
1023 }
1024
1025 //! @brief OP = = (LONG BITS, LONG BITS) BOOL.
1026
1027 void genie_eq_double_bits (NODE_T * p)
1028 {
1029 A68G_LONG_BITS i, j;
1030 POP_OBJECT (p, &j, A68G_LONG_BITS);
1031 POP_OBJECT (p, &i, A68G_LONG_BITS);
1032 BOOL_T u = HW (VALUE (&i)) == HW (VALUE (&j));
1033 BOOL_T v = LW (VALUE (&i)) == LW (VALUE (&j));
1034 PUSH_VALUE (p, (BOOL_T) (u & v ? A68G_TRUE : A68G_FALSE), A68G_BOOL);
1035 }
1036
1037 //! @brief OP ~= = (LONG BITS, LONG BITS) BOOL.
1038
1039 void genie_ne_double_bits (NODE_T * p)
1040 {
1041 A68G_LONG_BITS i, j;
1042 POP_OBJECT (p, &j, A68G_LONG_BITS); // (i ~= j) == ~ (i = j)
1043 POP_OBJECT (p, &i, A68G_LONG_BITS);
1044 BOOL_T u = HW (VALUE (&i)) == HW (VALUE (&j));
1045 BOOL_T v = LW (VALUE (&i)) == LW (VALUE (&j));
1046 PUSH_VALUE (p, (BOOL_T) (u & v ? A68G_FALSE : A68G_TRUE), A68G_BOOL);
1047 }
1048
1049 //! @brief OP <= = (LONG BITS, LONG BITS) BOOL
1050
1051 void genie_le_double_bits (NODE_T * p)
1052 {
1053 A68G_LONG_BITS i, j;
1054 POP_OBJECT (p, &j, A68G_LONG_BITS);
1055 POP_OBJECT (p, &i, A68G_LONG_BITS);
1056 BOOL_T u = (HW (VALUE (&i)) | HW (VALUE (&j))) == HW (VALUE (&j));
1057 BOOL_T v = (LW (VALUE (&i)) | LW (VALUE (&j))) == LW (VALUE (&j));
1058 PUSH_VALUE (p, (BOOL_T) (u & v ? A68G_TRUE : A68G_FALSE), A68G_BOOL);
1059 }
1060
1061 //! @brief OP > = (LONG BITS, LONG BITS) BOOL
1062
1063 void genie_gt_double_bits (NODE_T * p)
1064 {
1065 A68G_LONG_BITS i, j;
1066 POP_OBJECT (p, &j, A68G_LONG_BITS); // (i > j) == ! (i <= j)
1067 POP_OBJECT (p, &i, A68G_LONG_BITS);
1068 BOOL_T u = (HW (VALUE (&i)) | HW (VALUE (&j))) == HW (VALUE (&j));
1069 BOOL_T v = (LW (VALUE (&i)) | LW (VALUE (&j))) == LW (VALUE (&j));
1070 PUSH_VALUE (p, (BOOL_T) (u & v ? A68G_FALSE : A68G_TRUE), A68G_BOOL);
1071 }
1072
1073 //! @brief OP >= = (LONG BITS, LONG BITS) BOOL
1074
1075 void genie_ge_double_bits (NODE_T * p)
1076 {
1077 A68G_LONG_BITS i, j;
1078 POP_OBJECT (p, &j, A68G_LONG_BITS); // (i >= j) == (j <= i)
1079 POP_OBJECT (p, &i, A68G_LONG_BITS);
1080 BOOL_T u = (HW (VALUE (&i)) | HW (VALUE (&j))) == HW (VALUE (&i));
1081 BOOL_T v = (LW (VALUE (&i)) | LW (VALUE (&j))) == LW (VALUE (&i));
1082 PUSH_VALUE (p, (BOOL_T) (u & v ? A68G_TRUE : A68G_FALSE), A68G_BOOL);
1083 }
1084
1085 //! @brief OP < = (LONG BITS, LONG BITS) BOOL
1086
1087 void genie_lt_double_bits (NODE_T * p)
1088 {
1089 A68G_LONG_BITS i, j;
1090 POP_OBJECT (p, &j, A68G_LONG_BITS); // (i < j) == ! (i >= j)
1091 POP_OBJECT (p, &i, A68G_LONG_BITS);
1092 BOOL_T u = (HW (VALUE (&i)) | HW (VALUE (&j))) == HW (VALUE (&i));
1093 BOOL_T v = (LW (VALUE (&i)) | LW (VALUE (&j))) == LW (VALUE (&i));
1094 PUSH_VALUE (p, (BOOL_T) (u & v ? A68G_FALSE : A68G_TRUE), A68G_BOOL);
1095 }
1096
1097 //! @brief PROC long bits pack = ([] BOOL) BITS
1098
1099 void genie_double_bits_pack (NODE_T * p)
1100 {
1101 A68G_REF z;
1102 POP_REF (p, &z);
1103 CHECK_REF (p, z, M_ROW_BOOL);
1104 A68G_ARRAY *arr; A68G_TUPLE *tup;
1105 GET_DESCRIPTOR (arr, tup, &z);
1106 size_t size = ROW_SIZE (tup);
1107 PRELUDE_ERROR (size < 0 || size > A68G_LONG_BITS_WIDTH, p, ERROR_OUT_OF_BOUNDS, M_ROW_BOOL);
1108 DOUBLE_NUM_T w;
1109 set_lw (w, 0x0);
1110 if (ROW_SIZE (tup) > 0) {
1111 UNSIGNED_T bit = 0x0;
1112 BYTE_T *base = DEREF (BYTE_T, &ARRAY (arr));
1113 int n = 0;
1114 for (int k = UPB (tup); k >= LWB (tup); k--) {
1115 A68G_BOOL *boo = (A68G_BOOL *) & (base[INDEX_1_DIM (arr, tup, k)]);
1116 CHECK_INIT (p, INITIALISED (boo), M_BOOL);
1117 if (n == 0 || n == A68G_BITS_WIDTH) {
1118 bit = 0x1;
1119 }
1120 if (VALUE (boo)) {
1121 if (n < A68G_BITS_WIDTH) {
1122 LW (w) |= bit;
1123 } else {
1124 HW (w) |= bit;
1125 };
1126 }
1127 n++;
1128 bit <<= 1;
1129 }
1130 }
1131 PUSH_VALUE (p, w, A68G_LONG_BITS);
1132 }
1133
1134 //! @brief OP AND = (LONG BITS, LONG BITS) LONG BITS
1135
1136 void genie_and_double_bits (NODE_T * p)
1137 {
1138 A68G_LONG_BITS i, j;
1139 POP_OBJECT (p, &j, A68G_LONG_BITS);
1140 POP_OBJECT (p, &i, A68G_LONG_BITS);
1141 DOUBLE_NUM_T w;
1142 HW (w) = HW (VALUE (&i)) & HW (VALUE (&j));
1143 LW (w) = LW (VALUE (&i)) & LW (VALUE (&j));
1144 PUSH_VALUE (p, w, A68G_LONG_BITS);
1145 }
1146
1147 //! @brief OP OR = (LONG BITS, LONG BITS) LONG BITS
1148
1149 void genie_or_double_bits (NODE_T * p)
1150 {
1151 A68G_LONG_BITS i, j;
1152 POP_OBJECT (p, &j, A68G_LONG_BITS);
1153 POP_OBJECT (p, &i, A68G_LONG_BITS);
1154 DOUBLE_NUM_T w;
1155 HW (w) = HW (VALUE (&i)) | HW (VALUE (&j));
1156 LW (w) = LW (VALUE (&i)) | LW (VALUE (&j));
1157 PUSH_VALUE (p, w, A68G_LONG_BITS);
1158 }
1159
1160 //! @brief OP XOR = (LONG BITS, LONG BITS) LONG BITS
1161
1162 void genie_xor_double_bits (NODE_T * p)
1163 {
1164 A68G_LONG_BITS i, j;
1165 POP_OBJECT (p, &j, A68G_LONG_BITS);
1166 POP_OBJECT (p, &i, A68G_LONG_BITS);
1167 DOUBLE_NUM_T w;
1168 HW (w) = HW (VALUE (&i)) ^ HW (VALUE (&j));
1169 LW (w) = LW (VALUE (&i)) ^ LW (VALUE (&j));
1170 PUSH_VALUE (p, w, A68G_LONG_BITS);
1171 }
1172
1173 //! @brief OP + = (LONG BITS, LONG BITS) LONG BITS
1174
1175 void genie_add_double_bits (NODE_T * p)
1176 {
1177 A68G_LONG_BITS i, j;
1178 POP_OBJECT (p, &j, A68G_LONG_BITS);
1179 POP_OBJECT (p, &i, A68G_LONG_BITS);
1180 DOUBLE_NUM_T w;
1181 add_double (p, M_LONG_BITS, w, VALUE (&i), VALUE (&j));
1182 PUSH_VALUE (p, w, A68G_LONG_BITS);
1183 }
1184
1185 //! @brief OP - = (LONG BITS, LONG BITS) LONG BITS
1186
1187 void genie_sub_double_bits (NODE_T * p)
1188 {
1189 A68G_LONG_BITS i, j;
1190 POP_OBJECT (p, &j, A68G_LONG_BITS);
1191 POP_OBJECT (p, &i, A68G_LONG_BITS);
1192 DOUBLE_NUM_T w;
1193 sub_double (p, M_LONG_BITS, w, VALUE (&i), VALUE (&j));
1194 PUSH_VALUE (p, w, A68G_LONG_BITS);
1195 }
1196
1197 //! @brief OP * = (LONG BITS, LONG BITS) LONG BITS
1198
1199 void genie_times_double_bits (NODE_T * p)
1200 {
1201 A68G_LONG_BITS i, j;
1202 POP_OBJECT (p, &j, A68G_LONG_BITS);
1203 POP_OBJECT (p, &i, A68G_LONG_BITS);
1204 DOUBLE_NUM_T w = double_umul (p, M_LONG_BITS, VALUE (&i), VALUE (&j));
1205 PUSH_VALUE (p, w, A68G_LONG_BITS);
1206 }
1207
1208 //! @brief OP OVER = (LONG BITS, LONG BITS) LONG BITS
1209
1210 void genie_over_double_bits (NODE_T * p)
1211 {
1212 A68G_LONG_BITS i, j;
1213 POP_OBJECT (p, &j, A68G_LONG_BITS);
1214 POP_OBJECT (p, &i, A68G_LONG_BITS);
1215 DOUBLE_NUM_T w = double_udiv (p, M_LONG_BITS, VALUE (&i), VALUE (&j), 0);
1216 PUSH_VALUE (p, w, A68G_LONG_BITS);
1217 }
1218
1219 //! @brief OP MOD = (LONG BITS, LONG BITS) LONG BITS
1220
1221 void genie_mod_double_bits (NODE_T * p)
1222 {
1223 A68G_LONG_BITS i, j;
1224 DOUBLE_NUM_T w;
1225 POP_OBJECT (p, &j, A68G_LONG_BITS);
1226 POP_OBJECT (p, &i, A68G_LONG_BITS);
1227 w = double_udiv (p, M_LONG_BITS, VALUE (&i), VALUE (&j), 1);
1228 PUSH_VALUE (p, w, A68G_LONG_BITS);
1229 }
1230
1231 //! @brief OP ELEM = (INT, LONG BITS) BOOL
1232
1233 void genie_elem_double_bits (NODE_T * p)
1234 {
1235 A68G_LONG_BITS j; A68G_INT i;
1236 POP_OBJECT (p, &j, A68G_LONG_BITS);
1237 POP_OBJECT (p, &i, A68G_INT);
1238 int k = VALUE (&i);
1239 PRELUDE_ERROR (k < 1 || k > A68G_LONG_BITS_WIDTH, p, ERROR_OUT_OF_BOUNDS, M_INT);
1240 UNSIGNED_T mask = 0x1, *w;
1241 if (k <= A68G_BITS_WIDTH) {
1242 w = &(HW (VALUE (&j)));
1243 } else {
1244 w = &(LW (VALUE (&j)));
1245 k -= A68G_BITS_WIDTH;
1246 }
1247 for (int n = 0; n < A68G_BITS_WIDTH - k; n++) {
1248 mask = mask << 1;
1249 }
1250 PUSH_VALUE (p, (BOOL_T) ((*w & mask) ? A68G_TRUE : A68G_FALSE), A68G_BOOL);
1251 }
1252
1253 //! @brief OP SET = (INT, LONG BITS) LONG BITS
1254
1255 void genie_set_double_bits (NODE_T * p)
1256 {
1257 A68G_LONG_BITS j; A68G_INT i;
1258 POP_OBJECT (p, &j, A68G_LONG_BITS);
1259 POP_OBJECT (p, &i, A68G_INT);
1260 int k = VALUE (&i);
1261 PRELUDE_ERROR (k < 1 || k > A68G_LONG_BITS_WIDTH, p, ERROR_OUT_OF_BOUNDS, M_INT);
1262 UNSIGNED_T mask = 0x1, *w;
1263 if (k <= A68G_BITS_WIDTH) {
1264 w = &(HW (VALUE (&j)));
1265 } else {
1266 w = &(LW (VALUE (&j)));
1267 k -= A68G_BITS_WIDTH;
1268 }
1269 for (int n = 0; n < A68G_BITS_WIDTH - k; n++) {
1270 mask = mask << 1;
1271 }
1272 (*w) |= mask;
1273 PUSH_OBJECT (p, j, A68G_LONG_BITS);
1274 }
1275
1276 //! @brief OP CLEAR = (INT, LONG BITS) LONG BITS
1277
1278 void genie_clear_double_bits (NODE_T * p)
1279 {
1280 A68G_LONG_BITS j; A68G_INT i;
1281 POP_OBJECT (p, &j, A68G_LONG_BITS);
1282 POP_OBJECT (p, &i, A68G_INT);
1283 int k = VALUE (&i);
1284 PRELUDE_ERROR (k < 1 || k > A68G_LONG_BITS_WIDTH, p, ERROR_OUT_OF_BOUNDS, M_INT);
1285 UNSIGNED_T mask = 0x1, *w;
1286 if (k <= A68G_BITS_WIDTH) {
1287 w = &(HW (VALUE (&j)));
1288 } else {
1289 w = &(LW (VALUE (&j)));
1290 k -= A68G_BITS_WIDTH;
1291 }
1292 for (int n = 0; n < A68G_BITS_WIDTH - k; n++) {
1293 mask = mask << 1;
1294 }
1295 (*w) &= ~mask;
1296 PUSH_OBJECT (p, j, A68G_LONG_BITS);
1297 }
1298
1299 //! @brief OP SHL = (LONG BITS, INT) LONG BITS
1300
1301 void genie_shl_double_bits (NODE_T * p)
1302 {
1303 A68G_LONG_BITS i; A68G_INT j;
1304 POP_OBJECT (p, &j, A68G_INT);
1305 POP_OBJECT (p, &i, A68G_LONG_BITS);
1306 DOUBLE_NUM_T *w = &VALUE (&i);
1307 int k = VALUE (&j);
1308 if (VALUE (&j) >= 0) {
1309 for (int n = 0; n < k; n++) {
1310 UNSIGNED_T carry = ((LW (*w) & D_SIGN) ? 0x1 : 0x0);
1311 PRELUDE_ERROR (MODCHK (p, M_LONG_BITS, HW (*w) & D_SIGN), p, ERROR_MATH, M_LONG_BITS);
1312 HW (*w) = (HW (*w) << 1) | carry;
1313 LW (*w) = (LW (*w) << 1);
1314 }
1315 } else {
1316 k = -k;
1317 for (int n = 0; n < k; n++) {
1318 UNSIGNED_T carry = ((HW (*w) & 0x1) ? D_SIGN : 0x0);
1319 HW (*w) = (HW (*w) >> 1);
1320 LW (*w) = (LW (*w) >> 1) | carry;
1321 }
1322 }
1323 PUSH_OBJECT (p, i, A68G_LONG_BITS);
1324 }
1325
1326 //! @brief OP SHR = (LONG BITS, INT) LONG BITS
1327
1328 void genie_shr_double_bits (NODE_T * p)
1329 {
1330 A68G_INT *j;
1331 POP_OPERAND_ADDRESS (p, j, A68G_INT);
1332 VALUE (j) = -VALUE (j);
1333 genie_shl_double_bits (p); // Conform RR
1334 }
1335
1336 //! @brief OP ROL = (LONG BITS, INT) LONG BITS
1337
1338 void genie_rol_double_bits (NODE_T * p)
1339 {
1340 A68G_LONG_BITS i; A68G_INT j;
1341 DOUBLE_NUM_T *w = &VALUE (&i);
1342 POP_OBJECT (p, &j, A68G_INT);
1343 POP_OBJECT (p, &i, A68G_LONG_BITS);
1344 int k = VALUE (&j);
1345 if (k >= 0) {
1346 for (int n = 0; n < k; n++) {
1347 UNSIGNED_T carry = ((HW (*w) & D_SIGN) ? 0x1 : 0x0);
1348 UNSIGNED_T carry_between = ((LW (*w) & D_SIGN) ? 0x1 : 0x0);
1349 HW (*w) = (HW (*w) << 1) | carry_between;
1350 LW (*w) = (LW (*w) << 1) | carry;
1351 }
1352 } else {
1353 k = -k;
1354 for (int n = 0; n < k; n++) {
1355 UNSIGNED_T carry = ((LW (*w) & 0x1) ? D_SIGN : 0x0);
1356 UNSIGNED_T carry_between = ((HW (*w) & 0x1) ? D_SIGN : 0x0);
1357 HW (*w) = (HW (*w) >> 1) | carry;
1358 LW (*w) = (LW (*w) >> 1) | carry_between;
1359 }
1360 }
1361 PUSH_OBJECT (p, i, A68G_LONG_BITS);
1362 }
1363
1364 //! @brief OP ROR = (LONG BITS, INT) LONG BITS
1365
1366 void genie_ror_double_bits (NODE_T * p)
1367 {
1368 A68G_INT *j;
1369 POP_OPERAND_ADDRESS (p, j, A68G_INT);
1370 VALUE (j) = -VALUE (j);
1371 genie_rol_double_bits (p); // Conform RR
1372 }
1373
1374 //! @brief OP BIN = (LONG INT) LONG BITS
1375
1376 void genie_bin_double_int (NODE_T * p)
1377 {
1378 A68G_LONG_INT i;
1379 POP_OBJECT (p, &i, A68G_LONG_INT);
1380 // RR does not convert negative numbers
1381 if (D_NEG (VALUE (&i))) {
1382 errno = EDOM;
1383 diagnostic (A68G_RUNTIME_ERROR, p, ERROR_OUT_OF_BOUNDS, M_BITS);
1384 exit_genie (p, A68G_RUNTIME_ERROR);
1385 }
1386 PUSH_OBJECT (p, i, A68G_LONG_BITS);
1387 }
1388
1389 //! @brief OP +* = (LONG REAL, LONG REAL) LONG COMPLEX
1390
1391 void genie_i_double_compl (NODE_T * p)
1392 {
1393 (void) p;
1394 }
1395
1396 //! @brief OP SHORTEN = (LONG COMPLEX) COMPLEX
1397
1398 void genie_shorten_double_compl_to_complex (NODE_T * p)
1399 {
1400 A68G_LONG_REAL re, im;
1401 POP_OBJECT (p, &im, A68G_LONG_REAL);
1402 POP_OBJECT (p, &re, A68G_LONG_REAL);
1403 REAL_T w = VALUE (&re).f;
1404 PUSH_VALUE (p, w, A68G_REAL);
1405 w = VALUE (&im).f;
1406 PUSH_VALUE (p, w, A68G_REAL);
1407 }
1408
1409 //! @brief OP LENG = (LONG COMPLEX) LONG LONG COMPLEX
1410
1411 void genie_lengthen_double_compl_to_long_mp_complex (NODE_T * p)
1412 {
1413 int digs = DIGITS (M_LONG_LONG_REAL);
1414 A68G_LONG_REAL re, im;
1415 POP_OBJECT (p, &im, A68G_LONG_REAL);
1416 POP_OBJECT (p, &re, A68G_LONG_REAL);
1417 MP_T *z = nil_mp (p, digs);
1418 (void) double_to_mp (p, z, VALUE (&re).f, A68G_DOUBLE_DIG, A68G_FALSE, digs);
1419 MP_STATUS (z) = (MP_T) INIT_MASK;
1420 z = nil_mp (p, digs);
1421 (void) double_to_mp (p, z, VALUE (&im).f, A68G_DOUBLE_DIG, A68G_FALSE, digs);
1422 MP_STATUS (z) = (MP_T) INIT_MASK;
1423 }
1424
1425 //! @brief OP +* = (LONG INT, LONG INT) LONG COMPLEX
1426
1427 void genie_i_int_double_compl (NODE_T * p)
1428 {
1429 A68G_LONG_INT re, im;
1430 POP_OBJECT (p, &im, A68G_LONG_INT);
1431 POP_OBJECT (p, &re, A68G_LONG_INT);
1432 PUSH_VALUE (p, double_int_to_double (p, VALUE (&re)), A68G_LONG_REAL);
1433 PUSH_VALUE (p, double_int_to_double (p, VALUE (&im)), A68G_LONG_REAL);
1434 }
1435
1436 //! @brief OP RE = (LONG COMPLEX) LONG REAL
1437
1438 void genie_re_double_compl (NODE_T * p)
1439 {
1440 DECREMENT_STACK_POINTER (p, SIZE (M_LONG_REAL));
1441 }
1442
1443 //! @brief OP IM = (LONG COMPLEX) LONG REAL
1444
1445 void genie_im_double_compl (NODE_T * p)
1446 {
1447 A68G_LONG_REAL re, im;
1448 POP_OBJECT (p, &im, A68G_LONG_REAL);
1449 POP_OBJECT (p, &re, A68G_LONG_REAL);
1450 PUSH_OBJECT (p, im, A68G_LONG_REAL);
1451 }
1452
1453 //! @brief OP - = (LONG COMPLEX) LONG COMPLEX
1454
1455 void genie_minus_double_compl (NODE_T * p)
1456 {
1457 A68G_LONG_REAL re, im;
1458 POP_OBJECT (p, &im, A68G_LONG_REAL);
1459 POP_OBJECT (p, &re, A68G_LONG_REAL);
1460 VALUE (&re).f = -VALUE (&re).f;
1461 VALUE (&im).f = -VALUE (&im).f;
1462 PUSH_OBJECT (p, im, A68G_LONG_REAL);
1463 PUSH_OBJECT (p, re, A68G_LONG_REAL);
1464 }
1465
1466 //! @brief OP ABS = (LONG COMPLEX) LONG REAL
1467
1468 void genie_abs_double_compl (NODE_T * p)
1469 {
1470 A68G_LONG_REAL re, im;
1471 POP_LONG_COMPLEX (p, &re, &im);
1472 PUSH_VALUE (p, dble (a68g_hypot_double (VALUE (&re).f, VALUE (&im).f)), A68G_LONG_REAL);
1473 }
1474
1475 //! @brief OP ARG = (LONG COMPLEX) LONG REAL
1476
1477 void genie_arg_double_compl (NODE_T * p)
1478 {
1479 A68G_LONG_REAL re, im;
1480 POP_LONG_COMPLEX (p, &re, &im);
1481 PRELUDE_ERROR (VALUE (&re).f == 0.0q && VALUE (&im).f == 0.0q, p, ERROR_INVALID_ARGUMENT, M_LONG_COMPLEX);
1482 PUSH_VALUE (p, dble (atan2_double (VALUE (&im).f, VALUE (&re).f)), A68G_LONG_REAL);
1483 }
1484
1485 //! @brief OP CONJ = (LONG COMPLEX) LONG COMPLEX
1486
1487 void genie_conj_double_compl (NODE_T * p)
1488 {
1489 A68G_LONG_REAL im;
1490 POP_OBJECT (p, &im, A68G_LONG_REAL);
1491 VALUE (&im).f = -VALUE (&im).f;
1492 PUSH_OBJECT (p, im, A68G_LONG_REAL);
1493 }
1494
1495 //! @brief OP + = (COMPLEX, COMPLEX) COMPLEX
1496
1497 void genie_add_double_compl (NODE_T * p)
1498 {
1499 A68G_LONG_REAL re_x, im_x, re_y, im_y;
1500 POP_LONG_COMPLEX (p, &re_y, &im_y);
1501 POP_LONG_COMPLEX (p, &re_x, &im_x);
1502 VALUE (&re_x).f += VALUE (&re_y).f;
1503 VALUE (&im_x).f += VALUE (&im_y).f;
1504 CHECK_DOUBLE_COMPLEX (p, VALUE (&im_x).f, VALUE (&im_y).f);
1505 PUSH_OBJECT (p, re_x, A68G_LONG_REAL);
1506 PUSH_OBJECT (p, im_x, A68G_LONG_REAL);
1507 }
1508
1509 //! @brief OP - = (COMPLEX, COMPLEX) COMPLEX
1510
1511 void genie_sub_double_compl (NODE_T * p)
1512 {
1513 A68G_LONG_REAL re_x, im_x, re_y, im_y;
1514 POP_LONG_COMPLEX (p, &re_y, &im_y);
1515 POP_LONG_COMPLEX (p, &re_x, &im_x);
1516 VALUE (&re_x).f -= VALUE (&re_y).f;
1517 VALUE (&im_x).f -= VALUE (&im_y).f;
1518 CHECK_DOUBLE_COMPLEX (p, VALUE (&im_x).f, VALUE (&im_y).f);
1519 PUSH_OBJECT (p, re_x, A68G_LONG_REAL);
1520 PUSH_OBJECT (p, im_x, A68G_LONG_REAL);
1521 }
1522
1523 //! @brief OP * = (COMPLEX, COMPLEX) COMPLEX
1524
1525 void genie_mul_double_compl (NODE_T * p)
1526 {
1527 A68G_LONG_REAL re_x, im_x, re_y, im_y;
1528 POP_LONG_COMPLEX (p, &re_y, &im_y);
1529 POP_LONG_COMPLEX (p, &re_x, &im_x);
1530 DOUBLE_T re = VALUE (&re_x).f * VALUE (&re_y).f - VALUE (&im_x).f * VALUE (&im_y).f;
1531 DOUBLE_T im = VALUE (&im_x).f * VALUE (&re_y).f + VALUE (&re_x).f * VALUE (&im_y).f;
1532 CHECK_DOUBLE_COMPLEX (p, VALUE (&im_x).f, VALUE (&im_y).f);
1533 PUSH_VALUE (p, dble (re), A68G_LONG_REAL);
1534 PUSH_VALUE (p, dble (im), A68G_LONG_REAL);
1535 }
1536
1537 //! @brief OP / = (COMPLEX, COMPLEX) COMPLEX
1538
1539 void genie_div_double_compl (NODE_T * p)
1540 {
1541 A68G_LONG_REAL re_x, im_x, re_y, im_y;
1542 DOUBLE_T re = 0.0, im = 0.0;
1543 POP_LONG_COMPLEX (p, &re_y, &im_y);
1544 POP_LONG_COMPLEX (p, &re_x, &im_x);
1545 PRELUDE_ERROR (VALUE (&re_y).f == 0.0q && VALUE (&im_y).f == 0.0q, p, ERROR_DIVISION_BY_ZERO, M_LONG_COMPLEX);
1546 if (ABSQ (VALUE (&re_y).f) >= ABSQ (VALUE (&im_y).f)) {
1547 DOUBLE_T r = VALUE (&im_y).f / VALUE (&re_y).f, den = VALUE (&re_y).f + r * VALUE (&im_y).f;
1548 re = (VALUE (&re_x).f + r * VALUE (&im_x).f) / den;
1549 im = (VALUE (&im_x).f - r * VALUE (&re_x).f) / den;
1550 } else {
1551 DOUBLE_T r = VALUE (&re_y).f / VALUE (&im_y).f, den = VALUE (&im_y).f + r * VALUE (&re_y).f;
1552 re = (VALUE (&re_x).f * r + VALUE (&im_x).f) / den;
1553 im = (VALUE (&im_x).f * r - VALUE (&re_x).f) / den;
1554 }
1555 PUSH_VALUE (p, dble (re), A68G_LONG_REAL);
1556 PUSH_VALUE (p, dble (im), A68G_LONG_REAL);
1557 }
1558
1559 //! @brief OP ** = (LONG COMPLEX, INT) LONG COMPLEX
1560
1561 void genie_pow_double_compl_int (NODE_T * p)
1562 {
1563 A68G_LONG_REAL re_x, im_x;
1564 A68G_INT j;
1565 BOOL_T negative;
1566 POP_OBJECT (p, &j, A68G_INT);
1567 POP_LONG_COMPLEX (p, &re_x, &im_x);
1568 DOUBLE_T re_z = 1.0q, im_z = 0.0q;
1569 DOUBLE_T re_y = VALUE (&re_x).f, im_y = VALUE (&im_x).f;
1570 negative = (BOOL_T) (VALUE (&j) < 0);
1571 if (negative) {
1572 VALUE (&j) = -VALUE (&j);
1573 }
1574 INT_T expo = 1;
1575 while ((UNSIGNED_T) expo <= (UNSIGNED_T) (VALUE (&j))) {
1576 DOUBLE_T z;
1577 if (expo & VALUE (&j)) {
1578 z = re_z * re_y - im_z * im_y;
1579 im_z = re_z * im_y + im_z * re_y;
1580 re_z = z;
1581 }
1582 z = re_y * re_y - im_y * im_y;
1583 im_y = im_y * re_y + re_y * im_y;
1584 re_y = z;
1585 CHECK_DOUBLE_COMPLEX (p, re_y, im_y);
1586 CHECK_DOUBLE_COMPLEX (p, re_z, im_z);
1587 expo <<= 1;
1588 }
1589 if (negative) {
1590 PUSH_VALUE (p, dble (1.0q), A68G_LONG_REAL);
1591 PUSH_VALUE (p, dble (0.0q), A68G_LONG_REAL);
1592 PUSH_VALUE (p, dble (re_z), A68G_LONG_REAL);
1593 PUSH_VALUE (p, dble (im_z), A68G_LONG_REAL);
1594 genie_div_double_compl (p);
1595 } else {
1596 PUSH_VALUE (p, dble (re_z), A68G_LONG_REAL);
1597 PUSH_VALUE (p, dble (im_z), A68G_LONG_REAL);
1598 }
1599 }
1600
1601 //! @brief OP = = (COMPLEX, COMPLEX) BOOL
1602
1603 void genie_eq_double_compl (NODE_T * p)
1604 {
1605 A68G_LONG_REAL re_x, im_x, re_y, im_y;
1606 POP_LONG_COMPLEX (p, &re_y, &im_y);
1607 POP_LONG_COMPLEX (p, &re_x, &im_x);
1608 PUSH_VALUE (p, (BOOL_T) ((VALUE (&re_x).f == VALUE (&re_y).f) && (VALUE (&im_x).f == VALUE (&im_y).f)), A68G_BOOL);
1609 }
1610
1611 //! @brief OP /= = (COMPLEX, COMPLEX) BOOL
1612
1613 void genie_ne_double_compl (NODE_T * p)
1614 {
1615 A68G_LONG_REAL re_x, im_x, re_y, im_y;
1616 POP_LONG_COMPLEX (p, &re_y, &im_y);
1617 POP_LONG_COMPLEX (p, &re_x, &im_x);
1618 PUSH_VALUE (p, (BOOL_T) ! ((VALUE (&re_x).f == VALUE (&re_y).f) && (VALUE (&im_x).f == VALUE (&im_y).f)), A68G_BOOL);
1619 }
1620
1621 //! @brief OP +:= = (REF COMPLEX, COMPLEX) REF COMPLEX
1622
1623 void genie_plusab_double_compl (NODE_T * p)
1624 {
1625 genie_f_and_becomes (p, M_REF_LONG_COMPLEX, genie_add_double_compl);
1626 }
1627
1628 //! @brief OP -:= = (REF COMPLEX, COMPLEX) REF COMPLEX
1629
1630 void genie_minusab_double_compl (NODE_T * p)
1631 {
1632 genie_f_and_becomes (p, M_REF_LONG_COMPLEX, genie_sub_double_compl);
1633 }
1634
1635 //! @brief OP *:= = (REF COMPLEX, COMPLEX) REF COMPLEX
1636
1637 void genie_timesab_double_compl (NODE_T * p)
1638 {
1639 genie_f_and_becomes (p, M_REF_LONG_COMPLEX, genie_mul_double_compl);
1640 }
1641
1642 //! @brief OP /:= = (REF COMPLEX, COMPLEX) REF COMPLEX
1643
1644 void genie_divab_double_compl (NODE_T * p)
1645 {
1646 genie_f_and_becomes (p, M_REF_LONG_COMPLEX, genie_div_double_compl);
1647 }
1648
1649 //! @brief OP LENG = (COMPLEX) LONG COMPLEX
1650
1651 void genie_lengthen_complex_to_double_compl (NODE_T * p)
1652 {
1653 A68G_REAL i;
1654 POP_OBJECT (p, &i, A68G_REAL);
1655 genie_lengthen_real_to_double (p);
1656 PUSH_OBJECT (p, i, A68G_REAL);
1657 genie_lengthen_real_to_double (p);
1658 }
1659
1660 // Functions
1661
1662 //! @brief OP ROUND = (LONG REAL) LONG INT
1663
1664 void genie_round_double (NODE_T * p)
1665 {
1666 A68G_LONG_REAL x;
1667 POP_OBJECT (p, &x, A68G_LONG_REAL);
1668 DOUBLE_NUM_T u = VALUE (&x);
1669 if (u.f < 0.0q) {
1670 u.f = u.f - 0.5q;
1671 } else {
1672 u.f = u.f + 0.5q;
1673 }
1674 PUSH_VALUE (p, double_to_double_int (p, u), A68G_LONG_INT);
1675 }
1676
1677 //! @brief OP ENTIER = (LONG REAL) LONG INT
1678
1679 void genie_entier_double (NODE_T * p)
1680 {
1681 A68G_LONG_REAL x;
1682 POP_OBJECT (p, &x, A68G_LONG_REAL);
1683 DOUBLE_NUM_T u = VALUE (&x);
1684 u.f = floor_double (u.f);
1685 PUSH_VALUE (p, double_to_double_int (p, u), A68G_LONG_INT);
1686 }
1687
1688 //! @brief OP CEIL = (LONG REAL) LONG INT
1689
1690 void genie_ceil_double (NODE_T * p)
1691 {
1692 A68G_LONG_REAL x;
1693 POP_OBJECT (p, &x, A68G_LONG_REAL);
1694 DOUBLE_NUM_T u = VALUE (&x);
1695 u.f = ceil_double (u.f);
1696 PUSH_VALUE (p, double_to_double_int (p, u), A68G_LONG_INT);
1697 }
1698
1699 //! @brief OP TRUNC = (LONG REAL) LONG INT
1700
1701 void genie_trunc_double (NODE_T * p)
1702 {
1703 A68G_LONG_REAL x;
1704 POP_OBJECT (p, &x, A68G_LONG_REAL);
1705 DOUBLE_NUM_T u = VALUE (&x);
1706 PUSH_VALUE (p, double_to_double_int (p, u), A68G_LONG_INT);
1707 }
1708
1709 //! @brief OP FRAC = (LONG REAL) LONG REAL
1710
1711 void genie_frac_double (NODE_T * p)
1712 {
1713 A68G_LONG_REAL x;
1714 POP_OBJECT (p, &x, A68G_LONG_REAL);
1715 DOUBLE_NUM_T u = VALUE (&x), v, w;
1716 v.f = fabs_double (u.f);
1717 w.f = v.f - floor_double (v.f);
1718 if (u.f < 0.0q) {
1719 w.f = -w.f;
1720 }
1721 PUSH_VALUE (p, w, A68G_LONG_REAL);
1722 }
1723
1724 #define CD_FUNCTION(name, fun)\
1725 void name (NODE_T * p) {\
1726 A68G_LONG_REAL *x;\
1727 POP_OPERAND_ADDRESS (p, x, A68G_LONG_REAL);\
1728 errno = 0;\
1729 VALUE (x).f = fun (VALUE (x).f);\
1730 MATH_RTE (p, errno != 0, M_LONG_REAL, NO_TEXT);\
1731 }
1732
1733 CD_FUNCTION (genie_acos_double, acos_double);
1734 CD_FUNCTION (genie_acosh_double, acosh_double);
1735 CD_FUNCTION (genie_asinh_double, asinh_double);
1736 CD_FUNCTION (genie_atanh_double, atanh_double);
1737 CD_FUNCTION (genie_asin_double, asin_double);
1738 CD_FUNCTION (genie_atan_double, atan_double);
1739 CD_FUNCTION (genie_cosh_double, cosh_double);
1740 CD_FUNCTION (genie_cos_double, cos_double);
1741 CD_FUNCTION (genie_curt_double, cbrt_double);
1742 CD_FUNCTION (genie_exp_double, exp_double);
1743 CD_FUNCTION (genie_ln_double, log_double);
1744 CD_FUNCTION (genie_log_double, log10_double);
1745 CD_FUNCTION (genie_sinh_double, sinh_double);
1746 CD_FUNCTION (genie_sin_double, sin_double);
1747 CD_FUNCTION (genie_sqrt_double, sqrt_double);
1748 CD_FUNCTION (genie_tanh_double, tanh_double);
1749 CD_FUNCTION (genie_tan_double, tan_double);
1750 CD_FUNCTION (genie_erf_double, erf_double);
1751 CD_FUNCTION (genie_erfc_double, erfc_double);
1752 CD_FUNCTION (genie_lngamma_double, lgamma_double);
1753 CD_FUNCTION (genie_gamma_double, tgamma_double);
1754 CD_FUNCTION (genie_csc_double, a68g_csc_double);
1755 CD_FUNCTION (genie_cscdg_double, a68g_cscdg_double);
1756 CD_FUNCTION (genie_acsc_double, a68g_acsc_double);
1757 CD_FUNCTION (genie_acscdg_double, a68g_acscdg_double);
1758 CD_FUNCTION (genie_sec_double, a68g_sec_double);
1759 CD_FUNCTION (genie_secdg_double, a68g_secdg_double);
1760 CD_FUNCTION (genie_asec_double, a68g_asec_double);
1761 CD_FUNCTION (genie_asecdg_double, a68g_asecdg_double);
1762 CD_FUNCTION (genie_cot_double, a68g_cot_double);
1763 CD_FUNCTION (genie_acot_double, a68g_acot_double);
1764 CD_FUNCTION (genie_sindg_double, a68g_sindg_double);
1765 CD_FUNCTION (genie_cas_double, a68g_cas_double);
1766 CD_FUNCTION (genie_cosdg_double, a68g_cosdg_double);
1767 CD_FUNCTION (genie_tandg_double, a68g_tandg_double);
1768 CD_FUNCTION (genie_asindg_double, a68g_asindg_double);
1769 CD_FUNCTION (genie_acosdg_double, a68g_acosdg_double);
1770 CD_FUNCTION (genie_atandg_double, a68g_atandg_double);
1771 CD_FUNCTION (genie_cotdg_double, a68g_cotdg_double);
1772 CD_FUNCTION (genie_acotdg_double, a68g_acotdg_double);
1773 CD_FUNCTION (genie_sinpi_double, a68g_sinpi_double);
1774 CD_FUNCTION (genie_cospi_double, a68g_cospi_double);
1775 CD_FUNCTION (genie_tanpi_double, a68g_tanpi_double);
1776 CD_FUNCTION (genie_cotpi_double, a68g_cotpi_double);
1777
1778 //! @brief PROC long arctan2 = (LONG REAL) LONG REAL
1779
1780 void genie_atan2_double (NODE_T * p)
1781 {
1782 A68G_LONG_REAL x, y;
1783 POP_OBJECT (p, &y, A68G_LONG_REAL);
1784 POP_OBJECT (p, &x, A68G_LONG_REAL);
1785 errno = 0;
1786 PRELUDE_ERROR (VALUE (&x).f == 0.0q && VALUE (&y).f == 0.0q, p, ERROR_INVALID_ARGUMENT, M_LONG_REAL);
1787 VALUE (&x).f = a68g_atan2_real (VALUE (&y).f, VALUE (&x).f);
1788 PRELUDE_ERROR (errno != 0, p, ERROR_MATH_EXCEPTION, NO_TEXT);
1789 PUSH_OBJECT (p, x, A68G_LONG_REAL);
1790 }
1791
1792 //! @brief PROC long arctan2dg = (LONG REAL) LONG REAL
1793
1794 void genie_atan2dg_double (NODE_T * p)
1795 {
1796 A68G_LONG_REAL x, y;
1797 POP_OBJECT (p, &y, A68G_LONG_REAL);
1798 POP_OBJECT (p, &x, A68G_LONG_REAL);
1799 errno = 0;
1800 PRELUDE_ERROR (VALUE (&x).f == 0.0q && VALUE (&y).f == 0.0q, p, ERROR_INVALID_ARGUMENT, M_LONG_REAL);
1801 VALUE (&x).f = CONST_180_OVER_PI_Q * a68g_atan2_real (VALUE (&y).f, VALUE (&x).f);
1802 PRELUDE_ERROR (errno != 0, p, ERROR_MATH_EXCEPTION, NO_TEXT);
1803 PUSH_OBJECT (p, x, A68G_LONG_REAL);
1804 }
1805
1806 //! @brief PROC (LONG REAL) LONG REAL inverf
1807
1808 void genie_inverf_double (NODE_T * _p_)
1809 {
1810 A68G_LONG_REAL x;
1811 DOUBLE_T y, z;
1812 POP_OBJECT (_p_, &x, A68G_LONG_REAL);
1813 errno = 0;
1814 y = VALUE (&x).f;
1815 z = inverf_double (y);
1816 MATH_RTE (_p_, errno != 0, M_LONG_REAL, NO_TEXT);
1817 CHECK_DOUBLE_REAL (_p_, z);
1818 PUSH_VALUE (_p_, dble (z), A68G_LONG_REAL);
1819 }
1820
1821 //! @brief PROC (LONG REAL) LONG REAL inverfc
1822
1823 void genie_inverfc_double (NODE_T * p)
1824 {
1825 A68G_LONG_REAL *u;
1826 POP_OPERAND_ADDRESS (p, u, A68G_LONG_REAL);
1827 VALUE (u).f = 1.0q - (VALUE (u).f);
1828 genie_inverf_double (p);
1829 }
1830
1831 #define CD_C_FUNCTION(p, g)\
1832 A68G_LONG_REAL re, im;\
1833 DOUBLE_COMPLEX_T z;\
1834 POP_OBJECT (p, &im, A68G_LONG_REAL);\
1835 POP_OBJECT (p, &re, A68G_LONG_REAL);\
1836 errno = 0;\
1837 z = VALUE (&re).f + VALUE (&im).f * _Complex_I;\
1838 z = g (z);\
1839 PUSH_VALUE (p, dble ((DOUBLE_T) creal_double (z)), A68G_LONG_REAL);\
1840 PUSH_VALUE (p, dble ((DOUBLE_T) cimag_double (z)), A68G_LONG_REAL);\
1841 MATH_RTE (p, errno != 0, M_COMPLEX, NO_TEXT);
1842
1843 //! @brief PROC long csqrt = (LONG COMPLEX) LONG COMPLEX
1844
1845 void genie_sqrt_double_compl (NODE_T * p)
1846 {
1847 CD_C_FUNCTION (p, csqrt_double);
1848 }
1849
1850 //! @brief PROC long csin = (LONG COMPLEX) LONG COMPLEX
1851
1852 void genie_sin_double_compl (NODE_T * p)
1853 {
1854 CD_C_FUNCTION (p, csin_double);
1855 }
1856
1857 //! @brief PROC long ccos = (LONG COMPLEX) LONG COMPLEX
1858
1859 void genie_cos_double_compl (NODE_T * p)
1860 {
1861 CD_C_FUNCTION (p, ccos_double);
1862 }
1863
1864 //! @brief PROC long ctan = (LONG COMPLEX) LONG COMPLEX
1865
1866 void genie_tan_double_compl (NODE_T * p)
1867 {
1868 CD_C_FUNCTION (p, ctan_double);
1869 }
1870
1871 //! @brief PROC long casin = (LONG COMPLEX) LONG COMPLEX
1872
1873 void genie_asin_double_compl (NODE_T * p)
1874 {
1875 CD_C_FUNCTION (p, casin_double);
1876 }
1877
1878 //! @brief PROC long cacos = (LONG COMPLEX) LONG COMPLEX
1879
1880 void genie_acos_double_compl (NODE_T * p)
1881 {
1882 CD_C_FUNCTION (p, cacos_double);
1883 }
1884
1885 //! @brief PROC long catan = (LONG COMPLEX) LONG COMPLEX
1886
1887 void genie_atan_double_compl (NODE_T * p)
1888 {
1889 CD_C_FUNCTION (p, catan_double);
1890 }
1891
1892 //! @brief PROC long cexp = (LONG COMPLEX) LONG COMPLEX
1893
1894 void genie_exp_double_compl (NODE_T * p)
1895 {
1896 CD_C_FUNCTION (p, cexp_double);
1897 }
1898
1899 //! @brief PROC long cln = (LONG COMPLEX) LONG COMPLEX
1900
1901 void genie_ln_double_compl (NODE_T * p)
1902 {
1903 CD_C_FUNCTION (p, clog_double);
1904 }
1905
1906 //! @brief PROC long csinh = (LONG COMPLEX) LONG COMPLEX
1907
1908 void genie_sinh_double_compl (NODE_T * p)
1909 {
1910 CD_C_FUNCTION (p, csinh_double);
1911 }
1912
1913 //! @brief PROC long ccosh = (LONG COMPLEX) LONG COMPLEX
1914
1915 void genie_cosh_double_compl (NODE_T * p)
1916 {
1917 CD_C_FUNCTION (p, ccosh_double);
1918 }
1919
1920 //! @brief PROC long ctanh = (LONG COMPLEX) LONG COMPLEX
1921
1922 void genie_tanh_double_compl (NODE_T * p)
1923 {
1924 CD_C_FUNCTION (p, ctanh_double);
1925 }
1926
1927 //! @brief PROC long casinh = (LONG COMPLEX) LONG COMPLEX
1928
1929 void genie_asinh_double_compl (NODE_T * p)
1930 {
1931 CD_C_FUNCTION (p, casinh_double);
1932 }
1933
1934 //! @brief PROC long cacosh = (LONG COMPLEX) LONG COMPLEX
1935
1936 void genie_acosh_double_compl (NODE_T * p)
1937 {
1938 CD_C_FUNCTION (p, cacosh_double);
1939 }
1940
1941 //! @brief PROC long catanh = (LONG COMPLEX) LONG COMPLEX
1942
1943 void genie_atanh_double_compl (NODE_T * p)
1944 {
1945 CD_C_FUNCTION (p, catanh_double);
1946 }
1947
1948 //! @brief PROC next long random = LONG REAL
1949
1950 void genie_next_random_double (NODE_T * p)
1951 {
1952 // This is 'real width' digits only.
1953 genie_next_random (p);
1954 genie_lengthen_real_to_double (p);
1955 }
1956
1957 #define CALL(g, x, y) {\
1958 ADDR_T pop_sp = A68G_SP;\
1959 A68G_LONG_REAL *z = (A68G_LONG_REAL *) STACK_TOP;\
1960 DOUBLE_NUM_T _w_;\
1961 _w_.f = (x);\
1962 PUSH_VALUE (_p_, _w_, A68G_LONG_REAL);\
1963 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);\
1964 (y) = VALUE (z).f;\
1965 A68G_SP = pop_sp;\
1966 }
1967
1968 //! @brief Transform string into real-16.
1969
1970 DOUBLE_T string_to_double (char *s, char **end)
1971 {
1972 errno = 0;
1973 DOUBLE_T y[A68G_DOUBLE_DIG];
1974 for (int i = 0; i < A68G_DOUBLE_DIG; i++) {
1975 y[i] = 0.0q;
1976 }
1977 if (end != NO_REF) {
1978 (*end) = &(s[0]);
1979 }
1980 while (IS_SPACE (s[0])) {
1981 s++;
1982 }
1983 // Scan mantissa digits and put them into "y".
1984 DOUBLE_T W;
1985 if (s[0] == '-') {
1986 W = -1.0q;
1987 } else {
1988 W = 1.0q;
1989 }
1990 if (s[0] == '+' || s[0] == '-') {
1991 s++;
1992 }
1993 while (s[0] == '0') {
1994 s++;
1995 }
1996 int dot = -1, pos = 0, pow = 0;
1997 while (pow < A68G_DOUBLE_DIG && s[pos] != NULL_CHAR && (IS_DIGIT (s[pos]) || s[pos] == POINT_CHAR)) {
1998 if (s[pos] == POINT_CHAR) {
1999 dot = pos;
2000 } else {
2001 int val = (int) s[pos] - (int) '0';
2002 y[pow] = W * val;
2003 W /= 10.0q;
2004 pow++;
2005 }
2006 pos++;
2007 }
2008 while (IS_DIGIT (s[pos])) {
2009 pos++;
2010 }
2011 if (end != NO_REF) {
2012 (*end) = &(s[pos]);
2013 }
2014 // Sum from low to high to preserve precision.
2015 DOUBLE_T sum = 0.0q;
2016 for (int i = A68G_DOUBLE_DIG - 1; i >= 0; i--) {
2017 sum = sum + y[i];
2018 }
2019 // See if there is an exponent.
2020 int expo;
2021 if (s[pos] != NULL_CHAR && TO_UPPER (s[pos]) == TO_UPPER (EXPONENT_CHAR)) {
2022 expo = (int) strtol (&(s[++pos]), end, 10);
2023 } else {
2024 expo = 0;
2025 }
2026 // Standardise.
2027 if (dot >= 0) {
2028 expo += dot - 1;
2029 } else {
2030 expo += pow - 1;
2031 }
2032 while (sum != 0.0q && fabs_double (sum) < 1.0q) {
2033 sum *= 10.0q;
2034 expo -= 1;
2035 }
2036 if (errno == 0) {
2037 return sum * ten_up_double (expo);
2038 } else {
2039 return 0.0q;
2040 }
2041 }
2042
2043 void genie_beta_inc_cf_double (NODE_T * p)
2044 {
2045 A68G_LONG_REAL x, s, t;
2046 POP_OBJECT (p, &x, A68G_LONG_REAL);
2047 POP_OBJECT (p, &t, A68G_LONG_REAL);
2048 POP_OBJECT (p, &s, A68G_LONG_REAL);
2049 errno = 0;
2050 PUSH_VALUE (p, dble (a68g_beta_inc_double (VALUE (&s).f, VALUE (&t).f, VALUE (&x).f)), A68G_LONG_REAL);
2051 MATH_RTE (p, errno != 0, M_LONG_REAL, NO_TEXT);
2052 }
2053
2054 void genie_beta_double (NODE_T * p)
2055 {
2056 A68G_LONG_REAL a, b;
2057 POP_OBJECT (p, &b, A68G_LONG_REAL);
2058 POP_OBJECT (p, &a, A68G_LONG_REAL);
2059 errno = 0;
2060 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);
2061 MATH_RTE (p, errno != 0, M_LONG_REAL, NO_TEXT);
2062 }
2063
2064 void genie_ln_beta_double (NODE_T * p)
2065 {
2066 A68G_LONG_REAL a, b;
2067 POP_OBJECT (p, &b, A68G_LONG_REAL);
2068 POP_OBJECT (p, &a, A68G_LONG_REAL);
2069 errno = 0;
2070 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);
2071 MATH_RTE (p, errno != 0, M_LONG_REAL, NO_TEXT);
2072 }
2073
2074 // LONG REAL infinity
2075
2076 void genie_infinity_double (NODE_T * p)
2077 {
2078 PUSH_VALUE (p, dble (A68G_POSINF_DOUBLE), A68G_LONG_REAL);
2079 }
2080
2081 // LONG REAL minus infinity
2082
2083 void genie_minus_infinity_double (NODE_T * p)
2084 {
2085 PUSH_VALUE (p, dble (A68G_MININF_DOUBLE), A68G_LONG_REAL);
2086 }
2087
2088 // Range check.
2089
2090 BOOL_T a68g_finite_double (REAL_T x)
2091 {
2092 return !isinfq (x) && !isnanq (x);
2093 }
2094
2095 // Range check.
2096
2097 BOOL_T a68g_isnan_double (REAL_T x)
2098 {
2099 return isnanq (x);
2100 }
2101
2102 // Range check.
2103
2104 BOOL_T a68g_isinf_double (REAL_T x)
2105 {
2106 return isinfq (x);
2107 }
2108
2109 #endif
© 2002-2026 J.M. van der Veer
jmvdveer@algol68genie.nl
