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 int k = sign_double_int (double_ssub (p, VALUE (&i), VALUE (&j)));\
958 PUSH_VALUE (p, (BOOL_T) (k OP 0), A68_BOOL);\
959 }
960
961 A68_CMP_INT (genie_eq_double_int, ==);
962 A68_CMP_INT (genie_ne_double_int, !=);
963 A68_CMP_INT (genie_lt_double_int, <);
964 A68_CMP_INT (genie_gt_double_int, >);
965 A68_CMP_INT (genie_le_double_int, <=);
966 A68_CMP_INT (genie_ge_double_int, >=);
967
968 // OP (LONG REAL, LONG REAL) BOOL.
969 #define A68_CMP_REAL(n, OP)\
970 void n (NODE_T * p) {\
971 A68_LONG_REAL i, j;\
972 POP_OBJECT (p, &j, A68_LONG_REAL);\
973 POP_OBJECT (p, &i, A68_LONG_REAL);\
974 PUSH_VALUE (p, (BOOL_T) (VALUE (&i).f OP VALUE (&j).f), A68_BOOL);\
975 }
976
977 A68_CMP_REAL (genie_eq_double, ==);
978 A68_CMP_REAL (genie_ne_double, !=);
979 A68_CMP_REAL (genie_lt_double, <);
980 A68_CMP_REAL (genie_gt_double, >);
981 A68_CMP_REAL (genie_le_double, <=);
982 A68_CMP_REAL (genie_ge_double, >=);
983
984 //! @brief OP NOT = (LONG BITS) LONG BITS
985
986 void genie_not_double_bits (NODE_T * p)
987 {
988 A68_LONG_BITS i;
989 POP_OBJECT (p, &i, A68_LONG_BITS);
990 DOUBLE_NUM_T w;
991 HW (w) = ~HW (VALUE (&i));
992 LW (w) = ~LW (VALUE (&i));
993 PUSH_VALUE (p, w, A68_LONG_BITS);
994 }
995
996 //! @brief OP = = (LONG BITS, LONG BITS) BOOL.
997
998 void genie_eq_double_bits (NODE_T * p)
999 {
1000 A68_LONG_BITS i, j;
1001 POP_OBJECT (p, &j, A68_LONG_BITS);
1002 POP_OBJECT (p, &i, A68_LONG_BITS);
1003 BOOL_T u = HW (VALUE (&i)) == HW (VALUE (&j));
1004 BOOL_T v = LW (VALUE (&i)) == LW (VALUE (&j));
1005 PUSH_VALUE (p, (BOOL_T) (u & v ? A68_TRUE : A68_FALSE), A68_BOOL);
1006 }
1007
1008 //! @brief OP ~= = (LONG BITS, LONG BITS) BOOL.
1009
1010 void genie_ne_double_bits (NODE_T * p)
1011 {
1012 A68_LONG_BITS i, j;
1013 POP_OBJECT (p, &j, A68_LONG_BITS); // (i ~= j) == ~ (i = j)
1014 POP_OBJECT (p, &i, A68_LONG_BITS);
1015 BOOL_T u = HW (VALUE (&i)) == HW (VALUE (&j));
1016 BOOL_T v = LW (VALUE (&i)) == LW (VALUE (&j));
1017 PUSH_VALUE (p, (BOOL_T) (u & v ? A68_FALSE : A68_TRUE), A68_BOOL);
1018 }
1019
1020 //! @brief OP <= = (LONG BITS, LONG BITS) BOOL
1021
1022 void genie_le_double_bits (NODE_T * p)
1023 {
1024 A68_LONG_BITS i, j;
1025 POP_OBJECT (p, &j, A68_LONG_BITS);
1026 POP_OBJECT (p, &i, A68_LONG_BITS);
1027 BOOL_T u = (HW (VALUE (&i)) | HW (VALUE (&j))) == HW (VALUE (&j));
1028 BOOL_T v = (LW (VALUE (&i)) | LW (VALUE (&j))) == LW (VALUE (&j));
1029 PUSH_VALUE (p, (BOOL_T) (u & v ? A68_TRUE : A68_FALSE), A68_BOOL);
1030 }
1031
1032 //! @brief OP > = (LONG BITS, LONG BITS) BOOL
1033
1034 void genie_gt_double_bits (NODE_T * p)
1035 {
1036 A68_LONG_BITS i, j;
1037 POP_OBJECT (p, &j, A68_LONG_BITS); // (i > j) == ! (i <= j)
1038 POP_OBJECT (p, &i, A68_LONG_BITS);
1039 BOOL_T u = (HW (VALUE (&i)) | HW (VALUE (&j))) == HW (VALUE (&j));
1040 BOOL_T v = (LW (VALUE (&i)) | LW (VALUE (&j))) == LW (VALUE (&j));
1041 PUSH_VALUE (p, (BOOL_T) (u & v ? A68_FALSE : A68_TRUE), A68_BOOL);
1042 }
1043
1044 //! @brief OP >= = (LONG BITS, LONG BITS) BOOL
1045
1046 void genie_ge_double_bits (NODE_T * p)
1047 {
1048 A68_LONG_BITS i, j;
1049 POP_OBJECT (p, &j, A68_LONG_BITS); // (i >= j) == (j <= i)
1050 POP_OBJECT (p, &i, A68_LONG_BITS);
1051 BOOL_T u = (HW (VALUE (&i)) | HW (VALUE (&j))) == HW (VALUE (&i));
1052 BOOL_T v = (LW (VALUE (&i)) | LW (VALUE (&j))) == LW (VALUE (&i));
1053 PUSH_VALUE (p, (BOOL_T) (u & v ? A68_TRUE : A68_FALSE), A68_BOOL);
1054 }
1055
1056 //! @brief OP < = (LONG BITS, LONG BITS) BOOL
1057
1058 void genie_lt_double_bits (NODE_T * p)
1059 {
1060 A68_LONG_BITS i, j;
1061 POP_OBJECT (p, &j, A68_LONG_BITS); // (i < j) == ! (i >= j)
1062 POP_OBJECT (p, &i, A68_LONG_BITS);
1063 BOOL_T u = (HW (VALUE (&i)) | HW (VALUE (&j))) == HW (VALUE (&i));
1064 BOOL_T v = (LW (VALUE (&i)) | LW (VALUE (&j))) == LW (VALUE (&i));
1065 PUSH_VALUE (p, (BOOL_T) (u & v ? A68_FALSE : A68_TRUE), A68_BOOL);
1066 }
1067
1068 //! @brief PROC long bits pack = ([] BOOL) BITS
1069
1070 void genie_double_bits_pack (NODE_T * p)
1071 {
1072 A68_REF z;
1073 POP_REF (p, &z);
1074 CHECK_REF (p, z, M_ROW_BOOL);
1075 A68_ARRAY *arr; A68_TUPLE *tup;
1076 GET_DESCRIPTOR (arr, tup, &z);
1077 int size = ROW_SIZE (tup);
1078 PRELUDE_ERROR (size < 0 || size > A68_BITS_WIDTH, p, ERROR_OUT_OF_BOUNDS, M_ROW_BOOL);
1079 DOUBLE_NUM_T w;
1080 set_lw (w, 0x0);
1081 if (ROW_SIZE (tup) > 0) {
1082 UNSIGNED_T bit = 0x0;
1083 BYTE_T *base = DEREF (BYTE_T, &ARRAY (arr));
1084 int n = 0;
1085 for (int k = UPB (tup); k >= LWB (tup); k--) {
1086 A68_BOOL *boo = (A68_BOOL *) & (base[INDEX_1_DIM (arr, tup, k)]);
1087 CHECK_INIT (p, INITIALISED (boo), M_BOOL);
1088 if (n == 0 || n == A68_BITS_WIDTH) {
1089 bit = 0x1;
1090 }
1091 if (VALUE (boo)) {
1092 if (n > A68_BITS_WIDTH) {
1093 LW (w) |= bit;
1094 } else {
1095 HW (w) |= bit;
1096 };
1097 }
1098 n++;
1099 bit <<= 1;
1100 }
1101 }
1102 PUSH_VALUE (p, w, A68_LONG_BITS);
1103 }
1104
1105 //! @brief OP AND = (LONG BITS, LONG BITS) LONG BITS
1106
1107 void genie_and_double_bits (NODE_T * p)
1108 {
1109 A68_LONG_BITS i, j;
1110 POP_OBJECT (p, &j, A68_LONG_BITS);
1111 POP_OBJECT (p, &i, A68_LONG_BITS);
1112 DOUBLE_NUM_T w;
1113 HW (w) = HW (VALUE (&i)) & HW (VALUE (&j));
1114 LW (w) = LW (VALUE (&i)) & LW (VALUE (&j));
1115 PUSH_VALUE (p, w, A68_LONG_BITS);
1116 }
1117
1118 //! @brief OP OR = (LONG BITS, LONG BITS) LONG BITS
1119
1120 void genie_or_double_bits (NODE_T * p)
1121 {
1122 A68_LONG_BITS i, j;
1123 POP_OBJECT (p, &j, A68_LONG_BITS);
1124 POP_OBJECT (p, &i, A68_LONG_BITS);
1125 DOUBLE_NUM_T w;
1126 HW (w) = HW (VALUE (&i)) | HW (VALUE (&j));
1127 LW (w) = LW (VALUE (&i)) | LW (VALUE (&j));
1128 PUSH_VALUE (p, w, A68_LONG_BITS);
1129 }
1130
1131 //! @brief OP XOR = (LONG BITS, LONG BITS) LONG BITS
1132
1133 void genie_xor_double_bits (NODE_T * p)
1134 {
1135 A68_LONG_BITS i, j;
1136 POP_OBJECT (p, &j, A68_LONG_BITS);
1137 POP_OBJECT (p, &i, A68_LONG_BITS);
1138 DOUBLE_NUM_T w;
1139 HW (w) = HW (VALUE (&i)) ^ HW (VALUE (&j));
1140 LW (w) = LW (VALUE (&i)) ^ LW (VALUE (&j));
1141 PUSH_VALUE (p, w, A68_LONG_BITS);
1142 }
1143
1144 //! @brief OP + = (LONG BITS, LONG BITS) LONG BITS
1145
1146 void genie_add_double_bits (NODE_T * p)
1147 {
1148 A68_LONG_BITS i, j;
1149 POP_OBJECT (p, &j, A68_LONG_BITS);
1150 POP_OBJECT (p, &i, A68_LONG_BITS);
1151 DOUBLE_NUM_T w;
1152 add_double (p, M_LONG_BITS, w, VALUE (&i), VALUE (&j));
1153 PUSH_VALUE (p, w, A68_LONG_BITS);
1154 }
1155
1156 //! @brief OP - = (LONG BITS, LONG BITS) LONG BITS
1157
1158 void genie_sub_double_bits (NODE_T * p)
1159 {
1160 A68_LONG_BITS i, j;
1161 POP_OBJECT (p, &j, A68_LONG_BITS);
1162 POP_OBJECT (p, &i, A68_LONG_BITS);
1163 DOUBLE_NUM_T w;
1164 sub_double (p, M_LONG_BITS, w, VALUE (&i), VALUE (&j));
1165 PUSH_VALUE (p, w, A68_LONG_BITS);
1166 }
1167
1168 //! @brief OP * = (LONG BITS, LONG BITS) LONG BITS
1169
1170 void genie_times_double_bits (NODE_T * p)
1171 {
1172 A68_LONG_BITS i, j;
1173 POP_OBJECT (p, &j, A68_LONG_BITS);
1174 POP_OBJECT (p, &i, A68_LONG_BITS);
1175 DOUBLE_NUM_T w = double_umul (p, M_LONG_BITS, VALUE (&i), VALUE (&j));
1176 PUSH_VALUE (p, w, A68_LONG_BITS);
1177 }
1178
1179 //! @brief OP OVER = (LONG BITS, LONG BITS) LONG BITS
1180
1181 void genie_over_double_bits (NODE_T * p)
1182 {
1183 A68_LONG_BITS i, j;
1184 POP_OBJECT (p, &j, A68_LONG_BITS);
1185 POP_OBJECT (p, &i, A68_LONG_BITS);
1186 DOUBLE_NUM_T w = double_udiv (p, M_LONG_BITS, VALUE (&i), VALUE (&j), 0);
1187 PUSH_VALUE (p, w, A68_LONG_BITS);
1188 }
1189
1190 //! @brief OP MOD = (LONG BITS, LONG BITS) LONG BITS
1191
1192 void genie_mod_double_bits (NODE_T * p)
1193 {
1194 A68_LONG_BITS i, j;
1195 DOUBLE_NUM_T w;
1196 POP_OBJECT (p, &j, A68_LONG_BITS);
1197 POP_OBJECT (p, &i, A68_LONG_BITS);
1198 w = double_udiv (p, M_LONG_BITS, VALUE (&i), VALUE (&j), 1);
1199 PUSH_VALUE (p, w, A68_LONG_BITS);
1200 }
1201
1202 //! @brief OP ELEM = (INT, LONG BITS) BOOL
1203
1204 void genie_elem_double_bits (NODE_T * p)
1205 {
1206 A68_LONG_BITS j; A68_INT i;
1207 POP_OBJECT (p, &j, A68_LONG_BITS);
1208 POP_OBJECT (p, &i, A68_INT);
1209 int k = VALUE (&i);
1210 PRELUDE_ERROR (k < 1 || k > A68_LONG_BITS_WIDTH, p, ERROR_OUT_OF_BOUNDS, M_INT);
1211 UNSIGNED_T mask = 0x1, *w;
1212 if (k <= A68_BITS_WIDTH) {
1213 w = &(LW (VALUE (&j)));
1214 } else {
1215 w = &(HW (VALUE (&j)));
1216 }
1217 for (int n = 0; n < (A68_BITS_WIDTH - VALUE (&i)); n++) {
1218 mask = mask << 1;
1219 }
1220 PUSH_VALUE (p, (BOOL_T) ((*w & mask) ? A68_TRUE : A68_FALSE), A68_BOOL);
1221 }
1222
1223 //! @brief OP SET = (INT, LONG BITS) LONG BITS
1224
1225 void genie_set_double_bits (NODE_T * p)
1226 {
1227 A68_LONG_BITS j; A68_INT i;
1228 POP_OBJECT (p, &j, A68_LONG_BITS);
1229 POP_OBJECT (p, &i, A68_INT);
1230 int k = VALUE (&i);
1231 PRELUDE_ERROR (k < 1 || k > A68_LONG_BITS_WIDTH, p, ERROR_OUT_OF_BOUNDS, M_INT);
1232 UNSIGNED_T mask = 0x1, *w;
1233 if (k <= A68_BITS_WIDTH) {
1234 w = &(LW (VALUE (&j)));
1235 } else {
1236 w = &(HW (VALUE (&j)));
1237 }
1238 for (int n = 0; n < (A68_BITS_WIDTH - VALUE (&i)); n++) {
1239 mask = mask << 1;
1240 }
1241 (*w) |= mask;
1242 PUSH_OBJECT (p, j, A68_LONG_BITS);
1243 }
1244
1245 //! @brief OP CLEAR = (INT, LONG BITS) LONG BITS
1246
1247 void genie_clear_double_bits (NODE_T * p)
1248 {
1249 A68_LONG_BITS j; A68_INT i;
1250 POP_OBJECT (p, &j, A68_LONG_BITS);
1251 POP_OBJECT (p, &i, A68_INT);
1252 int k = VALUE (&i);
1253 PRELUDE_ERROR (k < 1 || k > A68_LONG_BITS_WIDTH, p, ERROR_OUT_OF_BOUNDS, M_INT);
1254 UNSIGNED_T mask = 0x1, *w;
1255 if (k <= A68_BITS_WIDTH) {
1256 w = &(LW (VALUE (&j)));
1257 } else {
1258 w = &(HW (VALUE (&j)));
1259 }
1260 for (int n = 0; n < (A68_BITS_WIDTH - VALUE (&i)); n++) {
1261 mask = mask << 1;
1262 }
1263 (*w) &= ~mask;
1264 PUSH_OBJECT (p, j, A68_LONG_BITS);
1265 }
1266
1267 //! @brief OP SHL = (LONG BITS, INT) LONG BITS
1268
1269 void genie_shl_double_bits (NODE_T * p)
1270 {
1271 A68_LONG_BITS i; A68_INT j;
1272 POP_OBJECT (p, &j, A68_INT);
1273 POP_OBJECT (p, &i, A68_LONG_BITS);
1274 DOUBLE_NUM_T *w = &VALUE (&i);
1275 int k = VALUE (&j);
1276 if (VALUE (&j) >= 0) {
1277 for (int n = 0; n < k; n++) {
1278 UNSIGNED_T carry = ((LW (*w) & D_SIGN) ? 0x1 : 0x0);
1279 PRELUDE_ERROR (MODCHK (p, M_LONG_BITS, HW (*w) | D_SIGN), p, ERROR_MATH, M_LONG_BITS);
1280 HW (*w) = (HW (*w) << 1) | carry;
1281 LW (*w) = (LW (*w) << 1);
1282 }
1283 } else {
1284 k = -k;
1285 for (int n = 0; n < k; n++) {
1286 UNSIGNED_T carry = ((HW (*w) & 0x1) ? D_SIGN : 0x0);
1287 HW (*w) = (HW (*w) >> 1);
1288 LW (*w) = (LW (*w) >> 1) | carry;
1289 }
1290 }
1291 PUSH_OBJECT (p, i, A68_LONG_BITS);
1292 }
1293
1294 //! @brief OP SHR = (LONG BITS, INT) LONG BITS
1295
1296 void genie_shr_double_bits (NODE_T * p)
1297 {
1298 A68_INT *j;
1299 POP_OPERAND_ADDRESS (p, j, A68_INT);
1300 VALUE (j) = -VALUE (j);
1301 genie_shl_double_bits (p); // Conform RR
1302 }
1303
1304 //! @brief OP ROL = (LONG BITS, INT) LONG BITS
1305
1306 void genie_rol_double_bits (NODE_T * p)
1307 {
1308 A68_LONG_BITS i; A68_INT j;
1309 DOUBLE_NUM_T *w = &VALUE (&i);
1310 POP_OBJECT (p, &j, A68_INT);
1311 POP_OBJECT (p, &i, A68_LONG_BITS);
1312 int k = VALUE (&j);
1313 if (k >= 0) {
1314 for (int n = 0; n < k; n++) {
1315 UNSIGNED_T carry = ((HW (*w) & D_SIGN) ? 0x1 : 0x0);
1316 UNSIGNED_T carry_between = ((LW (*w) & D_SIGN) ? 0x1 : 0x0);
1317 HW (*w) = (HW (*w) << 1) | carry_between;
1318 LW (*w) = (LW (*w) << 1) | carry;
1319 }
1320 } else {
1321 k = -k;
1322 for (int n = 0; n < k; n++) {
1323 UNSIGNED_T carry = ((LW (*w) & 0x1) ? D_SIGN : 0x0);
1324 UNSIGNED_T carry_between = ((HW (*w) & 0x1) ? D_SIGN : 0x0);
1325 HW (*w) = (HW (*w) >> 1) | carry;
1326 LW (*w) = (LW (*w) >> 1) | carry_between;
1327 }
1328 }
1329 PUSH_OBJECT (p, i, A68_LONG_BITS);
1330 }
1331
1332 //! @brief OP ROR = (LONG BITS, INT) LONG BITS
1333
1334 void genie_ror_double_bits (NODE_T * p)
1335 {
1336 A68_INT *j;
1337 POP_OPERAND_ADDRESS (p, j, A68_INT);
1338 VALUE (j) = -VALUE (j);
1339 genie_rol_double_bits (p); // Conform RR
1340 }
1341
1342 //! @brief OP BIN = (LONG INT) LONG BITS
1343
1344 void genie_bin_double_int (NODE_T * p)
1345 {
1346 A68_LONG_INT i;
1347 POP_OBJECT (p, &i, A68_LONG_INT);
1348 // RR does not convert negative numbers
1349 if (D_NEG (VALUE (&i))) {
1350 errno = EDOM;
1351 diagnostic (A68_RUNTIME_ERROR, p, ERROR_OUT_OF_BOUNDS, M_BITS);
1352 exit_genie (p, A68_RUNTIME_ERROR);
1353 }
1354 PUSH_OBJECT (p, i, A68_LONG_BITS);
1355 }
1356
1357 //! @brief OP +* = (LONG REAL, LONG REAL) LONG COMPLEX
1358
1359 void genie_i_double_compl (NODE_T * p)
1360 {
1361 (void) p;
1362 }
1363
1364 //! @brief OP SHORTEN = (LONG COMPLEX) COMPLEX
1365
1366 void genie_shorten_double_compl_to_complex (NODE_T * p)
1367 {
1368 A68_LONG_REAL re, im;
1369 POP_OBJECT (p, &im, A68_LONG_REAL);
1370 POP_OBJECT (p, &re, A68_LONG_REAL);
1371 REAL_T w = VALUE (&re).f;
1372 PUSH_VALUE (p, w, A68_REAL);
1373 w = VALUE (&im).f;
1374 PUSH_VALUE (p, w, A68_REAL);
1375 }
1376
1377 //! @brief OP LENG = (LONG COMPLEX) LONG LONG COMPLEX
1378
1379 void genie_lengthen_double_compl_to_long_mp_complex (NODE_T * p)
1380 {
1381 int digs = DIGITS (M_LONG_LONG_REAL);
1382 A68_LONG_REAL re, im;
1383 POP_OBJECT (p, &im, A68_LONG_REAL);
1384 POP_OBJECT (p, &re, A68_LONG_REAL);
1385 MP_T *z = nil_mp (p, digs);
1386 (void) double_to_mp (p, z, VALUE (&re).f, digs);
1387 MP_STATUS (z) = (MP_T) INIT_MASK;
1388 z = nil_mp (p, digs);
1389 (void) double_to_mp (p, z, VALUE (&im).f, digs);
1390 MP_STATUS (z) = (MP_T) INIT_MASK;
1391 }
1392
1393 //! @brief OP +* = (LONG INT, LONG INT) LONG COMPLEX
1394
1395 void genie_i_int_double_compl (NODE_T * p)
1396 {
1397 A68_LONG_INT re, im;
1398 POP_OBJECT (p, &im, A68_LONG_INT);
1399 POP_OBJECT (p, &re, A68_LONG_INT);
1400 PUSH_VALUE (p, double_int_to_double (p, VALUE (&re)), A68_LONG_REAL);
1401 PUSH_VALUE (p, double_int_to_double (p, VALUE (&im)), A68_LONG_REAL);
1402 }
1403
1404 //! @brief OP RE = (LONG COMPLEX) LONG REAL
1405
1406 void genie_re_double_compl (NODE_T * p)
1407 {
1408 DECREMENT_STACK_POINTER (p, SIZE (M_LONG_REAL));
1409 }
1410
1411 //! @brief OP IM = (LONG COMPLEX) LONG REAL
1412
1413 void genie_im_double_compl (NODE_T * p)
1414 {
1415 A68_LONG_REAL re, im;
1416 POP_OBJECT (p, &im, A68_LONG_REAL);
1417 POP_OBJECT (p, &re, A68_LONG_REAL);
1418 PUSH_OBJECT (p, im, A68_LONG_REAL);
1419 }
1420
1421 //! @brief OP - = (LONG COMPLEX) LONG COMPLEX
1422
1423 void genie_minus_double_compl (NODE_T * p)
1424 {
1425 A68_LONG_REAL re, im;
1426 POP_OBJECT (p, &im, A68_LONG_REAL);
1427 POP_OBJECT (p, &re, A68_LONG_REAL);
1428 VALUE (&re).f = -VALUE (&re).f;
1429 VALUE (&im).f = -VALUE (&im).f;
1430 PUSH_OBJECT (p, im, A68_LONG_REAL);
1431 PUSH_OBJECT (p, re, A68_LONG_REAL);
1432 }
1433
1434 //! @brief OP ABS = (LONG COMPLEX) LONG REAL
1435
1436 void genie_abs_double_compl (NODE_T * p)
1437 {
1438 A68_LONG_REAL re, im;
1439 POP_LONG_COMPLEX (p, &re, &im);
1440 PUSH_VALUE (p, dble (a68_hypot_double (VALUE (&re).f, VALUE (&im).f)), A68_LONG_REAL);
1441 }
1442
1443 //! @brief OP ARG = (LONG COMPLEX) LONG REAL
1444
1445 void genie_arg_double_compl (NODE_T * p)
1446 {
1447 A68_LONG_REAL re, im;
1448 POP_LONG_COMPLEX (p, &re, &im);
1449 PRELUDE_ERROR (VALUE (&re).f == 0.0q && VALUE (&im).f == 0.0q, p, ERROR_INVALID_ARGUMENT, M_LONG_COMPLEX);
1450 PUSH_VALUE (p, dble (atan2_double (VALUE (&im).f, VALUE (&re).f)), A68_LONG_REAL);
1451 }
1452
1453 //! @brief OP CONJ = (LONG COMPLEX) LONG COMPLEX
1454
1455 void genie_conj_double_compl (NODE_T * p)
1456 {
1457 A68_LONG_REAL im;
1458 POP_OBJECT (p, &im, A68_LONG_REAL);
1459 VALUE (&im).f = -VALUE (&im).f;
1460 PUSH_OBJECT (p, im, A68_LONG_REAL);
1461 }
1462
1463 //! @brief OP + = (COMPLEX, COMPLEX) COMPLEX
1464
1465 void genie_add_double_compl (NODE_T * p)
1466 {
1467 A68_LONG_REAL re_x, im_x, re_y, im_y;
1468 POP_LONG_COMPLEX (p, &re_y, &im_y);
1469 POP_LONG_COMPLEX (p, &re_x, &im_x);
1470 VALUE (&re_x).f += VALUE (&re_y).f;
1471 VALUE (&im_x).f += VALUE (&im_y).f;
1472 CHECK_DOUBLE_COMPLEX (p, VALUE (&im_x).f, VALUE (&im_y).f);
1473 PUSH_OBJECT (p, re_x, A68_LONG_REAL);
1474 PUSH_OBJECT (p, im_x, A68_LONG_REAL);
1475 }
1476
1477 //! @brief OP - = (COMPLEX, COMPLEX) COMPLEX
1478
1479 void genie_sub_double_compl (NODE_T * p)
1480 {
1481 A68_LONG_REAL re_x, im_x, re_y, im_y;
1482 POP_LONG_COMPLEX (p, &re_y, &im_y);
1483 POP_LONG_COMPLEX (p, &re_x, &im_x);
1484 VALUE (&re_x).f -= VALUE (&re_y).f;
1485 VALUE (&im_x).f -= VALUE (&im_y).f;
1486 CHECK_DOUBLE_COMPLEX (p, VALUE (&im_x).f, VALUE (&im_y).f);
1487 PUSH_OBJECT (p, re_x, A68_LONG_REAL);
1488 PUSH_OBJECT (p, im_x, A68_LONG_REAL);
1489 }
1490
1491 //! @brief OP * = (COMPLEX, COMPLEX) COMPLEX
1492
1493 void genie_mul_double_compl (NODE_T * p)
1494 {
1495 A68_LONG_REAL re_x, im_x, re_y, im_y;
1496 POP_LONG_COMPLEX (p, &re_y, &im_y);
1497 POP_LONG_COMPLEX (p, &re_x, &im_x);
1498 DOUBLE_T re = VALUE (&re_x).f * VALUE (&re_y).f - VALUE (&im_x).f * VALUE (&im_y).f;
1499 DOUBLE_T im = VALUE (&im_x).f * VALUE (&re_y).f + VALUE (&re_x).f * VALUE (&im_y).f;
1500 CHECK_DOUBLE_COMPLEX (p, VALUE (&im_x).f, VALUE (&im_y).f);
1501 PUSH_VALUE (p, dble (re), A68_LONG_REAL);
1502 PUSH_VALUE (p, dble (im), A68_LONG_REAL);
1503 }
1504
1505 //! @brief OP / = (COMPLEX, COMPLEX) COMPLEX
1506
1507 void genie_div_double_compl (NODE_T * p)
1508 {
1509 A68_LONG_REAL re_x, im_x, re_y, im_y;
1510 DOUBLE_T re = 0.0, im = 0.0;
1511 POP_LONG_COMPLEX (p, &re_y, &im_y);
1512 POP_LONG_COMPLEX (p, &re_x, &im_x);
1513 PRELUDE_ERROR (VALUE (&re_y).f == 0.0q && VALUE (&im_y).f == 0.0q, p, ERROR_DIVISION_BY_ZERO, M_LONG_COMPLEX);
1514 if (ABSQ (VALUE (&re_y).f) >= ABSQ (VALUE (&im_y).f)) {
1515 DOUBLE_T r = VALUE (&im_y).f / VALUE (&re_y).f, den = VALUE (&re_y).f + r * VALUE (&im_y).f;
1516 re = (VALUE (&re_x).f + r * VALUE (&im_x).f) / den;
1517 im = (VALUE (&im_x).f - r * VALUE (&re_x).f) / den;
1518 } else {
1519 DOUBLE_T r = VALUE (&re_y).f / VALUE (&im_y).f, den = VALUE (&im_y).f + r * VALUE (&re_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 }
1523 PUSH_VALUE (p, dble (re), A68_LONG_REAL);
1524 PUSH_VALUE (p, dble (im), A68_LONG_REAL);
1525 }
1526
1527 //! @brief OP ** = (LONG COMPLEX, INT) LONG COMPLEX
1528
1529 void genie_pow_double_compl_int (NODE_T * p)
1530 {
1531 A68_LONG_REAL re_x, im_x;
1532 A68_INT j;
1533 BOOL_T negative;
1534 POP_OBJECT (p, &j, A68_INT);
1535 POP_LONG_COMPLEX (p, &re_x, &im_x);
1536 DOUBLE_T re_z = 1.0q, im_z = 0.0q;
1537 DOUBLE_T re_y = VALUE (&re_x).f, im_y = VALUE (&im_x).f;
1538 negative = (BOOL_T) (VALUE (&j) < 0);
1539 if (negative) {
1540 VALUE (&j) = -VALUE (&j);
1541 }
1542 INT_T expo = 1;
1543 while ((UNSIGNED_T) expo <= (UNSIGNED_T) (VALUE (&j))) {
1544 DOUBLE_T z;
1545 if (expo & VALUE (&j)) {
1546 z = re_z * re_y - im_z * im_y;
1547 im_z = re_z * im_y + im_z * re_y;
1548 re_z = z;
1549 }
1550 z = re_y * re_y - im_y * im_y;
1551 im_y = im_y * re_y + re_y * im_y;
1552 re_y = z;
1553 CHECK_DOUBLE_COMPLEX (p, re_y, im_y);
1554 CHECK_DOUBLE_COMPLEX (p, re_z, im_z);
1555 expo <<= 1;
1556 }
1557 if (negative) {
1558 PUSH_VALUE (p, dble (1.0q), A68_LONG_REAL);
1559 PUSH_VALUE (p, dble (0.0q), A68_LONG_REAL);
1560 PUSH_VALUE (p, dble (re_z), A68_LONG_REAL);
1561 PUSH_VALUE (p, dble (im_z), A68_LONG_REAL);
1562 genie_div_double_compl (p);
1563 } else {
1564 PUSH_VALUE (p, dble (re_z), A68_LONG_REAL);
1565 PUSH_VALUE (p, dble (im_z), A68_LONG_REAL);
1566 }
1567 }
1568
1569 //! @brief OP = = (COMPLEX, COMPLEX) BOOL
1570
1571 void genie_eq_double_compl (NODE_T * p)
1572 {
1573 A68_LONG_REAL re_x, im_x, re_y, im_y;
1574 POP_LONG_COMPLEX (p, &re_y, &im_y);
1575 POP_LONG_COMPLEX (p, &re_x, &im_x);
1576 PUSH_VALUE (p, (BOOL_T) ((VALUE (&re_x).f == VALUE (&re_y).f) && (VALUE (&im_x).f == VALUE (&im_y).f)), A68_BOOL);
1577 }
1578
1579 //! @brief OP /= = (COMPLEX, COMPLEX) BOOL
1580
1581 void genie_ne_double_compl (NODE_T * p)
1582 {
1583 A68_LONG_REAL re_x, im_x, re_y, im_y;
1584 POP_LONG_COMPLEX (p, &re_y, &im_y);
1585 POP_LONG_COMPLEX (p, &re_x, &im_x);
1586 PUSH_VALUE (p, (BOOL_T) ! ((VALUE (&re_x).f == VALUE (&re_y).f) && (VALUE (&im_x).f == VALUE (&im_y).f)), A68_BOOL);
1587 }
1588
1589 //! @brief OP +:= = (REF COMPLEX, COMPLEX) REF COMPLEX
1590
1591 void genie_plusab_double_compl (NODE_T * p)
1592 {
1593 genie_f_and_becomes (p, M_REF_LONG_COMPLEX, genie_add_double_compl);
1594 }
1595
1596 //! @brief OP -:= = (REF COMPLEX, COMPLEX) REF COMPLEX
1597
1598 void genie_minusab_double_compl (NODE_T * p)
1599 {
1600 genie_f_and_becomes (p, M_REF_LONG_COMPLEX, genie_sub_double_compl);
1601 }
1602
1603 //! @brief OP *:= = (REF COMPLEX, COMPLEX) REF COMPLEX
1604
1605 void genie_timesab_double_compl (NODE_T * p)
1606 {
1607 genie_f_and_becomes (p, M_REF_LONG_COMPLEX, genie_mul_double_compl);
1608 }
1609
1610 //! @brief OP /:= = (REF COMPLEX, COMPLEX) REF COMPLEX
1611
1612 void genie_divab_double_compl (NODE_T * p)
1613 {
1614 genie_f_and_becomes (p, M_REF_LONG_COMPLEX, genie_div_double_compl);
1615 }
1616
1617 //! @brief OP LENG = (COMPLEX) LONG COMPLEX
1618
1619 void genie_lengthen_complex_to_double_compl (NODE_T * p)
1620 {
1621 A68_REAL i;
1622 POP_OBJECT (p, &i, A68_REAL);
1623 genie_lengthen_real_to_double (p);
1624 PUSH_OBJECT (p, i, A68_REAL);
1625 genie_lengthen_real_to_double (p);
1626 }
1627
1628 // Functions
1629
1630 #define CD_FUNCTION(name, fun)\
1631 void name (NODE_T * p) {\
1632 A68_LONG_REAL *x;\
1633 POP_OPERAND_ADDRESS (p, x, A68_LONG_REAL);\
1634 errno = 0;\
1635 VALUE (x).f = fun (VALUE (x).f);\
1636 MATH_RTE (p, errno != 0, M_LONG_REAL, NO_TEXT);\
1637 }
1638
1639 CD_FUNCTION (genie_acos_double, acos_double);
1640 CD_FUNCTION (genie_acosh_double, acosh_double);
1641 CD_FUNCTION (genie_asinh_double, asinh_double);
1642 CD_FUNCTION (genie_atanh_double, atanh_double);
1643 CD_FUNCTION (genie_asin_double, asin_double);
1644 CD_FUNCTION (genie_atan_double, atan_double);
1645 CD_FUNCTION (genie_cosh_double, cosh_double);
1646 CD_FUNCTION (genie_cos_double, cos_double);
1647 CD_FUNCTION (genie_curt_double, cbrt_double);
1648 CD_FUNCTION (genie_exp_double, exp_double);
1649 CD_FUNCTION (genie_ln_double, log_double);
1650 CD_FUNCTION (genie_log_double, log10_double);
1651 CD_FUNCTION (genie_sinh_double, sinh_double);
1652 CD_FUNCTION (genie_sin_double, sin_double);
1653 CD_FUNCTION (genie_sqrt_double, sqrt_double);
1654 CD_FUNCTION (genie_tanh_double, tanh_double);
1655 CD_FUNCTION (genie_tan_double, tan_double);
1656 CD_FUNCTION (genie_erf_double, erf_double);
1657 CD_FUNCTION (genie_erfc_double, erfc_double);
1658 CD_FUNCTION (genie_lngamma_double, lgamma_double);
1659 CD_FUNCTION (genie_gamma_double, tgamma_double);
1660 CD_FUNCTION (genie_csc_double, a68_csc_double);
1661 CD_FUNCTION (genie_cscdg_double, a68_cscdg_double);
1662 CD_FUNCTION (genie_acsc_double, a68_acsc_double);
1663 CD_FUNCTION (genie_acscdg_double, a68_acscdg_double);
1664 CD_FUNCTION (genie_sec_double, a68_sec_double);
1665 CD_FUNCTION (genie_secdg_double, a68_secdg_double);
1666 CD_FUNCTION (genie_asec_double, a68_asec_double);
1667 CD_FUNCTION (genie_asecdg_double, a68_asecdg_double);
1668 CD_FUNCTION (genie_cot_double, a68_cot_double);
1669 CD_FUNCTION (genie_acot_double, a68_acot_double);
1670 CD_FUNCTION (genie_sindg_double, a68_sindg_double);
1671 CD_FUNCTION (genie_cas_double, a68_cas_double);
1672 CD_FUNCTION (genie_cosdg_double, a68_cosdg_double);
1673 CD_FUNCTION (genie_tandg_double, a68_tandg_double);
1674 CD_FUNCTION (genie_asindg_double, a68_asindg_double);
1675 CD_FUNCTION (genie_acosdg_double, a68_acosdg_double);
1676 CD_FUNCTION (genie_atandg_double, a68_atandg_double);
1677 CD_FUNCTION (genie_cotdg_double, a68_cotdg_double);
1678 CD_FUNCTION (genie_acotdg_double, a68_acotdg_double);
1679 CD_FUNCTION (genie_sinpi_double, a68_sinpi_double);
1680 CD_FUNCTION (genie_cospi_double, a68_cospi_double);
1681 CD_FUNCTION (genie_tanpi_double, a68_tanpi_double);
1682 CD_FUNCTION (genie_cotpi_double, a68_cotpi_double);
1683
1684 //! @brief PROC long arctan2 = (LONG REAL) LONG REAL
1685
1686 void genie_atan2_double (NODE_T * p)
1687 {
1688 A68_LONG_REAL x, y;
1689 POP_OBJECT (p, &y, A68_LONG_REAL);
1690 POP_OBJECT (p, &x, A68_LONG_REAL);
1691 errno = 0;
1692 PRELUDE_ERROR (VALUE (&x).f == 0.0q && VALUE (&y).f == 0.0q, p, ERROR_INVALID_ARGUMENT, M_LONG_REAL);
1693 VALUE (&x).f = a68_atan2_real (VALUE (&y).f, VALUE (&x).f);
1694 PRELUDE_ERROR (errno != 0, p, ERROR_MATH_EXCEPTION, NO_TEXT);
1695 PUSH_OBJECT (p, x, A68_LONG_REAL);
1696 }
1697
1698 //! @brief PROC long arctan2dg = (LONG REAL) LONG REAL
1699
1700 void genie_atan2dg_double (NODE_T * p)
1701 {
1702 A68_LONG_REAL x, y;
1703 POP_OBJECT (p, &y, A68_LONG_REAL);
1704 POP_OBJECT (p, &x, A68_LONG_REAL);
1705 errno = 0;
1706 PRELUDE_ERROR (VALUE (&x).f == 0.0q && VALUE (&y).f == 0.0q, p, ERROR_INVALID_ARGUMENT, M_LONG_REAL);
1707 VALUE (&x).f = CONST_180_OVER_PI_Q * a68_atan2_real (VALUE (&y).f, VALUE (&x).f);
1708 PRELUDE_ERROR (errno != 0, p, ERROR_MATH_EXCEPTION, NO_TEXT);
1709 PUSH_OBJECT (p, x, A68_LONG_REAL);
1710 }
1711
1712 //! @brief PROC (LONG REAL) LONG REAL inverf
1713
1714 void genie_inverf_double (NODE_T * _p_)
1715 {
1716 A68_LONG_REAL x;
1717 DOUBLE_T y, z;
1718 POP_OBJECT (_p_, &x, A68_LONG_REAL);
1719 errno = 0;
1720 y = VALUE (&x).f;
1721 z = inverf_double (y);
1722 MATH_RTE (_p_, errno != 0, M_LONG_REAL, NO_TEXT);
1723 CHECK_DOUBLE_REAL (_p_, z);
1724 PUSH_VALUE (_p_, dble (z), A68_LONG_REAL);
1725 }
1726
1727 //! @brief PROC (LONG REAL) LONG REAL inverfc
1728
1729 void genie_inverfc_double (NODE_T * p)
1730 {
1731 A68_LONG_REAL *u;
1732 POP_OPERAND_ADDRESS (p, u, A68_LONG_REAL);
1733 VALUE (u).f = 1.0q - (VALUE (u).f);
1734 genie_inverf_double (p);
1735 }
1736
1737 #define CD_C_FUNCTION(p, g)\
1738 A68_LONG_REAL re, im;\
1739 DOUBLE_COMPLEX_T z;\
1740 POP_OBJECT (p, &im, A68_LONG_REAL);\
1741 POP_OBJECT (p, &re, A68_LONG_REAL);\
1742 errno = 0;\
1743 z = VALUE (&re).f + VALUE (&im).f * _Complex_I;\
1744 z = g (z);\
1745 PUSH_VALUE (p, dble ((DOUBLE_T) creal_double (z)), A68_LONG_REAL);\
1746 PUSH_VALUE (p, dble ((DOUBLE_T) cimag_double (z)), A68_LONG_REAL);\
1747 MATH_RTE (p, errno != 0, M_COMPLEX, NO_TEXT);
1748
1749 //! @brief PROC long csqrt = (LONG COMPLEX) LONG COMPLEX
1750
1751 void genie_sqrt_double_compl (NODE_T * p)
1752 {
1753 CD_C_FUNCTION (p, csqrt_double);
1754 }
1755
1756 //! @brief PROC long csin = (LONG COMPLEX) LONG COMPLEX
1757
1758 void genie_sin_double_compl (NODE_T * p)
1759 {
1760 CD_C_FUNCTION (p, csin_double);
1761 }
1762
1763 //! @brief PROC long ccos = (LONG COMPLEX) LONG COMPLEX
1764
1765 void genie_cos_double_compl (NODE_T * p)
1766 {
1767 CD_C_FUNCTION (p, ccos_double);
1768 }
1769
1770 //! @brief PROC long ctan = (LONG COMPLEX) LONG COMPLEX
1771
1772 void genie_tan_double_compl (NODE_T * p)
1773 {
1774 CD_C_FUNCTION (p, ctan_double);
1775 }
1776
1777 //! @brief PROC long casin = (LONG COMPLEX) LONG COMPLEX
1778
1779 void genie_asin_double_compl (NODE_T * p)
1780 {
1781 CD_C_FUNCTION (p, casin_double);
1782 }
1783
1784 //! @brief PROC long cacos = (LONG COMPLEX) LONG COMPLEX
1785
1786 void genie_acos_double_compl (NODE_T * p)
1787 {
1788 CD_C_FUNCTION (p, cacos_double);
1789 }
1790
1791 //! @brief PROC long catan = (LONG COMPLEX) LONG COMPLEX
1792
1793 void genie_atan_double_compl (NODE_T * p)
1794 {
1795 CD_C_FUNCTION (p, catan_double);
1796 }
1797
1798 //! @brief PROC long cexp = (LONG COMPLEX) LONG COMPLEX
1799
1800 void genie_exp_double_compl (NODE_T * p)
1801 {
1802 CD_C_FUNCTION (p, cexp_double);
1803 }
1804
1805 //! @brief PROC long cln = (LONG COMPLEX) LONG COMPLEX
1806
1807 void genie_ln_double_compl (NODE_T * p)
1808 {
1809 CD_C_FUNCTION (p, clog_double);
1810 }
1811
1812 //! @brief PROC long csinh = (LONG COMPLEX) LONG COMPLEX
1813
1814 void genie_sinh_double_compl (NODE_T * p)
1815 {
1816 CD_C_FUNCTION (p, csinh_double);
1817 }
1818
1819 //! @brief PROC long ccosh = (LONG COMPLEX) LONG COMPLEX
1820
1821 void genie_cosh_double_compl (NODE_T * p)
1822 {
1823 CD_C_FUNCTION (p, ccosh_double);
1824 }
1825
1826 //! @brief PROC long ctanh = (LONG COMPLEX) LONG COMPLEX
1827
1828 void genie_tanh_double_compl (NODE_T * p)
1829 {
1830 CD_C_FUNCTION (p, ctanh_double);
1831 }
1832
1833 //! @brief PROC long casinh = (LONG COMPLEX) LONG COMPLEX
1834
1835 void genie_asinh_double_compl (NODE_T * p)
1836 {
1837 CD_C_FUNCTION (p, casinh_double);
1838 }
1839
1840 //! @brief PROC long cacosh = (LONG COMPLEX) LONG COMPLEX
1841
1842 void genie_acosh_double_compl (NODE_T * p)
1843 {
1844 CD_C_FUNCTION (p, cacosh_double);
1845 }
1846
1847 //! @brief PROC long catanh = (LONG COMPLEX) LONG COMPLEX
1848
1849 void genie_atanh_double_compl (NODE_T * p)
1850 {
1851 CD_C_FUNCTION (p, catanh_double);
1852 }
1853
1854 //! @brief PROC next long random = LONG REAL
1855
1856 void genie_next_random_double (NODE_T * p)
1857 {
1858 // This is 'real width' digits only.
1859 genie_next_random (p);
1860 genie_lengthen_real_to_double (p);
1861 }
1862
1863 #define CALL(g, x, y) {\
1864 ADDR_T pop_sp = A68_SP;\
1865 A68_LONG_REAL *z = (A68_LONG_REAL *) STACK_TOP;\
1866 DOUBLE_NUM_T _w_;\
1867 _w_.f = (x);\
1868 PUSH_VALUE (_p_, _w_, A68_LONG_REAL);\
1869 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);\
1870 (y) = VALUE (z).f;\
1871 A68_SP = pop_sp;\
1872 }
1873
1874 //! @brief Transform string into real-16.
1875
1876 DOUBLE_T string_to_double (char *s, char **end)
1877 {
1878 errno = 0;
1879 DOUBLE_T y[A68_DOUBLE_DIG];
1880 for (int i = 0; i < A68_DOUBLE_DIG; i++) {
1881 y[i] = 0.0q;
1882 }
1883 while (IS_SPACE (s[0])) {
1884 s++;
1885 }
1886 // Scan mantissa digits and put them into "y".
1887 DOUBLE_T W;
1888 if (s[0] == '-') {
1889 W = -1.0q;
1890 } else {
1891 W = 1.0q;
1892 }
1893 if (s[0] == '+' || s[0] == '-') {
1894 s++;
1895 }
1896 while (s[0] == '0') {
1897 s++;
1898 }
1899 int dot = -1, pos = 0, pow = 0;
1900 while (pow < A68_DOUBLE_DIG && s[pos] != NULL_CHAR && (IS_DIGIT (s[pos]) || s[pos] == POINT_CHAR)) {
1901 if (s[pos] == POINT_CHAR) {
1902 dot = pos;
1903 } else {
1904 int val = (int) s[pos] - (int) '0';
1905 y[pow] = W * val;
1906 W /= 10.0q;
1907 pow++;
1908 }
1909 pos++;
1910 }
1911 (*end) = &(s[pos]);
1912 // Sum from low to high to preserve precision.
1913 DOUBLE_T sum = 0.0q;
1914 for (int i = A68_DOUBLE_DIG - 1; i >= 0; i--) {
1915 sum = sum + y[i];
1916 }
1917 // See if there is an exponent.
1918 int expo;
1919 if (s[pos] != NULL_CHAR && TO_UPPER (s[pos]) == TO_UPPER (EXPONENT_CHAR)) {
1920 expo = (int) strtol (&(s[++pos]), end, 10);
1921 } else {
1922 expo = 0;
1923 }
1924 // Standardise.
1925 if (dot >= 0) {
1926 expo += dot - 1;
1927 } else {
1928 expo += pow - 1;
1929 }
1930 while (sum != 0.0q && fabs_double (sum) < 1.0q) {
1931 sum *= 10.0q;
1932 expo -= 1;
1933 }
1934 if (errno == 0) {
1935 return sum * ten_up_double (expo);
1936 } else {
1937 return 0.0q;
1938 }
1939 }
1940
1941 void genie_beta_inc_cf_double (NODE_T * p)
1942 {
1943 A68_LONG_REAL x, s, t;
1944 POP_OBJECT (p, &x, A68_LONG_REAL);
1945 POP_OBJECT (p, &t, A68_LONG_REAL);
1946 POP_OBJECT (p, &s, A68_LONG_REAL);
1947 errno = 0;
1948 PUSH_VALUE (p, dble (a68_beta_inc_double (VALUE (&s).f, VALUE (&t).f, VALUE (&x).f)), A68_LONG_REAL);
1949 MATH_RTE (p, errno != 0, M_LONG_REAL, NO_TEXT);
1950 }
1951
1952 void genie_beta_double (NODE_T * p)
1953 {
1954 A68_LONG_REAL a, b;
1955 POP_OBJECT (p, &b, A68_LONG_REAL);
1956 POP_OBJECT (p, &a, A68_LONG_REAL);
1957 errno = 0;
1958 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);
1959 MATH_RTE (p, errno != 0, M_LONG_REAL, NO_TEXT);
1960 }
1961
1962 void genie_ln_beta_double (NODE_T * p)
1963 {
1964 A68_LONG_REAL a, b;
1965 POP_OBJECT (p, &b, A68_LONG_REAL);
1966 POP_OBJECT (p, &a, A68_LONG_REAL);
1967 errno = 0;
1968 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);
1969 MATH_RTE (p, errno != 0, M_LONG_REAL, NO_TEXT);
1970 }
1971
1972 // LONG REAL infinity
1973
1974 void genie_infinity_double (NODE_T * p)
1975 {
1976 PUSH_VALUE (p, dble (a68_posinf_real ()), A68_LONG_REAL);
1977 }
1978
1979 // LONG REAL minus infinity
1980
1981 void genie_minus_infinity_double (NODE_T * p)
1982 {
1983 PUSH_VALUE (p, dble (a68_neginf_double ()), A68_LONG_REAL);
1984 }
1985
1986 #endif
