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