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-2024 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 _aint (real_4 x)
95 {
96 return (x >= 0 ? floorf (x) : - floorf (- x));
97 }
98
99 real_4 _anint (real_4 x)
100 {
101 return (x >= 0 ? floorf (0.5 + x) : - floorf (0.5 - x));
102 }
103
104 real_4 _up_real_4 (real_4 x, int_4 n)
105 {
106 // Only positive n.
107 if (n < 0) {
108 return 1 / _up_real_4 (x, -n);
109 }
110 // Special cases.
111 if (x == 0 || x == 1) {
112 return x;
113 } else if (x == -1) {
114 return (n % 2 == 0 ? 1 : -1);
115 }
116 // General case.
117 unsigned bit = 1;
118 real_4 M = x, P = 1;
119 do {
120 if (n & bit) {
121 P *= M;
122 }
123 bit <<= 1;
124 if ((int_4) bit <= n) {
125 M *= M;
126 }
127 } while ((int_4) bit <= n);
128 return P;
129 }
130
131 real_4 cotanf (real_4 x)
132 {
133 return 1.0 / tanf (x);
134 }
135
136 // REAL*8
137
138 real_8 _dint (real_8 x)
139 {
140 return (x >= 0 ? floor (x) : - floor (- x));
141 }
142
143 real_8 _dnint (real_8 x)
144 {
145 return (x >= 0 ? floor (0.5 + x) : - floor (0.5 - x));
146 }
147
148 real_8 _up_real_8 (real_8 x, int_4 n)
149 {
150 // Only positive n.
151 if (n < 0) {
152 return 1 / _up_real_8 (x, -n);
153 }
154 // Special cases.
155 if (x == 0 || x == 1) {
156 return x;
157 } else if (x == -1) {
158 return (n % 2 == 0 ? 1 : -1);
159 }
160 // General case.
161 unsigned bit = 1;
162 real_8 M = x, P = 1;
163 do {
164 if (n & bit) {
165 P *= M;
166 }
167 bit <<= 1;
168 if ((int_4) bit <= n) {
169 M *= M;
170 }
171 } while ((int_4) bit <= n);
172 return P;
173 }
174
175 real_8 cotan (real_8 x)
176 {
177 return 1.0 / tan (x);
178 }
179
180 real_8 _vif_inverfc (real_8 y)
181 {
182 #define N_c_inverfc 34
183
184 static const real_8 c_inverfc[N_c_inverfc] = {
185 0.91646139826896400000,
186 0.48882664027310800000,
187 0.23172920032340500000,
188 0.12461045461371200000,
189 -0.07288467655856750000,
190 0.26999930867002900000,
191 0.15068904736022300000,
192 0.11606502534161400000,
193 0.49999930343979000000,
194 3.97886080735226000000,
195 0.00112648096188977922,
196 1.05739299623423047e-4,
197 0.00351287146129100025,
198 7.71708358954120939e-4,
199 0.00685649426074558612,
200 0.00339721910367775861,
201 0.01127491693325048700,
202 0.01185981170477711040,
203 0.01429619886978980180,
204 0.03464942077890999220,
205 0.00220995927012179067,
206 0.07434243572417848610,
207 0.10587217794159548800,
208 0.01472979383314851210,
209 0.31684763852013594400,
210 0.71365763586873036400,
211 1.05375024970847138000,
212 1.21448730779995237000,
213 1.16374581931560831000,
214 0.95646497474479900600,
215 0.68626594827409781600,
216 0.43439749233143011500,
217 0.24404451059319093500,
218 0.12078223763524522200
219 };
220
221 if (y < 0 || y > 2) {
222 errno = ERANGE;
223 return 0;
224 }
225 if (y == 0) {
226 return DBL_MAX;
227 } else if (y == 1) {
228 return 0;
229 } else if (y == 2) {
230 return -DBL_MAX;
231 } else {
232 // Next is based on code that originally contained following statement:
233 // Copyright (c) 1996 Takuya Ooura. You may use, copy, modify this
234 // code for any purpose and without fee.
235 real_8 s, t, u, v, x, z;
236 z = (y <= 1 ? y : 2 - y);
237 v = c_inverfc[0] - ln (z);
238 u = sqrt (v);
239 s = (ln (u) + c_inverfc[1]) / v;
240 t = 1 / (u + c_inverfc[2]);
241 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;
242 t = c_inverfc[9] / (x + c_inverfc[9]);
243 u = t - 0.5;
244 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];
245 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]);
246 x += s *(x *s + 1);
247 return (y <= 1 ? x : -x);
248 }
249 }
250
251 real_8 _vif_inverf (real_8 y)
252 {
253 return _vif_inverfc (1 - y);
254 }
255
256 void _merfi (real_8 *p, real_8 *x, int_4 *err)
257 {
258 (*x) = _vif_inverf (*p);
259 (*err) = errno;
260 }
261
262 // REAL*16
263
264 real_16 _qint (real_16 x)
265 {
266 return (x >= 0 ? floorq (x) : - floorq (- x));
267 }
268
269 real_16 _qext (real_8 x)
270 {
271 RECORD str;
272 _srecordf (str, "%.*le", DBL_DIG, x);
273 return _strtoquad (str, NULL);
274 }
275
276 real_16 _qmod (real_16 x, real_16 y)
277 {
278 real_16 q;
279 if ((q = x / y) >= 0) {
280 q = floorq (q);
281 } else {
282 q = -floorq (-q);
283 }
284 return (x - y *q);
285 }
286
287 real_16 _up_real_16 (real_16 x, int_4 n)
288 {
289 // Only positive n.
290 if (n < 0) {
291 return 1.0q / _up_real_16 (x, -n);
292 }
293 // Special cases.
294 if (x == 0.0q || x == 1.0q) {
295 return x;
296 } else if (x == -1.0q) {
297 return (n % 2 == 0 ? 1.0q : -1.0q);
298 }
299 // General case.
300 unsigned bit = 1;
301 real_16 M = x, P = 1.0q;
302 do {
303 if (n & bit) {
304 P *= M;
305 }
306 bit <<= 1;
307 if ((int_4) bit <= n) {
308 M *= M;
309 }
310 } while ((int_4) bit <= n);
311 return P;
312 }
313
314 real_16 cotanq (real_16 x)
315 {
316 return 1.0q / tanq (x);
317 }
318
319 // COMPLEX*8
320
321 complex_8 _cmplxd (complex_16 z)
322 {
323 return CMPLXF ((real_4) creal (z), (real_4) cimag (z));
324 }
325
326 complex_8 _up_complex_8 (complex_8 x, int_4 n)
327 {
328 // Only positive n.
329 if (n < 0) {
330 return 1 / _up_complex (x, -n);
331 }
332 // General case.
333 unsigned bit = 1;
334 complex_8 M = x, P = 1;
335 do {
336 if (n & bit) {
337 P *= M;
338 }
339 bit <<= 1;
340 if ((int_4) bit <= n) {
341 M *= M;
342 }
343 } while ((int_4) bit <= n);
344 return P;
345 }
346
347 // COMPLEX*16
348
349 complex_16 _dcmplxq (complex_32 z)
350 {
351 return CMPLX (crealq (z), cimagq (z));
352 }
353
354 complex_16 _up_complex (complex_16 x, int_4 n)
355 {
356 // Only positive n.
357 if (n < 0) {
358 return 1 / _up_complex (x, -n);
359 }
360 // General case.
361 unsigned bit = 1;
362 complex_16 M = x, P = 1;
363 do {
364 if (n & bit) {
365 P *= M;
366 }
367 bit <<= 1;
368 if ((int_4) bit <= n) {
369 M *= M;
370 }
371 } while ((int_4) bit <= n);
372 return P;
373 }
374
375 // COMPLEX*32
376
377 complex_32 _qcmplxd (complex_16 z)
378 {
379 return CMPLXQ (_qext (creal (z)), _qext (cimag (z)));
380 }
381
382 complex_32 _up_complex_32 (complex_32 x, int_4 n)
383 {
384 // Only positive n.
385 if (n < 0) {
386 return 1.0q / _up_complex_32 (x, -n);
387 }
388 // General case.
389 unsigned bit = 1;
390 complex_32 M = x, P = 1.0q;
391 do {
392 if (n & bit) {
393 P *= M;
394 }
395 bit <<= 1;
396 if ((int_4) bit <= n) {
397 M *= M;
398 }
399 } while ((int_4) bit <= n);
400 return P;
401 }
402
403 real_4 _zabs_8 (real_4 re, real_4 im)
404 {
405 return cabsf (CMPLXF (re, im));
406 }
407
408 real_8 _zabs_16 (real_8 re, real_8 im)
409 {
410 return cabs (CMPLX (re, im));
411 }
412
413 real_16 _zabs_32 (real_16 re, real_16 im)
414 {
415 return cabsq (CMPLXQ (re, im));
416 }
417
418 void _qhex (real_16 *u)
419 {
420 unt_2 *p = (unt_2 *) u;
421 fprintf (stdout, "{");
422 for (int_4 i = 0; i <= FLT128_LEN; i++) {
423 fprintf (stdout, "0x%04x", *(p++));
424 if (i < FLT128_LEN) {
425 fprintf (stdout, ", ");
426 }
427 }
428 fprintf (stdout, "}\n");
429 }
430
431 void _xhex (real_32 *u)
432 {
433 fprintf (stdout, "{");
434 for (int_4 i = 0; i <= FLT256_LEN; i++) {
435 fprintf (stdout, "0x%04x", u->value[i]);
436 if (i < FLT256_LEN) {
437 fprintf (stdout, ", ");
438 }
439 }
440 fprintf (stdout, "}\n");
441 }
442
443 int_4 _int4 (char *str)
444 {
445 int_4 sum = 0, fact = 1, len = strlen (str);
446 for (int_4 k = 0; k < 4 && k < len; k++) {
447 sum += fact * (int_4) str[k];
448 fact <<= 8;
449 }
450 return sum;
451 }
452
453 //
454
455 void _pi4 (real_4 * x)
456 {
457 *x = (real_4) M_PI;
458 }
459
460 void _pi8 (real_8 * x)
461 {
462 *x = M_PI;
463 }
464
465 void _pi16 (real_16 * x)
466 {
467 *x = M_PIq;
468 }
|
© 2002-2024 J.M. van der Veer (jmvdveer@xs4all.nl)
|
