1 /*
2 * go-quad.c: Extended precision routines.
3 *
4 * Authors
5 * Morten Welinder <terra@gnome.org>
6 *
7 * This library is free software; you can redistribute it and/or
8 * modify it under the terms of the GNU General Public License as
9 * published by the Free Software Foundation; either version 2 of the
10 * License, or (at your option) version 3.
11 *
12 * This library is distributed in the hope that it will be useful, but
13 * WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 * Library General Public License for more details.
16 *
17 * You should have received a copy of the GNU Library General Public
18 * License along with this library; if not, write to the Free Software
19 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
20 * USA.
21 *
22 * This follows "A Floating-Point Technique for Extending the Available
23 * Precision" by T. J. Dekker in _Numerische Mathematik_ 18.
24 * Springer Verlag 1971.
25 *
26 * Note: for this to work, the processor must evaluate with the right
27 * precision. For ix86 that means trouble as the default is to evaluate
28 * with long-double precision internally. We solve this by setting the
29 * relevant x87 control flag.
30 *
31 * Note: the compiler should not rearrange expressions.
32 */
33
34 #define GO_QUAD_IMPL
35
36 #include <goffice/goffice-config.h>
37 #include <goffice/goffice.h>
38 #include <math.h>
39
40 /* Normalize cpu id. */
41 #if !defined(i386) && (defined(__i386__) || defined(__i386))
42 #define i386 1
43 #endif
44
45 #ifndef DOUBLE
46
47 #define DEFINE_COMMON
48
49 #ifdef i386
50 #ifdef HAVE_FPU_CONTROL_H
51 #include <fpu_control.h>
52 #define USE_FPU_CONTROL
53 #elif defined(__GNUC__)
54 /* The next few lines from glibc licensed under lpgl 2.1 */
55 /* FPU control word bits. i387 version.
56 Copyright (C) 1993,1995-1998,2000,2001,2003 Free Software Foundation, Inc. */
57 #define _FPU_EXTENDED 0x300 /* libm requires double extended precision. */
58 #define _FPU_DOUBLE 0x200
59 #define _FPU_SINGLE 0x0
60 typedef unsigned int fpu_control_t __attribute__ ((__mode__ (__HI__)));
61 #define _FPU_GETCW(cw) __asm__ __volatile__ ("fnstcw %0" : "=m" (*&cw))
62 #define _FPU_SETCW(cw) __asm__ __volatile__ ("fldcw %0" : : "m" (*&cw))
63 #define USE_FPU_CONTROL
64 #endif
65 #endif
66
67 #define QUAD SUFFIX(GOQuad)
68
69 #define DOUBLE double
70 #define SUFFIX(_n) _n
71 #define DOUBLE_MANT_DIG DBL_MANT_DIG
72 #define DOUBLE_EPSILON DBL_EPSILON
73
74 #ifdef GOFFICE_WITH_LONG_DOUBLE
75 #include "go-quad.c"
76 #undef DEFINE_COMMON
77 #undef DOUBLE
78 #undef SUFFIX
79 #undef DOUBLE_MANT_DIG
80 #undef DOUBLE_EPSILON
81 #define DOUBLE long double
82 #define SUFFIX(_n) _n ## l
83 #define DOUBLE_MANT_DIG LDBL_MANT_DIG
84 #define DOUBLE_EPSILON LDBL_EPSILON
85 #endif
86
87 #endif
88
89 gboolean
SUFFIX(go_quad_functional)90 SUFFIX(go_quad_functional) (void)
91 {
92 if (FLT_RADIX != 2)
93 return FALSE;
94
95 #ifdef i386
96 if (sizeof (DOUBLE) != sizeof (double))
97 return TRUE;
98
99 #ifdef USE_FPU_CONTROL
100 return TRUE;
101 #else
102 return FALSE;
103 #endif
104 #else
105 return TRUE;
106 #endif
107 }
108
109 static guint SUFFIX(go_quad_depth) = 0;
110
111 static DOUBLE SUFFIX(CST);
112
113 #ifdef DEFINE_COMMON
114 /*
115 * Store constants in a way doesn't depend on the layout of DOUBLE. We use
116 * ~400 bits of data in the tables below -- that's way more than needed even
117 * for sunos' long double.
118 */
119
120 static const guint8 pi_hex_digits[] = {
121 0x03, 0x24, 0x3f, 0x6a, 0x88, 0x85, 0xa3, 0x08,
122 0xd3, 0x13, 0x19, 0x8a, 0x2e, 0x03, 0x70, 0x73,
123 0x44, 0xa4, 0x09, 0x38, 0x22, 0x29, 0x9f, 0x31,
124 0xd0, 0x08, 0x2e, 0xfa, 0x98, 0xec, 0x4e, 0x6c,
125 0x89, 0x45, 0x28, 0x21, 0xe6, 0x38, 0xd0, 0x13,
126 0x77, 0xbe, 0x54, 0x66, 0xcf, 0x34, 0xe9, 0x0c,
127 0x6c, 0xc0, 0xac
128 };
129
130 static const guint8 e_hex_digits[] = {
131 0x02, 0xb7, 0xe1, 0x51, 0x62, 0x8a, 0xed, 0x2a,
132 0x6a, 0xbf, 0x71, 0x58, 0x80, 0x9c, 0xf4, 0xf3,
133 0xc7, 0x62, 0xe7, 0x16, 0x0f, 0x38, 0xb4, 0xda,
134 0x56, 0xa7, 0x84, 0xd9, 0x04, 0x51, 0x90, 0xcf,
135 0xef, 0x32, 0x4e, 0x77, 0x38, 0x92, 0x6c, 0xfb,
136 0xe5, 0xf4, 0xbf, 0x8d, 0x8d, 0x8c, 0x31, 0xd7,
137 0x63, 0xda, 0x06
138 };
139
140 static const guint8 ln2_hex_digits[] = {
141 0xb1, 0x72, 0x17, 0xf7, 0xd1, 0xcf, 0x79, 0xab,
142 0xc9, 0xe3, 0xb3, 0x98, 0x03, 0xf2, 0xf6, 0xaf,
143 0x40, 0xf3, 0x43, 0x26, 0x72, 0x98, 0xb6, 0x2d,
144 0x8a, 0x0d, 0x17, 0x5b, 0x8b, 0xaa, 0xfa, 0x2b,
145 0xe7, 0xb8, 0x76, 0x20, 0x6d, 0xeb, 0xac, 0x98,
146 0x55, 0x95, 0x52, 0xfb, 0x4a, 0xfa, 0x1b, 0x10,
147 0xed, 0x2e
148 };
149
150 static const guint8 sqrt2_hex_digits[] = {
151 0x01, 0x6a, 0x09, 0xe6, 0x67, 0xf3, 0xbc, 0xc9,
152 0x08, 0xb2, 0xfb, 0x13, 0x66, 0xea, 0x95, 0x7d,
153 0x3e, 0x3a, 0xde, 0xc1, 0x75, 0x12, 0x77, 0x50,
154 0x99, 0xda, 0x2f, 0x59, 0x0b, 0x06, 0x67, 0x32,
155 0x2a, 0x95, 0xf9, 0x06, 0x08, 0x75, 0x71, 0x45,
156 0x87, 0x51, 0x63, 0xfc, 0xdf, 0xb9, 0x07, 0xb6,
157 0x72, 0x1e, 0xe9
158 };
159
160 static const guint8 euler_hex_digits[] = {
161 0x93, 0xc4, 0x67, 0xe3, 0x7d, 0xb0, 0xc7, 0xa4,
162 0xd1, 0xbe, 0x3f, 0x81, 0x01, 0x52, 0xcb, 0x56,
163 0xa1, 0xce, 0xcc, 0x3a, 0xf6, 0x5c, 0xc0, 0x19,
164 0x0c, 0x03, 0xdf, 0x34, 0x70, 0x9a, 0xff, 0xbd,
165 0x8e, 0x4b, 0x59, 0xfa, 0x03, 0xa9, 0xf0, 0xee,
166 0xd0, 0x64, 0x9c, 0xcb, 0x62, 0x10, 0x57, 0xd1,
167 0x10, 0x56
168 };
169 #endif
170
171 /**
172 * go_quad_start:
173 *
174 * Initializes #GOQuad arithmetic. Any use of #GOQuad must occur between calls
175 * to go_quad_start() and go_quad_end().
176 * Returns: (transfer full): a pointer to pass to go_quad_end() when done.
177 **/
178 /**
179 * go_quad_startl:
180 *
181 * Initializes #GOQuadl arithmetic. Any use of #GOQuadl must occur between calls
182 * to go_quad_startl() and go_quad_endl().
183 * Returns: (transfer full): a pointer to pass to go_quad_endl() when done.
184 **/
185 void *
SUFFIX(go_quad_start)186 SUFFIX(go_quad_start) (void)
187 {
188 void *res = NULL;
189 static gboolean first = TRUE;
190
191 if (SUFFIX(go_quad_depth)++ > 0)
192 return NULL;
193
194 if (!SUFFIX(go_quad_functional) () && first)
195 g_warning ("quad precision math may not be completely accurate.");
196
197 #ifdef i386
198 if (sizeof (DOUBLE) == sizeof (double)) {
199 #ifdef USE_FPU_CONTROL
200 fpu_control_t state, newstate;
201 fpu_control_t mask =
202 _FPU_EXTENDED | _FPU_DOUBLE | _FPU_SINGLE;
203
204 _FPU_GETCW (state);
205 res = g_memdup (&state, sizeof (state));
206
207 newstate = (state & ~mask) | _FPU_DOUBLE;
208 _FPU_SETCW (newstate);
209 #else
210 /* Hope for the best. */
211 #endif
212 }
213 #endif
214
215 if (first) {
216 first = FALSE;
217 SUFFIX(CST) = 1 + SUFFIX(ldexp) (1.0, (DOUBLE_MANT_DIG + 1) / 2);
218 SUFFIX(go_quad_constant8) (&SUFFIX(go_quad_pi),
219 pi_hex_digits,
220 G_N_ELEMENTS (pi_hex_digits),
221 256.0,
222 256.0);
223
224 SUFFIX(go_quad_constant8) (&SUFFIX(go_quad_2pi),
225 pi_hex_digits,
226 G_N_ELEMENTS (pi_hex_digits),
227 256.0,
228 512.0);
229
230 SUFFIX(go_quad_constant8) (&SUFFIX(go_quad_e),
231 e_hex_digits,
232 G_N_ELEMENTS (e_hex_digits),
233 256.0,
234 256.0);
235
236 SUFFIX(go_quad_constant8) (&SUFFIX(go_quad_ln2),
237 ln2_hex_digits,
238 G_N_ELEMENTS (ln2_hex_digits),
239 256.0,
240 1);
241
242 SUFFIX(go_quad_constant8) (&SUFFIX(go_quad_sqrt2),
243 sqrt2_hex_digits,
244 G_N_ELEMENTS (sqrt2_hex_digits),
245 256.0,
246 256.0);
247
248 SUFFIX(go_quad_constant8) (&SUFFIX(go_quad_euler),
249 euler_hex_digits,
250 G_N_ELEMENTS (euler_hex_digits),
251 256.0,
252 1);
253 }
254
255 return res;
256 }
257
258 /**
259 * go_quad_end:
260 * @state: state pointer from go_quad_start.
261 *
262 * This ends a section of quad precision arithmetic.
263 **/
264 /**
265 * go_quad_endl:
266 * @state: state pointer from go_quad_startl.
267 *
268 * This ends a section of quad precision arithmetic.
269 **/
270 void
SUFFIX(go_quad_end)271 SUFFIX(go_quad_end) (void *state)
272 {
273 SUFFIX(go_quad_depth)--;
274 if (!state)
275 return;
276
277 #ifdef i386
278 #ifdef USE_FPU_CONTROL
279 _FPU_SETCW (*(fpu_control_t*)state);
280 #endif
281 #endif
282
283 g_free (state);
284 }
285
286 const QUAD SUFFIX(go_quad_zero) = { 0, 0 };
287 const QUAD SUFFIX(go_quad_one) = { 1, 0 };
288 /*
289 * The following are non-const so we can initialize them. However,
290 * from other compilation units there are const. My reading of C99
291 * Section 6.2.7 says that is allowed.
292 */
293 QUAD SUFFIX(go_quad_pi);
294 QUAD SUFFIX(go_quad_2pi);
295 QUAD SUFFIX(go_quad_e);
296 QUAD SUFFIX(go_quad_ln2);
297 QUAD SUFFIX(go_quad_sqrt2);
298 QUAD SUFFIX(go_quad_euler);
299
300 /**
301 * go_quad_init:
302 * @res: (out): result location
303 * @h: a double precision value
304 *
305 * This stores the value @h in @res. As an exception, this may be called
306 * outside go_quad_start and go_quad_end sections.
307 **/
308 /**
309 * go_quad_initl:
310 * @res: (out): result location
311 * @h: a double precision value
312 *
313 * This stores the value @h in @res. As an exception, this may be called
314 * outside go_quad_startl and go_quad_endl sections.
315 **/
316 void
SUFFIX(go_quad_init)317 SUFFIX(go_quad_init) (QUAD *res, DOUBLE h)
318 {
319 res->h = h;
320 res->l = 0;
321 }
322
323 /**
324 * go_quad_value:
325 * @a: quad-precision value
326 *
327 * Returns: closest double precision value to @a. As an exception,
328 * this may be called outside go_quad_start and go_quad_end sections.
329 **/
330 /**
331 * go_quad_valuel:
332 * @a: quad-precision value
333 *
334 * Returns: closest double precision value to @a. As an exception,
335 * this may be called outside go_quad_startl and go_quad_endl sections.
336 **/
337 DOUBLE
SUFFIX(go_quad_value)338 SUFFIX(go_quad_value) (const QUAD *a)
339 {
340 return a->h + a->l;
341 }
342
343 /**
344 * go_quad_add:
345 * @res: (out): result location
346 * @a: quad-precision value
347 * @b: quad-precision value
348 *
349 * This function adds @a and @b, storing the result in @res.
350 **/
351 /**
352 * go_quad_addl:
353 * @res: (out): result location
354 * @a: quad-precision value
355 * @b: quad-precision value
356 *
357 * This function adds @a and @b, storing the result in @res.
358 **/
359 void
SUFFIX(go_quad_add)360 SUFFIX(go_quad_add) (QUAD *res, const QUAD *a, const QUAD *b)
361 {
362 DOUBLE r = a->h + b->h;
363 DOUBLE s = SUFFIX(fabs) (a->h) > SUFFIX(fabs) (b->h)
364 ? a->h - r + b->h + b->l + a->l
365 : b->h - r + a->h + a->l + b->l;
366 res->h = r + s;
367 res->l = r - res->h + s;
368
369 g_return_if_fail (SUFFIX(go_quad_depth) > 0);
370 }
371
372 /**
373 * go_quad_sub:
374 * @res: (out): result location
375 * @a: quad-precision value
376 * @b: quad-precision value
377 *
378 * This function subtracts @a and @b, storing the result in @res.
379 **/
380 /**
381 * go_quad_subl:
382 * @res: (out): result location
383 * @a: quad-precision value
384 * @b: quad-precision value
385 *
386 * This function subtracts @a and @b, storing the result in @res.
387 **/
388 void
SUFFIX(go_quad_sub)389 SUFFIX(go_quad_sub) (QUAD *res, const QUAD *a, const QUAD *b)
390 {
391 DOUBLE r = a->h - b->h;
392 DOUBLE s = SUFFIX(fabs) (a->h) > SUFFIX(fabs) (b->h)
393 ? +a->h - r - b->h - b->l + a->l
394 : -b->h - r + a->h + a->l - b->l;
395 res->h = r + s;
396 res->l = r - res->h + s;
397 }
398
399
400 #define SPLIT1(x,h,t) do { \
401 DOUBLE p = x * SUFFIX(CST); \
402 if (!SUFFIX(go_finite) (p) && SUFFIX(go_finite)(x)) { \
403 x *= DOUBLE_EPSILON; \
404 p = x * SUFFIX(CST); \
405 h = x - p + p; \
406 t = x - h; \
407 h *= (1 / DOUBLE_EPSILON); \
408 t *= (1 / DOUBLE_EPSILON); \
409 } else { \
410 h = x - p + p; \
411 t = x - h; \
412 } \
413 } while (0)
414
415 /**
416 * go_quad_mul12:
417 * @res: (out): result location
418 * @x: double precision value
419 * @y: double precision value
420 *
421 * This function multiplies @x and @y, storing the result in @res with full
422 * quad precision.
423 **/
424 /**
425 * go_quad_mul12l:
426 * @res: (out): result location
427 * @x: double precision value
428 * @y: double precision value
429 *
430 * This function multiplies @x and @y, storing the result in @res with full
431 * quad precision.
432 **/
433 void
SUFFIX(go_quad_mul12)434 SUFFIX(go_quad_mul12) (QUAD *res, DOUBLE x, DOUBLE y)
435 {
436 DOUBLE hx, tx, hy, ty, p, q;
437
438 SPLIT1 (x, hx, tx);
439 SPLIT1 (y, hy, ty);
440
441 p = hx * hy;
442 q = hx * ty + tx * hy;
443 res->h = p + q;
444 res->l = p - res->h + q + tx * ty;
445 }
446
447 #undef SPLIT1
448
449
450 /**
451 * go_quad_mul:
452 * @res: (out): result location
453 * @a: quad-precision value
454 * @b: quad-precision value
455 *
456 * This function multiplies @a and @b, storing the result in @res.
457 **/
458 /**
459 * go_quad_mull:
460 * @res: (out): result location
461 * @a: quad-precision value
462 * @b: quad-precision value
463 *
464 * This function multiplies @a and @b, storing the result in @res.
465 **/
466 void
SUFFIX(go_quad_mul)467 SUFFIX(go_quad_mul) (QUAD *res, const QUAD *a, const QUAD *b)
468 {
469 QUAD c;
470 SUFFIX(go_quad_mul12) (&c, a->h, b->h);
471 c.l = a->h * b->l + a->l * b->h + c.l;
472 res->h = c.h + c.l;
473 res->l = c.h - res->h + c.l;
474 }
475
476 /**
477 * go_quad_div:
478 * @res: (out): result location
479 * @a: quad-precision value
480 * @b: quad-precision value
481 *
482 * This function divides @a and @b, storing the result in @res.
483 **/
484 /**
485 * go_quad_divl:
486 * @res: (out): result location
487 * @a: quad-precision value
488 * @b: quad-precision value
489 *
490 * This function divides @a and @b, storing the result in @res.
491 **/
492 void
SUFFIX(go_quad_div)493 SUFFIX(go_quad_div) (QUAD *res, const QUAD *a, const QUAD *b)
494 {
495 QUAD c, u;
496 c.h = a->h / b->h;
497 SUFFIX(go_quad_mul12) (&u, c.h, b->h);
498 c.l = (a->h - u.h - u.l + a->l - c.h * b->l) / b->h;
499 res->h = c.h + c.l;
500 res->l = c.h - res->h + c.l;
501 }
502
503 /**
504 * go_quad_sqrt:
505 * @res: (out): result location
506 * @a: quad-precision value
507 *
508 * This function takes the square root of @a, storing the result in @res.
509 **/
510 /**
511 * go_quad_sqrtl:
512 * @res: (out): result location
513 * @a: quad-precision value
514 *
515 * This function takes the square root of @a, storing the result in @res.
516 **/
517 void
SUFFIX(go_quad_sqrt)518 SUFFIX(go_quad_sqrt) (QUAD *res, const QUAD *a)
519 {
520 if (a->h > 0) {
521 QUAD c, u;
522 c.h = SUFFIX(sqrt) (a->h);
523 SUFFIX(go_quad_mul12) (&u, c.h, c.h);
524 c.l = (a->h - u.h - u.l + a->l) * 0.5 / c.h;
525 res->h = c.h + c.l;
526 res->l = c.h - res->h + c.l;
527 } else
528 res->h = res->l = 0;
529 }
530
531 /**
532 * go_quad_floor:
533 * @res: (out): result location
534 * @a: quad-precision value
535 *
536 * This function takes the floor of @a, storing the result in @res.
537 **/
538 /**
539 * go_quad_floorl:
540 * @res: (out): result location
541 * @a: quad-precision value
542 *
543 * This function takes the floor of @a, storing the result in @res.
544 **/
545 void
SUFFIX(go_quad_floor)546 SUFFIX(go_quad_floor) (QUAD *res, const QUAD *a)
547 {
548 QUAD qh, ql, q, r;
549
550 SUFFIX(go_quad_init) (&qh, SUFFIX(floor)(a->h));
551 SUFFIX(go_quad_init) (&ql, SUFFIX(floor)(a->l));
552 SUFFIX(go_quad_add) (&q, &qh, &ql);
553
554 /* Due to dual floors, we might be off by one. */
555 SUFFIX(go_quad_sub) (&r, a, &q);
556 if (SUFFIX(go_quad_value) (&r) < 0)
557 SUFFIX(go_quad_sub) (res, &q, &SUFFIX(go_quad_one));
558 else {
559 SUFFIX(go_quad_sub) (&r, &r, &SUFFIX(go_quad_one));
560 if (SUFFIX(go_quad_value) (&r) < 0)
561 *res = q;
562 else
563 SUFFIX(go_quad_add) (res, &q, &SUFFIX(go_quad_one));
564 }
565 }
566
567 /**
568 * go_quad_dot_product:
569 * @res: (out): result location
570 * @a: (array length=n): vector of quad-precision values
571 * @b: (array length=n): vector of quad-precision values
572 * @n: length of vectors.
573 **/
574 /**
575 * go_quad_dot_productl:
576 * @res: (out): result location
577 * @a: (array length=n): vector of quad-precision values
578 * @b: (array length=n): vector of quad-precision values
579 * @n: length of vectors.
580 **/
581 void
SUFFIX(go_quad_dot_product)582 SUFFIX(go_quad_dot_product) (QUAD *res, const QUAD *a, const QUAD *b, int n)
583 {
584 int i;
585 SUFFIX(go_quad_init) (res, 0);
586 for (i = 0; i < n; i++) {
587 QUAD d;
588 SUFFIX(go_quad_mul) (&d, &a[i], &b[i]);
589 SUFFIX(go_quad_add) (res, res, &d);
590 }
591 }
592
593 /**
594 * go_quad_constant8:
595 * @res: (out): result location
596 * @data: (array length=n): vector of digits
597 * @base: base of vector's elements
598 * @n: length of digit vector.
599 * @scale: scaling value after interpreting digits
600 *
601 * This function interprets a vector of digits in a given base as a
602 * quad-precision value. It is mostly meant for internal use.
603 **/
604 /**
605 * go_quad_constant8l:
606 * @res: (out): result location
607 * @data: (array length=n): vector of digits
608 * @base: base of vector's elements
609 * @n: length of digit vector.
610 * @scale: scaling value after interpreting digits
611 *
612 * This function interprets a vector of digits in a given base as a
613 * quad-precision value. It is mostly meant for internal use.
614 **/
615 void
SUFFIX(go_quad_constant8)616 SUFFIX(go_quad_constant8) (QUAD *res, const guint8 *data, gsize n,
617 DOUBLE base, DOUBLE scale)
618 {
619 QUAD qbase, q;
620
621 *res = SUFFIX(go_quad_zero);
622 SUFFIX(go_quad_init) (&qbase, base);
623
624 while (n-- > 0) {
625 SUFFIX(go_quad_init) (&q, data[n]);
626 SUFFIX(go_quad_add) (res, res, &q);
627 SUFFIX(go_quad_div) (res, res, &qbase);
628 }
629
630 SUFFIX(go_quad_init) (&q, scale);
631 SUFFIX(go_quad_mul) (res, res, &q);
632 }
633
634 static void
SUFFIX(rescale2)635 SUFFIX(rescale2) (QUAD *x, DOUBLE *e)
636 {
637 int xe;
638
639 (void)SUFFIX(frexp) (SUFFIX(go_quad_value) (x), &xe);
640 if (xe != 0) {
641 QUAD qs;
642 SUFFIX(go_quad_init) (&qs, SUFFIX(ldexp) (1.0, -xe));
643 SUFFIX(go_quad_mul) (x, x, &qs);
644 *e += xe;
645 }
646 }
647
648
649 static void
SUFFIX(go_quad_pow_int)650 SUFFIX(go_quad_pow_int) (QUAD *res, DOUBLE *exp2, const QUAD *x, const QUAD *y)
651 {
652 QUAD xn;
653 DOUBLE dy;
654 DOUBLE xe = 0;
655
656 xn = *x;
657 *exp2 = 0;
658
659 dy = SUFFIX(go_quad_value) (y);
660 g_return_if_fail (dy >= 0);
661
662 *res = SUFFIX(go_quad_one);
663 SUFFIX(rescale2) (&xn, &xe);
664
665 while (dy > 0) {
666 if (SUFFIX(fmod) (dy, 2) > 0) {
667 SUFFIX(go_quad_mul) (res, res, &xn);
668 *exp2 += xe;
669 SUFFIX(rescale2) (res, exp2);
670 dy--;
671 if (dy == 0)
672 break;
673 }
674 dy /= 2;
675 SUFFIX(go_quad_mul) (&xn, &xn, &xn);
676 xe *= 2;
677 SUFFIX(rescale2) (&xn, &xe);
678 }
679 }
680
681 static void
SUFFIX(go_quad_sqrt1pm1)682 SUFFIX(go_quad_sqrt1pm1) (QUAD *res, const QUAD *a)
683 {
684 QUAD x0, x1, d;
685
686 SUFFIX(go_quad_add) (&x0, a, &SUFFIX(go_quad_one));
687 SUFFIX(go_quad_sqrt) (&x0, &x0);
688 SUFFIX(go_quad_sub) (&x0, &x0, &SUFFIX(go_quad_one));
689
690 /* Newton step. */
691 SUFFIX(go_quad_mul) (&x1, &x0, &x0);
692 SUFFIX(go_quad_add) (&x1, &x1, a);
693 SUFFIX(go_quad_add) (&d, &x0, &SUFFIX(go_quad_one));
694 SUFFIX(go_quad_add) (&d, &d, &d);
695 SUFFIX(go_quad_div) (&x1, &x1, &d);
696
697 *res = x1;
698 }
699
700 /*
701 * Compute pow(x,y) assuming y is in [0;1[. If r1m is true,
702 * compute 1-pow(x,y) instead.
703 */
704
705 static void
SUFFIX(go_quad_pow_frac)706 SUFFIX(go_quad_pow_frac) (QUAD *res, const QUAD *x, const QUAD *y,
707 gboolean r1m)
708 {
709 QUAD qx, qr, qy = *y;
710 DOUBLE dy;
711 gboolean x1m;
712
713 g_return_if_fail (SUFFIX(go_quad_value) (y) >= 0);
714 g_return_if_fail (SUFFIX(go_quad_value) (y) <= 1); /* =1 might happen */
715
716 /*
717 * "1m" mode refers to keeping 1-v instead of just v.
718 */
719 x1m = SUFFIX(fabs) (SUFFIX(go_quad_value) (x)) >= 0.5;
720 if (x1m) {
721 SUFFIX(go_quad_sub) (&qx, x, &SUFFIX(go_quad_one));
722 } else {
723 qx = *x;
724 }
725
726 qr = r1m ? SUFFIX(go_quad_zero) : SUFFIX(go_quad_one);
727
728 while ((dy = SUFFIX(go_quad_value) (&qy)) > 0) {
729 SUFFIX(go_quad_add) (&qy, &qy, &qy);
730 if (x1m)
731 SUFFIX(go_quad_sqrt1pm1) (&qx, &qx);
732 else {
733 SUFFIX(go_quad_sqrt) (&qx, &qx);
734 if (SUFFIX(go_quad_value) (&qx) >= 0.5) {
735 x1m = TRUE;
736 SUFFIX(go_quad_sub) (&qx, &qx, &SUFFIX(go_quad_one));
737 }
738 }
739 if (dy >= 0.5) {
740 QUAD qp;
741 SUFFIX(go_quad_sub) (&qy, &qy, &SUFFIX(go_quad_one));
742 SUFFIX(go_quad_mul) (&qp, &qx, &qr);
743 if (x1m && r1m) {
744 SUFFIX(go_quad_add) (&qr, &qr, &qp);
745 SUFFIX(go_quad_add) (&qr, &qr, &qx);
746 } else if (x1m) {
747 SUFFIX(go_quad_add) (&qr, &qr, &qp);
748 } else if (r1m) {
749 QUAD qy;
750 SUFFIX(go_quad_sub) (&qy, &qx, &SUFFIX(go_quad_one));
751 SUFFIX(go_quad_add) (&qr, &qp, &qy);
752 } else {
753 qr = qp;
754 }
755 }
756 }
757
758 *res = qr;
759 }
760
761 /*
762 * This computes pow(x,y) in the following way:
763 *
764 * 1. y is considered the sum of a number of one-bit values.
765 * 2. For each y bit in the integer part of y, the corresponding x^y
766 * is computed by repeated squaring.
767 * 3. For each y bit in the fractional part of y, the corresponding x^y
768 * is computed by repeating square rooting.
769 *
770 * Except that it's not quite as simple. Repeating square rooting of x
771 * will bring the value closer and closer to 1. That's no good, so we
772 * sometimes keep those values a 1-v instead of v. A square root in
773 * that representation is sqrt1pm1(v)=sqrt(1+v)-1.
774 *
775 * Finally we don't simply return one value. If the optional exp2 location
776 * is non-null, it is treated as a place to put a base 2 exponent with which
777 * to scale the returned result. This isn't much different from the exponent
778 * inside the floating point number, except that the range is *huge*. This
779 * function will happily compute exp(1e6):
780 * exp(1e+06) = 0.5143737638002868*2^1442696 [1.028747527600574*2^1442695]
781 * The bracked result is a reference value from Mathematica.
782 */
783
784 /**
785 * go_quad_pow:
786 * @res: (out): result location
787 * @exp2: (out): (allow-none): power-of-2 result scaling location
788 * @x: quad-precision value
789 * @y: quad-precision value
790 *
791 * This function computes @x to the power of @y, storing the result in @res.
792 * If the optional @exp2 is supplied, it is used to return a power of 2 by
793 * which the result should be scaled. This is useful to represent results
794 * much, much bigger than double precision can handle.
795 **/
796 /**
797 * go_quad_powl:
798 * @res: (out): result location
799 * @exp2: (out): (allow-none): power-of-2 result scaling location
800 * @x: quad-precision value
801 * @y: quad-precision value
802 *
803 * This function computes @x to the power of @y, storing the result in @res.
804 * If the optional @exp2 is supplied, it is used to return a power of 2 by
805 * which the result should be scaled. This is useful to represent results
806 * much, much bigger than double precision can handle.
807 **/
808 void
SUFFIX(go_quad_pow)809 SUFFIX(go_quad_pow) (QUAD *res, DOUBLE *exp2,
810 const QUAD *x, const QUAD *y)
811 {
812 DOUBLE dy, exp2ew;
813 QUAD qw, qf, qew, qef, qxm1;
814
815 dy = SUFFIX(go_quad_value) (y);
816
817 SUFFIX(go_quad_sub) (&qxm1, x, &SUFFIX(go_quad_one));
818 if (SUFFIX(go_quad_value) (&qxm1) == 0 || dy == 0) {
819 *res = SUFFIX(go_quad_one);
820 if (exp2) *exp2 = 0;
821 return;
822 }
823
824 SUFFIX(go_quad_floor) (&qw, y);
825 SUFFIX(go_quad_sub) (&qf, y, &qw);
826 if (SUFFIX(go_quad_value) (&qxm1) == 0 && dy > 0) {
827 gboolean wodd =
828 (SUFFIX(fmod)(SUFFIX(fabs)(qw.h),2) +
829 SUFFIX(fmod)(SUFFIX(fabs)(qw.l),2)) == 1;
830 if (SUFFIX(go_quad_value) (&qf) == 0 && wodd) {
831 /* 0 ^ (odd positive integer) */
832 *res = *x;
833 } else {
834 /* 0 ^ y, y positive, but not odd integer */
835 *res = SUFFIX(go_quad_zero);
836 }
837 if (exp2) *exp2 = 0;
838 return;
839 }
840
841 /* Lots of infinity/nan cases ignored */
842
843 if (dy < 0) {
844 QUAD my;
845 SUFFIX(go_quad_sub) (&my, &SUFFIX(go_quad_zero), y);
846 SUFFIX(go_quad_pow) (res, exp2, x, &my);
847 SUFFIX(go_quad_div) (res, &SUFFIX(go_quad_one), res);
848 if (exp2) *exp2 = 0 - *exp2;
849 return;
850 }
851
852 SUFFIX(go_quad_pow_int) (&qew, &exp2ew, x, &qw);
853 SUFFIX(go_quad_pow_frac) (&qef, x, &qf, FALSE);
854
855 SUFFIX(go_quad_mul) (res, &qew, &qef);
856 if (exp2)
857 *exp2 = exp2ew;
858 else {
859 QUAD qs;
860 int e = CLAMP (exp2ew, G_MININT, G_MAXINT);
861 SUFFIX(go_quad_init) (&qs, SUFFIX(ldexp)(1.0, e));
862 SUFFIX(go_quad_mul) (res, res, &qs);
863 }
864 }
865
866
867 /**
868 * go_quad_exp:
869 * @res: (out): result location
870 * @exp2: (out): (allow-none): power-of-2 result scaling location
871 * @a: quad-precision value
872 *
873 * This function computes the exponential function at @a, storing the result
874 * in @res. If the optional @exp2 is supplied, it is used to return a
875 * power of 2 by which the result should be scaled. This is useful to
876 * represent results much, much bigger than double precision can handle.
877 **/
878 /**
879 * go_quad_expl:
880 * @res: (out): result location
881 * @exp2: (out): (allow-none): power-of-2 result scaling location
882 * @a: quad-precision value
883 *
884 * This function computes the exponential function at @a, storing the result
885 * in @res. If the optional @exp2 is supplied, it is used to return a
886 * power of 2 by which the result should be scaled. This is useful to
887 * represent results much, much bigger than double precision can handle.
888 **/
889 void
SUFFIX(go_quad_exp)890 SUFFIX(go_quad_exp) (QUAD *res, DOUBLE *exp2, const QUAD *a)
891 {
892 SUFFIX(go_quad_pow) (res, exp2, &SUFFIX(go_quad_e), a);
893 }
894
895 /**
896 * go_quad_expm1:
897 * @res: (out): result location
898 * @a: quad-precision value
899 *
900 * This function computes the exponential function at @a with 1 subtracted,
901 * storing the difference in @res.
902 **/
903 /**
904 * go_quad_expm1l:
905 * @res: (out): result location
906 * @a: quad-precision value
907 *
908 * This function computes the exponential function at @a with 1 subtracted,
909 * storing the difference in @res.
910 **/
911 void
SUFFIX(go_quad_expm1)912 SUFFIX(go_quad_expm1) (QUAD *res, const QUAD *a)
913 {
914 DOUBLE da = SUFFIX(go_quad_value) (a);
915
916 if (!SUFFIX(go_finite) (da))
917 *res = *a;
918 else if (SUFFIX (fabs) (da) > 0.5) {
919 SUFFIX(go_quad_exp) (res, NULL, a);
920 SUFFIX(go_quad_sub) (res, res, &SUFFIX(go_quad_one));
921 } else if (da >= 0) {
922 SUFFIX(go_quad_pow_frac) (res, &SUFFIX(go_quad_e), a, TRUE);
923 } else {
924 QUAD ma, z, zp1;
925 SUFFIX(go_quad_sub) (&ma, &SUFFIX(go_quad_zero), a);
926 SUFFIX(go_quad_pow_frac) (&z, &SUFFIX(go_quad_e), &ma, TRUE);
927 SUFFIX(go_quad_add) (&zp1, &z, &SUFFIX(go_quad_one));
928 SUFFIX(go_quad_div) (res, &z, &zp1);
929 }
930 }
931
932 /**
933 * go_quad_log:
934 * @res: (out): result location
935 * @a: quad-precision value
936 *
937 * This function computes the natural logarithm at @a, storing the result
938 * in @res.
939 **/
940 /**
941 * go_quad_logl:
942 * @res: (out): result location
943 * @a: quad-precision value
944 *
945 * This function computes the natural logarithm at @a, storing the result
946 * in @res.
947 **/
948 void
SUFFIX(go_quad_log)949 SUFFIX(go_quad_log) (QUAD *res, const QUAD *a)
950 {
951 DOUBLE da = SUFFIX(go_quad_value) (a);
952
953 if (da == 0)
954 SUFFIX(go_quad_init) (res, SUFFIX(go_ninf));
955 else if (da < 0)
956 SUFFIX(go_quad_init) (res, SUFFIX(go_nan));
957 else if (!SUFFIX(go_finite) (da))
958 *res = *a;
959 else {
960 QUAD xi, yi, dx;
961 SUFFIX(go_quad_init) (&xi, SUFFIX(log) (da));
962
963 /* Newton step. */
964 SUFFIX(go_quad_exp) (&yi, NULL, &xi);
965 SUFFIX(go_quad_sub) (&dx, a, &yi);
966 SUFFIX(go_quad_div) (&dx, &dx, &yi);
967 SUFFIX(go_quad_add) (&xi, &xi, &dx);
968
969 *res = xi;
970 }
971 }
972
973 /**
974 * go_quad_hypot:
975 * @res: (out): result location
976 * @a: quad-precision value
977 * @b: quad-precision value
978 *
979 * This function computes the square root of @a^2 plugs @b^2, storing the
980 * result in @res.
981 **/
982 /**
983 * go_quad_hypotl:
984 * @res: (out): result location
985 * @a: quad-precision value
986 * @b: quad-precision value
987 *
988 * This function computes the square root of @a^2 plugs @b^2, storing the
989 * result in @res.
990 **/
991 void
SUFFIX(go_quad_hypot)992 SUFFIX(go_quad_hypot) (QUAD *res, const QUAD *a, const QUAD *b)
993 {
994 int e;
995 QUAD qa2, qb2, qn;
996
997 if (a->h == 0) {
998 res->h = SUFFIX(fabs)(b->h);
999 res->l = SUFFIX(fabs)(b->l);
1000 return;
1001 }
1002 if (b->h == 0) {
1003 res->h = SUFFIX(fabs)(a->h);
1004 res->l = SUFFIX(fabs)(a->l);
1005 return;
1006 }
1007
1008 /* Scale by power of 2 to protect against over- and underflow */
1009 (void)SUFFIX(frexp) (MAX (SUFFIX(fabs) (a->h), SUFFIX(fabs) (b->h)), &e);
1010
1011 qa2.h = SUFFIX(ldexp) (a->h, -e);
1012 qa2.l = SUFFIX(ldexp) (a->l, -e);
1013 SUFFIX(go_quad_mul) (&qa2, &qa2, &qa2);
1014
1015 qb2.h = SUFFIX(ldexp) (b->h, -e);
1016 qb2.l = SUFFIX(ldexp) (b->l, -e);
1017 SUFFIX(go_quad_mul) (&qb2, &qb2, &qb2);
1018
1019 SUFFIX(go_quad_add) (&qn, &qa2, &qb2);
1020 SUFFIX(go_quad_sqrt) (&qn, &qn);
1021 res->h = SUFFIX(ldexp) (qn.h, e);
1022 res->l = SUFFIX(ldexp) (qn.l, e);
1023 }
1024
1025 /* sqrt(1-a*a) helper */
1026 static void
SUFFIX(go_quad_ihypot)1027 SUFFIX(go_quad_ihypot) (QUAD *res, const QUAD *a)
1028 {
1029 QUAD qp;
1030
1031 SUFFIX(go_quad_mul) (&qp, a, a);
1032 SUFFIX(go_quad_sub) (&qp, &SUFFIX(go_quad_one), &qp);
1033 SUFFIX(go_quad_sqrt) (res, &qp);
1034 }
1035
1036
1037 #ifdef DEFINE_COMMON
1038 typedef enum {
1039 AGM_ARCSIN,
1040 AGM_ARCCOS,
1041 AGM_ARCTAN
1042 } AGM_Method;
1043 #endif
1044
1045 static void
SUFFIX(go_quad_agm_internal)1046 SUFFIX(go_quad_agm_internal) (QUAD *res, AGM_Method method, const QUAD *x)
1047 {
1048 QUAD g, gp, dk[20], dpk[20], qr, qrp, qnum;
1049 int n, k;
1050 gboolean converged = FALSE;
1051
1052 /*
1053 * This follows "An Algorithm for Computing Logarithms
1054 * and Arctangents" by B. C. Carlson in *Mathematics of
1055 * Computation*, Volume 26, Number 118, April 1972.
1056 *
1057 * If need be we can do log, arctanh, arcsinh, and
1058 * arccosh we the same code.
1059 */
1060
1061 qrp = SUFFIX(go_quad_zero);
1062
1063 switch (method) {
1064 case AGM_ARCSIN:
1065 g_return_if_fail (SUFFIX(fabs) (x->h) <= 1);
1066 SUFFIX(go_quad_ihypot) (&dpk[0], x);
1067 gp = SUFFIX(go_quad_one);
1068 qnum = *x;
1069 break;
1070 case AGM_ARCCOS:
1071 g_return_if_fail (SUFFIX(fabs) (x->h) <= 1);
1072 dpk[0] = *x;
1073 gp = SUFFIX(go_quad_one);
1074 SUFFIX(go_quad_ihypot) (&qnum, x);
1075 break;
1076 case AGM_ARCTAN:
1077 g_return_if_fail (SUFFIX(fabs) (x->h) <= 1);
1078 dpk[0] = SUFFIX(go_quad_one);
1079 SUFFIX(go_quad_hypot) (&gp, x, &SUFFIX(go_quad_one));
1080 qnum = *x;
1081 break;
1082 default:
1083 g_assert_not_reached ();
1084 }
1085
1086 for (n = 1; n < (int)G_N_ELEMENTS(dk); n++) {
1087 SUFFIX(go_quad_add) (&dk[0], &dpk[0], &gp);
1088 dk[0].h *= 0.5; dk[0].l *= 0.5;
1089
1090 SUFFIX(go_quad_mul) (&g, &dk[0], &gp);
1091 SUFFIX(go_quad_sqrt) (&g, &g);
1092
1093 for (k = 1; k <= n; k++) {
1094 QUAD f;
1095
1096 SUFFIX(go_quad_init) (&f, SUFFIX(ldexp) (1, -(2 * k)));
1097 SUFFIX(go_quad_mul) (&dk[k], &f, &dpk[k-1]);
1098 SUFFIX(go_quad_sub) (&dk[k], &dk[k-1], &dk[k]);
1099 SUFFIX(go_quad_sub) (&f, &SUFFIX(go_quad_one), &f);
1100 SUFFIX(go_quad_div) (&dk[k], &dk[k], &f);
1101 }
1102
1103 SUFFIX(go_quad_div) (&qr, &qnum, &dk[n]);
1104 SUFFIX(go_quad_sub) (&qrp, &qrp, &qr);
1105 if (SUFFIX(fabs)(qrp.h) <= SUFFIX(ldexp) (SUFFIX(fabs)(qr.h), -2 * (DOUBLE_MANT_DIG - 1))) {
1106 converged = TRUE;
1107 break;
1108 }
1109
1110 qrp = qr;
1111 gp = g;
1112 for (k = 0; k <= n; k++) dpk[k] = dk[k];
1113 }
1114
1115 if (!converged)
1116 g_warning ("go_quad_agm_internal(%.20g) failed to converge\n",
1117 (double)SUFFIX(go_quad_value) (x));
1118
1119 *res = qr;
1120 }
1121
1122 static gboolean
SUFFIX(go_quad_atan2_special)1123 SUFFIX(go_quad_atan2_special) (const QUAD *y, const QUAD *x, DOUBLE *f)
1124 {
1125 DOUBLE dy = SUFFIX(go_quad_value) (y);
1126 DOUBLE dx = SUFFIX(go_quad_value) (x);
1127
1128 if (dy == 0) {
1129 /* This assumes all zeros are +0. */
1130 *f = (dx >= 0 ? 0 : +1);
1131 return TRUE;
1132 }
1133
1134 if (dx == 0) {
1135 *f = (dy >= 0 ? 0.5 : -0.5);
1136 return TRUE;
1137 }
1138
1139 if (SUFFIX(fabs) (SUFFIX(fabs)(dx) - SUFFIX(fabs)(dy)) < 1e-10) {
1140 QUAD d;
1141 SUFFIX(go_quad_sub) (&d, x, y);
1142 if (d.h == 0) {
1143 *f = (dy >= 0 ? 0.25 : -0.75);
1144 return TRUE;
1145 }
1146 SUFFIX(go_quad_add) (&d, x, y);
1147 if (d.h == 0) {
1148 *f = (dy >= 0 ? +0.75 : -0.25);
1149 return TRUE;
1150 }
1151 }
1152
1153 return FALSE;
1154 }
1155
1156 /**
1157 * go_quad_atan2:
1158 * @res: (out): result location
1159 * @y: quad-precision value
1160 * @x: quad-precision value
1161 *
1162 * This function computes polar angle coordinate of the point (@x,@y), storing
1163 * the result in @res.
1164 **/
1165 /**
1166 * go_quad_atan2l:
1167 * @res: (out): result location
1168 * @y: quad-precision value
1169 * @x: quad-precision value
1170 *
1171 * This function computes polar angle coordinate of the point (@x,@y), storing
1172 * the result in @res.
1173 **/
1174 void
SUFFIX(go_quad_atan2)1175 SUFFIX(go_quad_atan2) (QUAD *res, const QUAD *y, const QUAD *x)
1176 {
1177 DOUBLE f;
1178 DOUBLE dy = SUFFIX(go_quad_value) (y);
1179 DOUBLE dx = SUFFIX(go_quad_value) (x);
1180 QUAD qr;
1181
1182 if (SUFFIX(go_quad_atan2_special) (y, x, &f)) {
1183 res->h = f * SUFFIX(go_quad_pi).h;
1184 res->l = f * SUFFIX(go_quad_pi).l;
1185 return;
1186 }
1187
1188 if (SUFFIX(fabs) (dx) >= SUFFIX(fabs) (dy)) {
1189 SUFFIX(go_quad_div) (&qr, y, x);
1190 SUFFIX(go_quad_agm_internal) (res, AGM_ARCTAN, &qr);
1191 } else {
1192 DOUBLE f;
1193 QUAD qa;
1194
1195 SUFFIX(go_quad_div) (&qr, x, y);
1196 SUFFIX(go_quad_agm_internal) (res, AGM_ARCTAN, &qr);
1197
1198 f = (qr.h >= 0) ? 0.5 : -0.5;
1199 qa.h = f * SUFFIX(go_quad_pi).h;
1200 qa.l = f * SUFFIX(go_quad_pi).l;
1201 SUFFIX(go_quad_sub) (res, &qa, res);
1202 }
1203
1204 if (dx < 0) {
1205 /* Correct for quadrants 2 and 3. */
1206 if (dy > 0)
1207 SUFFIX(go_quad_add) (res, res, &SUFFIX(go_quad_pi));
1208 else
1209 SUFFIX(go_quad_sub) (res, res, &SUFFIX(go_quad_pi));
1210 }
1211 }
1212
1213 /**
1214 * go_quad_atan2pi:
1215 * @res: (out): result location
1216 * @y: quad-precision value
1217 * @x: quad-precision value
1218 *
1219 * This function computes polar angle coordinate of the point (@x,@y) divided
1220 * by pi, storing the result in @res.
1221 **/
1222 /**
1223 * go_quad_atan2pil:
1224 * @res: (out): result location
1225 * @y: quad-precision value
1226 * @x: quad-precision value
1227 *
1228 * This function computes polar angle coordinate of the point (@x,@y) divided
1229 * by pi, storing the result in @res.
1230 **/
1231 void
SUFFIX(go_quad_atan2pi)1232 SUFFIX(go_quad_atan2pi) (QUAD *res, const QUAD *y, const QUAD *x)
1233 {
1234 DOUBLE f;
1235
1236 if (SUFFIX(go_quad_atan2_special) (y, x, &f)) {
1237 SUFFIX(go_quad_init) (res, f);
1238 return;
1239 }
1240
1241 /* Fallback. */
1242 SUFFIX(go_quad_atan2) (res, y, x);
1243 SUFFIX(go_quad_div) (res, res, &SUFFIX(go_quad_pi));
1244 }
1245
1246 static gboolean
SUFFIX(reduce_pi_half)1247 SUFFIX(reduce_pi_half) (QUAD *res, const QUAD *a, int *pk)
1248 {
1249 static QUAD pi_half;
1250 QUAD qa, qk, qh, qb;
1251 DOUBLE k;
1252 unsigned ui;
1253 static const DOUBLE pi_half_parts[] = {
1254 +0x1.921fb54442d18p+0,
1255 +0x1.1a62633145c04p-54,
1256 +0x1.707344a40938p-105,
1257 +0x1.114cf98e80414p-156,
1258 +0x1.bea63b139b224p-207,
1259 +0x1.14a08798e3404p-259,
1260 +0x1.bbdf2a33679a4p-312,
1261 +0x1.a431b302b0a6cp-363,
1262 +0x1.f25f14374fe1p-415,
1263 +0x1.ab6b6a8e122fp-466
1264 };
1265
1266 if (!SUFFIX(go_finite) (a->h))
1267 return TRUE;
1268
1269 if (SUFFIX(fabs) (a->h) > SUFFIX(ldexp) (1.0, DOUBLE_MANT_DIG)) {
1270 g_warning ("Reduced accuracy for very large trigonometric arguments");
1271 return TRUE;
1272 }
1273
1274 if (pi_half.h == 0) {
1275 pi_half.h = SUFFIX(go_quad_pi).h * 0.5;
1276 pi_half.l = SUFFIX(go_quad_pi).l * 0.5;
1277 }
1278
1279 SUFFIX(go_quad_div) (&qk, a, &pi_half);
1280 qh.h = 0.5; qh.l = 0;
1281 SUFFIX(go_quad_add) (&qk, &qk, &qh);
1282 SUFFIX(go_quad_floor) (&qk, &qk);
1283 k = SUFFIX(go_quad_value) (&qk);
1284 *pk = (int)(SUFFIX(fmod) (k, 4));
1285
1286 qa = *a;
1287 for (ui = 0; ui < G_N_ELEMENTS(pi_half_parts); ui++) {
1288 SUFFIX(go_quad_mul12) (&qb, pi_half_parts[ui], k);
1289 SUFFIX(go_quad_sub) (&qa, &qa, &qb);
1290 }
1291
1292 *res = qa;
1293
1294 return FALSE;
1295 }
1296
1297 static void
SUFFIX(reduce_half)1298 SUFFIX(reduce_half) (QUAD *res, const QUAD *a, int *pk)
1299 {
1300 static const QUAD half = { 0.5, 0 };
1301 int k = 0;
1302 QUAD qxr = *a;
1303
1304 if (a->h < 0) {
1305 QUAD aa;
1306 aa.h = -a->h; aa.l = -a->l;
1307 SUFFIX(reduce_half) (&qxr, &aa, &k);
1308 qxr.h = -qxr.h; qxr.l = -qxr.l;
1309 k = 4 - k;
1310 if (qxr.h <= -0.25 && qxr.l == 0) {
1311 SUFFIX(go_quad_add) (&qxr, &qxr, &half);
1312 k += 3;
1313 }
1314 } else {
1315 QUAD qdx;
1316
1317 /* Remove integer multiples of 2. */
1318 SUFFIX(go_quad_init) (&qdx, qxr.h - SUFFIX(fmod) (qxr.h, 2));
1319 SUFFIX(go_quad_sub) (&qxr, &qxr, &qdx);
1320
1321 /* Again, just in case it was really big. */
1322 SUFFIX(go_quad_init) (&qdx, qxr.h - SUFFIX(fmod) (qxr.h, 2));
1323 SUFFIX(go_quad_sub) (&qxr, &qxr, &qdx);
1324
1325 if (qxr.h >= 1) {
1326 SUFFIX(go_quad_sub) (&qxr, &qxr, &SUFFIX(go_quad_one));
1327 k += 2;
1328 }
1329 if (qxr.h >= 0.5) {
1330 SUFFIX(go_quad_sub) (&qxr, &qxr, &half);
1331 k++;
1332 }
1333 if (qxr.h > 0.25) {
1334 SUFFIX(go_quad_sub) (&qxr, &qxr, &half);
1335 k++;
1336 }
1337 }
1338
1339 *pk = (k & 3);
1340 *res = qxr;
1341 }
1342
1343 static void
SUFFIX(do_sin)1344 SUFFIX(do_sin) (QUAD *res, const QUAD *a, int k)
1345 {
1346 QUAD qr;
1347
1348 if (k & 1) {
1349 QUAD qn, qd, qq, aa;
1350
1351 aa.h = SUFFIX(fabs)(a->h);
1352 aa.l = SUFFIX(fabs)(a->l);
1353 SUFFIX(go_quad_init) (&qr, SUFFIX(cos) (aa.h));
1354
1355 /* Newton step */
1356 SUFFIX(go_quad_acos) (&qn, &qr);
1357 SUFFIX(go_quad_sub) (&qn, &qn, &aa);
1358 SUFFIX(go_quad_ihypot) (&qd, &qr);
1359 SUFFIX(go_quad_mul) (&qq, &qn, &qd);
1360 SUFFIX(go_quad_add) (&qr, &qr, &qq);
1361 } else {
1362 QUAD qn, qd, qq;
1363 SUFFIX(go_quad_init) (&qr, SUFFIX(sin) (a->h));
1364
1365 /* Newton step */
1366 SUFFIX(go_quad_asin) (&qn, &qr);
1367 SUFFIX(go_quad_sub) (&qn, &qn, a);
1368 SUFFIX(go_quad_ihypot) (&qd, &qr);
1369 SUFFIX(go_quad_mul) (&qq, &qn, &qd);
1370 SUFFIX(go_quad_sub) (&qr, &qr, &qq);
1371 }
1372
1373 if (k & 2) {
1374 qr.h = 0 - qr.h;
1375 qr.l = 0 - qr.l;
1376 }
1377
1378 *res = qr;
1379 }
1380
1381 static void
SUFFIX(do_sinpi)1382 SUFFIX(do_sinpi) (QUAD *res, const QUAD *a, int k)
1383 {
1384 QUAD qr;
1385
1386 if (a->h == 0)
1387 SUFFIX(go_quad_init) (&qr, k & 1);
1388 else if (a->h == 0.25 && a->l == 0)
1389 SUFFIX(go_quad_div) (&qr,
1390 &SUFFIX(go_quad_one),
1391 &SUFFIX(go_quad_sqrt2));
1392 else {
1393 QUAD api;
1394 SUFFIX(go_quad_mul) (&api, a, &SUFFIX(go_quad_pi));
1395 SUFFIX(do_sin) (&qr, &api, k & 1);
1396 }
1397
1398 if (k & 2) {
1399 qr.h = 0 - qr.h;
1400 qr.l = 0 - qr.l;
1401 }
1402
1403 *res = qr;
1404 }
1405
1406 /**
1407 * go_quad_sin:
1408 * @res: (out): result location
1409 * @a: quad-precision value
1410 *
1411 * This function computes the sine of @a, storing the result in @res.
1412 **/
1413 /**
1414 * go_quad_sinl:
1415 * @res: (out): result location
1416 * @a: quad-precision value
1417 *
1418 * This function computes the sine of @a, storing the result in @res.
1419 **/
1420 void
SUFFIX(go_quad_sin)1421 SUFFIX(go_quad_sin) (QUAD *res, const QUAD *a)
1422 {
1423 int k;
1424 QUAD a0;
1425
1426 if (SUFFIX(reduce_pi_half) (&a0, a, &k))
1427 SUFFIX(go_quad_init) (res, SUFFIX(sin) (a->h));
1428 else
1429 SUFFIX(do_sin) (res, &a0, k);
1430 }
1431
1432 /**
1433 * go_quad_sinpi:
1434 * @res: (out): result location
1435 * @a: quad-precision value
1436 *
1437 * This function computes the sine of @a times pi, storing the result in @res.
1438 * This is more accurate than actually doing the multiplication.
1439 **/
1440 /**
1441 * go_quad_sinpil:
1442 * @res: (out): result location
1443 * @a: quad-precision value
1444 *
1445 * This function computes the sine of @a times pi, storing the result in @res.
1446 * This is more accurate than actually doing the multiplication.
1447 **/
1448 void
SUFFIX(go_quad_sinpi)1449 SUFFIX(go_quad_sinpi) (QUAD *res, const QUAD *a)
1450 {
1451 int k;
1452 QUAD a0;
1453
1454 SUFFIX(reduce_half) (&a0, a, &k);
1455 SUFFIX(do_sinpi) (res, &a0, k);
1456 }
1457
1458 /**
1459 * go_quad_asin:
1460 * @res: (out): result location
1461 * @a: quad-precision value
1462 *
1463 * This function computes the arc sine of @a, storing the result in @res.
1464 **/
1465 /**
1466 * go_quad_asinl:
1467 * @res: (out): result location
1468 * @a: quad-precision value
1469 *
1470 * This function computes the arc sine of @a, storing the result in @res.
1471 **/
1472 void
SUFFIX(go_quad_asin)1473 SUFFIX(go_quad_asin) (QUAD *res, const QUAD *a)
1474 {
1475 QUAD aa, aam1;
1476
1477 aa.h = SUFFIX(fabs) (a->h);
1478 aa.l = SUFFIX(fabs) (a->l);
1479 SUFFIX(go_quad_sub) (&aam1, &aa, &SUFFIX(go_quad_one));
1480 if (aam1.h > 0) {
1481 SUFFIX(go_quad_init) (res, SUFFIX(go_nan));
1482 return;
1483 }
1484
1485 SUFFIX(go_quad_agm_internal) (res, AGM_ARCSIN, a);
1486 }
1487
1488 /**
1489 * go_quad_cos:
1490 * @res: (out): result location
1491 * @a: quad-precision value
1492 *
1493 * This function computes the cosine of @a, storing the result in @res.
1494 **/
1495 /**
1496 * go_quad_cosl:
1497 * @res: (out): result location
1498 * @a: quad-precision value
1499 *
1500 * This function computes the cosine of @a, storing the result in @res.
1501 **/
1502 void
SUFFIX(go_quad_cos)1503 SUFFIX(go_quad_cos) (QUAD *res, const QUAD *a)
1504 {
1505 int k;
1506 QUAD a0;
1507
1508 if (SUFFIX(reduce_pi_half) (&a0, a, &k))
1509 SUFFIX(go_quad_init) (res, SUFFIX(cos) (a->h));
1510 else
1511 SUFFIX(do_sin) (res, &a0, k + 1);
1512 }
1513
1514 /**
1515 * go_quad_cospi:
1516 * @res: (out): result location
1517 * @a: quad-precision value
1518 *
1519 * This function computes the cosine of @a times pi, storing the result in @res.
1520 * This is more accurate than actually doing the multiplication.
1521 **/
1522 /**
1523 * go_quad_cospil:
1524 * @res: (out): result location
1525 * @a: quad-precision value
1526 *
1527 * This function computes the cosine of @a times pi, storing the result in @res.
1528 * This is more accurate than actually doing the multiplication.
1529 **/
1530 void
SUFFIX(go_quad_cospi)1531 SUFFIX(go_quad_cospi) (QUAD *res, const QUAD *a)
1532 {
1533 int k;
1534 QUAD a0;
1535
1536 SUFFIX(reduce_half) (&a0, a, &k);
1537 SUFFIX(do_sinpi) (res, &a0, k + 1);
1538 }
1539
1540 /**
1541 * go_quad_acos:
1542 * @res: (out): result location
1543 * @a: quad-precision value
1544 *
1545 * This function computes the arc cosine of @a, storing the result in @res.
1546 **/
1547 /**
1548 * go_quad_acosl:
1549 * @res: (out): result location
1550 * @a: quad-precision value
1551 *
1552 * This function computes the arc cosine of @a, storing the result in @res.
1553 **/
1554 void
SUFFIX(go_quad_acos)1555 SUFFIX(go_quad_acos) (QUAD *res, const QUAD *a)
1556 {
1557 QUAD aa, aam1;
1558
1559 aa.h = SUFFIX(fabs) (a->h);
1560 aa.l = SUFFIX(fabs) (a->l);
1561 SUFFIX(go_quad_sub) (&aam1, &aa, &SUFFIX(go_quad_one));
1562 if (aam1.h > 0) {
1563 SUFFIX(go_quad_init) (res, SUFFIX(go_nan));
1564 return;
1565 }
1566
1567 SUFFIX(go_quad_agm_internal) (res, AGM_ARCCOS, &aa);
1568 if (a->h < 0)
1569 SUFFIX(go_quad_sub) (res, &SUFFIX(go_quad_pi), res);
1570 }
1571