quad.c
1 //! @file quad.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 library.
26
27 // A68G has an own LONG LONG REAL library with variable precision.
28 // This fixed-length library serves as more economical double precision
29 // to double real, for instance in computing the incomplete gamma function.
30 //
31 // This code is based on the HPA Library, a branch of the CCMath library.
32 // The CCMath library and derived HPA Library are free software under the terms
33 // of the GNU Lesser General Public License version 2.1 or any later version.
34 //
35 // Sources:
36 // [https://directory.fsf.org/wiki/HPAlib]
37 // [http://download-mirror.savannah.gnu.org/releases/hpalib/]
38 // [http://savannah.nongnu.org/projects/hpalib]
39 //
40 // Copyright 2000 Daniel A. Atkinson [DanAtk@aol.com]
41 // Copyright 2004 Ivano Primi [ivprimi@libero.it]
42 // Copyright 2023 Marcel van der Veer [algol68g@xs4all.nl] - A68G version.
43
44 // IEEE 754 floating point standard is assumed.
45 // A quad real number is represented as:
46 //
47 // Sign bit(s): 0 -> positive, 1 -> negative.
48 // Exponent(e): 15-bit biased integer (bias = 16383).
49 // Mantissa(m): 15 words of 16 bit length including leading 1.
50 //
51 // The range of representable numbers is then given by
52 //
53 // 2^16384 > x > 2^[-16383]
54 // 1.19*10^4932 > x > 1.68*10^-[4932]
55 //
56 // Special values of the exponent are:
57 //
58 // All 1 -> infinity (floating point overflow).
59 // All 0 -> number equals zero.
60 //
61 // Underflow in operations is handled by a flush to zero. Thus, a number
62 // with the exponent zero and nonzero mantissa is invalid (not-a-number).
63 // A complex number is a structure formed by two REAL*32 numbers.
64 //
65 // HPA cannot extend precision beyond the preset, hardcoded precision.
66 // Hence some math routines will not achieve full precision.
67
68 #include "a68g.h"
69
70 #if (A68_LEVEL >= 3)
71
72 #include "a68g-quad.h"
73
74 // Constants, extended with respect to original HPA lib.
75
76 const int_2 QUAD_REAL_BIAS = 16383;
77 const int_2 QUAD_REAL_DBL_BIAS = 15360;
78 const int_2 QUAD_REAL_DBL_LEX = 12;
79 const int_2 QUAD_REAL_DBL_MAX = 2047;
80 const int_2 QUAD_REAL_K_LIN = -8 * FLT256_LEN;
81 const int_2 QUAD_REAL_MAX_P = 16 * FLT256_LEN;
82 const QUAD_T QUAD_REAL_E2MAX = {{0x400c, 0xfffb}}; // +16382.75
83 const QUAD_T QUAD_REAL_E2MIN = {{0xc00c, 0xfffb}}; // -16382.75
84 const QUAD_T QUAD_REAL_EMAX = {{0x400c, 0xb16c}};
85 const QUAD_T QUAD_REAL_EMIN = {{0xc00c, 0xb16c}};
86 const QUAD_T QUAD_REAL_MINF = {{0xffff, 0x0}};
87 const QUAD_T QUAD_REAL_PINF = {{0x7fff, 0x0}};
88 const QUAD_T QUAD_REAL_VGV = {{0x4013, 0x8000}};
89 const QUAD_T QUAD_REAL_VSV = {{0x3ff2, 0x8000}};
90 const unt_2 QUAD_REAL_M_EXP = 0x7fff;
91 const unt_2 QUAD_REAL_M_SIGN = 0x8000;
92
93 const QUAD_T QUAD_REAL_ZERO =
94 {{0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000}};
95 const QUAD_T QUAD_REAL_TENTH =
96 {{0x3ffb, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0xcccc}};
97 const QUAD_T QUAD_REAL_HALF =
98 {{0x3ffe, 0x8000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000}};
99 const QUAD_T QUAD_REAL_ONE =
100 {{0x3fff, 0x8000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000}};
101 const QUAD_T QUAD_REAL_TWO =
102 {{0x4000, 0x8000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000}};
103 const QUAD_T QUAD_REAL_TEN =
104 {{0x4002, 0xa000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000}};
105 const QUAD_T QUAD_REAL_HUNDRED =
106 {{0x4005, 0xc800, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000}};
107 const QUAD_T QUAD_REAL_THOUSAND =
108 {{0x4008, 0xfa00, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000}};
109 const QUAD_T QUAD_REAL_PI4 =
110 {{0x3ffe, 0xc90f, 0xdaa2, 0x2168, 0xc234, 0xc4c6, 0x628b, 0x80dc, 0x1cd1, 0x2902, 0x4e08, 0x8a67, 0xcc74, 0x020b, 0xbea6, 0x3b14}};
111 const QUAD_T QUAD_REAL_PI2 =
112 {{0x3fff, 0xc90f, 0xdaa2, 0x2168, 0xc234, 0xc4c6, 0x628b, 0x80dc, 0x1cd1, 0x2902, 0x4e08, 0x8a67, 0xcc74, 0x020b, 0xbea6, 0x3b14}};
113 const QUAD_T QUAD_REAL_PI =
114 {{0x4000, 0xc90f, 0xdaa2, 0x2168, 0xc234, 0xc4c6, 0x628b, 0x80dc, 0x1cd1, 0x2902, 0x4e08, 0x8a67, 0xcc74, 0x020b, 0xbea6, 0x3b14}};
115 const QUAD_T QUAD_REAL_LN2 =
116 {{0x3ffe, 0xb172, 0x17f7, 0xd1cf, 0x79ab, 0xc9e3, 0xb398, 0x03f2, 0xf6af, 0x40f3, 0x4326, 0x7298, 0xb62d, 0x8a0d, 0x175b, 0x8bab}};
117 const QUAD_T QUAD_REAL_LN10 =
118 {{0x4000, 0x935d, 0x8ddd, 0xaaa8, 0xac16, 0xea56, 0xd62b, 0x82d3, 0x0a28, 0xe28f, 0xecf9, 0xda5d, 0xf90e, 0x83c6, 0x1e82, 0x01f0}};
119 const QUAD_T QUAD_REAL_SQRT2 =
120 {{0x3fff, 0xb504, 0xf333, 0xf9de, 0x6484, 0x597d, 0x89b3, 0x754a, 0xbe9f, 0x1d6f, 0x60ba, 0x893b, 0xa84c, 0xed17, 0xac85, 0x8334}};
121 const QUAD_T QUAD_REAL_LOG2_E =
122 {{0x3fff, 0xb8aa, 0x3b29, 0x5c17, 0xf0bb, 0xbe87, 0xfed0, 0x691d, 0x3e88, 0xeb57, 0x7aa8, 0xdd69, 0x5a58, 0x8b25, 0x166c, 0xd1a1}};
123 const QUAD_T QUAD_REAL_LOG2_10 =
124 {{0x4000, 0xd49a, 0x784b, 0xcd1b, 0x8afe, 0x492b, 0xf6ff, 0x4daf, 0xdb4c, 0xd96c, 0x55fe, 0x37b3, 0xad4e, 0x91b6, 0xac80, 0x82e8}};
125 const QUAD_T QUAD_REAL_LOG10_E =
126 {{0x3ffd, 0xde5b, 0xd8a9, 0x3728, 0x7195, 0x355b, 0xaaaf, 0xad33, 0xdc32, 0x3ee3, 0x4602, 0x45c9, 0xa202, 0x3a3f, 0x2d44, 0xf78f}};
127 const QUAD_T QUAD_REAL_RNDCORR =
128 {{0x3ffe, 0x8000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x01ae}};
129 const QUAD_T QUAD_REAL_FIXCORR =
130 {{0x3f17, 0xc000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000}};
131 const QUAD_T QUAD_REAL_NAN =
132 {{0x0000, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff}};
133
134 static const char *errmsg[] = {
135 "No error",
136 "Division by zero",
137 "Out of domain",
138 "Bad exponent",
139 "Floating point overflow",
140 "Invalid error code"
141 };
142
143 static QUAD_T c_tan (QUAD_T u);
144 static QUAD_T rred (QUAD_T u, int i, int *p);
145
146 //! @brief _pi32
147
148 void _pi32 (QUAD_T *x)
149 {
150 *x = QUAD_REAL_PI;
151 }
152
153 //! @brief sigerr_quad_real
154
155 int sigerr_quad_real (int errcond, int errcode, const char *where)
156 {
157 // errcode must come from the evaluation of an error condition.
158 // errcode, which should describe the type of the error,
159 // should always be one between QUAD_REAL_EDIV, QUAD_REAL_EDOM, QUAD_REAL_EBADEXP and QUAD_REAL_FPOFLOW.
160 if (errcond == 0) {
161 errcode = 0;
162 }
163 if (errcode < 0 || errcode > QUAD_REAL_NERR) {
164 errcode = QUAD_REAL_EINV;
165 }
166 if (errcode != 0) {
167 QUAD_RTE (where, errmsg[errcode]);
168 }
169 return errcode;
170 }
171
172 // Elementary stuff.
173
174 //! @brief neg_quad_real
175
176 inline QUAD_T neg_quad_real (QUAD_T s)
177 {
178 unt_2 *p = (unt_2 *) &QV (s);
179 *p ^= QUAD_REAL_M_SIGN;
180 return s;
181 }
182
183 //! @brief abs_quad_real
184
185 inline QUAD_T abs_quad_real (QUAD_T s)
186 {
187 unt_2 *p = (unt_2 *) &QV (s);
188 *p &= QUAD_REAL_M_EXP;
189 return s;
190 }
191
192 //! @brief getexp_quad_real
193
194 inline int getexp_quad_real (const QUAD_T * ps)
195 {
196 unt_2 *q = (unt_2 *) &(ps->value);
197 return (*q & QUAD_REAL_M_EXP) - QUAD_REAL_BIAS;
198 }
199
200 //! @brief getsgn_quad_real
201
202 inline int getsgn_quad_real (const QUAD_T * ps)
203 {
204 unt_2 *q = (unt_2 *) &(ps->value);
205 return (*q & QUAD_REAL_M_SIGN);
206 }
207
208 //! @brief real_cmp_quad_real
209
210 int real_cmp_quad_real (const QUAD_T * pa, const QUAD_T * pb)
211 {
212 unt_2 *p = (unt_2 *) &(pa->value), *q = (unt_2 *) &(pb->value);
213 unt_2 p0, q0;
214 int m;
215 for (m = 1; m <= FLT256_LEN && p[m] == 0; m++);
216 if (m > FLT256_LEN && (*p & QUAD_REAL_M_EXP) < QUAD_REAL_M_EXP) {
217 // *pa is actually zero
218 p0 = 0;
219 } else {
220 p0 = *p;
221 }
222 for (m = 1; m <= FLT256_LEN && q[m] == 0; m++);
223 if (m > FLT256_LEN && (*q & QUAD_REAL_M_EXP) < QUAD_REAL_M_EXP) {
224 // *pb is actually zero
225 q0 = 0;
226 } else {
227 q0 = *q;
228 }
229 unt_2 e = p0 & QUAD_REAL_M_SIGN, k = q0 & QUAD_REAL_M_SIGN;
230 if (e && !k) {
231 return -1;
232 } else if (!e && k) {
233 return 1;
234 } else { // *pa and *pb have the same sign
235 m = (e) ? -1 : 1;
236 e = p0 & QUAD_REAL_M_EXP;
237 k = q0 & QUAD_REAL_M_EXP;
238 if (e > k) {
239 return m;
240 } else if (e < k) {
241 return -m;
242 } else {
243 for (e = 0; *(++p) == *(++q) && e < FLT256_LEN; ++e);
244 if (e < FLT256_LEN) {
245 return (*p > *q ? m : -m);
246 } else {
247 return 0;
248 }
249 }
250 }
251 }
252
253 //! @brief real_2_quad_real
254
255 QUAD_T real_2_quad_real (QUAD_T s, int m)
256 {
257 unt_2 *p = (unt_2 *) &QV (s);
258 int e;
259 for (e = 1; e <= FLT256_LEN && p[e] == 0; e++);
260 if (e <= FLT256_LEN) {
261 e = *p & QUAD_REAL_M_EXP; // biased exponent
262 if (e + m < 0) {
263 return QUAD_REAL_ZERO;
264 } else if ((sigerr_quad_real (e + m >= QUAD_REAL_M_EXP, QUAD_REAL_FPOFLOW, NULL))) {
265 return ((s.value[0] & QUAD_REAL_M_SIGN) ? QUAD_REAL_MINF : QUAD_REAL_PINF);
266 } else {
267 *p += m;
268 return s;
269 }
270 } else { // s is zero or +-Inf
271 return s;
272 }
273 }
274
275 //! @brief sfmod_quad_real
276
277 QUAD_T sfmod_quad_real (QUAD_T s, int *p)
278 {
279 unt_2 *pa = (unt_2 *) &QV (s);
280 unt_2 *pb = pa + 1;
281 int_2 e, k;
282 e = (*pa & QUAD_REAL_M_EXP) - QUAD_REAL_BIAS;
283 if ((sigerr_quad_real (e >= 15, QUAD_REAL_FPOFLOW, NULL))) {
284 *p = -1;
285 return s;
286 } else if (e < 0) {
287 *p = 0;
288 return s;
289 }
290 *p = *pb >> (15 - e);
291 lshift_quad_real (++e, pb, FLT256_LEN);
292 *pa -= e;
293 for (e = 0; *pb == 0 && e < QUAD_REAL_MAX_P; ++pb, e += 16);
294 if (e == QUAD_REAL_MAX_P) {
295 return QUAD_REAL_ZERO;
296 }
297 for (k = 0; !((*pb << k) & QUAD_REAL_M_SIGN); ++k);
298 if ((k += e)) {
299 lshift_quad_real (k, pa + 1, FLT256_LEN);
300 *pa -= k;
301 }
302 return s;
303 }
304
305 //! @brief lshift_quad_real
306
307 void lshift_quad_real (int n, unt_2 *pm, int m)
308 {
309 unt_2 *pc = pm + m - 1;
310 if (n < 16 * m) {
311 unt_2 *pa = pm + n / 16;
312 m = n % 16;
313 n = 16 - m;
314 while (pa < pc) {
315 *pm = (*pa++) << m;
316 *pm++ |= *pa >> n;
317 }
318 *pm++ = *pa << m;
319 }
320 while (pm <= pc) {
321 *pm++ = 0;
322 }
323 }
324
325 //! @brief rshift_quad_real
326
327 void rshift_quad_real (int n, unt_2 *pm, int m)
328 {
329 unt_2 *pc = pm + m - 1;
330 if (n < 16 * m) {
331 unt_2 *pa = pc - n / 16;
332 m = n % 16;
333 n = 16 - m;
334 while (pa > pm) {
335 *pc = (*pa--) >> m;
336 *pc-- |= *pa << n;
337 }
338 *pc-- = *pa >> m;
339 }
340 while (pc >= pm) {
341 *pc-- = 0;
342 }
343 }
344
345 //! @brief nint_quad_real
346
347 QUAD_T nint_quad_real (QUAD_T x)
348 {
349 if (getsgn_quad_real (&x) == 0) {
350 return floor_quad_real (add_quad_real (x, QUAD_REAL_HALF, A68_FALSE));
351 } else {
352 return neg_quad_real (floor_quad_real (add_quad_real (QUAD_REAL_HALF, x, A68_TRUE)));
353 }
354 }
355
356 //! @brief aint_quad_real
357
358 QUAD_T aint_quad_real (QUAD_T x)
359 {
360 if (getsgn_quad_real (&x) == 0) {
361 return trunc_quad_real (x);
362 } else {
363 return neg_quad_real (trunc_quad_real (x));
364 }
365 }
366
367 //! @brief max_quad_real_2
368
369 QUAD_T max_quad_real_2 (QUAD_T a, QUAD_T b)
370 {
371 if (ge_quad_real (a, b)) {
372 return a;
373 } else {
374 return b;
375 }
376 }
377
378 //! @brief min_quad_real_2
379
380 QUAD_T min_quad_real_2 (QUAD_T a, QUAD_T b)
381 {
382 if (le_quad_real (a, b)) {
383 return a;
384 } else {
385 return b;
386 }
387 }
388
389 //! @brief mod_quad_real
390
391 QUAD_T mod_quad_real (QUAD_T a, QUAD_T b)
392 {
393 QUAD_T q = div_quad_real (a, b);
394 if (sgn_quad_real (&q) >= 0) {
395 q = floor_quad_real (q);
396 } else {
397 q = neg_quad_real (floor_quad_real (neg_quad_real (q)));
398 }
399 return add_quad_real (a, mul_quad_real (b, q), 1);
400 }
401
402 //! @brief add_quad_real
403
404 QUAD_T add_quad_real (QUAD_T s, QUAD_T t, int f)
405 {
406 REAL32 pe;
407 unt_2 *pc, *pf = pe;
408 unt_2 *pa = (unt_2 *) &QV (s);
409 unt_2 *pb = (unt_2 *) &QV (t);
410 int_2 e = *pa & QUAD_REAL_M_EXP;
411 int_2 k = *pb & QUAD_REAL_M_EXP;
412 if (f != 0) {
413 *pb ^= QUAD_REAL_M_SIGN;
414 }
415 unt_2 u = (*pb ^ *pa) & QUAD_REAL_M_SIGN;
416 f = 0;
417 if (e > k) {
418 if ((k = e - k) >= QUAD_REAL_MAX_P) {
419 return s;
420 }
421 rshift_quad_real (k, pb + 1, FLT256_LEN);
422 } else if (e < k) {
423 if ((e = k - e) >= QUAD_REAL_MAX_P) {
424 return t;
425 }
426 rshift_quad_real (e, pa + 1, FLT256_LEN);
427 e = k;
428 pc = pa;
429 pa = pb;
430 pb = pc;
431 } else if (u != 0) {
432 for (pc = pa, pf = pb; *(++pc) == *(++pf) && f < FLT256_LEN; ++f);
433 if (f >= FLT256_LEN) {
434 return QUAD_REAL_ZERO;
435 }
436 if (*pc < *pf) {
437 pc = pa;
438 pa = pb;
439 pb = pc;
440 }
441 pf = pe + f;
442 }
443 unt_2 h = *pa & QUAD_REAL_M_SIGN;
444 unt n = 0;
445 if (u != 0) {
446 for (pc = pb + FLT256_LEN; pc > pb; --pc) {
447 *pc = ~(*pc);
448 }
449 n = 1L;
450 }
451 for (pc = pe + FLT256_LEN, pa += FLT256_LEN, pb += FLT256_LEN; pc > pf;) {
452 n += *pa;
453 pa--;
454 n += *pb;
455 pb--;
456 *pc = n;
457 pc--;
458 n >>= 16;
459 }
460 if (u != 0) {
461 for (; *(++pc) == 0; ++f);
462 for (k = 0; !((*pc << k) & QUAD_REAL_M_SIGN); ++k);
463 if ((k += 16 * f)) {
464 if ((e -= k) <= 0) {
465 return QUAD_REAL_ZERO;
466 }
467 lshift_quad_real (k, pe + 1, FLT256_LEN);
468 }
469 } else {
470 if (n != 0) {
471 ++e;
472 if ((sigerr_quad_real (e == (short) QUAD_REAL_M_EXP, QUAD_REAL_FPOFLOW, NULL))) {
473 return (!h ? QUAD_REAL_PINF : QUAD_REAL_MINF);
474 }
475 ++pf;
476 rshift_quad_real (1, pf, FLT256_LEN);
477 *pf |= QUAD_REAL_M_SIGN;
478 }
479 }
480 *pe = e;
481 *pe |= h;
482 return *(QUAD_T *) pe;
483 }
484
485 //! @brief mul_quad_real
486
487 QUAD_T mul_quad_real (QUAD_T s, QUAD_T t)
488 {
489 unt_2 pe[FLT256_LEN + 2], *q0, *q1, h;
490 unt_2 *pa, *pb, *pc;
491 unt m, n, p;
492 int_2 e;
493 int_2 k;
494 q0 = (unt_2 *) &QV (s);
495 q1 = (unt_2 *) &QV (t);
496 e = (*q0 & QUAD_REAL_M_EXP) - QUAD_REAL_BIAS;
497 k = (*q1 & QUAD_REAL_M_EXP) + 1;
498 if ((sigerr_quad_real (e > (short) QUAD_REAL_M_EXP - k, QUAD_REAL_FPOFLOW, NULL))) {
499 return (((s.value[0] & QUAD_REAL_M_SIGN) ^ (t.value[0] & QUAD_REAL_M_SIGN)) ? QUAD_REAL_MINF : QUAD_REAL_PINF);
500 }
501 if ((e += k) <= 0) {
502 return QUAD_REAL_ZERO;
503 }
504 h = (*q0 ^ *q1) & QUAD_REAL_M_SIGN;
505 for (++q1, k = FLT256_LEN, p = n = 0L, pc = pe + FLT256_LEN + 1; k > 0; --k) {
506 for (pa = q0 + k, pb = q1; pa > q0;) {
507 m = *pa--;
508 m *= *pb++;
509 n += (m & 0xffffL);
510 p += (m >> 16);
511 }
512 *pc-- = n;
513 n = p + (n >> 16);
514 p = 0L;
515 }
516 *pc = n;
517 if (!(*pc & QUAD_REAL_M_SIGN)) {
518 --e;
519 if (e <= 0) {
520 return QUAD_REAL_ZERO;
521 }
522 lshift_quad_real (1, pc, FLT256_LEN + 1);
523 }
524 if ((sigerr_quad_real (e == (short) QUAD_REAL_M_EXP, QUAD_REAL_FPOFLOW, NULL))) {
525 return (!h ? QUAD_REAL_PINF : QUAD_REAL_MINF);
526 }
527 *pe = e;
528 *pe |= h;
529 return *(QUAD_T *) pe;
530 }
531
532 //! @brief div_quad_real
533
534 QUAD_T div_quad_real (QUAD_T s, QUAD_T t)
535 {
536 QUAD_T a;
537 unt_2 *pc, e, i;
538 pc = (unt_2 *) &QV (t);
539 e = *pc;
540 *pc = QUAD_REAL_BIAS;
541 if ((sigerr_quad_real (real_cmp_quad_real (&t, &QUAD_REAL_ZERO) == 0, QUAD_REAL_EDIV, "div_quad_real"))) {
542 return QUAD_REAL_ZERO;
543 } else {
544 a = real_to_quad_real (1 / quad_real_to_real (t));
545 *pc = e;
546 pc = (unt_2 *) &QV (a);
547 *pc += QUAD_REAL_BIAS - (e & QUAD_REAL_M_EXP);
548 *pc |= e & QUAD_REAL_M_SIGN;
549 for (i = 0; i < QUAD_REAL_ITT_DIV; ++i) {
550 a = mul_quad_real (a, add_quad_real (QUAD_REAL_TWO, mul_quad_real (a, t), 1));
551 }
552 return mul_quad_real (s, a);
553 }
554 }
555
556 //! @brief exp2_quad_real
557
558 QUAD_T exp2_quad_real (QUAD_T x)
559 {
560 QUAD_T s, d, f;
561 unt_2 *pf = (unt_2 *) &QV (x);
562 int m, k;
563 if (real_cmp_quad_real (&x, &QUAD_REAL_E2MIN) < 0) {
564 return QUAD_REAL_ZERO;
565 } else if ((sigerr_quad_real (real_cmp_quad_real (&x, &QUAD_REAL_E2MAX) > 0, QUAD_REAL_FPOFLOW, NULL))) {
566 return QUAD_REAL_PINF;
567 } else {
568 m = (*pf & QUAD_REAL_M_SIGN) ? 1 : 0;
569 x = sfmod_quad_real (x, &k);
570 if (m != 0) {
571 k *= -1;
572 }
573 // -QUAD_REAL_BIAS <= k <= +QUAD_REAL_BIAS
574 x = mul_quad_real (x, QUAD_REAL_LN2);
575 if (getexp_quad_real (&x) > -QUAD_REAL_BIAS) {
576 x = real_2_quad_real (x, -1);
577 s = mul_quad_real (x, x);
578 f = QUAD_REAL_ZERO;
579 for (d = int_to_quad_real (m = QUAD_REAL_MS_EXP); m > 1; m -= 2, d = int_to_quad_real (m)) {
580 f = div_quad_real (s, _add_quad_real_ (d, f));
581 }
582 f = div_quad_real (x, _add_quad_real_ (d, f));
583 f = div_quad_real (_add_quad_real_ (d, f), _sub_quad_real_ (d, f));
584 } else {
585 f = QUAD_REAL_ONE;
586 }
587 pf = (unt_2 *) &QV (f);
588 if (-k > *pf) {
589 return QUAD_REAL_ZERO;
590 } else {
591 *pf += k;
592 if ((sigerr_quad_real (*pf >= QUAD_REAL_M_EXP, QUAD_REAL_FPOFLOW, NULL))) {
593 return QUAD_REAL_PINF;
594 } else {
595 return f;
596 }
597 }
598 }
599 }
600
601 //! @brief exp_quad_real
602
603 QUAD_T exp_quad_real (QUAD_T z)
604 {
605 return exp2_quad_real (mul_quad_real (z, QUAD_REAL_LOG2_E));
606 }
607
608 //! @brief exp10_quad_real
609
610 QUAD_T exp10_quad_real (QUAD_T z)
611 {
612 return exp2_quad_real (mul_quad_real (z, QUAD_REAL_LOG2_10));
613 }
614
615 //! @brief fmod_quad_real
616
617 QUAD_T fmod_quad_real (QUAD_T s, QUAD_T t, QUAD_T * q)
618 {
619 if ((sigerr_quad_real (real_cmp_quad_real (&t, &QUAD_REAL_ZERO) == 0, QUAD_REAL_EDIV, "fmod_quad_real"))) {
620 return QUAD_REAL_ZERO;
621 } else {
622 unt_2 *p, mask = 0xffff;
623 int_2 e, i;
624 int u;
625 *q = div_quad_real (s, t);
626 p = (unt_2 *) &(q->value);
627 u = (*p & QUAD_REAL_M_SIGN) ? 0 : 1;
628 e = (*p &= QUAD_REAL_M_EXP); // biased exponent of *q
629 e = e < QUAD_REAL_BIAS ? 0 : e - QUAD_REAL_BIAS + 1;
630 for (i = 1; e / 16 > 0; i++, e -= 16);
631 if (i <= FLT256_LEN) {
632 // e = 0, ..., 15
633 mask <<= 16 - e;
634 p[i] &= mask;
635 for (i++; i <= FLT256_LEN; p[i] = 0, i++);
636 }
637 // Now *q == abs(quotient of (s/t))
638 return add_quad_real (s, mul_quad_real (t, *q), u);
639 }
640 }
641
642 //! @brief frexp_quad_real
643
644 QUAD_T frexp_quad_real (QUAD_T s, int *p)
645 {
646 unt_2 *ps = (unt_2 *) &QV (s), u;
647 *p = (*ps & QUAD_REAL_M_EXP) - QUAD_REAL_BIAS + 1;
648 u = *ps & QUAD_REAL_M_SIGN;
649 *ps = QUAD_REAL_BIAS - 1;
650 *ps |= u;
651 return s;
652 }
653
654 //! @brief nullify
655
656 static void nullify (int skip, unt_2 *p, int k)
657 {
658 // After skipping the first 'skip' bytes of the vector 'p',
659 // this function nullifies all the remaining ones. 'k' is
660 // the number of words forming the vector p.
661 // Warning: 'skip' must be positive !
662 int i;
663 unt_2 mask = 0xffff;
664 for (i = 0; skip / 16 > 0; i++, skip -= 16);
665 if (i < k) {
666 // skip = 0, ..., 15
667 mask <<= 16 - skip;
668 p[i] &= mask;
669 for (i++; i < k; p[i] = 0, i++);
670 }
671 }
672
673 //! @brief canonic_form
674
675 static void canonic_form (QUAD_T * px)
676 {
677 unt_2 *p, u;
678 int_2 e, i, j, skip;
679 p = (unt_2 *) &(px->value);
680 e = (*p & QUAD_REAL_M_EXP); // biased exponent of x
681 u = (*p & QUAD_REAL_M_SIGN); // sign of x
682 if (e < QUAD_REAL_BIAS - 1) {
683 return;
684 } else {
685 unt_2 mask = 0xffff;
686 // e >= QUAD_REAL_BIAS - 1
687 for (i = 1, skip = e + 1 - QUAD_REAL_BIAS; skip / 16 > 0; i++, skip -= 16);
688 if (i <= FLT256_LEN) {
689 // skip = 0, ..., 15
690 mask >>= skip;
691 if ((p[i] & mask) != mask) {
692 return;
693 } else {
694 for (j = i + 1; j <= FLT256_LEN && p[j] == 0xffff; j++);
695 if (j > FLT256_LEN) {
696 p[i] -= mask;
697 for (j = i + 1; j <= FLT256_LEN; p[j] = 0, j++);
698 if (!(p[1] & 0x8000)) {
699 p[1] = 0x8000;
700 *p = ++e;
701 *p |= u;
702 } else if (u != 0) {
703 *px = add_quad_real (*px, QUAD_REAL_ONE, 1);
704 } else {
705 *px = add_quad_real (*px, QUAD_REAL_ONE, 0);
706 }
707 }
708 }
709 }
710 }
711 }
712
713 //! @brief frac_quad_real
714
715 QUAD_T frac_quad_real (QUAD_T x)
716 {
717 // frac_quad_real(x) returns the fractional part of the number x.
718 // frac_quad_real(x) has the same sign as x.
719 unt_2 u, *p;
720 int_2 e;
721 int n;
722 canonic_form (&x);
723 p = (unt_2 *) &QV (x);
724 e = (*p & QUAD_REAL_M_EXP); // biased exponent of x
725 if (e < QUAD_REAL_BIAS) {
726 return x; // The integer part of x is zero
727 } else {
728 u = *p & QUAD_REAL_M_SIGN; // sign of x
729 n = e - QUAD_REAL_BIAS + 1;
730 lshift_quad_real (n, p + 1, FLT256_LEN);
731 e = QUAD_REAL_BIAS - 1;
732 // Now I have to take in account the rule
733 // of the leading one.
734 while (e > 0 && !(p[1] & QUAD_REAL_M_SIGN)) {
735 lshift_quad_real (1, p + 1, FLT256_LEN);
736 e -= 1;
737 }
738 // Now p+1 points to the fractionary part of x,
739 // u is its sign, e is its biased exponent.
740 p[0] = e;
741 p[0] |= u;
742 return *(QUAD_T *) p;
743 }
744 }
745
746 //! @brief trunc_quad_real
747
748 QUAD_T trunc_quad_real (QUAD_T x)
749 {
750 // trunc_quad_real(x) returns the integer part of the number x.
751 // trunc_quad_real(x) has the same sign as x.
752 unt_2 *p;
753 int_2 e;
754 canonic_form (&x);
755 p = (unt_2 *) &QV (x);
756 e = (*p & QUAD_REAL_M_EXP); // biased exponent of x
757 if (e < QUAD_REAL_BIAS) {
758 return QUAD_REAL_ZERO; // The integer part of x is zero
759 } else {
760 nullify (e - QUAD_REAL_BIAS + 1, p + 1, FLT256_LEN);
761 return *(QUAD_T *) p;
762 }
763 }
764
765 //! @brief round_quad_real
766
767 QUAD_T round_quad_real (QUAD_T x)
768 {
769 return trunc_quad_real (add_quad_real (x, QUAD_REAL_RNDCORR, x.value[0] & QUAD_REAL_M_SIGN));
770 }
771
772 //! @brief ceil_quad_real
773
774 QUAD_T ceil_quad_real (QUAD_T x)
775 {
776 unt_2 *ps = (unt_2 *) &QV (x);
777 if ((*ps & QUAD_REAL_M_SIGN)) {
778 return trunc_quad_real (x);
779 } else {
780 QUAD_T y = frac_quad_real (x);
781 // y has the same sign as x (see above).
782 return (real_cmp_quad_real (&y, &QUAD_REAL_ZERO) > 0 ? add_quad_real (trunc_quad_real (x), QUAD_REAL_ONE, 0) : x);
783 }
784 }
785
786 //! @brief floor_quad_real
787
788 QUAD_T floor_quad_real (QUAD_T x)
789 {
790 unt_2 *ps = (unt_2 *) &QV (x);
791 if ((*ps & QUAD_REAL_M_SIGN)) {
792 QUAD_T y = frac_quad_real (x);
793 // y has the same sign as x (see above).
794 return (real_cmp_quad_real (&y, &QUAD_REAL_ZERO) < 0 ? add_quad_real (trunc_quad_real (x), QUAD_REAL_ONE, 1) : x);
795 } else {
796 return trunc_quad_real (x);
797 }
798 }
799
800 //! @brief add_correction_quad_real
801
802 static void add_correction_quad_real (QUAD_T * px, int k)
803 {
804 int_2 e = (px->value[0] & QUAD_REAL_M_EXP) - QUAD_REAL_BIAS;
805 // e = (e > 0 ? e : 0);
806 *px = add_quad_real (*px, real_2_quad_real (QUAD_REAL_FIXCORR, e), k);
807 }
808
809 //! @brief fix_quad_real
810
811 QUAD_T fix_quad_real (QUAD_T x)
812 {
813 unt_2 *p;
814 int_2 e;
815 add_correction_quad_real (&x, x.value[0] & QUAD_REAL_M_SIGN);
816 p = (unt_2 *) &QV (x);
817 e = (*p & QUAD_REAL_M_EXP); // biased exponent of x
818 if (e < QUAD_REAL_BIAS) {
819 return QUAD_REAL_ZERO; // The integer part of x is zero
820 } else {
821 nullify (e - QUAD_REAL_BIAS + 1, p + 1, FLT256_LEN);
822 return *(QUAD_T *) p;
823 }
824 }
825
826 //! @brief tanh_quad_real
827
828 QUAD_T tanh_quad_real (QUAD_T z)
829 {
830 QUAD_T s, d, f;
831 int m, k;
832 if ((k = getexp_quad_real (&z)) > QUAD_REAL_K_TANH) {
833 if (getsgn_quad_real (&z)) {
834 return neg_quad_real (QUAD_REAL_ONE);
835 } else {
836 return QUAD_REAL_ONE;
837 }
838 }
839 if (k < QUAD_REAL_K_LIN) {
840 return z;
841 }
842 ++k;
843 if (k > 0) {
844 z = real_2_quad_real (z, -k);
845 }
846 s = mul_quad_real (z, z);
847 f = QUAD_REAL_ZERO;
848 for (d = int_to_quad_real (m = QUAD_REAL_MS_HYP); m > 1;) {
849 f = div_quad_real (s, _add_quad_real_ (d, f));
850 d = int_to_quad_real (m -= 2);
851 }
852 f = div_quad_real (z, _add_quad_real_ (d, f));
853 for (; k > 0; --k) {
854 f = div_quad_real (real_2_quad_real (f, 1), add_quad_real (d, mul_quad_real (f, f), 0));
855 }
856 return f;
857 }
858
859 //! @brief sinh_quad_real
860
861 QUAD_T sinh_quad_real (QUAD_T z)
862 {
863 int k;
864 if ((k = getexp_quad_real (&z)) < QUAD_REAL_K_LIN) {
865 return z;
866 } else if (k < 0) {
867 z = tanh_quad_real (real_2_quad_real (z, -1));
868 return div_quad_real (real_2_quad_real (z, 1), add_quad_real (QUAD_REAL_ONE, mul_quad_real (z, z), 1));
869 } else {
870 z = exp_quad_real (z);
871 return real_2_quad_real (add_quad_real (z, div_quad_real (QUAD_REAL_ONE, z), 1), -1);
872 }
873 }
874
875 //! @brief cosh_quad_real
876
877 QUAD_T cosh_quad_real (QUAD_T z)
878 {
879 if (getexp_quad_real (&z) < QUAD_REAL_K_LIN) {
880 return QUAD_REAL_ONE;
881 }
882 z = exp_quad_real (z);
883 return real_2_quad_real (add_quad_real (z, div_quad_real (QUAD_REAL_ONE, z), 0), -1);
884 }
885
886 //! @brief atanh_quad_real
887
888 QUAD_T atanh_quad_real (QUAD_T x)
889 {
890 QUAD_T y = x;
891 y.value[0] &= QUAD_REAL_M_EXP; // Now y == abs(x)
892 if ((sigerr_quad_real (real_cmp_quad_real (&y, &QUAD_REAL_ONE) >= 0, QUAD_REAL_EDOM, "xatanh"))) {
893 return ((x.value[0] & QUAD_REAL_M_SIGN) ? QUAD_REAL_MINF : QUAD_REAL_PINF);
894 } else {
895 y = div_quad_real (add_quad_real (QUAD_REAL_ONE, x, 0), add_quad_real (QUAD_REAL_ONE, x, 1));
896 return real_2_quad_real (log_quad_real (y), -1);
897 }
898 }
899
900 //! @brief asinh_quad_real
901
902 QUAD_T asinh_quad_real (QUAD_T x)
903 {
904 QUAD_T y = mul_quad_real (x, x);
905 y = sqrt_quad_real (add_quad_real (QUAD_REAL_ONE, y, 0));
906 if ((x.value[0] & QUAD_REAL_M_SIGN)) {
907 return neg_quad_real (log_quad_real (_sub_quad_real_ (y, x)));
908 } else {
909 return log_quad_real (_add_quad_real_ (x, y));
910 }
911 }
912
913 //! @brief acosh_quad_real
914
915 QUAD_T acosh_quad_real (QUAD_T x)
916 {
917 if ((sigerr_quad_real (real_cmp_quad_real (&x, &QUAD_REAL_ONE) < 0, QUAD_REAL_EDOM, "acosh_quad_real"))) {
918 return QUAD_REAL_ZERO;
919 } else {
920 QUAD_T y = mul_quad_real (x, x);
921 y = sqrt_quad_real (add_quad_real (y, QUAD_REAL_ONE, 1));
922 return log_quad_real (_add_quad_real_ (x, y));
923 }
924 }
925
926 //! @brief atan_quad_real
927
928 QUAD_T atan_quad_real (QUAD_T z)
929 {
930 QUAD_T s, f;
931 int k, m;
932 if ((k = getexp_quad_real (&z)) < QUAD_REAL_K_LIN) {
933 return z;
934 }
935 if (k >= 0) {
936 // k>=0 is equivalent to abs(z) >= 1.0
937 z = div_quad_real (QUAD_REAL_ONE, z);
938 m = 1;
939 } else {
940 m = 0;
941 }
942 f = real_to_quad_real (atan (quad_real_to_real (z)));
943 s = add_quad_real (QUAD_REAL_ONE, mul_quad_real (z, z), 0);
944 for (k = 0; k < QUAD_REAL_ITT_DIV; ++k) {
945 f = add_quad_real (f, div_quad_real (add_quad_real (z, tan_quad_real (f), 1), s), 0);
946 }
947 if (m != 0) {
948 if (getsgn_quad_real (&f)) {
949 return add_quad_real (neg_quad_real (QUAD_REAL_PI2), f, 1);
950 } else {
951 return add_quad_real (QUAD_REAL_PI2, f, 1);
952 }
953 } else {
954 return f;
955 }
956 }
957
958 //! @brief asin_quad_real
959
960 QUAD_T asin_quad_real (QUAD_T z)
961 {
962 QUAD_T u = z;
963 u.value[0] &= QUAD_REAL_M_EXP;
964 if ((sigerr_quad_real (real_cmp_quad_real (&u, &QUAD_REAL_ONE) > 0, QUAD_REAL_EDOM, "asin_quad_real"))) {
965 return ((getsgn_quad_real (&z)) ? neg_quad_real (QUAD_REAL_PI2) : QUAD_REAL_PI2);
966 } else {
967 if (getexp_quad_real (&z) < QUAD_REAL_K_LIN) {
968 return z;
969 }
970 u = sqrt_quad_real (add_quad_real (QUAD_REAL_ONE, mul_quad_real (z, z), 1));
971 if (getexp_quad_real (&u) == -QUAD_REAL_BIAS) {
972 return ((getsgn_quad_real (&z)) ? neg_quad_real (QUAD_REAL_PI2) : QUAD_REAL_PI2);
973 }
974 return atan_quad_real (div_quad_real (z, u));
975 }
976 }
977
978 //! @brief acos_quad_real
979
980 QUAD_T acos_quad_real (QUAD_T z)
981 {
982 QUAD_T u = z;
983 u.value[0] &= QUAD_REAL_M_EXP;
984 if ((sigerr_quad_real (real_cmp_quad_real (&u, &QUAD_REAL_ONE) > 0, QUAD_REAL_EDOM, "acos_quad_real"))) {
985 return ((getsgn_quad_real (&z)) ? QUAD_REAL_PI : QUAD_REAL_ZERO);
986 } else {
987 if (getexp_quad_real (&z) == -QUAD_REAL_BIAS) {
988 return QUAD_REAL_PI2;
989 }
990 u = sqrt_quad_real (add_quad_real (QUAD_REAL_ONE, mul_quad_real (z, z), 1));
991 u = atan_quad_real (div_quad_real (u, z));
992 if (getsgn_quad_real (&z)) {
993 return add_quad_real (QUAD_REAL_PI, u, 0);
994 } else {
995 return u;
996 }
997 }
998 }
999
1000 //! @brief atan2_quad_real
1001
1002 QUAD_T atan2_quad_real (QUAD_T y, QUAD_T x)
1003 {
1004 int rs, is;
1005 rs = sgn_quad_real (&x);
1006 is = sgn_quad_real (&y);
1007 if (rs > 0) {
1008 return atan_quad_real (div_quad_real (y, x));
1009 } else if (rs < 0) {
1010 x.value[0] ^= QUAD_REAL_M_SIGN;
1011 y.value[0] ^= QUAD_REAL_M_SIGN;
1012 if (is >= 0) {
1013 return add_quad_real (QUAD_REAL_PI, atan_quad_real (div_quad_real (y, x)), 0);
1014 } else {
1015 return add_quad_real (atan_quad_real (div_quad_real (y, x)), QUAD_REAL_PI, 1);
1016 }
1017 } else { // x is zero !
1018 if (!sigerr_quad_real (is == 0, QUAD_REAL_EDOM, "atan2_quad_real")) {
1019 return (is > 0 ? QUAD_REAL_PI2 : neg_quad_real (QUAD_REAL_PI2));
1020 } else {
1021 return QUAD_REAL_ZERO; // Dummy value :)
1022 }
1023 }
1024 }
1025
1026 //! @brief log_quad_real
1027
1028 QUAD_T log_quad_real (QUAD_T z)
1029 {
1030 QUAD_T f, h;
1031 int k, m;
1032 if ((sigerr_quad_real ((getsgn_quad_real (&z)) || getexp_quad_real (&z) == -QUAD_REAL_BIAS, QUAD_REAL_EDOM, "log_quad_real"))) {
1033 return QUAD_REAL_MINF;
1034 } else if (real_cmp_quad_real (&z, &QUAD_REAL_ONE) == 0) {
1035 return QUAD_REAL_ZERO;
1036 } else {
1037 z = frexp_quad_real (z, &m);
1038 z = mul_quad_real (z, QUAD_REAL_SQRT2);
1039 z = div_quad_real (add_quad_real (z, QUAD_REAL_ONE, 1), add_quad_real (z, QUAD_REAL_ONE, 0));
1040 h = real_2_quad_real (z, 1);
1041 z = mul_quad_real (z, z);
1042 for (f = h, k = 1; getexp_quad_real (&h) > -QUAD_REAL_MAX_P;) {
1043 h = mul_quad_real (h, z);
1044 f = add_quad_real (f, div_quad_real (h, int_to_quad_real (k += 2)), 0);
1045 }
1046 return add_quad_real (f, mul_quad_real (QUAD_REAL_LN2, real_to_quad_real (m - .5)), 0);
1047 }
1048 }
1049
1050 //! @brief log2_quad_real
1051
1052 QUAD_T log2_quad_real (QUAD_T z)
1053 {
1054 QUAD_T f, h;
1055 int k, m;
1056 if ((sigerr_quad_real ((getsgn_quad_real (&z)) || getexp_quad_real (&z) == -QUAD_REAL_BIAS, QUAD_REAL_EDOM, "log2_quad_real"))) {
1057 return QUAD_REAL_MINF;
1058 } else if (real_cmp_quad_real (&z, &QUAD_REAL_ONE) == 0) {
1059 return QUAD_REAL_ZERO;
1060 } else {
1061 z = frexp_quad_real (z, &m);
1062 z = mul_quad_real (z, QUAD_REAL_SQRT2);
1063 z = div_quad_real (add_quad_real (z, QUAD_REAL_ONE, 1), add_quad_real (z, QUAD_REAL_ONE, 0));
1064 h = real_2_quad_real (z, 1);
1065 z = mul_quad_real (z, z);
1066 for (f = h, k = 1; getexp_quad_real (&h) > -QUAD_REAL_MAX_P;) {
1067 h = mul_quad_real (h, z);
1068 f = add_quad_real (f, div_quad_real (h, int_to_quad_real (k += 2)), 0);
1069 }
1070 return add_quad_real (mul_quad_real (f, QUAD_REAL_LOG2_E), real_to_quad_real (m - .5), 0);
1071 }
1072 }
1073
1074 //! @brief log10_quad_real
1075
1076 QUAD_T log10_quad_real (QUAD_T z)
1077 {
1078 QUAD_T w = log_quad_real (z);
1079 if (real_cmp_quad_real (&w, &QUAD_REAL_MINF) <= 0) {
1080 return QUAD_REAL_MINF;
1081 } else {
1082 return mul_quad_real (w, QUAD_REAL_LOG10_E);
1083 }
1084 }
1085 //! @brief eq_quad_real
1086
1087 int eq_quad_real (QUAD_T x1, QUAD_T x2)
1088 {
1089 return (real_cmp_quad_real (&x1, &x2) == 0);
1090 }
1091
1092 //! @brief neq_quad_real
1093
1094 int neq_quad_real (QUAD_T x1, QUAD_T x2)
1095 {
1096 return (real_cmp_quad_real (&x1, &x2) != 0);
1097 }
1098
1099 //! @brief gt_quad_real
1100
1101 int gt_quad_real (QUAD_T x1, QUAD_T x2)
1102 {
1103 return (real_cmp_quad_real (&x1, &x2) > 0);
1104 }
1105
1106 //! @brief ge_quad_real
1107
1108 int ge_quad_real (QUAD_T x1, QUAD_T x2)
1109 {
1110 return (real_cmp_quad_real (&x1, &x2) >= 0);
1111 }
1112
1113 //! @brief lt_quad_real
1114
1115 int lt_quad_real (QUAD_T x1, QUAD_T x2)
1116 {
1117 return (real_cmp_quad_real (&x1, &x2) < 0);
1118 }
1119
1120 //! @brief le_quad_real
1121
1122 int le_quad_real (QUAD_T x1, QUAD_T x2)
1123 {
1124 return (real_cmp_quad_real (&x1, &x2) <= 0);
1125 }
1126 // isNaN_quad_real (&x) returns 1 if and only if x is not a valid number
1127
1128 int isNaN_quad_real (const QUAD_T * u)
1129 {
1130 unt_2 *p = (unt_2 *) &(u->value);
1131 if (*p != 0) {
1132 return 0;
1133 } else {
1134 int i;
1135 for (i = 1; i <= FLT256_LEN && p[i] == 0x0; i++);
1136 return (i <= FLT256_LEN ? 1 : 0);
1137 }
1138 }
1139
1140 //! @brief is0_quad_real
1141
1142 int is0_quad_real (const QUAD_T * u)
1143 {
1144 unt_2 *p = (unt_2 *) &(u->value);
1145 int m;
1146 for (m = 1; m <= FLT256_LEN && p[m] == 0; m++);
1147 return (m > FLT256_LEN && (*p & QUAD_REAL_M_EXP) < QUAD_REAL_M_EXP ? 1 : 0);
1148 }
1149
1150 //! @brief not0_quad_real
1151
1152 int not0_quad_real (const QUAD_T * u)
1153 {
1154 unt_2 *p = (unt_2 *) &(u->value);
1155 int m;
1156 for (m = 1; m <= FLT256_LEN && p[m] == 0; m++);
1157 return (m > FLT256_LEN && (*p & QUAD_REAL_M_EXP) < QUAD_REAL_M_EXP ? 0 : 1);
1158 }
1159
1160 //! @brief sgn_quad_real
1161
1162 int sgn_quad_real (const QUAD_T * u)
1163 {
1164 unt_2 *p = (unt_2 *) &(u->value);
1165 int m;
1166 for (m = 1; m <= FLT256_LEN && p[m] == 0; m++);
1167 if ((m > FLT256_LEN && (*p & QUAD_REAL_M_EXP) < QUAD_REAL_M_EXP) || !*p) {
1168 return 0;
1169 } else {
1170 return ((*p & QUAD_REAL_M_SIGN) ? -1 : 1);
1171 }
1172 }
1173
1174 //! @brief isPinf_quad_real
1175
1176 int isPinf_quad_real (const QUAD_T * u)
1177 {
1178 return (*u->value == QUAD_REAL_M_EXP ? 1 : 0);
1179 }
1180
1181 //! @brief isMinf_quad_real
1182
1183 int isMinf_quad_real (const QUAD_T * u)
1184 {
1185 return (*u->value == (QUAD_REAL_M_EXP | QUAD_REAL_M_SIGN) ? 1 : 0);
1186 }
1187
1188 //! @brief isordnumb_quad_real
1189
1190 int isordnumb_quad_real (const QUAD_T * u)
1191 {
1192 int isNaN, isfinite;
1193 unt_2 *p = (unt_2 *) &(u->value);
1194 if (*p != 0) {
1195 isNaN = 0;
1196 } else {
1197 int i;
1198 for (i = 1; i <= FLT256_LEN && p[i] == 0x0; i++);
1199 isNaN = i <= FLT256_LEN;
1200 }
1201 isfinite = (*p & QUAD_REAL_M_EXP) < QUAD_REAL_M_EXP;
1202 return (!isNaN && (isfinite) ? 1 : 0);
1203 }
1204
1205 //! @brief pwr_quad_real
1206
1207 QUAD_T pwr_quad_real (QUAD_T s, int n)
1208 {
1209 QUAD_T t;
1210 unsigned k, m;
1211 t = QUAD_REAL_ONE;
1212 if (n < 0) {
1213 m = -n;
1214 if ((sigerr_quad_real (real_cmp_quad_real (&s, &QUAD_REAL_ZERO) == 0, QUAD_REAL_EBADEXP, "pwr_quad_real"))) {
1215 return QUAD_REAL_ZERO;
1216 }
1217 s = div_quad_real (QUAD_REAL_ONE, s);
1218 } else {
1219 m = n;
1220 }
1221 if (m != 0) {
1222 k = 1;
1223 while (1) {
1224 if (k & m) {
1225 t = mul_quad_real (s, t);
1226 }
1227 if ((k <<= 1) <= m) {
1228 s = mul_quad_real (s, s);
1229 } else {
1230 break;
1231 }
1232 }
1233 } else {
1234 sigerr_quad_real (real_cmp_quad_real (&s, &QUAD_REAL_ZERO) == 0, QUAD_REAL_EBADEXP, "pwr_quad_real");
1235 }
1236 return t;
1237 }
1238
1239 //! @brief pow_quad_real
1240
1241 QUAD_T pow_quad_real (QUAD_T x, QUAD_T y)
1242 {
1243 if (sigerr_quad_real ((getsgn_quad_real (&x)) || getexp_quad_real (&x) == -QUAD_REAL_BIAS, QUAD_REAL_EDOM, "pow_quad_real")) {
1244 return QUAD_REAL_ZERO;
1245 } else {
1246 return exp2_quad_real (mul_quad_real (log2_quad_real (x), y));
1247 }
1248 }
1249
1250 //! @brief sqrt_quad_real
1251
1252 QUAD_T sqrt_quad_real (QUAD_T z)
1253 {
1254 if ((sigerr_quad_real ((getsgn_quad_real (&z)), QUAD_REAL_EDOM, "sqrt_quad_real"))) {
1255 return QUAD_REAL_ZERO;
1256 } else {
1257 unt_2 *pc = (unt_2 *) &QV (z);
1258 if (*pc == 0) {
1259 return QUAD_REAL_ZERO;
1260 }
1261 int_2 e = *pc - QUAD_REAL_BIAS;
1262 *pc = QUAD_REAL_BIAS + (e % 2);
1263 e /= 2;
1264 QUAD_T h, s = real_to_quad_real (sqrt (quad_real_to_real (z)));
1265 for (int_2 m = 0; m < QUAD_REAL_ITT_DIV; ++m) {
1266 h = div_quad_real (add_quad_real (z, mul_quad_real (s, s), 1), real_2_quad_real (s, 1));
1267 s = _add_quad_real_ (s, h);
1268 }
1269 pc = (unt_2 *) &QV (s);
1270 *pc += e;
1271 return s;
1272 }
1273 }
1274
1275 static int odd_quad_real (QUAD_T x)
1276 {
1277 unt_2 *p = (unt_2 *) &QV (x);
1278 int_2 e, i;
1279 e = (*p & QUAD_REAL_M_EXP) - QUAD_REAL_BIAS; // exponent of x
1280 if (e < 0) {
1281 return 0;
1282 } else {
1283 for (i = 1; e / 16 > 0; i++, e -= 16);
1284 // Now e = 0, ..., 15
1285 return (i <= FLT256_LEN ? p[i] & 0x8000 >> e : 0);
1286 }
1287 }
1288
1289 //! @brief tan_quad_real
1290
1291 QUAD_T tan_quad_real (QUAD_T z)
1292 {
1293 int k, m;
1294 z = rred (z, 't', &k);
1295 if ((sigerr_quad_real (real_cmp_quad_real (&z, &QUAD_REAL_PI2) >= 0, QUAD_REAL_EDOM, "tan_quad_real"))) {
1296 return (!k ? QUAD_REAL_PINF : QUAD_REAL_MINF);
1297 } else {
1298 if (real_cmp_quad_real (&z, &QUAD_REAL_PI4) == 1) {
1299 m = 1;
1300 z = add_quad_real (QUAD_REAL_PI2, z, 1);
1301 } else {
1302 m = 0;
1303 }
1304 if (k != 0) {
1305 z = neg_quad_real (c_tan (z));
1306 } else {
1307 z = c_tan (z);
1308 }
1309 if (m != 0) {
1310 return div_quad_real (QUAD_REAL_ONE, z);
1311 } else {
1312 return z;
1313 }
1314 }
1315 }
1316
1317 //! @brief cos_quad_real
1318
1319 QUAD_T cos_quad_real (QUAD_T z)
1320 {
1321 int k;
1322 z = rred (z, 'c', &k);
1323 if (getexp_quad_real (&z) < QUAD_REAL_K_LIN) {
1324 if (k != 0) {
1325 return neg_quad_real (QUAD_REAL_ONE);
1326 } else {
1327 return QUAD_REAL_ONE;
1328 }
1329 }
1330 z = c_tan (real_2_quad_real (z, -1));
1331 z = mul_quad_real (z, z);
1332 z = div_quad_real (add_quad_real (QUAD_REAL_ONE, z, 1), add_quad_real (QUAD_REAL_ONE, z, 0));
1333 if (k != 0) {
1334 return neg_quad_real (z);
1335 } else {
1336 return z;
1337 }
1338 }
1339
1340 //! @brief sin_quad_real
1341
1342 QUAD_T sin_quad_real (QUAD_T z)
1343 {
1344 int k;
1345 z = rred (z, 's', &k);
1346 if (getexp_quad_real (&z) >= QUAD_REAL_K_LIN) {
1347 z = c_tan (real_2_quad_real (z, -1));
1348 z = div_quad_real (real_2_quad_real (z, 1), add_quad_real (QUAD_REAL_ONE, mul_quad_real (z, z), 0));
1349 }
1350 if (k != 0) {
1351 return neg_quad_real (z);
1352 } else {
1353 return z;
1354 }
1355 }
1356
1357 static QUAD_T c_tan (QUAD_T z)
1358 {
1359 QUAD_T s, f, d;
1360 int m;
1361 unt_2 k;
1362 if (getexp_quad_real (&z) < QUAD_REAL_K_LIN)
1363 return z;
1364 s = neg_quad_real (mul_quad_real (z, z));
1365 for (k = 1; k <= FLT256_LEN && s.value[k] == 0; k++);
1366 if ((sigerr_quad_real (s.value[0] == 0xffff && k > FLT256_LEN, QUAD_REAL_FPOFLOW, NULL))) {
1367 return QUAD_REAL_ZERO;
1368 } else {
1369 f = QUAD_REAL_ZERO;
1370 for (d = int_to_quad_real (m = QUAD_REAL_MS_TRG); m > 1;) {
1371 f = div_quad_real (s, _add_quad_real_ (d, f));
1372 d = int_to_quad_real (m -= 2);
1373 }
1374 return div_quad_real (z, _add_quad_real_ (d, f));
1375 }
1376 }
1377
1378 static QUAD_T rred (QUAD_T z, int kf, int *ps)
1379 {
1380 QUAD_T is, q;
1381 if (getsgn_quad_real (&z)) {
1382 z = neg_quad_real (z);
1383 is = QUAD_REAL_ONE;
1384 } else {
1385 is = QUAD_REAL_ZERO;
1386 }
1387 z = fmod_quad_real (z, QUAD_REAL_PI, &q);
1388 if (kf == 't') {
1389 q = is;
1390 } else if (kf == 's') {
1391 q = _add_quad_real_ (q, is);
1392 }
1393 if (real_cmp_quad_real (&z, &QUAD_REAL_PI2) == 1) {
1394 z = add_quad_real (QUAD_REAL_PI, z, 1);
1395 if (kf == 'c' || kf == 't') {
1396 q = add_quad_real (q, QUAD_REAL_ONE, 0);
1397 }
1398 }
1399 *ps = (odd_quad_real (q)) ? 1 : 0;
1400 return z;
1401 }
1402
1403 // VIF additions (REAL*32)
1404
1405 //! @brief cotan_quad_real
1406
1407 QUAD_T cotan_quad_real (QUAD_T x)
1408 {
1409 return div_quad_real (QUAD_REAL_ONE, tan_quad_real (x));
1410 }
1411
1412 //! @brief acotan_quad_real
1413
1414 QUAD_T acotan_quad_real (QUAD_T x)
1415 {
1416 return atan_quad_real (div_quad_real (QUAD_REAL_ONE, x));
1417 }
1418
1419 //! @brief sgn_quad_real_2
1420
1421 QUAD_T sgn_quad_real_2 (QUAD_T a, QUAD_T b)
1422 {
1423 QUAD_T x = (getsgn_quad_real (&a) == 0 ? a : neg_quad_real (a));
1424 return (getsgn_quad_real (&b) == 0 ? x : neg_quad_real (x));
1425 }
1426
1427 #endif
