rts-math.c
1 //! @file rts-math.c
2 //! @author J. Marcel van der Veer
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 math operations.
25
26 #include <vif.h>
27
28 // INTEGER*4
29
30 int_4 _up_int_4 (int_4 m, int_4 n)
31 {
32 // Special cases.
33 if (m == 2 && n >= 0) {
34 return 1 << n;
35 } else if (m == 1) {
36 return m;
37 } else if (m == -1) {
38 return (n % 2 == 0 ? 1 : -1);
39 } else if (m == 0) {
40 return (n == 0 ? 1 : 0);
41 } else if (n < 0) {
42 return 0;
43 }
44 // General case.
45 unsigned bit = 1;
46 int_4 M = m, P = 1;
47 do {
48 if (n & bit) {
49 P *= M;
50 }
51 bit <<= 1;
52 if ((int_4) bit <= n) {
53 M *= M;
54 }
55 } while ((int_4) bit <= n);
56 return P;
57 }
58
59 // INTEGER*8
60
61 int_8 _up_int_8 (int_8 m, int_4 n)
62 {
63 // Special cases.
64 if (m == 2 && n >= 0) {
65 int_8 s = 1;
66 s <<= n;
67 return s;
68 } else if (m == 1) {
69 return m;
70 } else if (m == -1) {
71 return (n % 2 == 0 ? 1 : -1);
72 } else if (m == 0) {
73 return (n == 0 ? 1 : 0);
74 } else if (n < 0) {
75 return 0;
76 }
77 // General case.
78 unsigned bit = 1;
79 int_8 M = m, P = 1;
80 do {
81 if (n & bit) {
82 P *= M;
83 }
84 bit <<= 1;
85 if ((int_4) bit <= n) {
86 M *= M;
87 }
88 } while ((int_4) bit <= n);
89 return P;
90 }
91
92 // REAL*4
93
94 real_4 _up_real_4 (real_4 x, int_4 n)
95 {
96 // Only positive n.
97 if (n < 0) {
98 return 1 / _up_real_4 (x, -n);
99 }
100 // Special cases.
101 if (x == 0 || x == 1) {
102 return x;
103 } else if (x == -1) {
104 return (n % 2 == 0 ? 1 : -1);
105 }
106 // General case.
107 unsigned bit = 1;
108 real_4 M = x, P = 1;
109 do {
110 if (n & bit) {
111 P *= M;
112 }
113 bit <<= 1;
114 if ((int_4) bit <= n) {
115 M *= M;
116 }
117 } while ((int_4) bit <= n);
118 return P;
119 }
120
121 real_4 cotanf (real_4 x)
122 {
123 return 1.0 / tanf (x);
124 }
125
126 real_4 acotanf (real_4 x)
127 {
128 return atanf (1 / x);
129 }
130
131 // REAL*8
132
133 real_8 _up_real_8 (real_8 x, int_4 n)
134 {
135 // Only positive n.
136 if (n < 0) {
137 return 1 / _up_real_8 (x, -n);
138 }
139 // Special cases.
140 if (x == 0 || x == 1) {
141 return x;
142 } else if (x == -1) {
143 return (n % 2 == 0 ? 1 : -1);
144 }
145 // General case.
146 unsigned bit = 1;
147 real_8 M = x, P = 1;
148 do {
149 if (n & bit) {
150 P *= M;
151 }
152 bit <<= 1;
153 if ((int_4) bit <= n) {
154 M *= M;
155 }
156 } while ((int_4) bit <= n);
157 return P;
158 }
159
160 real_8 cotan (real_8 x)
161 {
162 return 1.0 / tan (x);
163 }
164
165 real_8 acotan (real_8 x)
166 {
167 return atan (1 / x);
168 }
169
170 real_8 _vif_inverfc (real_8 y)
171 {
172 #define N_c_inverfc 34
173
174 static const real_8 c_inverfc[N_c_inverfc] = {
175 0.91646139826896400000,
176 0.48882664027310800000,
177 0.23172920032340500000,
178 0.12461045461371200000,
179 -0.07288467655856750000,
180 0.26999930867002900000,
181 0.15068904736022300000,
182 0.11606502534161400000,
183 0.49999930343979000000,
184 3.97886080735226000000,
185 0.00112648096188977922,
186 1.05739299623423047e-4,
187 0.00351287146129100025,
188 7.71708358954120939e-4,
189 0.00685649426074558612,
190 0.00339721910367775861,
191 0.01127491693325048700,
192 0.01185981170477711040,
193 0.01429619886978980180,
194 0.03464942077890999220,
195 0.00220995927012179067,
196 0.07434243572417848610,
197 0.10587217794159548800,
198 0.01472979383314851210,
199 0.31684763852013594400,
200 0.71365763586873036400,
201 1.05375024970847138000,
202 1.21448730779995237000,
203 1.16374581931560831000,
204 0.95646497474479900600,
205 0.68626594827409781600,
206 0.43439749233143011500,
207 0.24404451059319093500,
208 0.12078223763524522200
209 };
210
211 if (y < 0 || y > 2) {
212 errno = ERANGE;
213 return 0;
214 }
215 if (y == 0) {
216 return DBL_MAX;
217 } else if (y == 1) {
218 return 0;
219 } else if (y == 2) {
220 return -DBL_MAX;
221 } else {
222 // Next is based on code that originally contained following statement:
223 // Copyright (c) 1996 Takuya Ooura. You may use, copy, modify this
224 // code for any purpose and without fee.
225 real_8 s, t, u, v, x, z;
226 z = (y <= 1 ? y : 2 - y);
227 v = c_inverfc[0] - ln (z);
228 u = sqrt (v);
229 s = (ln (u) + c_inverfc[1]) / v;
230 t = 1 / (u + c_inverfc[2]);
231 x = u *(1 - s *(s *c_inverfc[3] + 0.5)) - ((((c_inverfc[4] *t + c_inverfc[5]) *t + c_inverfc[6]) *t + c_inverfc[7]) *t + c_inverfc[8]) *t;
232 t = c_inverfc[9] / (x + c_inverfc[9]);
233 u = t - 0.5;
234 s = (((((((((c_inverfc[10] *u + c_inverfc[11]) *u - c_inverfc[12]) *u - c_inverfc[13]) *u + c_inverfc[14]) *u + c_inverfc[15]) *u - c_inverfc[16]) *u - c_inverfc[17]) *u + c_inverfc[18]) *u + c_inverfc[19]) *u + c_inverfc[20];
235 s = ((((((((((((s *u - c_inverfc[21]) *u - c_inverfc[22]) *u + c_inverfc[23]) *u + c_inverfc[24]) *u + c_inverfc[25]) *u + c_inverfc[26]) *u + c_inverfc[27]) *u + c_inverfc[28]) *u + c_inverfc[29]) *u + c_inverfc[30]) *u + c_inverfc[31]) *u + c_inverfc[32]) *t - z *exp (x *x - c_inverfc[33]);
236 x += s *(x *s + 1);
237 return (y <= 1 ? x : -x);
238 }
239 }
240
241 real_8 _vif_inverf (real_8 y)
242 {
243 return _vif_inverfc (1 - y);
244 }
245
246 void _merfi (real_8 *p, real_8 *x, int_4 *err)
247 {
248 (*x) = _vif_inverf (*p);
249 (*err) = errno;
250 }
251
252 // REAL*16
253
254 real_16 _qext (real_8 x)
255 {
256 /*
257 NEW_RECORD (str);
258 _srecordf (str, "%.*le", DBL_DIG, x);
259 return _strtoquad (str, NO_REF_TEXT);
260 */
261 return (real_16) x;
262 }
263
264 real_16 _qmod (real_16 x, real_16 y)
265 {
266 real_16 q;
267 if ((q = x / y) >= 0) {
268 q = floorq (q);
269 } else {
270 q = -floorq (-q);
271 }
272 return (x - y *q);
273 }
274
275 real_16 _up_real_16 (real_16 x, int_4 n)
276 {
277 // Only positive n.
278 if (n < 0) {
279 return 1.0q / _up_real_16 (x, -n);
280 }
281 // Special cases.
282 if (x == 0.0q || x == 1.0q) {
283 return x;
284 } else if (x == -1.0q) {
285 return (n % 2 == 0 ? 1.0q : -1.0q);
286 }
287 // General case.
288 unsigned bit = 1;
289 real_16 M = x, P = 1.0q;
290 do {
291 if (n & bit) {
292 P *= M;
293 }
294 bit <<= 1;
295 if ((int_4) bit <= n) {
296 M *= M;
297 }
298 } while ((int_4) bit <= n);
299 return P;
300 }
301
302 real_16 cotanq (real_16 x)
303 {
304 return 1.0q / tanq (x);
305 }
306
307 real_16 acotanq (real_16 x)
308 {
309 return atanq (1 / x);
310 }
311
312 // COMPLEX*8
313
314 complex_8 _cmplxd (complex_16 z)
315 {
316 return CMPLXF ((real_4) creal (z), (real_4) cimag (z));
317 }
318
319 complex_8 _up_complex_8 (complex_8 x, int_4 n)
320 {
321 // Only positive n.
322 if (n < 0) {
323 return 1 / _up_complex (x, -n);
324 }
325 // General case.
326 unsigned bit = 1;
327 complex_8 M = x, P = 1;
328 do {
329 if (n & bit) {
330 P *= M;
331 }
332 bit <<= 1;
333 if ((int_4) bit <= n) {
334 M *= M;
335 }
336 } while ((int_4) bit <= n);
337 return P;
338 }
339
340 // COMPLEX*16
341
342 complex_16 _dcmplxq (complex_32 z)
343 {
344 return CMPLX (crealq (z), cimagq (z));
345 }
346
347 complex_16 _up_complex (complex_16 x, int_4 n)
348 {
349 // Only positive n.
350 if (n < 0) {
351 return 1 / _up_complex (x, -n);
352 }
353 // General case.
354 unsigned bit = 1;
355 complex_16 M = x, P = 1;
356 do {
357 if (n & bit) {
358 P *= M;
359 }
360 bit <<= 1;
361 if ((int_4) bit <= n) {
362 M *= M;
363 }
364 } while ((int_4) bit <= n);
365 return P;
366 }
367
368 // COMPLEX*32
369
370 complex_32 _qcmplxd (complex_16 z)
371 {
372 return CMPLXQ (_qext (creal (z)), _qext (cimag (z)));
373 }
374
375 complex_32 _up_complex_32 (complex_32 x, int_4 n)
376 {
377 // Only positive n.
378 if (n < 0) {
379 return 1.0q / _up_complex_32 (x, -n);
380 }
381 // General case.
382 unsigned bit = 1;
383 complex_32 M = x, P = 1.0q;
384 do {
385 if (n & bit) {
386 P *= M;
387 }
388 bit <<= 1;
389 if ((int_4) bit <= n) {
390 M *= M;
391 }
392 } while ((int_4) bit <= n);
393 return P;
394 }
395
396 real_4 _zabs_8 (real_4 re, real_4 im)
397 {
398 return cabsf (CMPLXF (re, im));
399 }
400
401 real_8 _zabs_16 (real_8 re, real_8 im)
402 {
403 return cabs (CMPLX (re, im));
404 }
405
406 real_16 _zabs_32 (real_16 re, real_16 im)
407 {
408 return cabsq (CMPLXQ (re, im));
409 }
410
411 void _qhex (real_16 *u)
412 {
413 unt_2 *p = (unt_2 *) u;
414 fprintf (stdout, "{");
415 for (int_4 i = 0; i <= FLT128_LEN; i++) {
416 fprintf (stdout, "0x%04x", *(p++));
417 if (i < FLT128_LEN) {
418 fprintf (stdout, ", ");
419 }
420 }
421 fprintf (stdout, "}\n");
422 }
423
424 void _xhex (real_32 *u)
425 {
426 fprintf (stdout, "{");
427 for (int_4 i = 0; i <= FLT256_LEN; i++) {
428 fprintf (stdout, "0x%04x", u->value[i]);
429 if (i < FLT256_LEN) {
430 fprintf (stdout, ", ");
431 }
432 }
433 fprintf (stdout, "}\n");
434 }
435
436 //
437
438 void _pi4 (real_4 * x)
439 {
440 *x = (real_4) M_PI;
441 }
442
443 void _pi8 (real_8 * x)
444 {
445 *x = M_PI;
446 }
447
448 void _pi16 (real_16 * x)
449 {
450 *x = M_PIq;
451 }
|
© 2002-2025 J.M. van der Veer (jmvdveer@xs4all.nl)
|
