RFA: Add ! LDBL_EQ_DBL implementations of cosl, sinl and tanl
Gregory Pietsch
gpietsch@comcast.net
Tue Mar 3 16:24:00 GMT 2015
Here's my attempt at all of the functions in the background of this.
Feel free to criticize and carve it up like a Thanksgiving turkey. --
Gregory
On 3/3/2015 4:48 AM, Corinna Vinschen wrote:
> On Mar 2 13:18, Craig Howland wrote:
>> On 03/02/2015 12:52 PM, Nick Clifton wrote:
>>> Hi Jeff, Hi Corinna,
>>>
>>> Attached is a patch that implements cosl(), sinl() and tanl() when
>>> long double is not the same as double. In the course of implementing
>>> these functions I found that I also needed to implement the fabls()
>>> and floorl() functions as well as several internal library support
>>> routines.
>>>
>>> Tested in the usual way with an x86_64-linux-gnu toolchain with
>>> direct linking to the math functions.
>>>
>>> OK to apply ?
>>>
>>> Cheers
>>> Nick
>>>
>>> newlib/ChangeLog
>>> 2015-03-02 Nick Clifton <...>
>>>
>>> * libm/common/cosl.c (cosl): Add implementation for when long
>>> double is not the same as double.
>>> * libm/common/fabls.c (fabls): Likewise.
>>> * libm/common/floorl.c (floorl): Likewise.
>>> * libm/common/sinl.c (sinl): Likewise.
>>> * libm/common/tanl.c (tanl): Likewise.
>>> * libm/common/fdlibmh.h: Add prototypes for __ieee754_rem_pio2l,
>>> __kernel_cosl, __kernel_cosl, __kernel_sinl and __kernel_tanl.
>>> * libm/math/e_rem_pio2l.c: New file. Implements
>>> __ieee754_rem_pio2l.
>>> * libm/math/k_cosl.c: New file. Implements __kernel_cosl.
>>> * libm/math/k_sinl.c: New file. Implements __kernel_sinl.
>>> * libm/math/k_tanl.c: New file. Implements __kernel_tanl.
>>> * libm/math/Makefile.am (src): Add k_cosl.c, k_sinl.c, k_tanl.c
>>> and e_rem_pio2l.c.
>>> * libm/math/Makefile.in: Regenerate.
> First of all, thanks for this patch, Nick.
>
>> math.h needs to be changed, too, to take the new routines out of the long
>> double==double checks.
> That's right, a math.h patch is missing.
>
>> These routines are not fully C99/POSIX-compliant in terms of floating point
>> exceptions and error returns if IEC 60559 is being used. (See, for example,
>> http://pubs.opengroup.org/onlinepubs/9699919799/functions/cosl.html, under
>> the "MX" markings, for the error-return considerations.)
> Neither are all our non-long double functions. Have a look at
> libm/math/s_cos.c, for instance.
>
>> Do we want to start out with this consideration missing? I would expect
>> that the vast majority of people probably do not care about the floating
>> point exception aspect. Setting errno, however, is quite a different
>> question, for which I have no idea how to gauge interest.
> Setting errno should be easy to add:
>
> - cosl, sinl and tanl just need an extra check for Inf to set errno to
> EDOM.
>
> - tanl needs an extra check if the result of a non-Inf arg is Inf to
> return ERANGE.
>
> - fabsl doesn't define errors>
>
> - floorl doesn't produce errors (see Linux man page),
>
> Nick, any problem to add this?
>
> However, that also means our long double functions will behave
> differently than our non-long double functions. Is anybody willing to
> making a sweep through our non-long double functions to add errno
> setting where missing?
>
>
> Thanks,
> Corinna
>
-------------- next part --------------
/* fabs, fabsf, fabsl - absolute value
AUTHOR: Gregory Pietsch
SYNOPSIS:
#include <math.h>
double fabs(double x);
float fabsf(float x);
long double fabsl(long double x);
DESCRIPTION:
The functionality described on this reference page is aligned with the ISO
C standard. Any conflict between the requirements described here and the
ISO C standard is unintentional. This volume of POSIX.1-2008 defers to the
ISO C standard.
These functions shall compute the absolute value of their argument x, |x|.
RETURN VALUE:
Upon successful completion, these functions shall return the absolute
value of x.
If x is NaN, a NaN shall be returned.
If x is +/-0, +0 shall be returned.
If x is +/-Inf, +Inf shall be returned.
ERRORS:
No errors are defined.
*/
#include "xmath.h"
double (fabs) (double x)
{
if (fpclassify (x) != FP_NAN)
_Set_sign (&x, 0, sizeof (double));
return x;
}
float (fabsf) (float x)
{
if (fpclassify (x) != FP_NAN)
_Set_sign (&x, 0, sizeof (float));
return x;
}
long double (fabsl) (long double x)
{
if (fpclassify (x) != FP_NAN)
_Set_sign (&x, 0, sizeof (long double));
return x;
}
/* END OF FILE */
-------------- next part --------------
/* floor, floorf, floorl - floor function
AUTHOR: Gregory Pietsch
SYNOPSIS:
#include <math.h>
double floor(double x);
float floorf(float x);
long double floorl(long double x);
DESCRIPTION:
The functionality described on this reference page is aligned with the ISO
C standard. Any conflict between the requirements described here and the
ISO C standard is unintentional. This volume of POSIX.1-2008 defers to the
ISO C standard.
These functions shall compute the largest integral value not greater than x.
RETURN VALUE:
The result shall have the same sign as x.
Upon successful completion, floor(), floorf(), and floorl() shall return the
largest integral value not greater than x, expressed as a type double, float,
or long double, respectively.
If x is NaN, a NaN shall be returned.
If x is +/-0 or +/-Inf, x shall be returned.
ERRORS:
No errors are defined.
*/
#include "xmath.h"
double (floor) (double x)
{
double y;
switch (fpclassify (x))
{
case FP_NAN:
case FP_INFINITE:
case FP_ZERO:
return x;
default:
y = trunc (x);
return (y > x) ? y - 1.0 : y;
}
}
float (floorf) (float x)
{
float y;
switch (fpclassify (x))
{
case FP_NAN:
case FP_INFINITE:
case FP_ZERO:
return x;
default:
y = truncf (x);
return (y > x) ? y - 1.0F : y;
}
}
long double (floorl) (long double x)
{
long double y;
switch (fpclassify (x))
{
case FP_NAN:
case FP_INFINITE:
case FP_ZERO:
return x;
default:
y = truncl (x);
return (y > x) ? y - 1.0L : y;
}
}
/* END OF FILE */
-------------- next part --------------
/* matherr.c - the _Matherr function
AUTHOR: Gregory Pietsch
DESCRIPTION:
This function handles error conditions.
*/
#include "xmath.h"
void _Matherr(int err, int fmask)
{
if (math_errhandling & MATH_ERRNO)
errno = err;
if (math_errhandling & MATH_ERREXCEPT)
feraiseexcept (fmask);
}
/* END OF FILE */
-------------- next part --------------
/* setsign - set the sign bit
AUTHOR: Gregory Pietsch
*/
#include "xmath.h"
void
_Set_sign (void *x, int sign, size_t size)
{
unsigned char *u = x;
if (sign)
u[_Byte (0, size)] |= 0x80;
else
u[_Byte (0, size)] &= ~0x80;
}
/* END OF FILE */
-------------- next part --------------
/* trunc, truncf, truncl - round to truncated integer value
AUTHOR: Gregory Pietsch
SYNOPSIS:
#include <math.h>
double trunc(double x);
float truncf(float x);
long double truncl(long double x);
DESCRIPTION:
The functionality described on this reference page is aligned with the ISO
C standard. Any conflict between the requirements described here and the
ISO C standard is unintentional. This volume of POSIX.1-2008 defers to the
ISO C standard.
These functions shall round their argument to the integer value, in
floating format, nearest to but no larger in magnitude than the argument.
RETURN VALUE:
Upon successful completion, these functions shall return the truncated
integer value.
The result shall have the same sign as x.
If x is NaN, a NaN shall be returned.
If x is +/-0 or +/-Inf, x shall be returned.
ERRORS:
No errors are defined.
*/
#include "xmath.h"
static void
_Trunc (void *x, size_t size)
{
unsigned char *px = x;
size_t ebits, i;
unsigned short e, sign;
int p;
/* first, get biased exponent */
ebits = _Get_ebits (size);
e = px[_Byte (0, size)];
e <<= 8;
e |= px[_Byte (1, size)];
sign = e & 0x8000;
e ^= sign;
e >>= (15 - ebits);
p = (e - ~(~0 << (ebits - 1)));
if (p > 0 && p < ((size - 2) << 3) + (15 - ebits))
{
/* some int bits, but not all */
p = (size << 3) - 1 - ebits - p;
for (i = size - 1; p > 8; --i, p -= 8)
px[_Byte (i, size)] = 0;
px[_Byte (i, size)] &= (~0 << p);
}
else if (p < 0)
{
/* all frac bits */
for (i = 0; i < size; ++i)
px[i] = 0;
if (sign)
px[_Byte (0, size)] = 0x80;
}
}
double (trunc) (double x)
{
switch (fpclassify (x))
{
case FP_NAN:
case FP_INFINITE:
case FP_ZERO:
return x;
case FP_SUBNORMAL:
return signbit (x) ? -0.0 : +0.0;
default:
_Trunc (&x, sizeof (double));
return x;
}
}
float (truncf) (float x)
{
switch (fpclassify (x))
{
case FP_NAN:
case FP_INFINITE:
case FP_ZERO:
return x;
case FP_SUBNORMAL:
return signbit (x) ? -0.0F : +0.0F;
default:
_Trunc (&x, sizeof (float));
return x;
}
}
long double (truncl) (long double x)
{
switch (fpclassify (x))
{
case FP_NAN:
case FP_INFINITE:
case FP_ZERO:
return x;
case FP_SUBNORMAL:
return signbit (x) ? -0.0L : +0.0L;
default:
_Trunc (&x, sizeof (long double));
return x;
}
}
/* END OF FILE */
-------------- next part --------------
/* xmath.h - internal header for the math functions
AUTHOR: Gregory Pietsch
*/
#ifndef _XMATH
#define _XMATH
/* include files */
#include <errno.h>
#include <fenv.h>
//#include <float.h>
#include <limits.h>
#include <math.h>
#include <stddef.h>
/* macros */
#ifdef _BIG_ENDIAN
#define _Byte(_X,_S) (_X)
#define _B4(_0,_1,_2,_3) _0, _1, _2, _3
#define _B8(_0,_1,_2,_3,_4,_5,_6,_7) \
_0, _1, _2, _3, _4, _5, _6, _7
#define _B10(_0,_1,_2,_3,_4,_5,_6,_7,_8,_9)
_0, _1, _2, _3, _4, _5, _6, _7, _8, _9
#define _B16(_0,_1,_2,_3,_4,_5,_6,_7,\
_8,_9,_A,_B,_C,_D,_E,_F)\
_0, _1, _2, _3, _4, _5, _6, _7,\
_8, _9, _A, _B, _C, _D, _E, _F
#else
#define _Byte(_X,_S) ((_S)-1-(_X))
#define _B4(_0,_1,_2,_3) _3, _2, _1, _0
#define _B8(_0,_1,_2,_3,_4,_5,_6,_7) \
_7, _6, _5, _4, _3, _2, _1, _0
#define _B10(_0,_1,_2,_3,_4,_5,_6,_7,_8,_9) \
_9, _8, _7, _6, _5, _4, _3, _2, _1, _0
#define _B16(_0,_1,_2,_3,_4,_5,_6,_7,\
_8,_9,_A,_B,_C,_D,_E,_F)\
_F, _E, _D, _C, _B, _A, _9, _8,\
_7, _6, _5, _4, _3, _2, _1, _0
#endif
#define _M_E (2.7182818284590452353602874713526624977572)
#define _M_LOG2E (1.4426950408889634073599246810018921374266)
#define _M_LOG10E (0.4342944819032518276511289189166050822943)
#define _M_LN2 (0.6931471805599453094172321214581765680755)
#define _M_LN10 (2.3025850929940456840179914546843642076011)
#define _M_PI (3.1415926535897932384626433832795028841968)
#define _M_PI_2 (1.5707963267948966192313216916397514420984)
#define _M_PI_4 (0.7853981633974483096156608458198757210492)
#define _M_1_PI (0.3183098861837906715377675267450287240689)
#define _M_2_PI (0.6366197723675813430755350534900574481379)
#define _M_2_SQRTPI (1.1283791670955125738961589031215451716881)
#define _M_SQRT2 (1.4142135623730950488016887242096980785696)
#define _M_SQRT1_2 (0.7071067811865475244008443621048490392848)
#define _M_EF (2.7182818284590452353602874713526624977572F)
#define _M_LOG2EF (1.4426950408889634073599246810018921374266F)
#define _M_LOG10EF (0.4342944819032518276511289189166050822943F)
#define _M_LN2F (0.6931471805599453094172321214581765680755F)
#define _M_LN10F (2.3025850929940456840179914546843642076011F)
#define _M_PIF (3.1415926535897932384626433832795028841968F)
#define _M_PI_2F (1.5707963267948966192313216916397514420984F)
#define _M_PI_4F (0.7853981633974483096156608458198757210492F)
#define _M_1_PIF (0.3183098861837906715377675267450287240689F)
#define _M_2_PIF (0.6366197723675813430755350534900574481379F)
#define _M_2_SQRTPIF (1.1283791670955125738961589031215451716881F)
#define _M_SQRT2F (1.4142135623730950488016887242096980785696F)
#define _M_SQRT1_2F (0.7071067811865475244008443621048490392848F)
#define _M_EL (2.7182818284590452353602874713526624977572L)
#define _M_LOG2EL (1.4426950408889634073599246810018921374266L)
#define _M_LOG10EL (0.4342944819032518276511289189166050822943L)
#define _M_LN2L (0.6931471805599453094172321214581765680755L)
#define _M_LN10L (2.3025850929940456840179914546843642076011L)
#define _M_PIL (3.1415926535897932384626433832795028841968L)
#define _M_PI_2L (1.5707963267948966192313216916397514420984L)
#define _M_PI_4L (0.7853981633974483096156608458198757210492L)
#define _M_1_PIL (0.3183098861837906715377675267450287240689L)
#define _M_2_PIL (0.6366197723675813430755350534900574481379L)
#define _M_2_SQRTPIL (1.1283791670955125738961589031215451716881L)
#define _M_SQRT2L (1.4142135623730950488016887242096980785696L)
#define _M_SQRT1_2L (0.7071067811865475244008443621048490392848L)
/* functions */
#ifdef __cplusplus
extern "C"
{
#endif
size_t _Get_ebits (size_t);
void _Matherr (int, int);
void _Set_sign (void *, int, size_t);
#ifdef __cplusplus
};
#endif
#endif
/* END OF FILE */
-------------- next part --------------
/* fpclassify.c - classify real-floating type
AUTHOR: Gregory Pietsch
DESCRIPTION:
The functionality described on this reference page is aligned with the ISO
C Standard. Any conflict between the requirements described here and the
ISO C Standard is unintentional. This volume of POSIX.1-2008 defers to the
ISO C Standard.
The fpclassify() macro shall classify its argument value as NaN, infinite,
normal, subnormal, zero, or into another implementation-defined category.
First, an argument represented in a format wider than its semantic type is
converted to its semantic type. Then classification is based on the type of
the argument.
RETURN VALUE:
The fpclassify() macro shall return the value of the number classification
macro appropriate to the value of its argument.
ERRORS:
No errors are defined.
*/
#include "xmath.h"
static int
_Classify (unsigned char *x, size_t size)
{
unsigned short hdr, mask, f;
size_t ebits, i;
ebits = _Get_ebits(size);
hdr = x[_Byte (0, size)];
hdr <<= 8;
hdr |= x[_Byte (1, size)];
mask = ~(~0 << ebits) << (15 - ebits);
f = (hdr & ~(~0 << (15 - ebits)));
hdr &= mask;
for (i = 2; i < size; ++i)
f |= x[_Byte (i, size)];
if (hdr == 0)
return f ? FP_SUBNORMAL : FP_ZERO;
else if (hdr == mask)
return f ? FP_NAN : FP_INFINITE;
else
return FP_NORMAL;
}
int
_Fpclassify (double x)
{
return _Classify ((unsigned char *) (&x), sizeof (double));
}
int
_Fpclassifyf (float x)
{
return _Classify ((unsigned char *) (&x), sizeof (float));
}
int
_Fpclassifyl (long double x)
{
return _Classify ((unsigned char *) (&x), sizeof (long double));
}
/* END OF FILE */
-------------- next part --------------
/* getebits.c - returns the number of bits in the exponent
AUTHOR: Gregory Pietsch
*/
#include "xmath.h"
size_t
_Get_ebits (size_t s)
{
switch (s)
{
case 4:
return 8;
case 8:
return 11;
case 10:
case 16:
return 15;
default:
return (size_t) (-1);
}
}
/* END OF FILE */
More information about the Newlib
mailing list