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