Rich Felker | 007c03f | 2012-07-02 20:25:28 -0400 | [diff] [blame] | 1 | #define _GNU_SOURCE |
nsz | 0144b45 | 2012-03-15 08:17:28 +0100 | [diff] [blame] | 2 | #include "libm.h" |
| 3 | |
| 4 | #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 |
| 5 | void sincosl(long double x, long double *sin, long double *cos) |
| 6 | { |
Szabolcs Nagy | 73c870e | 2014-04-11 17:57:30 +0200 | [diff] [blame] | 7 | double sind, cosd; |
| 8 | sincos(x, &sind, &cosd); |
| 9 | *sin = sind; |
| 10 | *cos = cosd; |
nsz | 0144b45 | 2012-03-15 08:17:28 +0100 | [diff] [blame] | 11 | } |
| 12 | #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 |
nsz | 0144b45 | 2012-03-15 08:17:28 +0100 | [diff] [blame] | 13 | void sincosl(long double x, long double *sin, long double *cos) |
| 14 | { |
Szabolcs Nagy | ea9bb95 | 2013-09-03 18:50:58 +0000 | [diff] [blame] | 15 | union ldshape u = {x}; |
Szabolcs Nagy | bfda379 | 2013-05-18 14:40:22 +0000 | [diff] [blame] | 16 | unsigned n; |
nsz | 0144b45 | 2012-03-15 08:17:28 +0100 | [diff] [blame] | 17 | long double y[2], s, c; |
| 18 | |
Szabolcs Nagy | ea9bb95 | 2013-09-03 18:50:58 +0000 | [diff] [blame] | 19 | u.i.se &= 0x7fff; |
| 20 | if (u.i.se == 0x7fff) { |
nsz | 0144b45 | 2012-03-15 08:17:28 +0100 | [diff] [blame] | 21 | *sin = *cos = x - x; |
| 22 | return; |
| 23 | } |
Szabolcs Nagy | ea9bb95 | 2013-09-03 18:50:58 +0000 | [diff] [blame] | 24 | if (u.f < M_PI_4) { |
| 25 | if (u.i.se < 0x3fff - LDBL_MANT_DIG) { |
Szabolcs Nagy | bfda379 | 2013-05-18 14:40:22 +0000 | [diff] [blame] | 26 | /* raise underflow if subnormal */ |
Szabolcs Nagy | ea9bb95 | 2013-09-03 18:50:58 +0000 | [diff] [blame] | 27 | if (u.i.se == 0) FORCE_EVAL(x*0x1p-120f); |
Szabolcs Nagy | bfda379 | 2013-05-18 14:40:22 +0000 | [diff] [blame] | 28 | *sin = x; |
| 29 | /* raise inexact if x!=0 */ |
| 30 | *cos = 1.0 + x; |
| 31 | return; |
| 32 | } |
nsz | 0144b45 | 2012-03-15 08:17:28 +0100 | [diff] [blame] | 33 | *sin = __sinl(x, 0, 0); |
| 34 | *cos = __cosl(x, 0); |
| 35 | return; |
| 36 | } |
nsz | 0144b45 | 2012-03-15 08:17:28 +0100 | [diff] [blame] | 37 | n = __rem_pio2l(x, y); |
| 38 | s = __sinl(y[0], y[1], 1); |
| 39 | c = __cosl(y[0], y[1]); |
| 40 | switch (n & 3) { |
| 41 | case 0: |
| 42 | *sin = s; |
| 43 | *cos = c; |
| 44 | break; |
| 45 | case 1: |
| 46 | *sin = c; |
| 47 | *cos = -s; |
| 48 | break; |
| 49 | case 2: |
| 50 | *sin = -s; |
| 51 | *cos = -c; |
| 52 | break; |
| 53 | case 3: |
| 54 | default: |
| 55 | *sin = -c; |
| 56 | *cos = s; |
| 57 | break; |
| 58 | } |
| 59 | } |
| 60 | #endif |