rts-conversion.c
1 //! @file rts-conversion.c
2 //! @author J. Marcel van der Veer
3
4 //! @section Copyright
5 //!
6 //! This file is part of Algol68G - an Algol 68 compiler-interpreter.
7 //! Copyright 2001-2026 J. Marcel van der Veer [algol68g@algol68genie.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 //! Conversion related routines.
25
26 #include "a68g.h"
27 #include "a68g-conversion.h"
28 #include "a68g-genie.h"
29 #include "a68g-prelude.h"
30 #include "a68g-mp.h"
31 #include "a68g-double.h"
32 #include "a68g-transput.h"
33
34
35 //! @brief 10 ** expo
36
37 REAL_T ten_up (int expo)
38 {
39 // Values for 10 ^ (2 ^ n) for conversion purposes on IEEE 754 platforms.
40 static REAL_T pow_10[] = {
41 10.0, 100.0, 1.0e4, 1.0e8, 1.0e16, 1.0e32, 1.0e64, 1.0e128, 1.0e256
42 };
43 BOOL_T neg_expo = (BOOL_T) (expo < 0);
44 if (neg_expo) {
45 expo = -expo;
46 }
47 PRELUDE_ERROR (expo > MAX_REAL_EXPO, NO_NODE, ERROR_INVALID_VALUE, M_REAL);
48 // This appears sufficiently accurate.
49 REAL_T dbl_expo = 1.0;
50 for (REAL_T *dep = pow_10; expo != 0; expo >>= 1, dep++) {
51 if (expo & 0x1) {
52 dbl_expo *= *dep;
53 }
54 }
55 return neg_expo ? 1 / dbl_expo : dbl_expo;
56 }
57
58
59 //! @brief Transform string into real-8.
60
61 REAL_T string_to_real (NODE_T *p, char *s)
62 {
63 errno = 0;
64 #if defined (PRECISE_CONVERSIONS)
65 return (REAL_T) string_to_double (p, s);
66 #else
67 (void) p;
68 return strtod (s, NO_REF);
69 #endif
70 }
71
72 #if (A68G_LEVEL >= 3)
73
74
75 //! @brief 10 ** expo
76
77 DOUBLE_T ten_up_double (int expo)
78 {
79 static DOUBLE_T pow_10_double[] = {
80 10.0q, 100.0q, 1.0e4q, 1.0e8q, 1.0e16q,
81 1.0e32q, 1.0e64q, 1.0e128q, 1.0e256q, 1.0e512q,
82 1.0e1024q, 1.0e2048q, 1.0e4096q
83 };
84 // This appears sufficiently accurate.
85 if (expo == 0) {
86 return 1.0q;
87 }
88 BOOL_T neg_expo = (BOOL_T) (expo < 0);
89 if (neg_expo) {
90 expo = -expo;
91 }
92 if (expo > MAX_DOUBLE_EXPO) {
93 expo = 0;
94 errno = EDOM;
95 }
96 PRELUDE_ERROR (expo > MAX_DOUBLE_EXPO, NO_NODE, ERROR_INVALID_VALUE, M_LONG_REAL);
97 DOUBLE_T dbl_expo = 1.0q;
98 for (DOUBLE_T *dep = pow_10_double; expo != 0; expo >>= 1, dep++) {
99 if (expo & 0x1) {
100 dbl_expo *= *dep;
101 }
102 }
103 return neg_expo ? 1.0q / dbl_expo : dbl_expo;
104 }
105
106 DOUBLE_NUM_T double_int_to_double (NODE_T * p, DOUBLE_NUM_T z)
107 {
108 int neg = D_NEG (z);
109 if (neg) {
110 z = abs_double_int (z);
111 }
112 DOUBLE_NUM_T w, radix;
113 w.f = 0.0q;
114 set_lw (radix, RADIX);
115 DOUBLE_T weight = 1.0q;
116 while (!D_ZERO (z)) {
117 DOUBLE_NUM_T digit;
118 digit = double_udiv (p, M_LONG_INT, z, radix, 1);
119 w.f = w.f + LW (digit) * weight;
120 z = double_udiv (p, M_LONG_INT, z, radix, 0);
121 weight = weight * RADIX_Q;
122 }
123 if (neg) {
124 w.f = -w.f;
125 }
126 return w;
127 }
128
129 DOUBLE_NUM_T double_to_double_int (NODE_T * p, DOUBLE_NUM_T z)
130 {
131 // This routines looks a lot like "strtol".
132 BOOL_T negative = (BOOL_T) (z.f < 0);
133 z.f = fabs_double (trunc_double (z.f));
134 if (z.f > CONST_2_UP_112_Q) {
135 errno = EDOM;
136 PRELUDE_ERROR (errno != 0, p, ERROR_MATH, M_LONG_REAL);
137 }
138 DOUBLE_NUM_T sum, weight, radix;
139 set_lw (sum, 0);
140 set_lw (weight, 1);
141 set_lw (radix, RADIX);
142 while (z.f > 0) {
143 DOUBLE_NUM_T term, digit, quot, rest;
144 quot.f = trunc_double (z.f / RADIX_Q);
145 rest.f = z.f - quot.f * RADIX_Q;
146 z.f = quot.f;
147 set_lw (digit, (INT_T) (rest.f));
148 term = double_umul (p, M_LONG_INT, digit, weight);
149 sum = double_uadd (p, M_LONG_INT, sum, term);
150 if (z.f > 0.0q) {
151 weight = double_umul (p, M_LONG_INT, weight, radix);
152 }
153 }
154 if (negative) {
155 return neg_double_int (sum);
156 } else {
157 return sum;
158 }
159 }
160
161
162 //! @brief LONG BITS value of LONG BITS denotation
163
164 DOUBLE_NUM_T string_to_double_bits (NODE_T * p, char *s)
165 {
166 errno = 0;
167 char *radix = NO_TEXT;
168 int base = (int) a68g_strtou (s, &radix, 10);
169 if (base < 2 || base > 16) {
170 diagnostic (A68G_RUNTIME_ERROR, p, ERROR_INVALID_RADIX, base);
171 exit_genie (p, A68G_RUNTIME_ERROR);
172 }
173 DOUBLE_NUM_T z;
174 set_lw (z, 0x0);
175 if (radix != NO_TEXT && TO_UPPER (radix[0]) == TO_UPPER (RADIX_CHAR) && errno == 0) {
176 DOUBLE_NUM_T w;
177 char *q = radix;
178 while (q[0] != NULL_CHAR) {
179 q++;
180 }
181 set_lw (w, 1);
182 while ((--q) != radix) {
183 int digit = char_value (q[0]);
184 if (digit < 0 && digit >= base) {
185 diagnostic (A68G_RUNTIME_ERROR, p, ERROR_IN_DENOTATION, M_LONG_BITS);
186 exit_genie (p, A68G_RUNTIME_ERROR);
187 } else {
188 DOUBLE_NUM_T v;
189 set_lw (v, digit);
190 v = double_umul (p, M_LONG_INT, v, w);
191 z = double_uadd (p, M_LONG_INT, z, v);
192 set_lw (v, base);
193 w = double_umul (p, M_LONG_INT, w, v);
194 }
195 }
196 } else {
197 diagnostic (A68G_RUNTIME_ERROR, p, ERROR_IN_DENOTATION, M_LONG_BITS);
198 exit_genie (p, A68G_RUNTIME_ERROR);
199 }
200 return (z);
201 }
202
203
204 //! @brief Value of LONG INT denotation
205
206 BOOL_T string_to_double_int (NODE_T * p, A68G_LONG_INT * z, char *s)
207 {
208 while (IS_SPACE (s[0])) {
209 s++;
210 }
211 // Get the sign
212 int sign = (s[0] == '-' ? -1 : 1);
213 if (s[0] == '+' || s[0] == '-') {
214 s++;
215 }
216 int end = 0;
217 while (s[end] != '\0') {
218 end++;
219 }
220 DOUBLE_NUM_T weight, ten, sum;
221 set_lw (sum, 0);
222 set_lw (weight, 1);
223 set_lw (ten, 10);
224 for (int k = end - 1; k >= 0; k--) {
225 DOUBLE_NUM_T term;
226 int digit = s[k] - '0';
227 set_lw (term, digit);
228 term = double_umul (p, M_LONG_INT, term, weight);
229 sum = double_uadd (p, M_LONG_INT, sum, term);
230 weight = double_umul (p, M_LONG_INT, weight, ten);
231 }
232 if (sign == -1) {
233 HW (sum) = HW (sum) | D_SIGN;
234 }
235 VALUE (z) = sum;
236 STATUS (z) = INIT_MASK;
237 return A68G_TRUE;
238 }
239
240
241 //! @brief Transform string into real-16.
242
243 DOUBLE_T string_to_double (NODE_T *p, char *s)
244 {
245 errno = 0;
246 #if defined (PRECISE_CONVERSIONS)
247 ADDR_T pop_sp = A68G_SP;
248 int digs = DIGITS (M_LONG_LONG_REAL);
249 A68G_SP += SIZE_MP (digs); // In case z is TOS.
250 MP_T *z = nil_mp (p, digs);
251 string_to_mp (p, z, s, digs);
252 if (A68G_NAN_MP (z)) {
253 return nanq ("");
254 }
255 DOUBLE_T u = mp_to_double (p, z, digs);
256 A68G_SP = pop_sp;
257 CHECK_DOUBLE_REAL (p, u, M_LONG_REAL);
258 return u;
259 #else
260 errno = 0;
261 DOUBLE_T y[A68G_DOUBLE_DIG];
262 for (int i = 0; i < A68G_DOUBLE_DIG; i++) {
263 y[i] = 0.0q;
264 }
265 while (IS_SPACE (s[0])) {
266 s++;
267 }
268 // Scan mantissa digits and put them into "y".
269 DOUBLE_T W;
270 if (s[0] == '-') {
271 W = -1.0q;
272 } else {
273 W = 1.0q;
274 }
275 if (s[0] == '+' || s[0] == '-') {
276 s++;
277 }
278 while (s[0] == '0') {
279 s++;
280 }
281 int dot = -1, pos = 0, pow = 0;
282 while (pow < A68G_DOUBLE_DIG && s[pos] != NULL_CHAR && (IS_DIGIT (s[pos]) || s[pos] == POINT_CHAR)) {
283 if (s[pos] == POINT_CHAR) {
284 dot = pos;
285 } else {
286 int val = (int) s[pos] - (int) '0';
287 y[pow] = W * val;
288 W /= 10.0q;
289 pow++;
290 }
291 pos++;
292 }
293 while (IS_DIGIT (s[pos])) {
294 pos++;
295 }
296 // Sum from low to high to preserve precision.
297 DOUBLE_T sum = a68g_neumaier_sum_double (y, A68G_DOUBLE_DIG);
298 // See if there is an exponent.
299 int expo;
300 if (s[pos] != NULL_CHAR && TO_UPPER (s[pos]) == TO_UPPER (EXPONENT_CHAR)) {
301 expo = (int) strtol (&(s[++pos]), NO_REF, 10);
302 } else {
303 expo = 0;
304 }
305 // Standardise.
306 if (dot >= 0) {
307 expo += dot - 1;
308 } else {
309 expo += pow - 1;
310 }
311 while (sum != 0.0q && fabs_double (sum) < 1.0q) {
312 sum *= 10.0q;
313 expo -= 1;
314 }
315 if (errno == 0) {
316 DOUBLE_T w = sum * ten_up_double (expo);
317 return w;
318 } else {
319 return 0.0q;
320 }
321 #endif
322 }
323
324 #endif
325
326
327 //! @brief Set "z" to 10 ** "n".
328
329 MP_T *ten_up_mp (NODE_T * p, MP_T * z, int n, int digs)
330 {
331 #if (A68G_LEVEL >= 3)
332 static MP_T y[LOG_MP_RADIX] = { 1, 10, 100, 1000, 10000, 100000, 1000000, 10000000, 100000000 };
333 #else
334 static MP_T y[LOG_MP_RADIX] = { 1, 10, 100, 1000, 10000, 100000, 1000000 };
335 #endif
336 if (n >= 0) {
337 (void) set_mp (z, y[n % LOG_MP_RADIX], n / LOG_MP_RADIX, digs);
338 } else {
339 (void) set_mp (z, y[(LOG_MP_RADIX + n % LOG_MP_RADIX) % LOG_MP_RADIX], (n + 1) / LOG_MP_RADIX - 1, digs);
340 }
341 check_mp_exp (p, z);
342 return z;
343 }
344
345
346 //! @brief INT value of BITS denotation
347
348 UNSIGNED_T bits_to_int (NODE_T * p, char *str)
349 {
350 errno = 0;
351 char *radix = NO_TEXT, *end = NO_TEXT;
352 int base = (int) a68g_strtou (str, &radix, 10);
353 if (radix != NO_TEXT && TO_UPPER (radix[0]) == TO_UPPER (RADIX_CHAR) && errno == 0) {
354 if (base < 2 || base > 16) {
355 diagnostic (A68G_RUNTIME_ERROR, p, ERROR_INVALID_RADIX, base);
356 exit_genie (p, A68G_RUNTIME_ERROR);
357 }
358 UNSIGNED_T bits = a68g_strtou (&(radix[1]), &end, base);
359 if (end != NO_TEXT && end[0] == NULL_CHAR && errno == 0) {
360 return bits;
361 }
362 }
363 diagnostic (A68G_RUNTIME_ERROR, p, ERROR_IN_DENOTATION, M_BITS);
364 exit_genie (p, A68G_RUNTIME_ERROR);
365 return 0;
366 }
367
368
369 //! @brief Convert to other radix, binary up to hexadecimal.
370
371 BOOL_T convert_radix (NODE_T * p, UNSIGNED_T z, int radix, int width)
372 {
373 reset_transput_buffer (EDIT_BUFFER);
374 if (radix < 2 || radix > 16) {
375 radix = 16;
376 }
377 if (width > 0) {
378 while (width > 0) {
379 int digit = (int) (z % (UNSIGNED_T) radix);
380 plusto_transput_buffer (p, digchar (digit), EDIT_BUFFER);
381 width--;
382 z /= (UNSIGNED_T) radix;
383 }
384 return z == 0;
385 } else if (width == 0) {
386 do {
387 int digit = (int) (z % (UNSIGNED_T) radix);
388 plusto_transput_buffer (p, digchar (digit), EDIT_BUFFER);
389 z /= (UNSIGNED_T) radix;
390 } while (z > 0);
391 return A68G_TRUE;
392 } else {
393 return A68G_FALSE;
394 }
395 }
396
397 #if (A68G_LEVEL >= 3)
398
399
400 //! @brief Convert to other radix, binary up to hexadecimal.
401
402 BOOL_T convert_radix_double (NODE_T * p, DOUBLE_NUM_T z, int radix, int width)
403 {
404 if (radix < 2 || radix > 16) {
405 radix = 16;
406 }
407 DOUBLE_NUM_T w, rad;
408 set_lw (rad, radix);
409 reset_transput_buffer (EDIT_BUFFER);
410 if (width > 0) {
411 while (width > 0) {
412 w = double_udiv (p, M_LONG_INT, z, rad, 1);
413 plusto_transput_buffer (p, digchar (LW (w)), EDIT_BUFFER);
414 width--;
415 z = double_udiv (p, M_LONG_INT, z, rad, 0);
416 }
417 return D_ZERO (z);
418 } else if (width == 0) {
419 do {
420 w = double_udiv (p, M_LONG_INT, z, rad, 1);
421 plusto_transput_buffer (p, digchar (LW (w)), EDIT_BUFFER);
422 z = double_udiv (p, M_LONG_INT, z, rad, 0);
423 }
424 while (!D_ZERO (z));
425 return A68G_TRUE;
426 } else {
427 return A68G_FALSE;
428 }
429 }
430
431 #endif
432
433
434 //! @brief Convert string to multi-precision number.
435
436 MP_T *string_to_mp (NODE_T * p, MP_T * z, char *s, int digs)
437 {
438 ABEND (s == NO_TEXT, ERROR_INTERNAL_CONSISTENCY, NO_TEXT);
439 BOOL_T ok = A68G_TRUE;
440 errno = 0;
441 SET_MP_ZERO (z, digs);
442 while (IS_SPACE (s[0])) {
443 s++;
444 }
445 // Get the sign.
446 int sign = (s[0] == '-' ? -1 : 1);
447 if (s[0] == '+' || s[0] == '-') {
448 s++;
449 }
450 // Scan mantissa digs and put them into "z".
451 while (s[0] == '0') {
452 s++;
453 }
454 int i = 0, dig = 1;
455 INT_T sum = 0, dot = -1, one = -1, pow = 0, W = MP_RADIX / 10;
456 while (s[i] != NULL_CHAR && dig <= digs && (IS_DIGIT (s[i]) || s[i] == POINT_CHAR)) {
457 if (s[i] == POINT_CHAR) {
458 dot = i;
459 } else {
460 int value = (int) s[i] - (int) '0';
461 if (one < 0 && value > 0) {
462 one = pow;
463 }
464 sum += W * value;
465 if (one >= 0) {
466 W /= 10;
467 }
468 pow++;
469 if (W < 1) {
470 MP_DIGIT (z, dig++) = (MP_T) sum;
471 sum = 0;
472 W = MP_RADIX / 10;
473 }
474 }
475 i++;
476 }
477 // Store the last digs.
478 if (dig <= digs) {
479 MP_DIGIT (z, dig++) = (MP_T) sum;
480 }
481 // See if there is an exponent.
482 INT_T expo;
483 if (s[i] != NULL_CHAR && TO_UPPER (s[i]) == TO_UPPER (EXPONENT_CHAR)) {
484 char *end;
485 expo = (int) strtol (&(s[++i]), &end, 10);
486 ok = (BOOL_T) (end[0] == NULL_CHAR);
487 } else {
488 expo = 0;
489 ok = (BOOL_T) (s[i] == NULL_CHAR);
490 }
491 // Calculate effective exponent.
492 if (dot >= 0) {
493 if (one > dot) {
494 expo -= one - dot + 1;
495 } else {
496 expo += dot - 1;
497 }
498 } else {
499 expo += pow - 1;
500 }
501 (void) align_mp (z, &expo, digs);
502 MP_EXPONENT (z) = (MP_DIGIT (z, 1) == 0 ? 0 : (MP_T) expo);
503 MP_DIGIT (z, 1) *= sign;
504 check_mp_exp (p, z);
505 if (errno == 0 && ok) {
506 return z;
507 } else {
508 set_nan_mp (z);
509 return z;
510 }
511 }
512
513
514 //! @brief Convert integer to multi-precison number.
515
516 MP_T *int_to_mp (NODE_T * p, MP_T * z, INT_T k, int digs)
517 {
518 int sign_k = 1;
519 if (k < 0) {
520 k = -k;
521 sign_k = -1;
522 }
523 INT_T m = k, n = 0;
524 while ((m /= MP_RADIX) != 0) {
525 n++;
526 }
527 (void) set_mp (z, 0, n, digs);
528 for (int j = 1 + n; j >= 1; j--) {
529 MP_DIGIT (z, j) = (MP_T) (k % MP_RADIX);
530 k /= MP_RADIX;
531 }
532 MP_DIGIT (z, 1) = sign_k * MP_DIGIT (z, 1);
533 check_mp_exp (p, z);
534 return z;
535 }
536
537
538 //! @brief Convert unt to multi-precison number.
539
540 MP_T *unt_to_mp (NODE_T * p, MP_T * z, UNSIGNED_T k, int digs)
541 {
542 int m = k, n = 0;
543 while ((m /= MP_RADIX) != 0) {
544 n++;
545 }
546 (void) set_mp (z, 0, n, digs);
547 for (int j = 1 + n; j >= 1; j--) {
548 MP_DIGIT (z, j) = (MP_T) (k % MP_RADIX);
549 k /= MP_RADIX;
550 }
551 check_mp_exp (p, z);
552 return z;
553 }
554
555
556 //! @brief Convert multi-precision number to integer.
557
558 INT_T mp_to_int (NODE_T * p, MP_T * z, int digs)
559 {
560 // This routines looks a lot like "strtol".
561 INT_T expo = (int) MP_EXPONENT (z), sum = 0, weight = 1;
562 if (expo >= digs) {
563 diagnostic (A68G_RUNTIME_ERROR, p, ERROR_OUT_OF_BOUNDS, MOID (p));
564 exit_genie (p, A68G_RUNTIME_ERROR);
565 }
566 BOOL_T negative = (BOOL_T) (MP_DIGIT (z, 1) < 0);
567 if (negative) {
568 MP_DIGIT (z, 1) = -MP_DIGIT (z, 1);
569 }
570 for (int j = 1 + expo; j >= 1; j--) {
571 if ((MP_INT_T) MP_DIGIT (z, j) > A68G_MAX_INT / weight) {
572 diagnostic (A68G_RUNTIME_ERROR, p, ERROR_OUT_OF_BOUNDS, M_INT);
573 exit_genie (p, A68G_RUNTIME_ERROR);
574 }
575 INT_T term = (MP_INT_T) MP_DIGIT (z, j) * weight;
576 if (sum > A68G_MAX_INT - term) {
577 diagnostic (A68G_RUNTIME_ERROR, p, ERROR_OUT_OF_BOUNDS, M_INT);
578 exit_genie (p, A68G_RUNTIME_ERROR);
579 }
580 sum += term;
581 weight *= MP_RADIX;
582 }
583 return negative ? -sum : sum;
584 }
585
586
587 //! @brief Convert REAL_T to multi-precison number.
588
589 MP_T *real_to_mp (NODE_T * p, MP_T * z, REAL_T x, int digs)
590 {
591 #if (A68G_LEVEL >= 3)
592 return double_to_mp (p, z, (DOUBLE_T) x, digs);
593 #else
594 if (a68g_isnan_real (x)) {
595 set_nan_mp (z);
596 } else if (x == A68G_PLUS_INF_REAL) {
597 set_plus_inf_mp (z);
598 } else if (x == A68G_MINUS_INF_REAL) {
599 set_minus_inf_mp (z);
600 } else {
601 // Legacy code.
602 SET_MP_ZERO (z, digs);
603 if (a68g_isinf_real (x)) {
604 if (x == A68G_PLUS_INF_REAL) {
605 set_plus_inf_mp (z);
606 } else {
607 set_minus_inf_mp (z);
608 }
609 return z;
610 } else if (a68g_isnan_real (x)) {
611 set_nan_mp (z);
612 return z;
613 }
614 if (x == 0.0) {
615 return z;
616 }
617 // Small integers can be done better by int_to_mp.
618 if (ABS (x) < MP_RADIX && trunc (x) == x) {
619 return int_to_mp (p, z, (INT_T) trunc (x), digs);
620 }
621 int sign_x = SIGN (x);
622 // Scale to [0, 0.1>.
623 REAL_T a = ABS (x);
624 INT_T expo = (int) log10 (a);
625 a /= ten_up (expo);
626 expo--;
627 if (a >= 1) {
628 a /= 10;
629 expo++;
630 }
631 // Transport digs of x to the mantissa of z.
632 for (int j = 1, k = 0; k <= A68G_REAL_DIG && j <= digs; j++, k += LOG_MP_RADIX) {
633 DOUBLE_T dig;
634 a = modf (a * MP_RADIX, &dig);
635 MP_DIGIT (z, j) = (INT_T) dig;
636 }
637 // Return the result.
638 (void) align_mp (z, &expo, digs);
639 MP_EXPONENT (z) = (MP_T) expo;
640 MP_DIGIT (z, 1) *= sign_x;
641 check_mp_exp (p, z);
642 }
643 return z;
644 #endif
645 }
646
647
648 //! @brief Convert multi-precision number to real.
649
650 REAL_T mp_to_real (NODE_T * p, MP_T * z, int digs)
651 {
652 #if (A68G_LEVEL >= 3)
653 return (REAL_T) mp_to_double (p, z, digs);
654 #else
655 if (A68G_NAN_MP (z)) {
656 return A68G_NAN_REAL;
657 } else if (A68G_PLUS_INF_MP (z)) {
658 return A68G_PLUS_INF_REAL;
659 } else if (A68G_MINUS_INF_MP (z)) {
660 return A68G_MINUS_INF_REAL;
661 } else {
662 // Legacy code, like "strtod".
663 (void) p;
664 if (MP_EXPONENT (z) * (MP_T) LOG_MP_RADIX <= (MP_T) A68G_REAL_MIN_EXP) {
665 return 0;
666 } else {
667 REAL_T terms[MP_MAX_DIGITS];
668 REAL_T weight = 1.0;
669 int lim = A68G_MIN (digs, MP_MAX_DIGITS);
670 for (int k = 0; k < lim; k++) {
671 terms[k] = ABS (MP_DIGIT (z, k + 1)) * weight;
672 weight /= MP_RADIX;
673 }
674 REAL_T sum = a68g_neumaier_sum_real (terms, lim) * ten_up ((int) (MP_EXPONENT (z) * LOG_MP_RADIX));
675 CHECK_REAL (p, sum, M_REAL);
676 return MP_DIGIT (z, 1) >= 0 ? sum : -sum;
677 }
678 }
679 #endif
680 }
681
682 #if (A68G_LEVEL >= 3)
683
684
685 //! @brief Convert double to multi-precison number.
686
687 MP_T *double_to_mp (NODE_T * p, MP_T * z, DOUBLE_T u, int digs)
688 {
689 SET_MP_ZERO (z, digs);
690 // Infinity and beyond.
691 if (a68g_isinf_double (u)) {
692 if (u == A68G_PLUS_INF_DOUBLE) {
693 set_plus_inf_mp (z);
694 } else {
695 set_minus_inf_mp (z);
696 }
697 return z;
698 } else if (a68g_isnan_double (u)) {
699 set_nan_mp (z);
700 return z;
701 }
702 if (u == 0.0q) {
703 return z;
704 }
705 // Small integers can be done better by int_to_mp.
706 if (ABS (u) < MP_RADIX && trunc_double (u) == u) {
707 return int_to_mp (p, z, (int) trunc_double (u), digs);
708 }
709 // Conversion depends on availability of real*16.
710 #if defined (PRECISE_CONVERSIONS)
711 DOUBLE_NUM_T w = dble (u);
712 if (SUBNORMAL (w)) {
713 // Strict flush-to-zero policy.
714 return z;
715 }
716 // The IEEE 128-bit float is dissected as mantissa * 2^exponent.
717 // Conversion to MP number is done with the decimal LONG LONG library.
718 // This is sufficiently accurate.
719 int sign_u = SIGN (u);
720 // Scale to [0.5, 1>
721 int exp2;
722 (void) frexpq (u, &exp2);
723 // Convert mantissa to MP.
724 IEEE_754_MANTISSA (w);
725 double_int_to_mp (p, z, w, digs);
726 // Compute MP value by multiplying by 2 ^ exponent.
727 // Scale the 113-bit unsigned mantissa back to [0.5, 1> by
728 // subtracting 113 from the exponent.
729 ADDR_T pop_sp = A68G_SP;
730 A68G_SP += SIZE_MP (digs); // In case z is TOS.
731 ldexp_mp (p, z, z, exp2 - 113, digs);
732 A68G_SP = pop_sp;
733 #else
734 // Legacy approach.
735 int sign_u = SIGN (u);
736 // Scale to [0, 0.1>.
737 DOUBLE_T a = ABS (u);
738 INT_T expo = (int) log10_double (a);
739 a /= ten_up_double (expo);
740 expo--;
741 if (a >= 1.0q) {
742 a /= 10.0q;
743 expo++;
744 }
745 // Transport digits of u to the mantissa of z.
746 for (int j = 1, k = 0; k <= 1 + A68G_DOUBLE_DIG && j <= digs; j++, k += LOG_MP_RADIX) {
747 DOUBLE_T dig;
748 a = modfq (a * MP_RADIX_Q, &dig);
749 MP_DIGIT (z, j) = (INT_T) dig;
750 }
751 (void) align_mp (z, &expo, digs);
752 MP_EXPONENT (z) = (MP_T) expo;
753 #endif
754 // Check and return.
755 MP_DIGIT (z, 1) *= sign_u;
756 check_mp_exp (p, z);
757 return z;
758 }
759
760
761 //! @brief Convert multi-precision number to real.
762
763 DOUBLE_T mp_to_double (NODE_T * p, MP_T * z, int digs)
764 {
765 if (A68G_PLUS_INF_MP (z)) {
766 return A68G_PLUS_INF_DOUBLE;
767 } else if (A68G_MINUS_INF_MP (z)) {
768 return A68G_MINUS_INF_DOUBLE;
769 } else if (A68G_NAN_MP (z)) {
770 return A68G_NAN_DOUBLE;
771 }
772 if (MP_EXPONENT (z) * (MP_T) LOG_MP_RADIX <= (MP_T) A68G_DOUBLE_MIN_EXP) {
773 return 0.0q;
774 } else {
775 // Compute the DOUBLE value of |z|.
776 int lim = A68G_MIN (digs, MP_MAX_DIGITS);
777 DOUBLE_T terms[MP_MAX_DIGITS], f = 1.0q;
778 for (int k = 0; k < lim; k++) {
779 terms[k] = ABS (MP_DIGIT (z, k + 1)) * f;
780 f /= MP_RADIX_Q;
781 }
782 DOUBLE_T u = a68g_neumaier_sum_double (terms, lim) * ten_up_double ((int) (MP_EXPONENT (z) * LOG_MP_RADIX));
783 #if defined (PRECISE_CONVERSIONS)
784 // We have a good approximation of the 128-bit real value, in 'u'.
785 // Conversion error is expected one ulp, that we correct here.
786 // This way, double_to_mp and mp_to_double are effectively 'mirrors' and
787 // the relation 'x = SHORTEN LENG x' holds for LONG REAL x.
788 #define H_MASK 0x0000ffffffffffffULL
789 #define L_MASK 0xffffffffffffffffULL
790 DOUBLE_NUM_T v = dble (u);
791 if (SUBNORMAL (v)) {
792 // Strict flush-to-zero policy.
793 return 0.0q;
794 }
795 // Rounding a converted number of form (2^k - 1 ulp), k<0.
796 // Doing this here guarantees correct exponent later.
797 if ((HW (v) & H_MASK) == H_MASK && (LW (v) == L_MASK)) {
798 LW (v) = 0;
799 HW (v)++;
800 }
801 // Scale to [0.5, 1>
802 int exp2;
803 (void) frexpq (u, &exp2);
804 // Convert u 113-bit mantissa to MP unsigned integer.
805 // Keep the exponent.
806 UNSIGNED_T expo = HW (v) & DOUBLE_EXPONENT_MASK;
807 IEEE_754_MANTISSA (v);
808 ADDR_T pop_sp = A68G_SP;
809 A68G_SP += SIZE_MP (digs); // In case z is TOS.
810 MP_T *mant = nil_mp (p, digs);
811 ldexp_mp (p, mant, z, 113 - exp2, digs);
812 round_mp (p, mant, mant, digs);
813 DOUBLE_NUM_T w = mp_to_double_int (p, mant, digs);
814 A68G_SP = pop_sp;
815 IEEE_754_MANTISSA (w);
816 // Restore the exponent.
817 HW (w) = (HW (w) & ~DOUBLE_EXPONENT_MASK) | expo;
818 u = w.f;
819 #undef H_MASK
820 #undef L_MASK
821 #endif
822 CHECK_DOUBLE_REAL (p, u, M_LONG_REAL);
823 return (MP_DIGIT (z, 1) >= 0 ? u : -u);
824 }
825 }
826
827 MP_T *double_int_to_mp (NODE_T * p, MP_T * z, DOUBLE_NUM_T k, int digs)
828 {
829 int negative = D_NEG (k);
830 if (negative) {
831 k = neg_double_int (k);
832 }
833 DOUBLE_NUM_T radix;
834 set_lw (radix, MP_RADIX);
835 DOUBLE_NUM_T k2 = k;
836 int n = 0;
837 do {
838 k2 = double_udiv (p, M_LONG_INT, k2, radix, 0);
839 if (!D_ZERO (k2)) {
840 n++;
841 }
842 }
843 while (!D_ZERO (k2));
844 SET_MP_ZERO (z, digs);
845 MP_EXPONENT (z) = (MP_T) n;
846 for (int j = 1 + n; j >= 1; j--) {
847 DOUBLE_NUM_T term = double_udiv (p, M_LONG_INT, k, radix, 1);
848 MP_DIGIT (z, j) = (MP_T) LW (term);
849 k = double_udiv (p, M_LONG_INT, k, radix, 0);
850 }
851 MP_DIGIT (z, 1) = (negative ? -MP_DIGIT (z, 1) : MP_DIGIT (z, 1));
852 check_mp_exp (p, z);
853 return z;
854 }
855
856
857 //! @brief Convert multi-precision number to integer.
858
859 DOUBLE_NUM_T mp_to_double_int (NODE_T * p, MP_T * z, int digs)
860 {
861 // This routines looks a lot like "strtol".
862 int expo = (int) MP_EXPONENT (z);
863 DOUBLE_NUM_T sum, weight;
864 set_lw (sum, 0);
865 set_lw (weight, 1);
866 BOOL_T negative;
867 if (expo >= digs) {
868 diagnostic (A68G_RUNTIME_ERROR, p, ERROR_OUT_OF_BOUNDS, MOID (p));
869 exit_genie (p, A68G_RUNTIME_ERROR);
870 }
871 negative = (BOOL_T) (MP_DIGIT (z, 1) < 0);
872 if (negative) {
873 MP_DIGIT (z, 1) = -MP_DIGIT (z, 1);
874 }
875 for (int j = 1 + expo; j >= 1; j--) {
876 DOUBLE_NUM_T term, digit, radix;
877 set_lw (digit, (MP_INT_T) MP_DIGIT (z, j));
878 term = double_umul (p, M_LONG_INT, digit, weight);
879 sum = double_uadd (p, M_LONG_INT, sum, term);
880 set_lw (radix, MP_RADIX);
881 weight = double_umul (p, M_LONG_INT, weight, radix);
882 }
883 if (negative) {
884 return neg_double_int (sum);
885 } else {
886 return sum;
887 }
888 }
889
890 #endif
891
892
893 //! @brief Convert string to required mode and store.
894
895 BOOL_T genie_string_to_value_internal (NODE_T * p, MOID_T * m, char *a, BYTE_T * item)
896 {
897 errno = 0;
898 // strto.. does not mind empty strings.
899 if (strlen (a) == 0) {
900 return A68G_FALSE;
901 }
902 if (m == M_INT) {
903 A68G_INT *z = (A68G_INT *) item;
904 char *end;
905 VALUE (z) = (INT_T) a68g_strtoi (a, &end, 10);
906 if (end[0] == NULL_CHAR && errno == 0) {
907 STATUS (z) = INIT_MASK;
908 return A68G_TRUE;
909 } else {
910 return A68G_FALSE;
911 }
912 }
913 if (m == M_REAL) {
914 A68G_REAL *z = (A68G_REAL *) item;
915 VALUE (z) = string_to_real (p, a);
916 if (errno == 0) {
917 STATUS (z) = INIT_MASK;
918 return A68G_TRUE;
919 } else {
920 return A68G_FALSE;
921 }
922 }
923 #if (A68G_LEVEL >= 3)
924 if (m == M_LONG_INT) {
925 A68G_LONG_INT *z = (A68G_LONG_INT *) item;
926 if (string_to_double_int (p, z, a) == A68G_FALSE) {
927 return A68G_FALSE;
928 }
929 STATUS (z) = INIT_MASK;
930 return A68G_TRUE;
931 }
932 if (m == M_LONG_REAL) {
933 A68G_LONG_REAL *z = (A68G_LONG_REAL *) item;
934 VALUE (z).f = string_to_double (p, a);
935 PRELUDE_ERROR (errno != 0, p, ERROR_MATH, M_LONG_REAL);
936 if (errno == 0) {
937 STATUS (z) = INIT_MASK;
938 return A68G_TRUE;
939 } else {
940 return A68G_FALSE;
941 }
942 }
943 if (m == M_LONG_BITS) {
944 A68G_LONG_BITS *z = (A68G_LONG_BITS *) item;
945 int ret = A68G_TRUE;
946 DOUBLE_NUM_T b;
947 set_lw (b, 0x0);
948 if (a[0] == FLIP_CHAR || a[0] == FLOP_CHAR) {
949 // [] BOOL denotation is "TTFFFFTFT ...".
950 if (strlen (a) > (size_t) A68G_LONG_BITS_WIDTH) {
951 errno = ERANGE;
952 ret = A68G_FALSE;
953 } else {
954 int n = 1;
955 UNSIGNED_T k = 0x1;
956 for (INT_T j = (INT_T) strlen (a) - 1; j >= 0; j--) {
957 if (a[j] == FLIP_CHAR) {
958 if (n <= A68G_LONG_BITS_WIDTH / 2) {
959 LW (b) |= k;
960 } else {
961 HW (b) |= k;
962 }
963 } else if (a[j] != FLOP_CHAR) {
964 ret = A68G_FALSE;
965 }
966 k <<= 1;
967 }
968 }
969 VALUE (z) = b;
970 } else {
971 // BITS denotation.
972 VALUE (z) = string_to_double_bits (p, a);
973 }
974 return ret;
975 }
976 #else
977 if (m == M_LONG_BITS || m == M_LONG_LONG_BITS) {
978 int digits = DIGITS (m);
979 int status = A68G_TRUE;
980 ADDR_T pop_sp = A68G_SP;
981 MP_T *z = (MP_T *) item;
982 if (a[0] == FLIP_CHAR || a[0] == FLOP_CHAR) {
983 // [] BOOL denotation is "TTFFFFTFT ...".
984 if (strlen (a) > (size_t) A68G_BITS_WIDTH) {
985 errno = ERANGE;
986 status = A68G_FALSE;
987 } else {
988 MP_T *w = lit_mp (p, 1, 0, digits);
989 SET_MP_ZERO (z, digits);
990 for (INT_T j = (INT_T) strlen (a) - 1; j >= 0; j--) {
991 if (a[j] == FLIP_CHAR) {
992 (void) add_mp (p, z, z, w, digits);
993 } else if (a[j] != FLOP_CHAR) {
994 status = A68G_FALSE;
995 }
996 (void) mul_mp_digit (p, w, w, (MP_T) 2, digits);
997 }
998 }
999 } else {
1000 // BITS denotation is also allowed.
1001 mp_strtou (p, z, a, m);
1002 }
1003 A68G_SP = pop_sp;
1004 if (errno != 0 || status == A68G_FALSE) {
1005 return A68G_FALSE;
1006 }
1007 SET_INIT_MP (z);
1008 return A68G_TRUE;
1009 }
1010 #endif
1011 if (m == M_LONG_INT || m == M_LONG_LONG_INT) {
1012 int digits = DIGITS (m);
1013 MP_T *z = (MP_T *) item;
1014 string_to_mp (p, z, a, digits);
1015 if (A68G_NAN_MP (z)) {
1016 return A68G_FALSE;
1017 }
1018 if (!check_mp_int (z, m)) {
1019 errno = ERANGE;
1020 return A68G_FALSE;
1021 }
1022 SET_INIT_MP (z);
1023 return A68G_TRUE;
1024 }
1025 if (m == M_LONG_REAL || m == M_LONG_LONG_REAL) {
1026 int digits = DIGITS (m);
1027 MP_T *z = (MP_T *) item;
1028 string_to_mp (p, z, a, digits);
1029 if (A68G_NAN_MP (z)) {
1030 return A68G_FALSE;
1031 }
1032 SET_INIT_MP (z);
1033 return A68G_TRUE;
1034 }
1035 if (m == M_BOOL) {
1036 A68G_BOOL *z = (A68G_BOOL *) item;
1037 char q = a[0], flip = FLIP_CHAR, flop = FLOP_CHAR;
1038 if (q == flip || q == flop) {
1039 VALUE (z) = (BOOL_T) (q == flip);
1040 STATUS (z) = INIT_MASK;
1041 return A68G_TRUE;
1042 } else {
1043 return A68G_FALSE;
1044 }
1045 }
1046 if (m == M_BITS) {
1047 A68G_BITS *z = (A68G_BITS *) item;
1048 int status = A68G_TRUE;
1049 if (a[0] == FLIP_CHAR || a[0] == FLOP_CHAR) {
1050 // [] BOOL denotation is "TTFFFFTFT ...".
1051 if (strlen (a) > (size_t) A68G_BITS_WIDTH) {
1052 errno = ERANGE;
1053 status = A68G_FALSE;
1054 } else {
1055 UNSIGNED_T k = 0x1;
1056 VALUE (z) = 0;
1057 for (INT_T j = (INT_T) strlen (a) - 1; j >= 0; j--) {
1058 if (a[j] == FLIP_CHAR) {
1059 VALUE (z) += k;
1060 } else if (a[j] != FLOP_CHAR) {
1061 status = A68G_FALSE;
1062 }
1063 k <<= 1;
1064 }
1065 }
1066 } else {
1067 // BITS denotation is also allowed.
1068 VALUE (z) = bits_to_int (p, a);
1069 }
1070 if (errno != 0 || status == A68G_FALSE) {
1071 return A68G_FALSE;
1072 }
1073 STATUS (z) = INIT_MASK;
1074 return A68G_TRUE;
1075 }
1076 return A68G_FALSE;
1077 }
© J.M. van der Veer • jmvdveer@algol68genie.nl
