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-2023 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 k -= A68_BITS_WIDTH;
1217 }
1218 for (int n = 0; n < (A68_BITS_WIDTH - VALUE (&i)); n++) {
1219 mask = mask << 1;
1220 }
1221 PUSH_VALUE (p, (BOOL_T) ((*w & mask) ? A68_TRUE : A68_FALSE), A68_BOOL);
1222 }
1223
1224 //! @brief OP SET = (INT, LONG BITS) LONG BITS
1225
1226 void genie_set_double_bits (NODE_T * p)
1227 {
1228 A68_LONG_BITS j; A68_INT i;
1229 POP_OBJECT (p, &j, A68_LONG_BITS);
1230 POP_OBJECT (p, &i, A68_INT);
1231 int k = VALUE (&i);
1232 PRELUDE_ERROR (k < 1 || k > A68_LONG_BITS_WIDTH, p, ERROR_OUT_OF_BOUNDS, M_INT);
1233 UNSIGNED_T mask = 0x1, *w;
1234 if (k <= A68_BITS_WIDTH) {
1235 w = &(LW (VALUE (&j)));
1236 } else {
1237 w = &(HW (VALUE (&j)));
1238 k -= A68_BITS_WIDTH;
1239 }
1240 for (int n = 0; n < (A68_BITS_WIDTH - VALUE (&i)); n++) {
1241 mask = mask << 1;
1242 }
1243 (*w) |= mask;
1244 PUSH_OBJECT (p, j, A68_LONG_BITS);
1245 }
1246
1247 //! @brief OP CLEAR = (INT, LONG BITS) LONG BITS
1248
1249 void genie_clear_double_bits (NODE_T * p)
1250 {
1251 A68_LONG_BITS j; A68_INT i;
1252 POP_OBJECT (p, &j, A68_LONG_BITS);
1253 POP_OBJECT (p, &i, A68_INT);
1254 int k = VALUE (&i);
1255 PRELUDE_ERROR (k < 1 || k > A68_LONG_BITS_WIDTH, p, ERROR_OUT_OF_BOUNDS, M_INT);
1256 UNSIGNED_T mask = 0x1, *w;
1257 if (k <= A68_BITS_WIDTH) {
1258 w = &(LW (VALUE (&j)));
1259 } else {
1260 w = &(HW (VALUE (&j)));
1261 k -= A68_BITS_WIDTH;
1262 }
1263 for (int n = 0; n < (A68_BITS_WIDTH - VALUE (&i)); n++) {
1264 mask = mask << 1;
1265 }
1266 (*w) &= ~mask;
1267 PUSH_OBJECT (p, j, A68_LONG_BITS);
1268 }
1269
1270 //! @brief OP SHL = (LONG BITS, INT) LONG BITS
1271
1272 void genie_shl_double_bits (NODE_T * p)
1273 {
1274 A68_LONG_BITS i; A68_INT j;
1275 POP_OBJECT (p, &j, A68_INT);
1276 POP_OBJECT (p, &i, A68_LONG_BITS);
1277 DOUBLE_NUM_T *w = &VALUE (&i);
1278 int k = VALUE (&j);
1279 if (VALUE (&j) >= 0) {
1280 for (int n = 0; n < k; n++) {
1281 UNSIGNED_T carry = ((LW (*w) & D_SIGN) ? 0x1 : 0x0);
1282 PRELUDE_ERROR (MODCHK (p, M_LONG_BITS, HW (*w) | D_SIGN), p, ERROR_MATH, M_LONG_BITS);
1283 HW (*w) = (HW (*w) << 1) | carry;
1284 LW (*w) = (LW (*w) << 1);
1285 }
1286 } else {
1287 k = -k;
1288 for (int n = 0; n < k; n++) {
1289 UNSIGNED_T carry = ((HW (*w) & 0x1) ? D_SIGN : 0x0);
1290 HW (*w) = (HW (*w) >> 1);
1291 LW (*w) = (LW (*w) >> 1) | carry;
1292 }
1293 }
1294 PUSH_OBJECT (p, i, A68_LONG_BITS);
1295 }
1296
1297 //! @brief OP SHR = (LONG BITS, INT) LONG BITS
1298
1299 void genie_shr_double_bits (NODE_T * p)
1300 {
1301 A68_INT *j;
1302 POP_OPERAND_ADDRESS (p, j, A68_INT);
1303 VALUE (j) = -VALUE (j);
1304 genie_shl_double_bits (p); // Conform RR
1305 }
1306
1307 //! @brief OP ROL = (LONG BITS, INT) LONG BITS
1308
1309 void genie_rol_double_bits (NODE_T * p)
1310 {
1311 A68_LONG_BITS i; A68_INT j;
1312 DOUBLE_NUM_T *w = &VALUE (&i);
1313 POP_OBJECT (p, &j, A68_INT);
1314 POP_OBJECT (p, &i, A68_LONG_BITS);
1315 int k = VALUE (&j);
1316 if (k >= 0) {
1317 for (int n = 0; n < k; n++) {
1318 UNSIGNED_T carry = ((HW (*w) & D_SIGN) ? 0x1 : 0x0);
1319 UNSIGNED_T carry_between = ((LW (*w) & D_SIGN) ? 0x1 : 0x0);
1320 HW (*w) = (HW (*w) << 1) | carry_between;
1321 LW (*w) = (LW (*w) << 1) | carry;
1322 }
1323 } else {
1324 k = -k;
1325 for (int n = 0; n < k; n++) {
1326 UNSIGNED_T carry = ((LW (*w) & 0x1) ? D_SIGN : 0x0);
1327 UNSIGNED_T carry_between = ((HW (*w) & 0x1) ? D_SIGN : 0x0);
1328 HW (*w) = (HW (*w) >> 1) | carry;
1329 LW (*w) = (LW (*w) >> 1) | carry_between;
1330 }
1331 }
1332 PUSH_OBJECT (p, i, A68_LONG_BITS);
1333 }
1334
1335 //! @brief OP ROR = (LONG BITS, INT) LONG BITS
1336
1337 void genie_ror_double_bits (NODE_T * p)
1338 {
1339 A68_INT *j;
1340 POP_OPERAND_ADDRESS (p, j, A68_INT);
1341 VALUE (j) = -VALUE (j);
1342 genie_rol_double_bits (p); // Conform RR
1343 }
1344
1345 //! @brief OP BIN = (LONG INT) LONG BITS
1346
1347 void genie_bin_double_int (NODE_T * p)
1348 {
1349 A68_LONG_INT i;
1350 POP_OBJECT (p, &i, A68_LONG_INT);
1351 // RR does not convert negative numbers
1352 if (D_NEG (VALUE (&i))) {
1353 errno = EDOM;
1354 diagnostic (A68_RUNTIME_ERROR, p, ERROR_OUT_OF_BOUNDS, M_BITS);
1355 exit_genie (p, A68_RUNTIME_ERROR);
1356 }
1357 PUSH_OBJECT (p, i, A68_LONG_BITS);
1358 }
1359
1360 //! @brief OP +* = (LONG REAL, LONG REAL) LONG COMPLEX
1361
1362 void genie_i_double_compl (NODE_T * p)
1363 {
1364 (void) p;
1365 }
1366
1367 //! @brief OP SHORTEN = (LONG COMPLEX) COMPLEX
1368
1369 void genie_shorten_double_compl_to_complex (NODE_T * p)
1370 {
1371 A68_LONG_REAL re, im;
1372 POP_OBJECT (p, &im, A68_LONG_REAL);
1373 POP_OBJECT (p, &re, A68_LONG_REAL);
1374 REAL_T w = VALUE (&re).f;
1375 PUSH_VALUE (p, w, A68_REAL);
1376 w = VALUE (&im).f;
1377 PUSH_VALUE (p, w, A68_REAL);
1378 }
1379
1380 //! @brief OP LENG = (LONG COMPLEX) LONG LONG COMPLEX
1381
1382 void genie_lengthen_double_compl_to_long_mp_complex (NODE_T * p)
1383 {
1384 int digs = DIGITS (M_LONG_LONG_REAL);
1385 A68_LONG_REAL re, im;
1386 POP_OBJECT (p, &im, A68_LONG_REAL);
1387 POP_OBJECT (p, &re, A68_LONG_REAL);
1388 MP_T *z = nil_mp (p, digs);
1389 (void) double_to_mp (p, z, VALUE (&re).f, digs);
1390 MP_STATUS (z) = (MP_T) INIT_MASK;
1391 z = nil_mp (p, digs);
1392 (void) double_to_mp (p, z, VALUE (&im).f, digs);
1393 MP_STATUS (z) = (MP_T) INIT_MASK;
1394 }
1395
1396 //! @brief OP +* = (LONG INT, LONG INT) LONG COMPLEX
1397
1398 void genie_i_int_double_compl (NODE_T * p)
1399 {
1400 A68_LONG_INT re, im;
1401 POP_OBJECT (p, &im, A68_LONG_INT);
1402 POP_OBJECT (p, &re, A68_LONG_INT);
1403 PUSH_VALUE (p, double_int_to_double (p, VALUE (&re)), A68_LONG_REAL);
1404 PUSH_VALUE (p, double_int_to_double (p, VALUE (&im)), A68_LONG_REAL);
1405 }
1406
1407 //! @brief OP RE = (LONG COMPLEX) LONG REAL
1408
1409 void genie_re_double_compl (NODE_T * p)
1410 {
1411 DECREMENT_STACK_POINTER (p, SIZE (M_LONG_REAL));
1412 }
1413
1414 //! @brief OP IM = (LONG COMPLEX) LONG REAL
1415
1416 void genie_im_double_compl (NODE_T * p)
1417 {
1418 A68_LONG_REAL re, im;
1419 POP_OBJECT (p, &im, A68_LONG_REAL);
1420 POP_OBJECT (p, &re, A68_LONG_REAL);
1421 PUSH_OBJECT (p, im, A68_LONG_REAL);
1422 }
1423
1424 //! @brief OP - = (LONG COMPLEX) LONG COMPLEX
1425
1426 void genie_minus_double_compl (NODE_T * p)
1427 {
1428 A68_LONG_REAL re, im;
1429 POP_OBJECT (p, &im, A68_LONG_REAL);
1430 POP_OBJECT (p, &re, A68_LONG_REAL);
1431 VALUE (&re).f = -VALUE (&re).f;
1432 VALUE (&im).f = -VALUE (&im).f;
1433 PUSH_OBJECT (p, im, A68_LONG_REAL);
1434 PUSH_OBJECT (p, re, A68_LONG_REAL);
1435 }
1436
1437 //! @brief OP ABS = (LONG COMPLEX) LONG REAL
1438
1439 void genie_abs_double_compl (NODE_T * p)
1440 {
1441 A68_LONG_REAL re, im;
1442 POP_LONG_COMPLEX (p, &re, &im);
1443 PUSH_VALUE (p, dble (a68_hypot_double (VALUE (&re).f, VALUE (&im).f)), A68_LONG_REAL);
1444 }
1445
1446 //! @brief OP ARG = (LONG COMPLEX) LONG REAL
1447
1448 void genie_arg_double_compl (NODE_T * p)
1449 {
1450 A68_LONG_REAL re, im;
1451 POP_LONG_COMPLEX (p, &re, &im);
1452 PRELUDE_ERROR (VALUE (&re).f == 0.0q && VALUE (&im).f == 0.0q, p, ERROR_INVALID_ARGUMENT, M_LONG_COMPLEX);
1453 PUSH_VALUE (p, dble (atan2_double (VALUE (&im).f, VALUE (&re).f)), A68_LONG_REAL);
1454 }
1455
1456 //! @brief OP CONJ = (LONG COMPLEX) LONG COMPLEX
1457
1458 void genie_conj_double_compl (NODE_T * p)
1459 {
1460 A68_LONG_REAL im;
1461 POP_OBJECT (p, &im, A68_LONG_REAL);
1462 VALUE (&im).f = -VALUE (&im).f;
1463 PUSH_OBJECT (p, im, A68_LONG_REAL);
1464 }
1465
1466 //! @brief OP + = (COMPLEX, COMPLEX) COMPLEX
1467
1468 void genie_add_double_compl (NODE_T * p)
1469 {
1470 A68_LONG_REAL re_x, im_x, re_y, im_y;
1471 POP_LONG_COMPLEX (p, &re_y, &im_y);
1472 POP_LONG_COMPLEX (p, &re_x, &im_x);
1473 VALUE (&re_x).f += VALUE (&re_y).f;
1474 VALUE (&im_x).f += VALUE (&im_y).f;
1475 CHECK_DOUBLE_COMPLEX (p, VALUE (&im_x).f, VALUE (&im_y).f);
1476 PUSH_OBJECT (p, re_x, A68_LONG_REAL);
1477 PUSH_OBJECT (p, im_x, A68_LONG_REAL);
1478 }
1479
1480 //! @brief OP - = (COMPLEX, COMPLEX) COMPLEX
1481
1482 void genie_sub_double_compl (NODE_T * p)
1483 {
1484 A68_LONG_REAL re_x, im_x, re_y, im_y;
1485 POP_LONG_COMPLEX (p, &re_y, &im_y);
1486 POP_LONG_COMPLEX (p, &re_x, &im_x);
1487 VALUE (&re_x).f -= VALUE (&re_y).f;
1488 VALUE (&im_x).f -= VALUE (&im_y).f;
1489 CHECK_DOUBLE_COMPLEX (p, VALUE (&im_x).f, VALUE (&im_y).f);
1490 PUSH_OBJECT (p, re_x, A68_LONG_REAL);
1491 PUSH_OBJECT (p, im_x, A68_LONG_REAL);
1492 }
1493
1494 //! @brief OP * = (COMPLEX, COMPLEX) COMPLEX
1495
1496 void genie_mul_double_compl (NODE_T * p)
1497 {
1498 A68_LONG_REAL re_x, im_x, re_y, im_y;
1499 POP_LONG_COMPLEX (p, &re_y, &im_y);
1500 POP_LONG_COMPLEX (p, &re_x, &im_x);
1501 DOUBLE_T re = VALUE (&re_x).f * VALUE (&re_y).f - VALUE (&im_x).f * VALUE (&im_y).f;
1502 DOUBLE_T im = VALUE (&im_x).f * VALUE (&re_y).f + VALUE (&re_x).f * VALUE (&im_y).f;
1503 CHECK_DOUBLE_COMPLEX (p, VALUE (&im_x).f, VALUE (&im_y).f);
1504 PUSH_VALUE (p, dble (re), A68_LONG_REAL);
1505 PUSH_VALUE (p, dble (im), A68_LONG_REAL);
1506 }
1507
1508 //! @brief OP / = (COMPLEX, COMPLEX) COMPLEX
1509
1510 void genie_div_double_compl (NODE_T * p)
1511 {
1512 A68_LONG_REAL re_x, im_x, re_y, im_y;
1513 DOUBLE_T re = 0.0, im = 0.0;
1514 POP_LONG_COMPLEX (p, &re_y, &im_y);
1515 POP_LONG_COMPLEX (p, &re_x, &im_x);
1516 PRELUDE_ERROR (VALUE (&re_y).f == 0.0q && VALUE (&im_y).f == 0.0q, p, ERROR_DIVISION_BY_ZERO, M_LONG_COMPLEX);
1517 if (ABSQ (VALUE (&re_y).f) >= ABSQ (VALUE (&im_y).f)) {
1518 DOUBLE_T r = VALUE (&im_y).f / VALUE (&re_y).f, den = VALUE (&re_y).f + r * VALUE (&im_y).f;
1519 re = (VALUE (&re_x).f + r * VALUE (&im_x).f) / den;
1520 im = (VALUE (&im_x).f - r * VALUE (&re_x).f) / den;
1521 } else {
1522 DOUBLE_T r = VALUE (&re_y).f / VALUE (&im_y).f, den = VALUE (&im_y).f + r * VALUE (&re_y).f;
1523 re = (VALUE (&re_x).f * r + VALUE (&im_x).f) / den;
1524 im = (VALUE (&im_x).f * r - VALUE (&re_x).f) / den;
1525 }
1526 PUSH_VALUE (p, dble (re), A68_LONG_REAL);
1527 PUSH_VALUE (p, dble (im), A68_LONG_REAL);
1528 }
1529
1530 //! @brief OP ** = (LONG COMPLEX, INT) LONG COMPLEX
1531
1532 void genie_pow_double_compl_int (NODE_T * p)
1533 {
1534 A68_LONG_REAL re_x, im_x;
1535 A68_INT j;
1536 BOOL_T negative;
1537 POP_OBJECT (p, &j, A68_INT);
1538 POP_LONG_COMPLEX (p, &re_x, &im_x);
1539 DOUBLE_T re_z = 1.0q, im_z = 0.0q;
1540 DOUBLE_T re_y = VALUE (&re_x).f, im_y = VALUE (&im_x).f;
1541 negative = (BOOL_T) (VALUE (&j) < 0);
1542 if (negative) {
1543 VALUE (&j) = -VALUE (&j);
1544 }
1545 INT_T expo = 1;
1546 while ((UNSIGNED_T) expo <= (UNSIGNED_T) (VALUE (&j))) {
1547 DOUBLE_T z;
1548 if (expo & VALUE (&j)) {
1549 z = re_z * re_y - im_z * im_y;
1550 im_z = re_z * im_y + im_z * re_y;
1551 re_z = z;
1552 }
1553 z = re_y * re_y - im_y * im_y;
1554 im_y = im_y * re_y + re_y * im_y;
1555 re_y = z;
1556 CHECK_DOUBLE_COMPLEX (p, re_y, im_y);
1557 CHECK_DOUBLE_COMPLEX (p, re_z, im_z);
1558 expo <<= 1;
1559 }
1560 if (negative) {
1561 PUSH_VALUE (p, dble (1.0q), A68_LONG_REAL);
1562 PUSH_VALUE (p, dble (0.0q), A68_LONG_REAL);
1563 PUSH_VALUE (p, dble (re_z), A68_LONG_REAL);
1564 PUSH_VALUE (p, dble (im_z), A68_LONG_REAL);
1565 genie_div_double_compl (p);
1566 } else {
1567 PUSH_VALUE (p, dble (re_z), A68_LONG_REAL);
1568 PUSH_VALUE (p, dble (im_z), A68_LONG_REAL);
1569 }
1570 }
1571
1572 //! @brief OP = = (COMPLEX, COMPLEX) BOOL
1573
1574 void genie_eq_double_compl (NODE_T * p)
1575 {
1576 A68_LONG_REAL re_x, im_x, re_y, im_y;
1577 POP_LONG_COMPLEX (p, &re_y, &im_y);
1578 POP_LONG_COMPLEX (p, &re_x, &im_x);
1579 PUSH_VALUE (p, (BOOL_T) ((VALUE (&re_x).f == VALUE (&re_y).f) && (VALUE (&im_x).f == VALUE (&im_y).f)), A68_BOOL);
1580 }
1581
1582 //! @brief OP /= = (COMPLEX, COMPLEX) BOOL
1583
1584 void genie_ne_double_compl (NODE_T * p)
1585 {
1586 A68_LONG_REAL re_x, im_x, re_y, im_y;
1587 POP_LONG_COMPLEX (p, &re_y, &im_y);
1588 POP_LONG_COMPLEX (p, &re_x, &im_x);
1589 PUSH_VALUE (p, (BOOL_T) ! ((VALUE (&re_x).f == VALUE (&re_y).f) && (VALUE (&im_x).f == VALUE (&im_y).f)), A68_BOOL);
1590 }
1591
1592 //! @brief OP +:= = (REF COMPLEX, COMPLEX) REF COMPLEX
1593
1594 void genie_plusab_double_compl (NODE_T * p)
1595 {
1596 genie_f_and_becomes (p, M_REF_LONG_COMPLEX, genie_add_double_compl);
1597 }
1598
1599 //! @brief OP -:= = (REF COMPLEX, COMPLEX) REF COMPLEX
1600
1601 void genie_minusab_double_compl (NODE_T * p)
1602 {
1603 genie_f_and_becomes (p, M_REF_LONG_COMPLEX, genie_sub_double_compl);
1604 }
1605
1606 //! @brief OP *:= = (REF COMPLEX, COMPLEX) REF COMPLEX
1607
1608 void genie_timesab_double_compl (NODE_T * p)
1609 {
1610 genie_f_and_becomes (p, M_REF_LONG_COMPLEX, genie_mul_double_compl);
1611 }
1612
1613 //! @brief OP /:= = (REF COMPLEX, COMPLEX) REF COMPLEX
1614
1615 void genie_divab_double_compl (NODE_T * p)
1616 {
1617 genie_f_and_becomes (p, M_REF_LONG_COMPLEX, genie_div_double_compl);
1618 }
1619
1620 //! @brief OP LENG = (COMPLEX) LONG COMPLEX
1621
1622 void genie_lengthen_complex_to_double_compl (NODE_T * p)
1623 {
1624 A68_REAL i;
1625 POP_OBJECT (p, &i, A68_REAL);
1626 genie_lengthen_real_to_double (p);
1627 PUSH_OBJECT (p, i, A68_REAL);
1628 genie_lengthen_real_to_double (p);
1629 }
1630
1631 // Functions
1632
1633 #define CD_FUNCTION(name, fun)\
1634 void name (NODE_T * p) {\
1635 A68_LONG_REAL *x;\
1636 POP_OPERAND_ADDRESS (p, x, A68_LONG_REAL);\
1637 errno = 0;\
1638 VALUE (x).f = fun (VALUE (x).f);\
1639 MATH_RTE (p, errno != 0, M_LONG_REAL, NO_TEXT);\
1640 }
1641
1642 CD_FUNCTION (genie_acos_double, acos_double);
1643 CD_FUNCTION (genie_acosh_double, acosh_double);
1644 CD_FUNCTION (genie_asinh_double, asinh_double);
1645 CD_FUNCTION (genie_atanh_double, atanh_double);
1646 CD_FUNCTION (genie_asin_double, asin_double);
1647 CD_FUNCTION (genie_atan_double, atan_double);
1648 CD_FUNCTION (genie_cosh_double, cosh_double);
1649 CD_FUNCTION (genie_cos_double, cos_double);
1650 CD_FUNCTION (genie_curt_double, cbrt_double);
1651 CD_FUNCTION (genie_exp_double, exp_double);
1652 CD_FUNCTION (genie_ln_double, log_double);
1653 CD_FUNCTION (genie_log_double, log10_double);
1654 CD_FUNCTION (genie_sinh_double, sinh_double);
1655 CD_FUNCTION (genie_sin_double, sin_double);
1656 CD_FUNCTION (genie_sqrt_double, sqrt_double);
1657 CD_FUNCTION (genie_tanh_double, tanh_double);
1658 CD_FUNCTION (genie_tan_double, tan_double);
1659 CD_FUNCTION (genie_erf_double, erf_double);
1660 CD_FUNCTION (genie_erfc_double, erfc_double);
1661 CD_FUNCTION (genie_lngamma_double, lgamma_double);
1662 CD_FUNCTION (genie_gamma_double, tgamma_double);
1663 CD_FUNCTION (genie_csc_double, a68_csc_double);
1664 CD_FUNCTION (genie_cscdg_double, a68_cscdg_double);
1665 CD_FUNCTION (genie_acsc_double, a68_acsc_double);
1666 CD_FUNCTION (genie_acscdg_double, a68_acscdg_double);
1667 CD_FUNCTION (genie_sec_double, a68_sec_double);
1668 CD_FUNCTION (genie_secdg_double, a68_secdg_double);
1669 CD_FUNCTION (genie_asec_double, a68_asec_double);
1670 CD_FUNCTION (genie_asecdg_double, a68_asecdg_double);
1671 CD_FUNCTION (genie_cot_double, a68_cot_double);
1672 CD_FUNCTION (genie_acot_double, a68_acot_double);
1673 CD_FUNCTION (genie_sindg_double, a68_sindg_double);
1674 CD_FUNCTION (genie_cas_double, a68_cas_double);
1675 CD_FUNCTION (genie_cosdg_double, a68_cosdg_double);
1676 CD_FUNCTION (genie_tandg_double, a68_tandg_double);
1677 CD_FUNCTION (genie_asindg_double, a68_asindg_double);
1678 CD_FUNCTION (genie_acosdg_double, a68_acosdg_double);
1679 CD_FUNCTION (genie_atandg_double, a68_atandg_double);
1680 CD_FUNCTION (genie_cotdg_double, a68_cotdg_double);
1681 CD_FUNCTION (genie_acotdg_double, a68_acotdg_double);
1682 CD_FUNCTION (genie_sinpi_double, a68_sinpi_double);
1683 CD_FUNCTION (genie_cospi_double, a68_cospi_double);
1684 CD_FUNCTION (genie_tanpi_double, a68_tanpi_double);
1685 CD_FUNCTION (genie_cotpi_double, a68_cotpi_double);
1686
1687 //! @brief PROC long arctan2 = (LONG REAL) LONG REAL
1688
1689 void genie_atan2_double (NODE_T * p)
1690 {
1691 A68_LONG_REAL x, y;
1692 POP_OBJECT (p, &y, A68_LONG_REAL);
1693 POP_OBJECT (p, &x, A68_LONG_REAL);
1694 errno = 0;
1695 PRELUDE_ERROR (VALUE (&x).f == 0.0q && VALUE (&y).f == 0.0q, p, ERROR_INVALID_ARGUMENT, M_LONG_REAL);
1696 VALUE (&x).f = a68_atan2_real (VALUE (&y).f, VALUE (&x).f);
1697 PRELUDE_ERROR (errno != 0, p, ERROR_MATH_EXCEPTION, NO_TEXT);
1698 PUSH_OBJECT (p, x, A68_LONG_REAL);
1699 }
1700
1701 //! @brief PROC long arctan2dg = (LONG REAL) LONG REAL
1702
1703 void genie_atan2dg_double (NODE_T * p)
1704 {
1705 A68_LONG_REAL x, y;
1706 POP_OBJECT (p, &y, A68_LONG_REAL);
1707 POP_OBJECT (p, &x, A68_LONG_REAL);
1708 errno = 0;
1709 PRELUDE_ERROR (VALUE (&x).f == 0.0q && VALUE (&y).f == 0.0q, p, ERROR_INVALID_ARGUMENT, M_LONG_REAL);
1710 VALUE (&x).f = CONST_180_OVER_PI_Q * a68_atan2_real (VALUE (&y).f, VALUE (&x).f);
1711 PRELUDE_ERROR (errno != 0, p, ERROR_MATH_EXCEPTION, NO_TEXT);
1712 PUSH_OBJECT (p, x, A68_LONG_REAL);
1713 }
1714
1715 //! @brief PROC (LONG REAL) LONG REAL inverf
1716
1717 void genie_inverf_double (NODE_T * _p_)
1718 {
1719 A68_LONG_REAL x;
1720 DOUBLE_T y, z;
1721 POP_OBJECT (_p_, &x, A68_LONG_REAL);
1722 errno = 0;
1723 y = VALUE (&x).f;
1724 z = inverf_double (y);
1725 MATH_RTE (_p_, errno != 0, M_LONG_REAL, NO_TEXT);
1726 CHECK_DOUBLE_REAL (_p_, z);
1727 PUSH_VALUE (_p_, dble (z), A68_LONG_REAL);
1728 }
1729
1730 //! @brief PROC (LONG REAL) LONG REAL inverfc
1731
1732 void genie_inverfc_double (NODE_T * p)
1733 {
1734 A68_LONG_REAL *u;
1735 POP_OPERAND_ADDRESS (p, u, A68_LONG_REAL);
1736 VALUE (u).f = 1.0q - (VALUE (u).f);
1737 genie_inverf_double (p);
1738 }
1739
1740 #define CD_C_FUNCTION(p, g)\
1741 A68_LONG_REAL re, im;\
1742 DOUBLE_COMPLEX_T z;\
1743 POP_OBJECT (p, &im, A68_LONG_REAL);\
1744 POP_OBJECT (p, &re, A68_LONG_REAL);\
1745 errno = 0;\
1746 z = VALUE (&re).f + VALUE (&im).f * _Complex_I;\
1747 z = g (z);\
1748 PUSH_VALUE (p, dble ((DOUBLE_T) creal_double (z)), A68_LONG_REAL);\
1749 PUSH_VALUE (p, dble ((DOUBLE_T) cimag_double (z)), A68_LONG_REAL);\
1750 MATH_RTE (p, errno != 0, M_COMPLEX, NO_TEXT);
1751
1752 //! @brief PROC long csqrt = (LONG COMPLEX) LONG COMPLEX
1753
1754 void genie_sqrt_double_compl (NODE_T * p)
1755 {
1756 CD_C_FUNCTION (p, csqrt_double);
1757 }
1758
1759 //! @brief PROC long csin = (LONG COMPLEX) LONG COMPLEX
1760
1761 void genie_sin_double_compl (NODE_T * p)
1762 {
1763 CD_C_FUNCTION (p, csin_double);
1764 }
1765
1766 //! @brief PROC long ccos = (LONG COMPLEX) LONG COMPLEX
1767
1768 void genie_cos_double_compl (NODE_T * p)
1769 {
1770 CD_C_FUNCTION (p, ccos_double);
1771 }
1772
1773 //! @brief PROC long ctan = (LONG COMPLEX) LONG COMPLEX
1774
1775 void genie_tan_double_compl (NODE_T * p)
1776 {
1777 CD_C_FUNCTION (p, ctan_double);
1778 }
1779
1780 //! @brief PROC long casin = (LONG COMPLEX) LONG COMPLEX
1781
1782 void genie_asin_double_compl (NODE_T * p)
1783 {
1784 CD_C_FUNCTION (p, casin_double);
1785 }
1786
1787 //! @brief PROC long cacos = (LONG COMPLEX) LONG COMPLEX
1788
1789 void genie_acos_double_compl (NODE_T * p)
1790 {
1791 CD_C_FUNCTION (p, cacos_double);
1792 }
1793
1794 //! @brief PROC long catan = (LONG COMPLEX) LONG COMPLEX
1795
1796 void genie_atan_double_compl (NODE_T * p)
1797 {
1798 CD_C_FUNCTION (p, catan_double);
1799 }
1800
1801 //! @brief PROC long cexp = (LONG COMPLEX) LONG COMPLEX
1802
1803 void genie_exp_double_compl (NODE_T * p)
1804 {
1805 CD_C_FUNCTION (p, cexp_double);
1806 }
1807
1808 //! @brief PROC long cln = (LONG COMPLEX) LONG COMPLEX
1809
1810 void genie_ln_double_compl (NODE_T * p)
1811 {
1812 CD_C_FUNCTION (p, clog_double);
1813 }
1814
1815 //! @brief PROC long csinh = (LONG COMPLEX) LONG COMPLEX
1816
1817 void genie_sinh_double_compl (NODE_T * p)
1818 {
1819 CD_C_FUNCTION (p, csinh_double);
1820 }
1821
1822 //! @brief PROC long ccosh = (LONG COMPLEX) LONG COMPLEX
1823
1824 void genie_cosh_double_compl (NODE_T * p)
1825 {
1826 CD_C_FUNCTION (p, ccosh_double);
1827 }
1828
1829 //! @brief PROC long ctanh = (LONG COMPLEX) LONG COMPLEX
1830
1831 void genie_tanh_double_compl (NODE_T * p)
1832 {
1833 CD_C_FUNCTION (p, ctanh_double);
1834 }
1835
1836 //! @brief PROC long casinh = (LONG COMPLEX) LONG COMPLEX
1837
1838 void genie_asinh_double_compl (NODE_T * p)
1839 {
1840 CD_C_FUNCTION (p, casinh_double);
1841 }
1842
1843 //! @brief PROC long cacosh = (LONG COMPLEX) LONG COMPLEX
1844
1845 void genie_acosh_double_compl (NODE_T * p)
1846 {
1847 CD_C_FUNCTION (p, cacosh_double);
1848 }
1849
1850 //! @brief PROC long catanh = (LONG COMPLEX) LONG COMPLEX
1851
1852 void genie_atanh_double_compl (NODE_T * p)
1853 {
1854 CD_C_FUNCTION (p, catanh_double);
1855 }
1856
1857 //! @brief PROC next long random = LONG REAL
1858
1859 void genie_next_random_double (NODE_T * p)
1860 {
1861 // This is 'real width' digits only.
1862 genie_next_random (p);
1863 genie_lengthen_real_to_double (p);
1864 }
1865
1866 #define CALL(g, x, y) {\
1867 ADDR_T pop_sp = A68_SP;\
1868 A68_LONG_REAL *z = (A68_LONG_REAL *) STACK_TOP;\
1869 DOUBLE_NUM_T _w_;\
1870 _w_.f = (x);\
1871 PUSH_VALUE (_p_, _w_, A68_LONG_REAL);\
1872 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);\
1873 (y) = VALUE (z).f;\
1874 A68_SP = pop_sp;\
1875 }
1876
1877 //! @brief Transform string into real-16.
1878
1879 DOUBLE_T string_to_double (char *s, char **end)
1880 {
1881 errno = 0;
1882 DOUBLE_T y[A68_DOUBLE_DIG];
1883 for (int i = 0; i < A68_DOUBLE_DIG; i++) {
1884 y[i] = 0.0q;
1885 }
1886 while (IS_SPACE (s[0])) {
1887 s++;
1888 }
1889 // Scan mantissa digits and put them into "y".
1890 DOUBLE_T W;
1891 if (s[0] == '-') {
1892 W = -1.0q;
1893 } else {
1894 W = 1.0q;
1895 }
1896 if (s[0] == '+' || s[0] == '-') {
1897 s++;
1898 }
1899 while (s[0] == '0') {
1900 s++;
1901 }
1902 int dot = -1, pos = 0, pow = 0;
1903 while (pow < A68_DOUBLE_DIG && s[pos] != NULL_CHAR && (IS_DIGIT (s[pos]) || s[pos] == POINT_CHAR)) {
1904 if (s[pos] == POINT_CHAR) {
1905 dot = pos;
1906 } else {
1907 int val = (int) s[pos] - (int) '0';
1908 y[pow] = W * val;
1909 W /= 10.0q;
1910 pow++;
1911 }
1912 pos++;
1913 }
1914 (*end) = &(s[pos]);
1915 // Sum from low to high to preserve precision.
1916 DOUBLE_T sum = 0.0q;
1917 for (int i = A68_DOUBLE_DIG - 1; i >= 0; i--) {
1918 sum = sum + y[i];
1919 }
1920 // See if there is an exponent.
1921 int expo;
1922 if (s[pos] != NULL_CHAR && TO_UPPER (s[pos]) == TO_UPPER (EXPONENT_CHAR)) {
1923 expo = (int) strtol (&(s[++pos]), end, 10);
1924 } else {
1925 expo = 0;
1926 }
1927 // Standardise.
1928 if (dot >= 0) {
1929 expo += dot - 1;
1930 } else {
1931 expo += pow - 1;
1932 }
1933 while (sum != 0.0q && fabs_double (sum) < 1.0q) {
1934 sum *= 10.0q;
1935 expo -= 1;
1936 }
1937 if (errno == 0) {
1938 return sum * ten_up_double (expo);
1939 } else {
1940 return 0.0q;
1941 }
1942 }
1943
1944 void genie_beta_inc_cf_double (NODE_T * p)
1945 {
1946 A68_LONG_REAL x, s, t;
1947 POP_OBJECT (p, &x, A68_LONG_REAL);
1948 POP_OBJECT (p, &t, A68_LONG_REAL);
1949 POP_OBJECT (p, &s, A68_LONG_REAL);
1950 errno = 0;
1951 PUSH_VALUE (p, dble (a68_beta_inc_double (VALUE (&s).f, VALUE (&t).f, VALUE (&x).f)), A68_LONG_REAL);
1952 MATH_RTE (p, errno != 0, M_LONG_REAL, NO_TEXT);
1953 }
1954
1955 void genie_beta_double (NODE_T * p)
1956 {
1957 A68_LONG_REAL a, b;
1958 POP_OBJECT (p, &b, A68_LONG_REAL);
1959 POP_OBJECT (p, &a, A68_LONG_REAL);
1960 errno = 0;
1961 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);
1962 MATH_RTE (p, errno != 0, M_LONG_REAL, NO_TEXT);
1963 }
1964
1965 void genie_ln_beta_double (NODE_T * p)
1966 {
1967 A68_LONG_REAL a, b;
1968 POP_OBJECT (p, &b, A68_LONG_REAL);
1969 POP_OBJECT (p, &a, A68_LONG_REAL);
1970 errno = 0;
1971 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);
1972 MATH_RTE (p, errno != 0, M_LONG_REAL, NO_TEXT);
1973 }
1974
1975 // LONG REAL infinity
1976
1977 void genie_infinity_double (NODE_T * p)
1978 {
1979 PUSH_VALUE (p, dble (a68_posinf_real ()), A68_LONG_REAL);
1980 }
1981
1982 // LONG REAL minus infinity
1983
1984 void genie_minus_infinity_double (NODE_T * p)
1985 {
1986 PUSH_VALUE (p, dble (a68_neginf_double ()), A68_LONG_REAL);
1987 }
1988
1989 #endif
