d/freebsd-dev
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

sinpi[fl] etc: Fix the ld128 implementations

Browse Source 
Mark Murray graciously provided access to an aarch64 system
to test the ld128 implementations.  This patch address
* Misuses of copysignl() in sinpil() and tanpil().
* Redo the splitting of argument 'x' into an integer part and
  remainder.  The remainder must satify 0 <= r < 1.
* Update the reduction of the integer part to something that can
  easily be seen as even or odd, e.g., sin(pi*x) = (-1)^n*sin(pi*r)
  with n <= 2^112 and we an reduce n by subtracting integer powers
  of 2.
* In s_cospil.c, fix typos where 'x' is used where 'ax', the
  remainder, is required.
* In tanpil(), fix the use of an uninitialized variable, ax = fabsl(ax),
  ax should be x in fabsl().

One item of note, in the limited tested on aarch64, the max ULP
for sinpil() and cospil() were less than 1.1 ULP, which is higher
that the desired max ULP less than 1.  This was traced to the
kernel for cosl() in the fundamental interval [0,pi/4].
The coefficients in the minmax polynomial likely need refinement.

PR:	218514
MFC after:	1 week
This commit is contained in:
Steve Kargl 2021-11-01 00:26:20 +02:00 committed by Konstantin Belousov
parent d5d2ce1c85
commit 4f889260c3
3 changed files with 46 additions and 37 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

31
lib/msun/ld128/s_cospil.c
View File

@ -1,5 +1,5 @@
/*-
* Copyright (c) 2017 Steven G. Kargl
* Copyright (c) 2017-2021 Steven G. Kargl
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
@ -26,9 +26,6 @@
/*
* See ../src/s_cospi.c for implementation details.
*
* FIXME: This has not been compiled nor has it been tested for accuracy.
* FIXME: This should use bit twiddling.
*/
#include "math.h"
@ -54,7 +51,7 @@ cospil(long double x)
ax = fabsl(x);
if (ax < 1) {
if (ax <= 1) {
if (ax < 0.25) {
if (ax < 0x1p-60) {
if ((int)x == 0)
@ -75,15 +72,21 @@ cospil(long double x)
}
if (ax < 0x1p112) {
xf = floorl(ax);
ax -= xf;
if (x < 0.5) {
if (x < 0.25)
/* 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) {
if (ax < 0.25)
c = ax == 0 ? 1 : __kernel_cospil(ax);
else
c = __kernel_sinpil(0.5 - ax);
} else {
if (x < 0.75) {
if (ax < 0.75) {
if (ax == 0.5)
return (0);
c = -__kernel_sinpil(ax - 0.5);
@ -91,10 +94,10 @@ cospil(long double x)
c = -__kernel_cospil(1 - ax);
}
if (xf > 0x1p50)
xf -= 0x1p50;
if (xf > 0x1p30)
xf -= 0x1p30;
if (xf > 0x1p64)
xf -= 0x1p64;
if (xf > 0x1p32)
xf -= 0x1p32;
ix = (uint32_t)xf;
return (ix & 1 ? -c : c);
}

29
lib/msun/ld128/s_sinpil.c
View File

@ -1,5 +1,5 @@
/*-
* Copyright (c) 2017 Steven G. Kargl
* Copyright (c) 2017-2021 Steven G. Kargl
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
@ -26,9 +26,6 @@
/*
* See ../src/s_sinpi.c for implementation details.
*
* FIXME: This has not been compiled nor has it been tested for accuracy.
* FIXME: This should use bit twiddling.
*/
#include "math.h"
@ -68,7 +65,7 @@ sinpil(long double x)
}
s = __kernel_sinpil(ax);
return (copysignl(s, x));
return (x < 0 ? -s : s);
}
if (ax < 0.5)
@ -77,12 +74,18 @@ sinpil(long double x)
s = __kernel_cospil(ax - 0.5);
else
s = __kernel_sinpil(1 - ax);
return (copysignl(s, x));
return (x < 0 ? -s : s);
}
if (ax < 0x1p112) {
xf = floorl(ax);
ax -= xf;
/* 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) {
s = 0;
} else {
@ -98,14 +101,14 @@ sinpil(long double x)
s = __kernel_sinpil(1 - ax);
}
if (xf > 0x1p50)
xf -= 0x1p50;
if (xf > 0x1p30)
xf -= 0x1p30;
if (xf > 0x1p64)
xf -= 0x1p64;
if (xf > 0x1p32)
xf -= 0x1p32;
ix = (uint32_t)xf;
if (ix & 1) s = -s;
}
return (copysignl(s, x));
return (x < 0 ? -s : s);
}
if (isinf(x) || isnan(x))

23
lib/msun/ld128/s_tanpil.c
View File

@ -1,5 +1,5 @@
/*-
* Copyright (c) 2017 Steven G. Kargl
* Copyright (c) 2017-2021 Steven G. Kargl
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
@ -26,9 +26,6 @@
/*
* See ../src/s_tanpi.c for implementation details.
*
* FIXME: This has not been compiled nor has it been tested for accuracy.
* FIXME: This should use bit twiddling.
*/
#include "math.h"
@ -75,7 +72,7 @@ tanpil(long double x)
long double ax, hi, lo, xf, t;
uint32_t ix;
ax = fabsl(ax);
ax = fabsl(x);
if (ax < 1) {
if (ax < 0.5) {
@ -94,19 +91,25 @@ tanpil(long double x)
return ((ax - ax) / (ax - ax));
else
t = -__kernel_tanpil(1 - ax);
return (copysignl(t, x));
return (x < 0 ? -t : t);
}
if (ix < 0x1p112) {
xf = floorl(ax);
ax -= xf;
/* 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));
else
t = -__kernel_tanpil(1 - ax);
return (copysignl(t, x));
return (x < 0 ? -t : t);
}
/* x = +-inf or nan. */
@ -114,7 +117,7 @@ tanpil(long double x)
return (vzero / vzero);
/*
* |x| >= 0x1p53 is always an integer, so return +-0.
* |x| >= 0x1p112 is always an integer, so return +-0.
*/
return (copysignl(0, x));
}