rts-real32.h
1 //! @file rts-real32.h
2 //! @author (see below)
3 //
4 //! @section Copyright
5 //
6 // This file is part of VIF - vintage FORTRAN compiler.
7 // Copyright 2020-2025 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 //! Runtime support for REAL*32 and COMPLEX*64.
25
26 // The code is based on the HPA Library, available from:
27 // <http://download-mirror.savannah.gnu.org/releases/hpalib/>
28 //
29 // Copyright (C) 2000 Daniel A. Atkinson <DanAtk@aol.com>
30 // Copyright (C) 2004 Ivano Primi <ivprimi@libero.it>
31 // Copyright (C) 2022 Marcel van der Veer <algol68g@xs4all.nl> - VIF version.
32 //
33 // The HPA Library is free software; you can redistribute it and/or
34 // modify it under the terms of the GNU Lesser General Public
35 // License as published by the Free Software Foundation; either
36 // version 2.1 of the License, or (at your option) any later version.
37 //
38 // The HPA Library is distributed in the hope that it will be useful,
39 // but WITHOUT ANY WARRANTY; without even the implied warranty of
40 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
41 // Lesser General Public License for more details.
42 //
43 // You should have received a copy of the GNU Lesser General Public
44 // License along with the HPA Library; if not, write to the Free
45 // Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA
46 // 02110-1301 USA.
47
48 #include <vif.h>
49
50 #define CX1I_CHAR 'i'
51 #define CX1I_STR "i"
52 #define CXDEF_LDEL '('
53 #define CXDEF_RDEL ')'
54 #define CX_EMPTY_SEP " "
55 #define CX_SEPARATOR ", "
56 #define CX_SEP_L (strlen (CX_SEPARATOR))
57 #define XDEF_LIM 6
58 #define XFMT_ALT 2
59 #define XFMT_RAW 1
60 #define XFMT_STD 0
61 #define XOUT_FIXED 0
62 #define XOUT_SCIENTIFIC 1
63
64 #define HPA_VERSION "1.7 VIF"
65 #define XERR_DFL 1
66 #define XLITTLE_ENDIAN 1
67 #define XMAX_10EX 4931
68 #define XMAX_DEGREE 50
69 #define XULONG_BITSIZE 32
70
71 #define cxconvert cxconv
72 #define CXCONV(x) (complex_64){x, X_0}
73 #define cxdiff cxsub
74 #define CXIM(z) (z).im
75 #define cxipow cxpwr
76 #define CXRE(z) (z).re
77 #define CXSWAP(z) (complex_64){(z).im, (z).re}
78
79 struct xoutflags
80 {
81 int_2 fmt, notat, sf, mfwd, lim;
82 signed char padding, ldel, rdel;
83 };
84
85 #if defined (XERR_IGN)
86 #define xsigerr(errcond, errcode, where) 0
87 #else
88 #define XENONE 0
89 #define XEDIV 1
90 #define XEDOM 2
91 #define XEBADEXP 3
92 #define XFPOFLOW 4
93 #define XNERR 4
94 #define XEINV (XNERR + 1)
95 int_4 xsigerr (int_4 errcond, int_4 errcode, const char *where);
96 #endif
97
98 #if (FLT256_LEN == 15)
99 static const int_4 xItt_div = 3;
100 static const int_4 xK_tanh = 6;
101 static const int_4 xMS_exp = 39;
102 static const int_4 xMS_hyp = 45;
103 static const int_4 xMS_trg = 55;
104 #else
105 #error invalid real*32 length
106 #endif
107
108 static const int_2 xD_bias = 15360;
109 static const int_2 xD_lex = 12;
110 static const int_2 xD_max = 2047;
111 static const int_2 X_EXPO_BIAS = 16383;
112 static const int_2 xF_bias = 16256;
113 static const int_2 xF_lex = 9;
114 static const int_2 xF_max = 255;
115 static const int_2 xK_lin = -8 * FLT256_LEN;
116 static const int_2 xMax_p = 16 * FLT256_LEN;
117 static const real_32 xE2max = {{0x400c, 0xfffb}}; // +16382.75
118 static const real_32 xE2min = {{0xc00c, 0xfffb}}; // -16382.75
119 static const real_32 xEmax = {{0x400c, 0xb16c}};
120 static const real_32 xEmin = {{0xc00c, 0xb16c}};
121 static const real_32 X_MINUS_INF = {{0xffff, 0x0}};
122 static const real_32 X_PLUS_INF = {{0x7fff, 0x0}};
123 static const real_32 X_VGV = {{0x4013, 0x8000}};
124 static const real_32 X_VSV = {{0x3ff2, 0x8000}};
125 static const unt_2 X_EXPO_MASK = 0x7fff;
126 static const unt_2 X_SIGN_MASK = 0x8000;
127
128 // List of constants, extended with respect to original HPA lib.
129
130 static const real_32 FLT256_MIN = // Minimum representable number
131 {{0x0001, 0x8000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000}};
132 static const real_32 FLT256_MAX = // Maximum representable number
133 {{0x7ffe, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff}};
134 static const real_32 FLT256_EPSILON = // Minimum positive number such that 1 + FLT256_EPSILON != 1.
135 {{0x3f11, 0x8000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000}};
136 static const real_32 FLT256_EPSILON_HALF = // Maximum positive number such that 1 + FLT256_EPSILON != 1.
137 {{0x3f10, 0x8000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000}};
138 static const real_32 X_LOG10_2 =
139 {{0x3ffd, 0x9a20, 0x9a84, 0xfbcf, 0xf798, 0x8f89, 0x59ac, 0x0b7c, 0x9178, 0x26ad, 0x30c5, 0x43d1, 0xf349, 0x8a5e, 0x6f26, 0xb7d3}};
140
141 static const real_32 X_MINUS_1 = // -1
142 {{0xbfff, 0x8000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000}};
143 static const real_32 X_0 =
144 {{0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000}};
145 static const real_32 X_1_OVER_10 =
146 {{0x3ffb, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0xcccc}};
147 static const real_32 X_1_OVER_2 =
148 {{0x3ffe, 0x8000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000}};
149 static const real_32 X_1 =
150 {{0x3fff, 0x8000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000}};
151 static const real_32 X_2 =
152 {{0x4000, 0x8000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000}};
153 static const real_32 X_3 =
154 {{0x4000, 0xc000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000}};
155 static const real_32 X_4 =
156 {{0x4001, 0x8000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000}} ;
157 static const real_32 X_5 =
158 {{0x4001, 0xa000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000}} ;
159 static const real_32 X_6 =
160 {{0x4001, 0xc000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000}} ;
161 static const real_32 X_7 =
162 {{0x4001, 0xe000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000}} ;
163 static const real_32 X_8 =
164 {{0x4002, 0x8000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000}} ;
165 static const real_32 X_9 =
166 {{0x4002, 0x9000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000}} ;
167 static const real_32 X_10 =
168 {{0x4002, 0xa000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000}};
169 static const real_32 X_100 =
170 {{0x4005, 0xc800, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000}};
171 static const real_32 X_1000=
172 {{0x4008, 0xfa00, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000}};
173
174 static const real_32 X_PI_OVER_4 =
175 {{0x3ffe, 0xc90f, 0xdaa2, 0x2168, 0xc234, 0xc4c6, 0x628b, 0x80dc, 0x1cd1, 0x2902, 0x4e08, 0x8a67, 0xcc74, 0x020b, 0xbea6, 0x3b14}};
176 static const real_32 X_PI_OVER_2 =
177 {{0x3fff, 0xc90f, 0xdaa2, 0x2168, 0xc234, 0xc4c6, 0x628b, 0x80dc, 0x1cd1, 0x2902, 0x4e08, 0x8a67, 0xcc74, 0x020b, 0xbea6, 0x3b14}};
178 static const real_32 X_PI =
179 {{0x4000, 0xc90f, 0xdaa2, 0x2168, 0xc234, 0xc4c6, 0x628b, 0x80dc, 0x1cd1, 0x2902, 0x4e08, 0x8a67, 0xcc74, 0x020b, 0xbea6, 0x3b14}};
180 static const real_32 X_LN_2 =
181 {{0x3ffe, 0xb172, 0x17f7, 0xd1cf, 0x79ab, 0xc9e3, 0xb398, 0x03f2, 0xf6af, 0x40f3, 0x4326, 0x7298, 0xb62d, 0x8a0d, 0x175b, 0x8bab}};
182 static const real_32 X_SQRT_2 =
183 {{0x3fff, 0xb504, 0xf333, 0xf9de, 0x6484, 0x597d, 0x89b3, 0x754a, 0xbe9f, 0x1d6f, 0x60ba, 0x893b, 0xa84c, 0xed17, 0xac85, 0x8334}};
184 static const real_32 X_LOG2_E =
185 {{0x3fff, 0xb8aa, 0x3b29, 0x5c17, 0xf0bb, 0xbe87, 0xfed0, 0x691d, 0x3e88, 0xeb57, 0x7aa8, 0xdd69, 0x5a58, 0x8b25, 0x166c, 0xd1a1}};
186 static const real_32 X_LOG2_10 =
187 {{0x4000, 0xd49a, 0x784b, 0xcd1b, 0x8afe, 0x492b, 0xf6ff, 0x4daf, 0xdb4c, 0xd96c, 0x55fe, 0x37b3, 0xad4e, 0x91b6, 0xac80, 0x82e8}};
188 static const real_32 X_LOG10_E =
189 {{0x3ffd, 0xde5b, 0xd8a9, 0x3728, 0x7195, 0x355b, 0xaaaf, 0xad33, 0xdc32, 0x3ee3, 0x4602, 0x45c9, 0xa202, 0x3a3f, 0x2d44, 0xf78f}};
190 static const real_32 X_RND_CORR =
191 {{0x3ffe, 0x8000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x01ae}};
192 static const real_32 X_FIX_CORR =
193 {{0x3f17, 0xc000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000}};
194 static const real_32 X_NAN =
195 {{0x0000, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff}};
196
197 static const complex_64 X_0_0I = {
198 {{0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000}},
199 {{0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000}}
200 };
201
202 static const complex_64 X_1_0I = {
203 {{0x3fff, 0x8000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000}},
204 {{0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000}}
205 };
206
207 static const complex_64 X_0_1I = {
208 {{0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000}},
209 {{0x3fff, 0x8000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000}}
210 };
211
212 extern complex_64 atocx (const char *);
213 extern complex_64 cxacos (complex_64);
214 extern complex_64 cxacosh (complex_64);
215 extern complex_64 cxadd (complex_64, complex_64, int_4);
216 extern complex_64 cxasin (complex_64);
217 extern complex_64 cxasinh (complex_64);
218 extern complex_64 cxatan (complex_64);
219 extern complex_64 cxatanh (complex_64);
220 extern complex_64 cxceil (complex_64);
221 extern complex_64 cxconj (complex_64);
222 extern complex_64 cxcos (complex_64);
223 extern complex_64 cxcosh (complex_64);
224 extern complex_64 cxdiv (complex_64, complex_64);
225 extern complex_64 cxdrot (complex_64);
226 extern complex_64 cxexp10 (complex_64);
227 extern complex_64 cxexp2 (complex_64);
228 extern complex_64 cxexp (complex_64);
229 extern complex_64 cxfix (complex_64);
230 extern complex_64 cxfloor (complex_64);
231 extern complex_64 cxfrac (complex_64);
232 extern complex_64 cxgdiv (complex_64, complex_64);
233 extern complex_64 cxgmod (complex_64, complex_64);
234 extern complex_64 cxidiv (complex_64, complex_64);
235 extern complex_64 cxinv (complex_64);
236 extern complex_64 cxlog10 (complex_64);
237 extern complex_64 cxlog2 (complex_64);
238 extern complex_64 cxlog (complex_64);
239 extern complex_64 cxlog_sqrt (complex_64);
240 extern complex_64 cxmod (complex_64, complex_64);
241 extern complex_64 cxmul (complex_64, complex_64);
242 extern complex_64 cxneg (complex_64);
243 extern complex_64 cxpow (complex_64, complex_64);
244 extern complex_64 cxpwr (complex_64, int_4);
245 extern complex_64 cxreset (real_32, real_32);
246 extern complex_64 cxrmul (real_32, complex_64);
247 extern complex_64 cxroot (complex_64, int_4, int_4);
248 extern complex_64 cxround (complex_64);
249 extern complex_64 cxrrot (complex_64);
250 extern complex_64 cxsin (complex_64);
251 extern complex_64 cxsinh (complex_64);
252 extern complex_64 cxsqr (complex_64);
253 extern complex_64 cxsqrt (complex_64);
254 extern complex_64 cxsub (complex_64, complex_64);
255 extern complex_64 cxsum (complex_64, complex_64);
256 extern complex_64 cxswap (complex_64);
257 extern complex_64 cxtan (complex_64);
258 extern complex_64 cxtanh (complex_64);
259 extern complex_64 cxtrunc (complex_64);
260 extern complex_64 dctocx (real_8, real_8);
261 extern complex_64 fctocx (real_4, real_4);
262 extern complex_64 ictocx (int_4, int_4);
263 extern complex_64 uctocx (unt_4, unt_4);
264
265 extern int_4 xgetexp (const real_32 *);
266 extern int_4 xgetsgn (const real_32 *);
267 extern int_4 xreal_cmp (const real_32 *, const real_32 *);
268
269 extern logical_4 cxeq (complex_64, complex_64);
270 extern logical_4 cxis0 (const complex_64*);
271 extern logical_4 cxneq (complex_64, complex_64);
272 extern logical_4 cxnot0 (const complex_64*);
273 extern logical_4 cxrec (complex_64, complex_64 *);
274 extern logical_4 xisordnumb (const real_32 *);
275
276 extern real_32 cxabs (complex_64);
277 extern real_32 cxarg (complex_64);
278 extern real_32 xacosh (real_32);
279 extern real_32 xacos (real_32);
280 extern real_32 xacotan (real_32);
281 extern real_32 _xaint_4 (real_32);
282 extern real_32 xasinh (real_32);
283 extern real_32 xasin (real_32);
284 extern real_32 xatan2 (real_32, real_32);
285 extern real_32 xatanh (real_32);
286 extern real_32 xatan (real_32);
287 extern real_32 xceil (real_32);
288 extern real_32* xchcof (int_4, real_32 (*) (real_32));
289 extern real_32 xcosh (real_32);
290 extern real_32 xcos (real_32);
291 extern real_32 xcotan (real_32);
292 extern real_32 xevtch (real_32, real_32 *, int_4);
293 extern real_32 xexp10 (real_32);
294 extern real_32 xexp2 (real_32);
295 extern real_32 xexp (real_32);
296 extern real_32 xfix (real_32);
297 extern real_32 xfloor (real_32);
298 extern real_32 xfmod (real_32, real_32, real_32 *);
299 extern real_32 xfrexp (real_32, int_4 *);
300 extern real_32 xlog2 (real_32);
301 extern real_32 xlog (real_32);
302 extern real_32 _xmax (real_32, real_32);
303 extern real_32 _xmin (real_32, real_32);
304 extern real_32 _xnint_4 (real_32);
305 extern real_32 xpow (real_32, real_32);
306 extern real_32 xreal_2 (real_32, int_4);
307 extern real_32 _xsgn (real_32, real_32);
308 extern real_32 xsinh (real_32);
309 extern real_32 xsin (real_32);
310 extern real_32 xsqrt (real_32);
311 extern real_32 xtanh (real_32);
312 extern real_32 xtan (real_32);
313 extern void cxtodc (const complex_64 *, real_8 *, real_8 *);
314 extern void cxtofc (const complex_64 *, real_4 *, real_4 *);
315 extern void _xdump (real_32 *);
316 extern void xlshift (int_4, unt_2 *, int_4);
317 extern void xrshift (int_4, unt_2 *, int_4);
318
319 static real_32 c_tan (real_32 u);
320 static real_32 rred (real_32 u, int_4 i, int_4 *p);
|
© 2002-2025 J.M. van der Veer (jmvdveer@xs4all.nl)
|
