mp-quad.c
1 //! @file mp-quad.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 //! Fixed-length LONG LONG REAL and COMPLEX 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-mp.h"
33 #include "a68g-double.h"
34 #include "a68g-quad.h"
35
36 //! @brief Convert mp to quad real.
37
38 QUAD_T mp_to_quad_real (NODE_T * p, MP_T * z, int digits)
39 {
40 // This routine looks a lot like "strtod".
41 (void) p;
42 QUAD_T u = QUAD_REAL_ZERO;
43 if (MP_EXPONENT (z) * (MP_T) LOG_MP_RADIX > (MP_T) REAL_MIN_10_EXP) {
44 BOOL_T neg = MP_DIGIT (z, 1) < 0;
45 MP_DIGIT (z, 1) = ABS (MP_DIGIT (z, 1));
46 int expo = (int) (MP_EXPONENT (z) * LOG_MP_RADIX);
47 QUAD_T w = ten_up_quad_real (expo), r = double_real_to_quad_real (MP_RADIX);
48 for (int j = 1; j <= digits; j++) {
49 QUAD_T term = double_real_to_quad_real (MP_DIGIT (z, j));
50 term = mul_quad_real (term, w);
51 u = _add_quad_real_ (u, term);
52 w = div_quad_real (w, r);
53 }
54 if (neg) {
55 u = neg_quad_real (u);
56 }
57 }
58 return u;
59 }
60
61 //! @brief Convert quad real to mp number.
62
63 void quad_real_to_mp (NODE_T * p, MP_T *z, QUAD_T x, int digits)
64 {
65 SET_MP_ZERO (z, digits);
66 if (is0_quad_real (&x)) {
67 return;
68 }
69 int sign_x = getsgn_quad_real (&x);
70 // Scale to [0, 0.1>.
71 // a = ABS (x);
72 QUAD_T u = abs_quad_real (x);
73 // expo = (int) log10q (a);
74 QUAD_T v = log10_quad_real (u);
75 INT_T expo = (INT_T) quad_real_to_double_real (v);
76 // v /= ten_up_double (expo);
77 v = ten_up_quad_real (expo);
78 u = div_quad_real (u, v);
79 expo--;
80 if (real_cmp_quad_real (&u, &QUAD_REAL_ONE) >= 0) {
81 u = div_quad_real (u, QUAD_REAL_TEN);
82 expo++;
83 }
84 // Transport digits of x to the mantissa of z.
85 INT_T sum = 0, W = (MP_RADIX / 10); int j, k;
86 for (k = 0, j = 1; j <= digits && k < QUAD_DIGITS; k++) {
87 QUAD_T t = mul_quad_real (u, QUAD_REAL_TEN);
88 v = floor_quad_real (t);
89 u = frac_quad_real (t);
90 sum += (INT_T) (W * quad_real_to_double_real (v));
91 W /= 10;
92 if (W < 1) {
93 MP_DIGIT (z, j++) = (MP_T) sum;
94 sum = 0;
95 W = (MP_RADIX / 10);
96 }
97 }
98 // Store the last digits.
99 if (j <= digits) {
100 MP_DIGIT (z, j) = (MP_T) sum;
101 }
102 (void) align_mp (z, &expo, digits);
103 MP_EXPONENT (z) = (MP_T) expo;
104 if (sign_x != 0) {
105 MP_DIGIT (z, 1) = -MP_DIGIT (z, 1);
106 }
107 check_mp_exp (p, z);
108 }
109
110 //! @brief PROC quad mp = (LONG LONG REAL) LONG LONG REAL
111
112 void genie_quad_mp (NODE_T * p)
113 {
114 MOID_T *mode = MOID (p);
115 int digits = DIGITS (mode);
116 int size = SIZE (mode);
117 MP_T *z = (MP_T *) STACK_OFFSET (-size);
118 QUAD_T u = mp_to_quad_real (p, z, digits);
119 quad_real_to_mp (p, z, u, digits);
120 }
121
122 #endif
