quad-conversion.c
1 //! @file quad-conversion.c
2 //!
3 //! @author (see below)
4 //!
5 //! @section Copyright
6 //!
7 //! This file is part of Algol68G - an Algol 68 compiler-interpreter.
8 //! Copyright 2001-2023 J. Marcel van der Veer [algol68g@xs4all.nl].
9 //!
10 //! @section License
11 //!
12 //! This program is free software; you can redistribute it and/or modify it
13 //! under the terms of the GNU General Public License as published by the
14 //! Free Software Foundation; either version 3 of the License, or
15 //! (at your option) any later version.
16 //!
17 //! This program is distributed in the hope that it will be useful, but
18 //! WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
19 //! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
20 //! more details. You should have received a copy of the GNU General Public
21 //! License along with this program. If not, see [http://www.gnu.org/licenses/].
22
23 //! @section Synopsis
24 //!
25 //! Fixed precision LONG LONG REAL/COMPLEX conversions.
26
27 // This code is based in part on the HPA Library, a branch of the CCMath library.
28 // The CCMath library and derived HPA Library are free software under the terms
29 // of the GNU Lesser General Public License version 2.1 or any later version.
30 //
31 // Sources:
32 // [https://directory.fsf.org/wiki/HPAlib]
33 // [http://download-mirror.savannah.gnu.org/releases/hpalib/]
34 // [http://savannah.nongnu.org/projects/hpalib]
35 //
36 // Copyright 2000 Daniel A. Atkinson [DanAtk@aol.com]
37 // Copyright 2004 Ivano Primi [ivprimi@libero.it]
38 // Copyright 2023 Marcel van der Veer [algol68g@xs4all.nl] - A68G version.
39
40 #include "a68g.h"
41
42 #if (A68_LEVEL >= 3)
43
44 #include "a68g-quad.h"
45
46 // Conversion.
47
48 //! @brief int_to_quad_real
49
50 QUAD_T int_to_quad_real (int n)
51 {
52 REAL32 pe;
53 unt_2 *pc, e;
54 unt k, h;
55 memset (pe, 0, sizeof (pe));
56 k = ABS (n);
57 pc = (unt_2 *) &k;
58 if (n == 0) {
59 return *(QUAD_T *) pe;
60 }
61
62 #if (__BYTE_ORDER == __LITTLE_ENDIAN)
63 pe[1] = *(pc + 1);
64 pe[2] = *pc;
65 #else
66 pe[1] = *pc;
67 pe[2] = *(pc + 1);
68 #endif
69
70 for (e = 0, h = 1; h <= k && e < ((8 * sizeof (unt)) - 1); h <<= 1, ++e);
71 if (h <= k) {
72 e += 1;
73 }
74 *pe = QUAD_REAL_BIAS + e - 1;
75 if (n < 0) {
76 *pe |= QUAD_REAL_M_SIGN;
77 }
78 lshift_quad_real ((8 * sizeof (unt)) - e, pe + 1, FLT256_LEN);
79 return *(QUAD_T *) pe;
80 }
81
82 //! @brief quad_real_to_real
83
84 REAL_T quad_real_to_real (QUAD_T s)
85 {
86 union {unt_2 pe[4]; REAL_T d;} v;
87 unt_2 *pc, u;
88 int_2 i, e;
89 pc = (unt_2 *) &QV (s);
90 u = *pc & QUAD_REAL_M_SIGN;
91 e = (*pc & QUAD_REAL_M_EXP) - QUAD_REAL_DBL_BIAS;
92 if (e >= QUAD_REAL_DBL_MAX) {
93 return (!u ? DBL_MAX : -DBL_MAX);
94 }
95 if (e <= 0) {
96 return 0.;
97 }
98 for (i = 0; i < 4; v.pe[i] = *++pc, i++);
99 v.pe[0] &= QUAD_REAL_M_EXP;
100 rshift_quad_real (QUAD_REAL_DBL_LEX - 1, v.pe, 4);
101 v.pe[0] |= (e << (16 - QUAD_REAL_DBL_LEX));
102 v.pe[0] |= u;
103 #if (__BYTE_ORDER == __LITTLE_ENDIAN)
104 u = v.pe[3];
105 v.pe[3] = v.pe[0];
106 v.pe[0] = u;
107 u = v.pe[2];
108 v.pe[2] = v.pe[1];
109 v.pe[1] = u;
110 #endif
111 return v.d;
112 }
113
114 //! @brief real_to_quad_real
115
116 QUAD_T real_to_quad_real (REAL_T y)
117 {
118 REAL32 pe;
119 unt_2 *pc, u;
120 int_2 i, e;
121 if (y < DBL_MIN && y > -DBL_MIN) {
122 return QUAD_REAL_ZERO;
123 }
124 memset (pe, 0, sizeof (pe));
125 pc = (unt_2 *) &y;
126 #if (__BYTE_ORDER == __LITTLE_ENDIAN)
127 pc += 3;
128 #endif
129 u = *pc & QUAD_REAL_M_SIGN;
130 e = QUAD_REAL_DBL_BIAS + ((*pc & QUAD_REAL_M_EXP) >> (16 - QUAD_REAL_DBL_LEX));
131 #if (__BYTE_ORDER == __LITTLE_ENDIAN)
132 for (i = 1; i < 5; pe[i] = *pc--, i++);
133 #else
134 for (i = 1; i < 5; pe[i] = *pc++, i++);
135 #endif
136 pc = pe + 1;
137 lshift_quad_real (QUAD_REAL_DBL_LEX - 1, pc, 4);
138 *pc |= QUAD_REAL_M_SIGN;
139 *pe = e;
140 *pe |= u;
141 return *(QUAD_T *) pe;
142 }
143
144 //! @brief double_real_to_quad_real
145
146 QUAD_T double_real_to_quad_real (DOUBLE_T x)
147 {
148 if (x == 0.0q) {
149 return QUAD_REAL_ZERO;
150 }
151 int sinf = isinfq (x);
152 if (sinf == 1) {
153 return QUAD_REAL_PINF;
154 } else if (sinf == -1) {
155 return QUAD_REAL_MINF;
156 }
157 if (isnanq (x)) {
158 return QUAD_REAL_NAN;
159 }
160 QUAD_T z = QUAD_REAL_ZERO;
161 unt_2 *y = (unt_2 *) &x;
162 for (unt i = 0; i < 8; i++) {
163 #if (__BYTE_ORDER == __LITTLE_ENDIAN)
164 QV (z)[i] = y[7 - i];
165 #else
166 QV (z)[i] = y[i];
167 #endif
168 }
169 // double_real skips the default first bit, HPA lib does not.
170 unt_2 cy = 0x8000;
171 for (unt i = 1; i < 8; i++) {
172 if (QV (z)[i] & 0x1) {
173 QV (z)[i] = (QV (z)[i] >> 1) | cy;
174 cy = 0x8000;
175 } else {
176 QV (z)[i] = (QV (z)[i] >> 1) | cy;
177 cy = 0x0;
178 }
179 }
180 QV (z)[8] |= cy;
181 return z;
182 }
183
184 DOUBLE_T quad_real_to_double_real (QUAD_T x)
185 {
186 REAL16 y; REAL32 z;
187 int i;
188 memcpy (z, x.value, sizeof (QUAD_T));
189 // Catch NaN, +-Inf is handled correctly.
190 if (isNaN_quad_real (&x)) {
191 z[0] = 0x7fff;
192 z[1] = 0xffff;
193 }
194 // double_real skips the default first bit, HPA lib does not.
195 unt_2 cy = (z[8] & 0x8000 ? 0x1 : 0x0);
196 for (i = 7; i > 0; i--) {
197 if (z[i] & 0x8000) {
198 z[i] = (z[i] << 1) | cy;
199 cy = 0x1;
200 } else {
201 z[i] = (z[i] << 1) | cy;
202 cy = 0x0;
203 }
204 }
205 for (i = 0; i < 8; i++) {
206 #if (__BYTE_ORDER == __LITTLE_ENDIAN)
207 y[i] = z[7 - i];
208 #else
209 y[i] = z[i];
210 #endif
211 }
212 // Avoid 'dereferencing type-punned pointer will break strict-aliasing rules'
213 DOUBLE_T u;
214 memcpy (&u, &y[0], sizeof (DOUBLE_T));
215 return u;
216 }
217
218 //! @brief string_to_quad_real
219
220 QUAD_T string_to_quad_real (const char *q, char **endptr)
221 {
222 #define QUAD_UPB (QUAD_REAL_MAX_10EX + 100)
223 QUAD_T s, f;
224 REAL32 pc;
225 unt_2 *pn, *pf, *pa, *pb;
226 unt_2 sfg, ibex, fbex;
227 unt n;
228 int j;
229 int idex = 0, fdex = 0, c, m;
230 int_2 noip, nofp;
231 const char *ptr;
232 pn = (unt_2 *) &QV (s);
233 pf = (unt_2 *) &QV (f);
234 for (j = 0; j <= FLT256_LEN; pn[j] = pf[j] = pc[j] = 0, ++j);
235 sfg = 0;
236 m = FLT256_LEN + 1;
237 if (endptr != NULL) {
238 *endptr = (char *) q;
239 }
240 // Skip leading white space
241 while ((isspace (*q))) {
242 q++;
243 }
244 // Sign
245 if (*q == '+') {
246 ++q;
247 } else if (*q == '-') {
248 sfg = 0x8000;
249 ++q;
250 }
251 // Integer part
252 for (ptr = q; (c = *q - '0') >= 0 && c <= 9; ++q) {
253 if (pn[0]) {
254 ++idex;
255 } else {
256 lshift_quad_real (1, pn, m);
257 for (j = 0; j < m; ++j) {
258 pc[j] = pn[j];
259 }
260 lshift_quad_real (2, pn, m);
261 for (n = (unt) c, pa = pn + FLT256_LEN, pb = pc + FLT256_LEN; pa >= pn; pa--, pb--) {
262 n += *pa + *pb;
263 *pa = n;
264 n >>= 16;
265 }
266 }
267 }
268 for (j = 0; j < m && pn[j] == 0; ++j);
269 if (j == m) {
270 ibex = 0;
271 } else {
272 ibex = QUAD_REAL_BIAS + QUAD_REAL_MAX_P - 1;
273 if (j != 0) {
274 j <<= 4;
275 ibex -= j;
276 lshift_quad_real (j, pn, m);
277 }
278 while (pn[0]) {
279 rshift_quad_real (1, pn, m);
280 ++ibex;
281 }
282 pn[0] = ibex | sfg;
283 }
284 noip = ptr == q;
285 // End Integer part
286 if (*q == '.') {
287 // Fractionary part
288 for (j = 0; j <= FLT256_LEN; ++j) {
289 pc[j] = 0;
290 }
291 for (ptr = ++q; (c = *q - '0') >= 0 && c <= 9 && pf[0] == 0; --fdex, ++q) {
292 lshift_quad_real (1, pf, m);
293 for (j = 0; j < m; ++j) {
294 pc[j] = pf[j];
295 }
296 lshift_quad_real (2, pf, m);
297 for (n = (unt) c, pa = pf + FLT256_LEN, pb = pc + FLT256_LEN; pa >= pf; pa--, pb--) {
298 n += *pa + *pb;
299 *pa = n;
300 n >>= 16;
301 }
302 }
303 for (j = 0; j < m && pf[j] == 0; ++j);
304 if (j == m) {
305 fbex = 0;
306 } else {
307 fbex = QUAD_REAL_BIAS + QUAD_REAL_MAX_P - 1;
308 if (j != 0) {
309 j <<= 4;
310 fbex -= j;
311 lshift_quad_real (j, pf, m);
312 }
313 while (pf[0]) {
314 rshift_quad_real (1, pf, m);
315 ++fbex;
316 }
317 pf[0] = fbex | sfg;
318 }
319 nofp = ptr == q;
320 } else {
321 nofp = 1;
322 }
323 if (noip && nofp) {
324 return QUAD_REAL_NAN;
325 } else {
326 // Ignore digits beyond precision.
327 for (; *q >= '0' && *q <= '9'; q++);
328 if (endptr != NULL) {
329 *endptr = (char *) q;
330 }
331 }
332 // Exponent
333 if (tolower (*q) == 'e' || tolower (*q) == 'q') {
334 ++q;
335 sfg = 0;
336 if (*q == '+') {
337 ++q;
338 } else if (*q == '-') {
339 sfg = 1;
340 ++q;
341 }
342 for (ptr = q, j = 0; (c = *q - '0') >= 0 && c <= 9 && j <= QUAD_UPB; ++q) {
343 j <<= 1;
344 m = j;
345 j <<= 2;
346 j += c + m;
347 }
348 if (ptr != q && (endptr != 0)) {
349 *endptr = (char *) q;
350 }
351 if (sfg != 0) {
352 j = -j;
353 }
354 idex += j;
355 fdex += j;
356 }
357 // Remark: s and f have the same sign (see above).
358 if (idex > QUAD_REAL_MAX_10EX || fdex > QUAD_REAL_MAX_10EX) {
359 return ((s.value[0] & QUAD_REAL_M_SIGN) ? QUAD_REAL_MINF : QUAD_REAL_PINF);
360 } else {
361 if (idex != 0) {
362 s = mul_quad_real (s, ten_up_quad_real (idex));
363 }
364 if (fdex != 0) {
365 f = mul_quad_real (f, ten_up_quad_real (fdex));
366 }
367 return _add_quad_real_ (s, f);
368 }
369 #undef QUAD_UPB
370 }
371
372 //! @brief atox
373
374 QUAD_T atox (const char *q)
375 {
376 return string_to_quad_real (q, NULL);
377 }
378
379 // _scale_ is set by nP in formats.
380
381 int _scale_ = 1;
382
383 //! @brief srecordf
384
385 int srecordf (char *s, const char *format, ...)
386 {
387 size_t N = BUFFER_SIZE;
388 va_list ap;
389 va_start (ap, format);
390 int vsnprintf (char *, size_t, const char *, va_list);
391 // Print in new string, just in case 's' is also an argument!
392 int M = N + 16; // A bit longer so we trap too long prints.
393 char *t = malloc (M);
394 int Np = vsnprintf (t, M, format, ap);
395 va_end (ap);
396 if (Np >= N) {
397 free (t);
398 QUAD_RTE ("srecordf", "overflow");
399 } else {
400 strcpy (s, t);
401 free (t);
402 }
403 return Np;
404 }
405
406 //! @brief strlcat
407
408 char *strlcat (char *dst, char *src, int len)
409 {
410 // Fortran string lengths are len+1, last one is for null.
411 int avail = len - (int) strlen (dst);
412 if (avail <= (int) strlen (src)) {
413 QUAD_RTE ("_strlcat", "overflow");
414 } else {
415 strncat (dst, src, avail);
416 }
417 dst[len] = '\0';
418 return dst;
419 }
420
421 //! @brief plusab
422
423 static char *plusab (char *buf, char c)
424 {
425 char z[2];
426 z[0] = c;
427 z[1] = '\0';
428 strlcat (buf, z, BUFFER_SIZE);
429 return buf;
430 }
431
432 //! @brief plusto
433
434 static char *plusto (char ch, char *buf)
435 {
436 memmove (&buf[1], &buf[0], strlen(buf) + 1);
437 buf[0] = ch;
438 return buf;
439 }
440
441 //! @brief leading_spaces
442
443 static char *leading_spaces (char *buf, int width)
444 {
445 if (width > 0) {
446 int j = ABS (width) - (int) strlen (buf);
447 while (--j >= 0) {
448 (void) plusto (' ', buf);
449 }
450 }
451 return buf;
452 }
453
454 //! @brief error_chars
455
456 static char *error_chars (char *buf, int n)
457 {
458 int k = (n != 0 ? ABS (n) : 1);
459 buf[k] = '\0';
460 while (--k >= 0) {
461 buf[k] = ERROR_CHAR;
462 }
463 return buf;
464 }
465
466 //! @brief strputc_quad_real
467
468 void strputc_quad_real (char c, char *buffer)
469 {
470 // Assumes that buffer is initially filled with \0.
471 char *ptr;
472 for (ptr = buffer; *ptr != '\0'; ptr++);
473 *ptr = c;
474 }
475
476 //! @brief sprintfmt_quad_real
477
478 void sprintfmt_quad_real (char *buffer, const char *fmt, ...)
479 {
480 BUFFER ibuff;
481 BUFCLR (ibuff);
482 va_list ap;
483 va_start (ap, fmt);
484 vsprintf (ibuff, fmt, ap);
485 va_end (ap);
486 strcat (buffer, ibuff);
487 }
488
489 static int special_value (char *s, QUAD_T u, int sign)
490 {
491 if ((isPinf_quad_real (&u))) {
492 if (sign != 0) {
493 *s++ = '+';
494 }
495 strcpy (s, "Inf");
496 return 1;
497 } else if ((isMinf_quad_real (&u))) {
498 strcpy (s, "-Inf");
499 return 1;
500 } else if ((isNaN_quad_real (&u))) {
501 if (sign != 0) {
502 *s++ = '\?';
503 }
504 strcpy (s, "NaN");
505 return 1;
506 } else {
507 return 0;
508 }
509 }
510
511 //! @brief subfixed_quad_real
512
513 char *subfixed_quad_real (char *buffer, QUAD_T u, int sign, int digs)
514 {
515 BUFCLR (buffer);
516 if ((special_value (buffer, u, sign))) {
517 return buffer;
518 }
519 // Added for VIF to be compatible with fortran F format.
520 QUAD_T t = abs_quad_real (u);
521 while (ge_quad_real (t, QUAD_REAL_TEN)) {
522 t = div_quad_real (t, QUAD_REAL_TEN);
523 digs++;
524 }
525 // -- end of addition.
526 digs = MIN (ABS (digs), FLT256_DIG);
527 // Put sign and take abs value.
528 unt_2 *pa = (unt_2 *) &QV (u);
529 if ((*pa & QUAD_REAL_M_SIGN)) {
530 *pa ^= QUAD_REAL_M_SIGN;
531 strputc_quad_real ('-', buffer);
532 } else if (sign) {
533 strputc_quad_real ('+', buffer);
534 }
535 // Zero as special case.
536 int k;
537 if ((is0_quad_real (&u))) {
538 sprintfmt_quad_real (buffer, "0.");
539 for (k = 0; k < digs; ++k) {
540 strputc_quad_real ('0', buffer);
541 }
542 return buffer;
543 }
544 // Reduce argument.
545 int mag = ((*pa & QUAD_REAL_M_EXP) - QUAD_REAL_BIAS);
546 mag = (int) ((REAL_T) (mag + 1) * M_LOG10_2);
547 if ((mag != 0)) {
548 u = mul_quad_real (u, pwr_quad_real (QUAD_REAL_TEN, -mag));
549 }
550 while ((*pa & QUAD_REAL_M_EXP) < QUAD_REAL_BIAS) {
551 --mag;
552 u = mul_quad_real (u, QUAD_REAL_TEN);
553 }
554 if (mag < 0) {
555 digs = MAX (0, digs + mag);
556 }
557 // We store digits as pseudo-chars.
558 BUFFER q;
559 BUFCLR (q);
560 char *p = q;
561 BUFCLR (q);
562 //
563 int dig;
564 for (*p = 0, k = 0; k <= digs; ++k) {
565 u = sfmod_quad_real (u, &dig);
566 *(++p) = (char) dig;
567 if (*pa == 0) {
568 break;
569 }
570 u = mul_quad_real (QUAD_REAL_TEN, u);
571 }
572 // Round fraction
573 if (*pa != 0) {
574 u = sfmod_quad_real (u, &dig);
575 if (dig >= 5) {
576 ++(*p);
577 }
578 while (*p == 10) {
579 *p = 0;
580 ++(*--p);
581 }
582 }
583 p = q;
584 if (*p == 0) {
585 ++p;
586 } else {
587 ++mag;
588 }
589 // Now sprintf
590 if (mag > QUAD_REAL_MAX_10EX) {
591 sprintfmt_quad_real (buffer, "Inf");
592 return buffer;
593 }
594 //
595 if (mag >= 0) {
596 for (k = 0; k <= mag; k++) {
597 if (k <= digs) {
598 strputc_quad_real ('0' + p[k], buffer);
599 } else {
600 strputc_quad_real ('0', buffer);
601 }
602 }
603 if (k <= digs) {
604 strputc_quad_real ('.', buffer);
605 for (; k <= digs; k++) {
606 strputc_quad_real ('0' + p[k], buffer);
607 }
608 }
609 } else {
610 sprintfmt_quad_real (buffer, "0.");
611 for (k = 1; k < ABS (mag) && digs >= 0; k++) {
612 strputc_quad_real ('0', buffer);
613 digs--;
614 }
615 for (k = 0; k <= digs; ++k) {
616 strputc_quad_real ('0' + *p++, buffer);
617 }
618 }
619 return buffer;
620 }
621
622 //! @brief fixed_quad_real
623
624 char *fixed_quad_real (char *buf, QUAD_T x, int width, int digs, int precision)
625 {
626 width = ABS (width);
627 digs = MIN (abs (digs), precision);
628 subfixed_quad_real (buf, x, 0, digs);
629 if (width > 0 && strlen (buf) > width) {
630 return error_chars (buf, width);
631 } else {
632 return leading_spaces (buf, width);
633 }
634 }
635
636 //! @brief float_quad_real
637
638 char *float_quad_real (char *buf, QUAD_T z, int width, int digs, int expos, int mult, int precision, char sym)
639 {
640 buf[0] = '\0';
641 width = ABS (width);
642 digs = MIN (abs (digs), precision);
643 expos = ABS (expos);
644 if (expos > 5) {
645 return error_chars (buf, width);
646 }
647 // Scientific notation mult = 1, Engineering notation mult = 3
648 mult = MAX (1, mult);
649 // Default _scale_ is 1.
650 int q = 1;
651 char *max = "1";
652 QUAD_T x = abs_quad_real (z), lwb, upb;
653 //
654 if (_scale_ < 0 || _scale_ > 3) {
655 _scale_ = 1;
656 }
657 if (mult == 1) {
658 if (_scale_ == 0) {
659 lwb = QUAD_REAL_TENTH;
660 upb = QUAD_REAL_ONE;
661 q = 1;
662 max = "0.1";
663 } else if (_scale_ == 1) {
664 lwb = QUAD_REAL_ONE;
665 upb = QUAD_REAL_TEN;
666 q = 0;
667 max = "1";
668 } else if (_scale_ == 2) {
669 lwb = QUAD_REAL_TEN;
670 upb = QUAD_REAL_HUNDRED;
671 q = -1;
672 max = "10";
673 } else if (_scale_ == 3) {
674 lwb = QUAD_REAL_HUNDRED;
675 upb = QUAD_REAL_THOUSAND;
676 max = "100";
677 q = -2;
678 }
679 }
680 // Standardize.
681 int p = 0;
682 if (not0_quad_real (&x)) {
683 p = (int) round (quad_real_to_real (log10_quad_real (abs_quad_real(x)))) + q;
684 x = div_quad_real (x, ten_up_quad_real (p));
685 if (le_quad_real (x, lwb)) {
686 x = mul_quad_real (x, QUAD_REAL_TEN);
687 p--;
688 }
689 if (ge_quad_real (x, upb)) {
690 x = div_quad_real (x, QUAD_REAL_TEN);
691 p++;
692 }
693 while (p % mult != 0) {
694 x = mul_quad_real (x, QUAD_REAL_TEN);
695 p--;
696 }
697 }
698 // Form number.
699 BUFFER mant;
700 subfixed_quad_real (mant, x, 0, digs);
701 // Correction of rounding issue by which |mant| equals UPB.
702 if (strchr (mant, '*') == NULL && ge_quad_real (abs_quad_real (string_to_quad_real (mant, NULL)), upb)) {
703 srecordf (mant, "%s", max);
704 if (digs > 0) {
705 plusab (mant, '.');
706 for (int k = 0; k < digs; k++) {
707 plusab (mant, '0');
708 }
709 }
710 p++;
711 }
712 //
713 BUFFER fmt;
714 if (sgn_quad_real (&z) >= 0) {
715 srecordf (fmt, " %%s%c%%+0%dd", sym, expos);
716 } else {
717 srecordf (fmt, "-%%s%c%%+0%dd", sym, expos);
718 }
719 srecordf (buf, fmt, mant, p);
720 if (width > 0 && (strchr (buf, '*') != NULL || strlen (buf) > width)) {
721 if (digs > 0) {
722 return float_quad_real (buf, z, width, digs - 1, expos, mult, precision, sym);
723 } else {
724 return error_chars (buf, width);
725 }
726 } else {
727 return leading_spaces (buf, width);
728 }
729 }
730
731 #endif
