system/freebsd-src
system/freebsd-src
1
0
Fork 0
mirror of https://github.com/freebsd/freebsd-src synced 2024-07-21 18:27:22 +00:00
Code Issues Actions 2 Packages Projects Releases Wiki Activity

Fixes for bugs in sinpi/cospi/tanpi

Browse source 
patch to fix half-cycle trigonometric functions

Paul Zimmermann, a MPFR developer, contacted me about large errors in
the half-cycle trigonometric functions.  I've have investigated these
issues and developed the attached patch. The float, double, and ld80
(long double) changes have been tested.

Caveat emptor: The ld128 changes have not been compiled.  The ld128
changes have not been tested.  I do not have access to a system that
uses ld128 floating point.

Here is an itemized list of changes:

* lib/msun/src/math_private.h:
  . Add fast floor macros to compute the integer part of |x| for
    0 <= |x| 01xp(N-1), where N is the precision of the type of x.
    These macros are used in the half-cycle trigonometric functions
    (e.g., sinpi(x)).
  . The FFLOOR80 macros is used with the Intel 80-bit extended double
    functions.  This macors corrects an off-by-one error, which led to
    enormous error for |x| > 0x1p32.

* lib/msun/src/s_cospif.c:
* lib/msun/src/s_cospi.c:
* lib/msun/ld80/s_cospil.c:
  . Update Copyright years.
  . Use FFLOOR*() macro to get integer part of |x|.
  . Correct handle the range 0x1p(N-1) <= |x| < 0x1pN.  Here, one needs
    to determine if the integral value of |x| is even or odd to choose
    +1 or -1.  If |x| >= 0x1pN, always return +1.

* lib/msun/src/s_sinpif.c:
* lib/msun/src/s_sinpi.c:
* lib/msun/ld80/s_sinpil.c:
  . Update Copyright years.
  . Use FFLOOR*() macro to get integer part of |x|.

* lib/msun/src/s_tanpif.c:
* lib/msun/src/s_tanpi.c:
* lib/msun/ld80/s_tanpil.c:
  . Update Copyright years.
  . For +-0.5, return +-inf.  Previously, tanpi[fl]() returned an NaN.
  . Use FFLOOR*() to get integer part of |x|.  Need to determine if the
    integer part is even or odd.  This is used to set +-0 for |x|
integral
    and +-inf for (n+1/2).
  . For 0x1p(N-1) <= |x| < 0x1pN need to determine if x is an even or
odd
    integer to select +0 or -0.  For |x| >= 0x1pN, it is always an even
    integer, select 0.
  . Note, tanpi[fl](x) is an odd function, so one needs to consider
    tanpi[fl](-|x|) = - tanpi[fl](|x|).

* lib/msun/ld128/s_cospil.c:
* lib/msun/ld128/s_sinpil.c:
* lib/msun/ld128/s_tanpil.c:
  . Update Copyright years.
  . These routines use an FFLOOR128 macros, which likely should be
    replaced by a bit twiddling algorithm.
  . The same considerations above are applied to 0x1p112 <= |x| <
0x1p113,
    and |x| >= 0x1p113 cases.
  . Note, even and odd determination used fmodl(x,2.), which is likely
    slow.

PR:	272742
MFC after:	1 week
This commit is contained in:
Steve Kargl 2023-08-01 01:34:48 +03:00 committed by Konstantin Belousov
parent c66a499e03
commit 2d3b0a687b
13 changed files with 183 additions and 193 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

50
lib/msun/ld128/s_cospil.c
View file

@ -1,5 +1,5 @@
/*-
* Copyright (c) 2017-2021 Steven G. Kargl
* Copyright (c) 2017-2023 Steven G. Kargl
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
@ -28,6 +28,7 @@
* See ../src/s_cospi.c for implementation details.
*/
#include "fpmath.h"
#include "math.h"
#include "math_private.h"
@ -46,8 +47,7 @@ volatile static const double vzero = 0;
long double
cospil(long double x)
{
long double ax, c, xf;
uint32_t ix;
long double ai, ar, ax, c;
ax = fabsl(x);
@ -72,41 +72,37 @@ cospil(long double x)
}
if (ax < 0x1p112) {
/* Split x = n + r with 0 <= r < 1. */
xf = (ax + 0x1p112L) - 0x1p112L; /* Integer part */
ax -= xf; /* Remainder */
if (ax < 0) {
ax += 1;
xf -= 1;
}
/* Split ax = ai + ar with 0 <= ar < 1. */
FFLOORL128(ax, ai, ar);
if (ax < 0.5) {
if (ax < 0.25)
c = ax == 0 ? 1 : __kernel_cospil(ax);
if (ar < 0.5) {
if (ar < 0.25)
c = ar == 0 ? 1 : __kernel_cospil(ar);
else
c = __kernel_sinpil(0.5 - ax);
c = __kernel_sinpil(0.5 - ar);
} else {
if (ax < 0.75) {
if (ax == 0.5)
if (ar < 0.75) {
if (ar == 0.5)
return (0);
c = -__kernel_sinpil(ax - 0.5);
c = -__kernel_sinpil(ar - 0.5);
} else
c = -__kernel_cospil(1 - ax);
c = -__kernel_cospil(1 - ar);
}
if (xf > 0x1p64)
xf -= 0x1p64;
if (xf > 0x1p32)
xf -= 0x1p32;
ix = (uint32_t)xf;
return (ix & 1 ? -c : c);
return (fmodl(ai, 2.L) == 0 ? c : -c);
}
if (isinf(x) || isnan(x))
return (vzero / vzero);
/*
* |x| >= 0x1p112 is always an even integer, so return 1.
* For |x| >= 0x1p113, it is always an even integer, so return 1.
*/
return (1);
if (ax >= 0x1p113)
return (1);
/*
* For 0x1p112 <= |x| < 0x1p113 need to determine if x is an even
* or odd integer to return 1 or -1.
*/
return (fmodl(ax, 2.L) == 0 ? 1 : -1);
}

38
lib/msun/ld128/s_sinpil.c
View file

@ -1,5 +1,5 @@
/*-
* Copyright (c) 2017-2021 Steven G. Kargl
* Copyright (c) 2017-2023 Steven G. Kargl
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
@ -28,6 +28,7 @@
* See ../src/s_sinpi.c for implementation details.
*/
#include "fpmath.h"
#include "math.h"
#include "math_private.h"
@ -46,8 +47,7 @@ volatile static const double vzero = 0;
long double
sinpil(long double x)
{
long double ax, hi, lo, s, xf, xhi, xlo;
uint32_t ix;
long double ai, ar, ax, hi, lo, s, xhi, xlo;
ax = fabsl(x);
@ -78,35 +78,25 @@ sinpil(long double x)
}
if (ax < 0x1p112) {
/* Split x = n + r with 0 <= r < 1. */
xf = (ax + 0x1p112L) - 0x1p112L; /* Integer part */
ax -= xf; /* Remainder */
if (ax < 0) {
ax += 1;
xf -= 1;
}
/* Split ax = ai + ar with 0 <= ar < 1. */
FFLOORL128(ax, ai, ar);
if (ax == 0) {
if (ar == 0) {
s = 0;
} else {
if (ax < 0.5) {
if (ax <= 0.25)
s = __kernel_sinpil(ax);
if (ar < 0.5) {
if (ar <= 0.25)
s = __kernel_sinpil(ar);
else
s = __kernel_cospil(0.5 - ax);
s = __kernel_cospil(0.5 - ar);
} else {
if (ax < 0.75)
s = __kernel_cospil(ax - 0.5);
if (ar < 0.75)
s = __kernel_cospil(ar - 0.5);
else
s = __kernel_sinpil(1 - ax);
s = __kernel_sinpil(1 - ar);
}
if (xf > 0x1p64)
xf -= 0x1p64;
if (xf > 0x1p32)
xf -= 0x1p32;
ix = (uint32_t)xf;
if (ix & 1) s = -s;
s = fmodl(ai, 2.L) == 0 ? s : -s;
}
return (x < 0 ? -s : s);
}

39
lib/msun/ld128/s_tanpil.c
View file

@ -1,5 +1,5 @@
/*-
* Copyright (c) 2017-2021 Steven G. Kargl
* Copyright (c) 2017-2023 Steven G. Kargl
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
@ -28,6 +28,7 @@
* See ../src/s_tanpi.c for implementation details.
*/
#include "fpmath.h"
#include "math.h"
#include "math_private.h"
@ -69,8 +70,8 @@ volatile static const double vzero = 0;
long double
tanpil(long double x)
{
long double ax, hi, lo, xf, t;
uint32_t ix;
long double ai, ar, ax, hi, lo, t;
double odd;
ax = fabsl(x);
@ -88,27 +89,22 @@ tanpil(long double x)
}
t = __kernel_tanpil(ax);
} else if (ax == 0.5)
return ((ax - ax) / (ax - ax));
t = 1 / vzero;
else
t = -__kernel_tanpil(1 - ax);
return (x < 0 ? -t : t);
}
if (ix < 0x1p112) {
/* Split x = n + r with 0 <= r < 1. */
xf = (ax + 0x1p112L) - 0x1p112L; /* Integer part */
ax -= xf; /* Remainder */
if (ax < 0) {
ax += 1;
xf -= 1;
}
if (ax < 0.5)
t = ax == 0 ? 0 : __kernel_tanpil(ax);
else if (ax == 0.5)
return ((ax - ax) / (ax - ax));
if (ax < 0x1p112) {
/* Split ax = ai + ar with 0 <= ar < 1. */
FFLOORL128(ax, ai, ar);
odd = fmodl(ai, 2.L) == 0 ? 1 : -1;
if (ar < 0.5)
t = ar == 0 ? copysign(0., odd) : __kernel_tanpil(ar);
else if (ar == 0.5)
t = odd / vzero;
else
t = -__kernel_tanpil(1 - ax);
t = -__kernel_tanpil(1 - ar);
return (x < 0 ? -t : t);
}
@ -117,7 +113,10 @@ tanpil(long double x)
return (vzero / vzero);
/*
* |x| >= 0x1p112 is always an integer, so return +-0.
* For 0x1p112 <= |x| < 0x1p113 need to determine if x is an even
* or odd integer to set t = +0 or -0.
* For |x| >= 0x1p113, it is always an even integer, so t = 0.
*/
return (copysignl(0, x));
t = fmodl(ax,2.L) == 0 ? 0 : copysign(0., -1.);
return (copysignl(t, x));
}

22
lib/msun/ld80/s_cospil.c
View file

@ -1,5 +1,5 @@
/*-
* Copyright (c) 2017 Steven G. Kargl
* Copyright (c) 2017, 2023 Steven G. Kargl
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
@ -80,18 +80,8 @@ cospil(long double x)
RETURNI(c);
}
if (ix < 0x403e) { /* 1 <= |x| < 0x1p63 */
/* Determine integer part of ax. */
j0 = ix - 0x3fff + 1;
if (j0 < 32) {
lx = (lx >> 32) << 32;
lx &= ~(((lx << 32)-1) >> j0);
} else {
m = (uint64_t)-1 >> (j0 + 1);
if (lx & m) lx &= ~m;
}
INSERT_LDBL80_WORDS(x, ix, lx);
if (ix < 0x403e) { /* 1 <= |x| < 0x1p63 */
FFLOORL80(x, j0, ix, lx); /* Integer part of ax. */
ax -= x;
EXTRACT_LDBL80_WORDS(ix, lx, ax);
@ -123,7 +113,9 @@ cospil(long double x)
RETURNI(vzero / vzero);
/*
* |x| >= 0x1p63 is always an even integer, so return 1.
* For 0x1p63 <= |x| < 0x1p64 need to determine if x is an even
* or odd integer to return t = +1 or -1.
* For |x| >= 0x1p64, it is always an even integer, so t = 1.
*/
RETURNI(1);
RETURNI(ix >= 0x403f ? 1 : ((lx & 1) ? -1 : 1));
}

16
lib/msun/ld80/s_sinpil.c
View file

@ -1,5 +1,5 @@
/*-
* Copyright (c) 2017 Steven G. Kargl
* Copyright (c) 2017, 2023 Steven G. Kargl
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
@ -88,18 +88,8 @@ sinpil(long double x)
RETURNI((hx & 0x8000) ? -s : s);
}
if (ix < 0x403e) { /* 1 <= |x| < 0x1p63 */
/* Determine integer part of ax. */
j0 = ix - 0x3fff + 1;
if (j0 < 32) {
lx = (lx >> 32) << 32;
lx &= ~(((lx << 32)-1) >> j0);
} else {
m = (uint64_t)-1 >> (j0 + 1);
if (lx & m) lx &= ~m;
}
INSERT_LDBL80_WORDS(x, ix, lx);
if (ix < 0x403e) { /* 1 <= |x| < 0x1p63 */
FFLOORL80(x, j0, ix, lx); /* Integer part of ax. */
ax -= x;
EXTRACT_LDBL80_WORDS(ix, lx, ax);

34
lib/msun/ld80/s_tanpil.c
View file

@ -1,5 +1,5 @@
/*-
* Copyright (c) 2017 Steven G. Kargl
* Copyright (c) 2017, 2023 Steven G. Kargl
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
@ -72,7 +72,7 @@ volatile static const double vzero = 0;
long double
tanpil(long double x)
{
long double ax, hi, lo, t;
long double ax, hi, lo, odd, t;
uint64_t lx, m;
uint32_t j0;
uint16_t hx, ix;
@ -98,31 +98,22 @@ tanpil(long double x)
}
t = __kernel_tanpil(ax);
} else if (ax == 0.5)
RETURNI((ax - ax) / (ax - ax));
t = 1 / vzero;
else
t = -__kernel_tanpil(1 - ax);
RETURNI((hx & 0x8000) ? -t : t);
}
if (ix < 0x403e) { /* 1 <= |x| < 0x1p63 */
/* Determine integer part of ax. */
j0 = ix - 0x3fff + 1;
if (j0 < 32) {
lx = (lx >> 32) << 32;
lx &= ~(((lx << 32)-1) >> j0);
} else {
m = (uint64_t)-1 >> (j0 + 1);
if (lx & m) lx &= ~m;
}
INSERT_LDBL80_WORDS(x, ix, lx);
if (ix < 0x403e) { /* 1 <= |x| < 0x1p63 */
FFLOORL80(x, j0, ix, lx); /* Integer part of ax. */
odd = (uint64_t)x & 1 ? -1 : 1;
ax -= x;
EXTRACT_LDBL80_WORDS(ix, lx, ax);
if (ix < 0x3ffe) /* |x| < 0.5 */
t = ax == 0 ? 0 : __kernel_tanpil(ax);
else if (ax == 0.5)
RETURNI((ax - ax) / (ax - ax));
t = ix == 0 ? copysignl(0, odd) : __kernel_tanpil(ax);
else if (ax == 0.5L)
t = odd / vzero;
else
t = -__kernel_tanpil(1 - ax);
RETURNI((hx & 0x8000) ? -t : t);
@ -133,7 +124,10 @@ tanpil(long double x)
RETURNI(vzero / vzero);
/*
* |x| >= 0x1p63 is always an integer, so return +-0.
* For 0x1p63 <= |x| < 0x1p64 need to determine if x is an even
* or odd integer to set t = +0 or -0.
* For |x| >= 0x1p64, it is always an even integer, so t = 0.
*/
RETURNI(copysignl(0, x));
t = ix >= 0x403f ? 0 : (copysignl(0, (lx & 1) ? -1 : 1));
RETURNI((hx & 0x8000) ? -t : t);
}

53
lib/msun/src/math_private.h
View file

@ -688,6 +688,59 @@ irintl(long double x)
}
#endif
/*
* The following are fast floor macros for 0 <= |x| < 0x1p(N-1), where
* N is the precision of the type of x. These macros are used in the
* half-cycle trignometric functions (e.g., sinpi(x)).
*/
#define FFLOORF(x, j0, ix) do { \
(j0) = (((ix) >> 23) & 0xff) - 0x7f; \
(ix) &= ~(0x007fffff >> (j0)); \
SET_FLOAT_WORD((x), (ix)); \
} while (0)
#define FFLOOR(x, j0, ix, lx) do { \
(j0) = (((ix) >> 20) & 0x7ff) - 0x3ff; \
if ((j0) < 20) { \
(ix) &= ~(0x000fffff >> (j0)); \
(lx) = 0; \
} else { \
(lx) &= ~((uint32_t)0xffffffff >> ((j0) - 20)); \
} \
INSERT_WORDS((x), (ix), (lx)); \
} while (0)
#define FFLOORL80(x, j0, ix, lx) do { \
j0 = ix - 0x3fff + 1; \
if ((j0) < 32) { \
(lx) = ((lx) >> 32) << 32; \
(lx) &= ~((((lx) << 32)-1) >> (j0)); \
} else { \
uint64_t _m; \
_m = (uint64_t)-1 >> (j0); \
if ((lx) & _m) (lx) &= ~_m; \
} \
INSERT_LDBL80_WORDS((x), (ix), (lx)); \
} while (0)
#define FFLOORL128(x, ai, ar) do { \
union IEEEl2bits u; \
uint64_t m; \
int e; \
u.e = (x); \
e = u.bits.exp - 16383; \
if (e < 48) { \
m = ((1llu << 49) - 1) >> (e + 1); \
u.bits.manh &= ~m; \
u.bits.manl = 0; \
} else { \
m = (uint64_t)-1 >> (e - 48); \
u.bits.manl &= ~m; \
} \
(ai) = u.e; \
(ar) = (x) - (ai); \
} while (0)
#ifdef DEBUG
#if defined(__amd64__) || defined(__i386__)
#define breakpoint() asm("int $3")

20
lib/msun/src/s_cospi.c
View file

@ -1,5 +1,5 @@
/*-
* Copyright (c) 2017 Steven G. Kargl
* Copyright (c) 2017, 2023 Steven G. Kargl
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
@ -104,20 +104,10 @@ cospi(double x)
}
if (ix < 0x43300000) { /* 1 <= |x| < 0x1p52 */
/* Determine integer part of ax. */
j0 = ((ix >> 20) & 0x7ff) - 0x3ff;
if (j0 < 20) {
ix &= ~(0x000fffff >> j0);
lx = 0;
} else {
lx &= ~((uint32_t)0xffffffff >> (j0 - 20));
}
INSERT_WORDS(x, ix, lx);
FFLOOR(x, j0, ix, lx); /* Integer part of ax. */
ax -= x;
EXTRACT_WORDS(ix, lx, ax);
if (ix < 0x3fe00000) { /* |x| < 0.5 */
if (ix < 0x3fd00000) /* |x| < 0.25 */
c = ix == 0 ? 1 : __kernel_cospi(ax);
@ -143,9 +133,11 @@ cospi(double x)
return (vzero / vzero);
/*
* |x| >= 0x1p52 is always an even integer, so return 1.
* For 0x1p52 <= |x| < 0x1p53 need to determine if x is an even
* or odd integer to return +1 or -1.
* For |x| >= 0x1p53, it is always an even integer, so return 1.
*/
return (1);
return (ix < 0x43400000 ? ((lx & 1) ? -1 : 1) : 1);
}
#if LDBL_MANT_DIG == 53

16
lib/msun/src/s_cospif.c
View file

@ -1,5 +1,5 @@
/*-
* Copyright (c) 2017 Steven G. Kargl
* Copyright (c) 2017,2023 Steven G. Kargl
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
@ -71,12 +71,8 @@ cospif(float x)
return (c);
}
if (ix < 0x4b000000) { /* 1 <= |x| < 0x1p23 */
/* Determine integer part of ax. */
j0 = ((ix >> 23) & 0xff) - 0x7f;
ix &= ~(0x007fffff >> j0);
SET_FLOAT_WORD(x, ix);
if (ix < 0x4b000000) { /* 1 <= |x| < 0x1p23 */
FFLOORF(x, j0, ix); /* Integer part of ax. */
ax -= x;
GET_FLOAT_WORD(ix, ax);
@ -103,7 +99,9 @@ cospif(float x)
return (vzero / vzero);
/*
* |x| >= 0x1p23 is always an even integer, so return 1.
* For 0x1p23 <= |x| < 0x1p24 need to determine if x is an even
* or odd integer to return +1 or -1.
* For |x| >= 0x1p24, it is always an even integer, so return 1.
*/
return (1);
return (ix < 0x4b800000 ? ((ix & 1) ? -1 : 1) : 1);
}

13
lib/msun/src/s_sinpi.c
View file

@ -1,5 +1,5 @@
/*-
* Copyright (c) 2017 Steven G. Kargl
* Copyright (c) 2017, 2023 Steven G. Kargl
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
@ -118,16 +118,7 @@ sinpi(double x)
}
if (ix < 0x43300000) { /* 1 <= |x| < 0x1p52 */
/* Determine integer part of ax. */
j0 = ((ix >> 20) & 0x7ff) - 0x3ff;
if (j0 < 20) {
ix &= ~(0x000fffff >> j0);
lx = 0;
} else {
lx &= ~((uint32_t)0xffffffff >> (j0 - 20));
}
INSERT_WORDS(x, ix, lx);
FFLOOR(x, j0, ix, lx); /* Integer part of ax. */
ax -= x;
EXTRACT_WORDS(ix, lx, ax);

10
lib/msun/src/s_sinpif.c
View file

@ -1,5 +1,5 @@
/*-
* Copyright (c) 2017 Steven G. Kargl
* Copyright (c) 2017,2023 Steven G. Kargl
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
@ -81,12 +81,8 @@ sinpif(float x)
return ((hx & 0x80000000) ? -s : s);
}
if (ix < 0x4b000000) { /* 1 <= |x| < 0x1p23 */
/* Determine integer part of ax. */
j0 = ((ix >> 23) & 0xff) - 0x7f;
ix &= ~(0x007fffff >> j0);
SET_FLOAT_WORD(x, ix);
if (ix < 0x4b000000) { /* 1 <= |x| < 0x1p23 */
FFLOORF(x, j0, ix); /* Integer part of ax. */
ax -= x;
GET_FLOAT_WORD(ix, ax);

41
lib/msun/src/s_tanpi.c
View file

@ -1,5 +1,5 @@
/*-
* Copyright (c) 2017 Steven G. Kargl
* Copyright (c) 2017, 2023 Steven G. Kargl
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
@ -56,11 +56,15 @@
* 5. Special cases:
*
* tanpi(+-0) = +-0
* tanpi(+-n) = +-0, for positive integers n.
* tanpi(n) = +0 for positive even and negative odd integer n.
* tanpi(n) = -0 for positive odd and negative even integer n.
* tanpi(+-n+1/4) = +-1, for positive integers n.
* tanpi(+-n+1/2) = NaN, for positive integers n.
* tanpi(+-inf) = NaN. Raises the "invalid" floating-point exception.
* tanpi(nan) = NaN. Raises the "invalid" floating-point exception.
* tanpi(n+1/2) = +inf and raises the FE_DIVBYZERO exception for
* even integers n.
* tanpi(n+1/2) = -inf and raises the FE_DIVBYZERO exception for
* odd integers n.
* tanpi(+-inf) = NaN and raises the FE_INVALID exception.
* tanpi(nan) = NaN and raises the FE_INVALID exception.
*/
#include <float.h>
@ -106,7 +110,7 @@ volatile static const double vzero = 0;
double
tanpi(double x)
{
double ax, hi, lo, t;
double ax, hi, lo, odd, t;
uint32_t hx, ix, j0, lx;
EXTRACT_WORDS(hx, lx, x);
@ -132,30 +136,22 @@ tanpi(double x)
}
t = __kernel_tanpi(ax);
} else if (ax == 0.5)
return ((ax - ax) / (ax - ax));
t = 1 / vzero;
else
t = - __kernel_tanpi(1 - ax);
return ((hx & 0x80000000) ? -t : t);
}
if (ix < 0x43300000) { /* 1 <= |x| < 0x1p52 */
/* Determine integer part of ax. */
j0 = ((ix >> 20) & 0x7ff) - 0x3ff;
if (j0 < 20) {
ix &= ~(0x000fffff >> j0);
lx = 0;
} else {
lx &= ~(((uint32_t)(0xffffffff)) >> (j0 - 20));
}
INSERT_WORDS(x,ix,lx);
FFLOOR(x, j0, ix, lx); /* Integer part of ax. */
odd = (uint64_t)x & 1 ? -1 : 1;
ax -= x;
EXTRACT_WORDS(ix, lx, ax);
if (ix < 0x3fe00000) /* |x| < 0.5 */
t = ax == 0 ? 0 : __kernel_tanpi(ax);
t = ix == 0 ? copysign(0, odd) : __kernel_tanpi(ax);
else if (ax == 0.5)
return ((ax - ax) / (ax - ax));
t = odd / vzero;
else
t = - __kernel_tanpi(1 - ax);
@ -167,9 +163,12 @@ tanpi(double x)
return (vzero / vzero);
/*
* |x| >= 0x1p52 is always an integer, so return +-0.
* For 0x1p52 <= |x| < 0x1p53 need to determine if x is an even
* or odd integer to set t = +0 or -0.
* For |x| >= 0x1p54, it is always an even integer, so t = 0.
*/
return (copysign(0, x));
t = ix >= 0x43400000 ? 0 : (copysign(0, (lx & 1) ? -1 : 1));
return ((hx & 0x80000000) ? -t : t);
}
#if LDBL_MANT_DIG == 53

24
lib/msun/src/s_tanpif.c
View file

@ -1,5 +1,5 @@
/*-
* Copyright (c) 2017 Steven G. Kargl
* Copyright (c) 2017,2023 Steven G. Kargl
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
@ -58,7 +58,7 @@ volatile static const float vzero = 0;
float
tanpif(float x)
{
float ax, hi, lo, t;
float ax, hi, lo, odd, t;
uint32_t hx, ix, j0;
GET_FLOAT_WORD(hx, x);
@ -79,25 +79,22 @@ tanpif(float x)
}
t = __kernel_tanpif(ax);
} else if (ix == 0x3f000000)
return ((ax - ax) / (ax - ax));
t = 1 / vzero;
else
t = - __kernel_tanpif(1 - ax);
return ((hx & 0x80000000) ? -t : t);
}
if (ix < 0x4b000000) { /* 1 <= |x| < 0x1p23 */
/* Determine integer part of ax. */
j0 = ((ix >> 23) & 0xff) - 0x7f;
ix &= ~(0x007fffff >> j0);
SET_FLOAT_WORD(x, ix);
FFLOORF(x, j0, ix); /* Integer part of ax. */
odd = (uint32_t)x & 1 ? -1 : 1;
ax -= x;
GET_FLOAT_WORD(ix, ax);
if (ix < 0x3f000000) /* |x| < 0.5 */
t = ix == 0 ? 0 : __kernel_tanpif(ax);
t = ix == 0 ? copysignf(0, odd) : __kernel_tanpif(ax);
else if (ix == 0x3f000000)
return ((ax - ax) / (ax - ax));
t = odd / vzero;
else
t = - __kernel_tanpif(1 - ax);
return ((hx & 0x80000000) ? -t : t);
@ -108,7 +105,10 @@ tanpif(float x)
return (vzero / vzero);
/*
* |x| >= 0x1p23 is always an integer, so return +-0.
* For 0x1p23 <= |x| < 0x1p24 need to determine if x is an even
* or odd integer to set t = +0 or -0.
* For |x| >= 0x1p24, it is always an even integer, so t = 0.
*/
return (copysignf(0, x));
t = ix >= 0x4b800000 ? 0 : (copysignf(0, (ix & 1) ? -1 : 1));
return ((hx & 0x80000000) ? -t : t);
}