Lines 1-5
Link Here
|
1 |
/*- |
1 |
/*- |
2 |
* Copyright (c) 2017 Steven G. Kargl |
2 |
* Copyright (c) 2017, 2023 Steven G. Kargl |
3 |
* All rights reserved. |
3 |
* All rights reserved. |
4 |
* |
4 |
* |
5 |
* Redistribution and use in source and binary forms, with or without |
5 |
* Redistribution and use in source and binary forms, with or without |
Lines 56-66
Link Here
|
56 |
* 5. Special cases: |
56 |
* 5. Special cases: |
57 |
* |
57 |
* |
58 |
* tanpi(+-0) = +-0 |
58 |
* tanpi(+-0) = +-0 |
59 |
* tanpi(+-n) = +-0, for positive integers n. |
59 |
* tanpi(n) = +0 for positive even and negative odd integer n. |
|
|
60 |
* tanpi(n) = -0 for positive odd and negative even integer n. |
60 |
* tanpi(+-n+1/4) = +-1, for positive integers n. |
61 |
* tanpi(+-n+1/4) = +-1, for positive integers n. |
61 |
* tanpi(+-n+1/2) = NaN, for positive integers n. |
62 |
* tanpi(n+1/2) = +inf and raises the FE_DIVBYZERO exception for |
62 |
* tanpi(+-inf) = NaN. Raises the "invalid" floating-point exception. |
63 |
* even integers n. |
63 |
* tanpi(nan) = NaN. Raises the "invalid" floating-point exception. |
64 |
* tanpi(n+1/2) = -inf and raises the FE_DIVBYZERO exception for |
|
|
65 |
* odd integers n. |
66 |
* tanpi(+-inf) = NaN and raises the FE_INVALID exception. |
67 |
* tanpi(nan) = NaN and raises the FE_INVALID exception. |
64 |
*/ |
68 |
*/ |
65 |
|
69 |
|
66 |
#include <float.h> |
70 |
#include <float.h> |
Lines 106-112
volatile static const double vzero = 0;
Link Here
|
106 |
double |
110 |
double |
107 |
tanpi(double x) |
111 |
tanpi(double x) |
108 |
{ |
112 |
{ |
109 |
double ax, hi, lo, t; |
113 |
double ax, hi, lo, odd, t; |
110 |
uint32_t hx, ix, j0, lx; |
114 |
uint32_t hx, ix, j0, lx; |
111 |
|
115 |
|
112 |
EXTRACT_WORDS(hx, lx, x); |
116 |
EXTRACT_WORDS(hx, lx, x); |
Lines 132-161
tanpi(double x)
Link Here
|
132 |
} |
136 |
} |
133 |
t = __kernel_tanpi(ax); |
137 |
t = __kernel_tanpi(ax); |
134 |
} else if (ax == 0.5) |
138 |
} else if (ax == 0.5) |
135 |
return ((ax - ax) / (ax - ax)); |
139 |
t = 1 / vzero; |
136 |
else |
140 |
else |
137 |
t = - __kernel_tanpi(1 - ax); |
141 |
t = - __kernel_tanpi(1 - ax); |
138 |
return ((hx & 0x80000000) ? -t : t); |
142 |
return ((hx & 0x80000000) ? -t : t); |
139 |
} |
143 |
} |
140 |
|
144 |
|
141 |
if (ix < 0x43300000) { /* 1 <= |x| < 0x1p52 */ |
145 |
if (ix < 0x43300000) { /* 1 <= |x| < 0x1p52 */ |
142 |
/* Determine integer part of ax. */ |
146 |
FFLOOR(x, j0, ix, lx); /* Integer part of ax. */ |
143 |
j0 = ((ix >> 20) & 0x7ff) - 0x3ff; |
147 |
odd = (uint64_t)x & 1 ? -1 : 1; |
144 |
if (j0 < 20) { |
|
|
145 |
ix &= ~(0x000fffff >> j0); |
146 |
lx = 0; |
147 |
} else { |
148 |
lx &= ~(((uint32_t)(0xffffffff)) >> (j0 - 20)); |
149 |
} |
150 |
INSERT_WORDS(x,ix,lx); |
151 |
|
152 |
ax -= x; |
148 |
ax -= x; |
153 |
EXTRACT_WORDS(ix, lx, ax); |
149 |
EXTRACT_WORDS(ix, lx, ax); |
154 |
|
150 |
|
155 |
if (ix < 0x3fe00000) /* |x| < 0.5 */ |
151 |
if (ix < 0x3fe00000) /* |x| < 0.5 */ |
156 |
t = ax == 0 ? 0 : __kernel_tanpi(ax); |
152 |
t = ix == 0 ? copysign(0, odd) : __kernel_tanpi(ax); |
157 |
else if (ax == 0.5) |
153 |
else if (ax == 0.5) |
158 |
return ((ax - ax) / (ax - ax)); |
154 |
t = odd / vzero; |
159 |
else |
155 |
else |
160 |
t = - __kernel_tanpi(1 - ax); |
156 |
t = - __kernel_tanpi(1 - ax); |
161 |
|
157 |
|
Lines 167-175
tanpi(double x)
Link Here
|
167 |
return (vzero / vzero); |
163 |
return (vzero / vzero); |
168 |
|
164 |
|
169 |
/* |
165 |
/* |
170 |
* |x| >= 0x1p52 is always an integer, so return +-0. |
166 |
* For 0x1p52 <= |x| < 0x1p53 need to determine if x is an even |
|
|
167 |
* or odd integer to set t = +0 or -0. |
168 |
* For |x| >= 0x1p54, it is always an even integer, so t = 0. |
171 |
*/ |
169 |
*/ |
172 |
return (copysign(0, x)); |
170 |
t = ix >= 0x43400000 ? 0 : (copysign(0, (lx & 1) ? -1 : 1)); |
|
|
171 |
return ((hx & 0x80000000) ? -t : t); |
173 |
} |
172 |
} |
174 |
|
173 |
|
175 |
#if LDBL_MANT_DIG == 53 |
174 |
#if LDBL_MANT_DIG == 53 |