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