diff options
Diffstat (limited to 'libm')
| -rw-r--r-- | libm/Makefile.in | 17 | ||||
| -rw-r--r-- | libm/fp_private.h | 90 | ||||
| -rw-r--r-- | libm/fpmacros.c | 303 | ||||
| -rw-r--r-- | libm/powerpc/classic/s_ceil.c | 113 | ||||
| -rw-r--r-- | libm/powerpc/classic/s_copysign.c | 59 | ||||
| -rw-r--r-- | libm/powerpc/classic/s_floor.c | 113 | ||||
| -rw-r--r-- | libm/powerpc/classic/s_frexp.c | 65 | ||||
| -rw-r--r-- | libm/powerpc/classic/s_ldexp.c | 49 | ||||
| -rw-r--r-- | libm/powerpc/classic/s_logb.c | 107 | ||||
| -rw-r--r-- | libm/powerpc/classic/s_modf.c | 344 | ||||
| -rw-r--r-- | libm/powerpc/classic/s_nearbyint.c | 38 | ||||
| -rw-r--r-- | libm/powerpc/classic/s_rint.c | 159 | ||||
| -rw-r--r-- | libm/powerpc/classic/s_round.c | 115 | ||||
| -rw-r--r-- | libm/powerpc/classic/s_trunc.c | 89 | ||||
| -rw-r--r-- | libm/powerpc/classic/w_scalb.c | 94 | ||||
| -rw-r--r-- | libm/s_fpclassify.c | 44 | ||||
| -rw-r--r-- | libm/s_fpclassifyf.c | 44 | ||||
| -rw-r--r-- | libm/s_isinf.c | 34 | ||||
| -rw-r--r-- | libm/s_isinff.c | 30 | ||||
| -rw-r--r-- | libm/s_isnan.c | 45 | ||||
| -rw-r--r-- | libm/s_isnanf.c | 43 | ||||
| -rw-r--r-- | libm/s_signbit.c | 34 | ||||
| -rw-r--r-- | libm/s_signbitf.c | 34 |
23 files changed, 313 insertions, 1750 deletions
diff --git a/libm/Makefile.in b/libm/Makefile.in index f3768d38e..a94989dcf 100644 --- a/libm/Makefile.in +++ b/libm/Makefile.in @@ -34,15 +34,10 @@ libm_OUT:=$(top_builddir)libm # Fix builds for powerpc as there are different cores in this # section now.` -ifeq ($(TARGET_ARCH),powerpc) -ifeq ($(CONFIG_E500),y) +ifeq ($(TARGET_ARCH)-$(CONFIG_E500),powerpc-y) libm_ARCH_DIR:=$(libm_DIR)/$(TARGET_ARCH)/e500 libm_ARCH_OUT:=$(libm_OUT)/$(TARGET_ARCH)/e500 else -libm_ARCH_DIR:=$(libm_DIR)/$(TARGET_ARCH)/classic -libm_ARCH_OUT:=$(libm_OUT)/$(TARGET_ARCH)/classic -endif -else libm_ARCH_DIR:=$(libm_DIR)/$(TARGET_ARCH) libm_ARCH_OUT:=$(libm_OUT)/$(TARGET_ARCH) endif @@ -74,7 +69,9 @@ libm_CSRC := \ w_cabs.c w_cosh.c w_drem.c w_exp.c w_fmod.c w_gamma.c w_gamma_r.c \ w_hypot.c w_j0.c w_j1.c w_jn.c w_lgamma.c w_lgamma_r.c \ w_log.c w_log10.c w_pow.c w_remainder.c w_scalb.c w_sinh.c \ - w_sqrt.c fpmacros.c nan.c carg.c s_llrint.c + w_sqrt.c nan.c carg.c s_llrint.c \ + s_fpclassify.c s_fpclassifyf.c s_signbit.c s_signbitf.c \ + s_isnan.c s_isnanf.c s_isnf.c s_isinff.c FL_MOBJ := \ acosf.o acoshf.o asinf.o asinhf.o atan2f.o atanf.o atanhf.o cbrtf.o \ ceilf.o copysignf.o cosf.o coshf.o erfcf.o erff.o exp2f.o expf.o \ @@ -103,13 +100,9 @@ endif ifeq ($(UCLIBC_HAS_FPU),y) ifeq ($(DO_C99_MATH),y) ifneq ($(strip $(libm_ARCH_OBJS)),) -ifeq ($(TARGET_ARCH),powerpc) -ifeq ($(CONFIG_E500),y) +ifeq ($(TARGET_ARCH)-$(CONFIG_E500),powerpc-y) CFLAGS-libm/$(TARGET_ARCH)/e500/ := $(CFLAGS-libm) else -CFLAGS-libm/$(TARGET_ARCH)/classic/ := $(CFLAGS-libm) -endif -else CFLAGS-libm/$(TARGET_ARCH)/ := $(CFLAGS-libm) endif diff --git a/libm/fp_private.h b/libm/fp_private.h deleted file mode 100644 index 0ddb616c4..000000000 --- a/libm/fp_private.h +++ /dev/null @@ -1,90 +0,0 @@ -/******************************************************************************* -* * -* File fp_private.h, * -* All pack 4 dependencies for the MathLib elems plus some defines used * -* throughout MathLib. * -* * -* Copyright © 1991 Apple Computer, Inc. All rights reserved. * -* * -* Written by Ali Sazegari, started on October 1991, * -* * -* W A R N I N G: This routine expects a 64 bit double model. * -* * -*******************************************************************************/ - -#define NoException 0 - -/******************************************************************************* -* Values of constants. * -*******************************************************************************/ - -//#define SgnMask 0x8000 -#define dSgnMask 0x80000000 -#define sSgnMask 0x7FFFFFFF - -//#define ExpMask 0x7FFF -#define dExpMask 0x7FF00000 -#define sExpMask 0xFF000000 - - /* according to rounding BIG & SMALL are: */ -#define BIG 1.1e+300 /* used to deliver ±° or largest number, */ -#define SMALL 1.1e-300 /* used to deliver ±0 or smallest number. */ -#define InfExp 0x7FF -#define dMaxExp 0x7FF00000 - -#define MaxExpP1 1024 -#define MaxExp 1023 - -#define DenormLimit -52 - -//#define ManMask 0x80000000 -#define dManMask 0x00080000 - -//#define IsItDenorm 0x80000000 -#define dIsItDenorm 0x00080000 - -//#define xIsItSNaN 0x40000000 -#define dIsItSNaN 0x00080000 - -#define dHighMan 0x000FFFFF -#define dFirstBitSet 0x00080000 -#define BIAS 0x3FF - -//#define GetSign 0x8000 -#define dGetSign 0x80000000 -#define sGetSign 0x80000000 - -//#define Infinity(x) ( x.hex.exponent & ExpMask ) == ExpMask -#define dInfinity(x) ( x.hex.high & dExpMask ) == dExpMask -#define sInfinity(x) ( ( x.hexsgl << 1 ) & sExpMask ) == sExpMask - -//#define Exponent(x) x.hex.exponent & ExpMask -#define dExponent(x) x.hex.high & dExpMask -#define sExponent(x) ( ( x.hexsgl << 1 ) & sExpMask ) - -#define sZero(x) ( x.hexsgl & sSgnMask ) == 0 -//#define Sign(x) ( x.hex.exponent & SgnMask ) == SgnMask - -/******************************************************************************* -* Types used in the auxiliary functions. * -*******************************************************************************/ - -#include <stdint.h> -#include <endian.h> - -typedef struct /* Hex representation of a double. */ - { -#if (__BYTE_ORDER == __BIG_ENDIAN) - uint32_t high; - uint32_t low; -#else - uint32_t low; - uint32_t high; -#endif - } dHexParts; - -typedef union - { - unsigned char byties[8]; - double dbl; - } DblInHex; diff --git a/libm/fpmacros.c b/libm/fpmacros.c deleted file mode 100644 index 0a079c016..000000000 --- a/libm/fpmacros.c +++ /dev/null @@ -1,303 +0,0 @@ -/*********************************************************************** -** File: fpmacros.c -** -** Contains: C source code for implementations of floating-point -** functions which involve float format numbers, as -** defined in header <fp.h>. In particular, this file -** contains implementations of functions -** __fpclassify(d,f), __isnormal(d,f), __isfinite(d,f), -** __isnan(d,f), and __signbit(d,f). This file targets -** PowerPC platforms. -** -** Written by: Robert A. Murley, Ali Sazegari -** -** Copyright: c 2001 by Apple Computer, Inc., all rights reserved -** -** Change History (most recent first): -** -** 07 Jul 01 ram First created from fpfloatfunc.c, fp.c, -** classify.c and sign.c in MathLib v3 Mac OS9. -** -***********************************************************************/ - -#include <features.h> -#include <sys/types.h> -#include <math.h> -#include "fp_private.h" - -#define SIGN_MASK 0x80000000 -#define NSIGN_MASK 0x7fffffff -#define FEXP_MASK 0x7f800000 -#define FFRAC_MASK 0x007fffff - -/*********************************************************************** - int __fpclassifyf(float x) returns the classification code of the - argument x, as defined in <fp.h>. - - Exceptions: INVALID signaled if x is a signaling NaN; in this case, - the FP_QNAN code is returned. - - Calls: none -***********************************************************************/ - -libm_hidden_proto(__fpclassifyf) -int __fpclassifyf ( float x ) -{ - unsigned int iexp; - - union { - u_int32_t lval; - float fval; - } z; - - z.fval = x; - iexp = z.lval & FEXP_MASK; /* isolate float exponent */ - - if (iexp == FEXP_MASK) { /* NaN or INF case */ - if ((z.lval & 0x007fffff) == 0) - return FP_INFINITE; - return FP_NAN; - } - - if (iexp != 0) /* normal float */ - return FP_NORMAL; - - if (x == 0.0) - return FP_ZERO; /* zero */ - else - return FP_SUBNORMAL; /* must be subnormal */ -} -libm_hidden_def(__fpclassifyf) - - -/*********************************************************************** - Function __fpclassify, - Implementation of classify of a double number for the PowerPC. - - Exceptions: INVALID signaled if x is a signaling NaN; in this case, - the FP_QNAN code is returned. - - Calls: none -***********************************************************************/ - -libm_hidden_proto(__fpclassify) -int __fpclassify ( double arg ) -{ - register unsigned int exponent; - union - { - dHexParts hex; - double dbl; - } x; - - x.dbl = arg; - - exponent = x.hex.high & dExpMask; - if ( exponent == dExpMask ) - { - if ( ( ( x.hex.high & dHighMan ) | x.hex.low ) == 0 ) - return FP_INFINITE; - else - return FP_NAN; - } - else if ( exponent != 0) - return FP_NORMAL; - else { - if ( arg == 0.0 ) - return FP_ZERO; - else - return FP_SUBNORMAL; - } -} -libm_hidden_def(__fpclassify) - - -/*********************************************************************** - int __isnormalf(float x) returns nonzero if and only if x is a - normalized float number and zero otherwise. - - Exceptions: INVALID is raised if x is a signaling NaN; in this case, - zero is returned. - - Calls: none -***********************************************************************/ - -int __isnormalf ( float x ); -int __isnormalf ( float x ) -{ - unsigned int iexp; - union { - u_int32_t lval; - float fval; - } z; - - z.fval = x; - iexp = z.lval & FEXP_MASK; /* isolate float exponent */ - return ((iexp != FEXP_MASK) && (iexp != 0)); -} - - -int __isnormal ( double x ); -int __isnormal ( double x ) -{ - return ( __fpclassify ( x ) == FP_NORMAL ); -} - - -/*********************************************************************** - int __isfinitef(float x) returns nonzero if and only if x is a - finite (normal, subnormal, or zero) float number and zero otherwise. - - Exceptions: INVALID is raised if x is a signaling NaN; in this case, - zero is returned. - - Calls: none -***********************************************************************/ - -int __finitef ( float x ) -{ - union { - u_int32_t lval; - float fval; - } z; - - z.fval = x; - return ((z.lval & FEXP_MASK) != FEXP_MASK); -} -strong_alias(__finitef,finitef) - -#if 0 /* use __finite in s_finite.c */ -int __finite ( double x ) -{ - return ( __fpclassify ( x ) >= FP_ZERO ); -} -strong_alias(__finite,finite) -#endif - - -/*********************************************************************** - int __signbitf(float x) returns nonzero if and only if the sign - bit of x is set and zero otherwise. - - Exceptions: INVALID is raised if x is a signaling NaN. - - Calls: none -***********************************************************************/ - -libm_hidden_proto(__signbitf) -int __signbitf ( float x ) -{ - union { - u_int32_t lval; - float fval; - } z; - - z.fval = x; - return ((z.lval & SIGN_MASK) != 0); -} -libm_hidden_def(__signbitf) - - -/*********************************************************************** - Function sign of a double. - Implementation of sign bit for the PowerPC. - - Calls: none -***********************************************************************/ - -libm_hidden_proto(__signbit) -int __signbit ( double arg ) -{ - union - { - dHexParts hex; - double dbl; - } x; - int sign; - - x.dbl = arg; - sign = ( ( x.hex.high & dSgnMask ) == dSgnMask ) ? 1 : 0; - return sign; -} -libm_hidden_def(__signbit) - - -/*********************************************************************** -* int __isinff(float x) returns -1 if value represents negative -* infinity, 1 if value represents positive infinity, -* and 0 otherwise. -* -* Calls: __signbit -* +***********************************************************************/ -int __isinff ( float x ) -{ - int class = __fpclassifyf(x); - if ( class == FP_INFINITE ) { - return ( (__signbitf(x)) ? -1 : 1); - } - return 0; -} -strong_alias(__isinff,isinff) - -int __isinf ( double x ) -{ - int class = __fpclassify(x); - if ( class == FP_INFINITE ) { - return ( (__signbit(x)) ? -1 : 1); - } - return 0; -} -strong_alias(__isinf,isinf) - -#if 0 -int __isinfl ( long double x ) -{ - int class = __fpclassify(x); - if ( class == FP_INFINITE ) { - return ( (__signbit(x)) ? -1 : 1); - } - return 0; -} -strong_alias(__isinfl,isinfl) -#endif - -/*********************************************************************** - int __isnanf(float x) returns nonzero if and only if x is a - NaN and zero otherwise. - - Exceptions: INVALID is raised if x is a signaling NaN; in this case, - nonzero is returned. - - Calls: none -***********************************************************************/ - -int __isnanf ( float x ) -{ - union { - u_int32_t lval; - float fval; - } z; - - z.fval = x; - return (((z.lval&FEXP_MASK) == FEXP_MASK) && ((z.lval&FFRAC_MASK) != 0)); -} -strong_alias(__isnanf,isnanf) - -libm_hidden_proto(__isnan) -int __isnan ( double x ) -{ - int class = __fpclassify(x); - return ( class == FP_NAN ); -} -libm_hidden_def(__isnan) -strong_alias(__isnan,isnan) - -#if 0 -int __isnanl ( long double x ) -{ - int class = __fpclassify(x); - return ( class == FP_NAN ); -} -strong_alias(__isnanl,isnanl) -#endif - diff --git a/libm/powerpc/classic/s_ceil.c b/libm/powerpc/classic/s_ceil.c deleted file mode 100644 index ee4ceb5fc..000000000 --- a/libm/powerpc/classic/s_ceil.c +++ /dev/null @@ -1,113 +0,0 @@ -/******************************************************************************* -* * -* File ceilfloor.c, * -* Function ceil(x) and floor(x), * -* Implementation of ceil and floor for the PowerPC. * -* * -* Copyright © 1991 Apple Computer, Inc. All rights reserved. * -* * -* Written by Ali Sazegari, started on November 1991, * -* * -* based on math.h, library code for Macintoshes with a 68881/68882 * -* by Jim Thomas. * -* * -* W A R N I N G: This routine expects a 64 bit double model. * -* * -* December 03 1992: first rs6000 port. * -* July 14 1993: comment changes and addition of #pragma fenv_access. * -* May 06 1997: port of the ibm/taligent ceil and floor routines. * -* April 11 2001: first port to os x using gcc. * -* June 13 2001: replaced __setflm with in-line assembly * -* * -*******************************************************************************/ - -#include <math.h> -#include <endian.h> - -static const double twoTo52 = 4503599627370496.0; -static const unsigned long signMask = 0x80000000ul; - -typedef union - { - struct { -#if (__BYTE_ORDER == __BIG_ENDIAN) - unsigned long int hi; - unsigned long int lo; -#else - unsigned long int lo; - unsigned long int hi; -#endif - } words; - double dbl; - } DblInHex; - -/******************************************************************************* -* Functions needed for the computation. * -*******************************************************************************/ - -/******************************************************************************* -* Ceil(x) returns the smallest integer not less than x. * -*******************************************************************************/ - -libm_hidden_proto(ceil) -double ceil ( double x ) - { - DblInHex xInHex,OldEnvironment; - register double y; - register unsigned long int xhi; - register int target; - - xInHex.dbl = x; - xhi = xInHex.words.hi & 0x7fffffffUL; // xhi is the high half of |x| - target = ( xInHex.words.hi < signMask ); - - if ( xhi < 0x43300000ul ) -/******************************************************************************* -* Is |x| < 2.0^52? * -*******************************************************************************/ - { - if ( xhi < 0x3ff00000ul ) -/******************************************************************************* -* Is |x| < 1.0? * -*******************************************************************************/ - { - if ( ( xhi | xInHex.words.lo ) == 0ul ) // zero x is exact case - return ( x ); - else - { // inexact case - __asm__ ("mffs %0" : "=f" (OldEnvironment.dbl)); - OldEnvironment.words.lo |= 0x02000000ul; - __asm__ ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl )); - if ( target ) - return ( 1.0 ); - else - return ( -0.0 ); - } - } -/******************************************************************************* -* Is 1.0 < |x| < 2.0^52? * -*******************************************************************************/ - if ( target ) - { - y = ( x + twoTo52 ) - twoTo52; // round at binary pt. - if ( y < x ) - return ( y + 1.0 ); - else - return ( y ); - } - - else - { - y = ( x - twoTo52 ) + twoTo52; // round at binary pt. - if ( y < x ) - return ( y + 1.0 ); - else - return ( y ); - } - } -/******************************************************************************* -* |x| >= 2.0^52 or x is a NaN. * -*******************************************************************************/ - return ( x ); - } -libm_hidden_def(ceil) diff --git a/libm/powerpc/classic/s_copysign.c b/libm/powerpc/classic/s_copysign.c deleted file mode 100644 index c6f1307a3..000000000 --- a/libm/powerpc/classic/s_copysign.c +++ /dev/null @@ -1,59 +0,0 @@ -/******************************************************************************* -* * -* File sign.c, * -* Functions copysign and __signbitd. * -* For PowerPC based machines. * -* * -* Copyright © 1991, 2001 Apple Computer, Inc. All rights reserved. * -* * -* Written by Ali Sazegari, started on June 1991. * -* * -* August 26 1991: no CFront Version 1.1d17 warnings. * -* September 06 1991: passes the test suite with invalid raised on * -* signaling nans. sane rom code behaves the same. * -* September 24 1992: took the ̉#include support.hÓ out. * -* Dcember 02 1992: PowerPC port. * -* July 20 1994: __fabs added * -* July 21 1994: deleted unnecessary functions: neg, COPYSIGNnew, * -* and SIGNNUMnew. * -* April 11 2001: first port to os x using gcc. * -* removed fabs and deffered to gcc for direct * -* instruction generation. * -* * -*******************************************************************************/ - -#include <math.h> -#include "../../fp_private.h" - -/******************************************************************************* -* * -* Function copysign. * -* Implementation of copysign for the PowerPC. * -* * -******************************************************************************** -* Note: The order of the operands in this function is reversed from that * -* suggested in the IEEE standard 754. * -*******************************************************************************/ - -libm_hidden_proto(copysign) -double copysign ( double arg2, double arg1 ) - { - union - { - dHexParts hex; - double dbl; - } x, y; - -/******************************************************************************* -* No need to flush NaNs out. * -*******************************************************************************/ - - x.dbl = arg1; - y.dbl = arg2; - - y.hex.high = y.hex.high & 0x7FFFFFFF; - y.hex.high = ( y.hex.high | ( x.hex.high & dSgnMask ) ); - - return y.dbl; - } -libm_hidden_def(copysign) diff --git a/libm/powerpc/classic/s_floor.c b/libm/powerpc/classic/s_floor.c deleted file mode 100644 index 2cd720eff..000000000 --- a/libm/powerpc/classic/s_floor.c +++ /dev/null @@ -1,113 +0,0 @@ -/******************************************************************************* -* * -* File ceilfloor.c, * -* Function ceil(x) and floor(x), * -* Implementation of ceil and floor for the PowerPC. * -* * -* Copyright © 1991 Apple Computer, Inc. All rights reserved. * -* * -* Written by Ali Sazegari, started on November 1991, * -* * -* based on math.h, library code for Macintoshes with a 68881/68882 * -* by Jim Thomas. * -* * -* W A R N I N G: This routine expects a 64 bit double model. * -* * -* December 03 1992: first rs6000 port. * -* July 14 1993: comment changes and addition of #pragma fenv_access. * -* May 06 1997: port of the ibm/taligent ceil and floor routines. * -* April 11 2001: first port to os x using gcc. * -* June 13 2001: replaced __setflm with in-line assembly * -* * -*******************************************************************************/ - -#include <math.h> -#include <endian.h> - -static const double twoTo52 = 4503599627370496.0; -static const unsigned long signMask = 0x80000000ul; - -typedef union - { - struct { -#if (__BYTE_ORDER == __BIG_ENDIAN) - unsigned long int hi; - unsigned long int lo; -#else - unsigned long int lo; - unsigned long int hi; -#endif - } words; - double dbl; - } DblInHex; - -/******************************************************************************* -* Functions needed for the computation. * -*******************************************************************************/ - -/******************************************************************************* -* Floor(x) returns the largest integer not greater than x. * -*******************************************************************************/ - -libm_hidden_proto(floor) -double floor ( double x ) - { - DblInHex xInHex,OldEnvironment; - register double y; - register unsigned long int xhi; - register long int target; - - xInHex.dbl = x; - xhi = xInHex.words.hi & 0x7fffffffUL; // xhi is the high half of |x| - target = ( xInHex.words.hi < signMask ); - - if ( xhi < 0x43300000ul ) -/******************************************************************************* -* Is |x| < 2.0^52? * -*******************************************************************************/ - { - if ( xhi < 0x3ff00000ul ) -/******************************************************************************* -* Is |x| < 1.0? * -*******************************************************************************/ - { - if ( ( xhi | xInHex.words.lo ) == 0ul ) // zero x is exact case - return ( x ); - else - { // inexact case - __asm__ ("mffs %0" : "=f" (OldEnvironment.dbl)); - OldEnvironment.words.lo |= 0x02000000ul; - __asm__ ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl )); - if ( target ) - return ( 0.0 ); - else - return ( -1.0 ); - } - } -/******************************************************************************* -* Is 1.0 < |x| < 2.0^52? * -*******************************************************************************/ - if ( target ) - { - y = ( x + twoTo52 ) - twoTo52; // round at binary pt. - if ( y > x ) - return ( y - 1.0 ); - else - return ( y ); - } - - else - { - y = ( x - twoTo52 ) + twoTo52; // round at binary pt. - if ( y > x ) - return ( y - 1.0 ); - else - return ( y ); - } - } -/******************************************************************************* -* |x| >= 2.0^52 or x is a NaN. * -*******************************************************************************/ - return ( x ); - } -libm_hidden_def(floor) diff --git a/libm/powerpc/classic/s_frexp.c b/libm/powerpc/classic/s_frexp.c deleted file mode 100644 index 001aaf708..000000000 --- a/libm/powerpc/classic/s_frexp.c +++ /dev/null @@ -1,65 +0,0 @@ -/******************************************************************************* -* * -* File frexpldexp.c, * -* Functions frexp(x) and ldexp(x), * -* Implementation of frexp and ldexp functions for the PowerPC. * -* * -* Copyright © 1991 Apple Computer, Inc. All rights reserved. * -* * -* Written by Ali Sazegari, started on January 1991, * -* * -* W A R N I N G: This routine expects a 64 bit double model. * -* * -* December03 1992: first rs6000 implementation. * -* October 05 1993: added special cases for NaN and ° in frexp. * -* May 27 1997: improved the performance of frexp by eliminating the * -* switch statement. * -* June 13 2001: (ram) rewrote frexp to eliminate calls to scalb and * -* logb. * -* * -*******************************************************************************/ - -#include <limits.h> -#include <math.h> -#include <endian.h> - -static const double two54 = 1.80143985094819840000e+16; /* 0x43500000, 0x00000000 */ - -typedef union - { - struct { -#if (__BYTE_ORDER == __BIG_ENDIAN) - unsigned long int hi; - unsigned long int lo; -#else - unsigned long int lo; - unsigned long int hi; -#endif - } words; - double dbl; - } DblInHex; - -libm_hidden_proto(frexp) -double frexp ( double value, int *eptr ) - { - DblInHex argument; - unsigned long int valueHead; - - argument.dbl = value; - valueHead = argument.words.hi & 0x7fffffffUL; // valueHead <- |x| - - *eptr = 0; - if ( valueHead >= 0x7ff00000 || ( valueHead | argument.words.lo ) == 0 ) - return value; // 0, inf, or NaN - - if ( valueHead < 0x00100000 ) - { // denorm - argument.dbl = two54 * value; - valueHead = argument.words.hi &0x7fffffff; - *eptr = -54; - } - *eptr += ( valueHead >> 20 ) - 1022; - argument.words.hi = ( argument.words.hi & 0x800fffff ) | 0x3fe00000; - return argument.dbl; - } -libm_hidden_def(frexp) diff --git a/libm/powerpc/classic/s_ldexp.c b/libm/powerpc/classic/s_ldexp.c deleted file mode 100644 index 10100d7c2..000000000 --- a/libm/powerpc/classic/s_ldexp.c +++ /dev/null @@ -1,49 +0,0 @@ -/******************************************************************************* -* * -* File frexpldexp.c, * -* Functions frexp(x) and ldexp(x), * -* Implementation of frexp and ldexp functions for the PowerPC. * -* * -* Copyright © 1991 Apple Computer, Inc. All rights reserved. * -* * -* Written by Ali Sazegari, started on January 1991, * -* * -* W A R N I N G: This routine expects a 64 bit double model. * -* * -* December03 1992: first rs6000 implementation. * -* October 05 1993: added special cases for NaN and ° in frexp. * -* May 27 1997: improved the performance of frexp by eliminating the * -* switch statement. * -* June 13 2001: (ram) rewrote frexp to eliminate calls to scalb and * -* logb. * -* * -*******************************************************************************/ - -#include <limits.h> -#include <math.h> -#include <endian.h> - -typedef union - { - struct { -#if (__BYTE_ORDER == __BIG_ENDIAN) - unsigned long int hi; - unsigned long int lo; -#else - unsigned long int lo; - unsigned long int hi; -#endif - } words; - double dbl; - } DblInHex; - -libm_hidden_proto(ldexp) -double ldexp ( double value, int exp ) - { - if ( exp > SHRT_MAX ) - exp = SHRT_MAX; - else if ( exp < -SHRT_MAX ) - exp = -SHRT_MAX; - return scalb ( value, exp ); - } -libm_hidden_def(ldexp) diff --git a/libm/powerpc/classic/s_logb.c b/libm/powerpc/classic/s_logb.c deleted file mode 100644 index 81daa412e..000000000 --- a/libm/powerpc/classic/s_logb.c +++ /dev/null @@ -1,107 +0,0 @@ -/******************************************************************************* -* * -* File logb.c, * -* Functions logb. * -* Implementation of logb for the PowerPC. * -* * -* Copyright © 1991 Apple Computer, Inc. All rights reserved. * -* * -* Written by Ali Sazegari, started on June 1991, * -* * -* August 26 1991: removed CFront Version 1.1d17 warnings. * -* August 27 1991: no errors reported by the test suite. * -* November 11 1991: changed CLASSEXTENDED to the macro CLASSIFY and * -* + or - infinity to constants. * -* November 18 1991: changed the macro CLASSIFY to CLASSEXTENDEDint to * -* improve performance. * -* February 07 1992: changed bit operations to macros ( object size is * -* unchanged ). * -* September24 1992: took the "#include support.h" out. * -* December 03 1992: first rs/6000 port. * -* August 30 1992: set the divide by zero for the zero argument case. * -* October 05 1993: corrected the environment. * -* October 17 1994: replaced all environmental functions with __setflm. * -* May 28 1997: made speed improvements. * -* April 30 2001: forst mac os x port using gcc. * -* * -******************************************************************************** -* The C math library offers a similar function called "frexp". It is * -* different in details from logb, but similar in spirit. This current * -* implementation of logb follows the recommendation in IEEE Standard 854 * -* which is different in its handling of denormalized numbers from the IEEE * -* Standard 754. * -*******************************************************************************/ - -#include <math.h> -#include <endian.h> - -typedef union - { - struct { -#if (__BYTE_ORDER == __BIG_ENDIAN) - unsigned long int hi; - unsigned long int lo; -#else - unsigned long int lo; - unsigned long int hi; -#endif - } words; - double dbl; - } DblInHex; - -static const double twoTo52 = 4.50359962737049600e15; // 0x1p52 -static const double klTod = 4503601774854144.0; // 0x1.000008p52 -static const unsigned long int signMask = 0x80000000ul; -static const DblInHex minusInf = {{ 0xFFF00000, 0x00000000 }}; - - -/******************************************************************************* -******************************************************************************** -* L O G B * -******************************************************************************** -*******************************************************************************/ - -libm_hidden_proto(logb) -double logb ( double x ) - { - DblInHex xInHex; - long int shiftedExp; - - xInHex.dbl = x; - shiftedExp = ( xInHex.words.hi & 0x7ff00000UL ) >> 20; - - if ( shiftedExp == 2047 ) - { // NaN or INF - if ( ( ( xInHex.words.hi & signMask ) == 0 ) || ( x != x ) ) - return x; // NaN or +INF return x - else - return -x; // -INF returns +INF - } - - if ( shiftedExp != 0 ) // normal number - shiftedExp -= 1023; // unbias exponent - - else if ( x == 0.0 ) - { // zero - xInHex.words.hi = 0x0UL; // return -infinity - return ( minusInf.dbl ); - } - - else - { // subnormal number - xInHex.dbl *= twoTo52; // scale up - shiftedExp = ( xInHex.words.hi & 0x7ff00000UL ) >> 20; - shiftedExp -= 1075; // unbias exponent - } - - if ( shiftedExp == 0 ) // zero result - return ( 0.0 ); - - else - { // nonzero result - xInHex.dbl = klTod; - xInHex.words.lo += shiftedExp; - return ( xInHex.dbl - klTod ); - } - } -libm_hidden_def(logb) diff --git a/libm/powerpc/classic/s_modf.c b/libm/powerpc/classic/s_modf.c deleted file mode 100644 index c2221bc0b..000000000 --- a/libm/powerpc/classic/s_modf.c +++ /dev/null @@ -1,344 +0,0 @@ -/******************************************************************************* -** File: rndint.c -** -** Contains: C source code for implementations of floating-point -** functions which round to integral value or format, as -** defined in header <fp.h>. In particular, this file -** contains implementations of functions rinttol, roundtol, -** modf and modfl. This file targets PowrPC or Power platforms. -** -** Written by: A. Sazegari, Apple AltiVec Group -** Created originally by Jon Okada, Apple Numerics Group -** -** Copyright: © 1992-2001 by Apple Computer, Inc., all rights reserved -** -** Change History (most recent first): -** -** 13 Jul 01 ram replaced --setflm calls with inline assembly -** 03 Mar 01 ali first port to os x using gcc, added the crucial __setflm -** definition. -** 1. removed double_t, put in double for now. -** 2. removed iclass from nearbyint. -** 3. removed wrong comments intrunc. -** 4. -** 13 May 97 ali made performance improvements in rint, rinttol, roundtol -** and trunc by folding some of the taligent ideas into this -** implementation. nearbyint is faster than the one in taligent, -** rint is more elegant, but slower by %30 than the taligent one. -** 09 Apr 97 ali deleted modfl and deferred to AuxiliaryDD.c -** 15 Sep 94 ali Major overhaul and performance improvements of all functions. -** 20 Jul 94 PAF New faster version -** 16 Jul 93 ali Added the modfl function. -** 18 Feb 93 ali Changed the return value of fenv functions -** feclearexcept and feraiseexcept to their new -** NCEG X3J11.1/93-001 definitions. -** 16 Dec 92 JPO Removed __itrunc implementation to a -** separate file. -** 15 Dec 92 JPO Added __itrunc implementation and modified -** rinttol to include conversion from double -** to long int format. Modified roundtol to -** call __itrunc. -** 10 Dec 92 JPO Added modf (double) implementation. -** 04 Dec 92 JPO First created. -** -*******************************************************************************/ - -#include <limits.h> -#include <math.h> -#include <endian.h> - -#define SET_INVALID 0x01000000UL - -typedef union - { - struct { -#if (__BYTE_ORDER == __BIG_ENDIAN) - unsigned long int hi; - unsigned long int lo; -#else - unsigned long int lo; - unsigned long int hi; -#endif - } words; - double dbl; - } DblInHex; - -static const unsigned long int signMask = 0x80000000ul; -static const double twoTo52 = 4503599627370496.0; -static const double doubleToLong = 4503603922337792.0; // 2^52 -static const DblInHex TOWARDZERO = {{ 0x00000000, 0x00000001 }}; - - -/******************************************************************************* -* * -* The function rinttol converts its double argument to integral value * -* according to the current rounding direction and returns the result in * -* long int format. This conversion signals invalid if the argument is a * -* NaN or the rounded intermediate result is out of range of the * -* destination long int format, and it delivers an unspecified result in * -* this case. This function signals inexact if the rounded result is * -* within range of the long int format but unequal to the operand. * -* * -*******************************************************************************/ - -long int rinttol ( double x ) - { - register double y; - DblInHex argument, OldEnvironment; - unsigned long int xHead; - register long int target; - - argument.dbl = x; - target = ( argument.words.hi < signMask ); // flag positive sign - xHead = argument.words.hi & 0x7ffffffful; // high 32 bits of x - - if ( target ) -/******************************************************************************* -* Sign of x is positive. * -*******************************************************************************/ - { - if ( xHead < 0x41dffffful ) - { // x is safely in long range - y = ( x + twoTo52 ) - twoTo52; // round at binary point - argument.dbl = y + doubleToLong; // force result into argument.words.lo - return ( ( long ) argument.words.lo ); - } - - __asm__ ("mffs %0" : "=f" (OldEnvironment.dbl)); // get environment - - if ( xHead > 0x41dffffful ) - { // x is safely out of long range - OldEnvironment.words.lo |= SET_INVALID; - __asm__ ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl )); - return ( LONG_MAX ); - } - -/******************************************************************************* -* x > 0.0 and may or may not be out of range of long. * -*******************************************************************************/ - - y = ( x + twoTo52 ) - twoTo52; // do rounding - if ( y > ( double ) LONG_MAX ) - { // out of range of long - OldEnvironment.words.lo |= SET_INVALID; - __asm__ ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl )); - return ( LONG_MAX ); - } - argument.dbl = y + doubleToLong; // in range - return ( ( long ) argument.words.lo ); // return result & flags - } - -/******************************************************************************* -* Sign of x is negative. * -*******************************************************************************/ - if ( xHead < 0x41e00000ul ) - { // x is safely in long range - y = ( x - twoTo52 ) + twoTo52; - argument.dbl = y + doubleToLong; - return ( ( long ) argument.words.lo ); - } - - __asm__ ("mffs %0" : "=f" (OldEnvironment.dbl)); // get environment - - if ( xHead > 0x41e00000ul ) - { // x is safely out of long range - OldEnvironment.words.lo |= SET_INVALID; - __asm__ ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl )); - return ( LONG_MIN ); - } - -/******************************************************************************* -* x < 0.0 and may or may not be out of range of long. * -*******************************************************************************/ - - y = ( x - twoTo52 ) + twoTo52; // do rounding - if ( y < ( double ) LONG_MIN ) - { // out of range of long - OldEnvironment.words.lo |= SET_INVALID; - __asm__ ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl )); - return ( LONG_MIN ); - } - argument.dbl = y + doubleToLong; // in range - return ( ( long ) argument.words.lo ); // return result & flags - } - -/******************************************************************************* -* * -* The function roundtol converts its double argument to integral format * -* according to the "add half to the magnitude and chop" rounding mode of * -* Pascal's Round function and FORTRAN's NINT function. This conversion * -* signals invalid if the argument is a NaN or the rounded intermediate * -* result is out of range of the destination long int format, and it * -* delivers an unspecified result in this case. This function signals * -* inexact if the rounded result is within range of the long int format but * -* unequal to the operand. * -* * -*******************************************************************************/ - -long int roundtol ( double x ) - { - register double y, z; - DblInHex argument, OldEnvironment; - register unsigned long int xhi; - register long int target; - const DblInHex kTZ = {{ 0x0, 0x1 }}; - const DblInHex kUP = {{ 0x0, 0x2 }}; - - argument.dbl = x; - xhi = argument.words.hi & 0x7ffffffful; // high 32 bits of x - target = ( argument.words.hi < signMask ); // flag positive sign - - if ( xhi > 0x41e00000ul ) -/******************************************************************************* -* Is x is out of long range or NaN? * -*******************************************************************************/ - { - __asm__ ("mffs %0" : "=f" (OldEnvironment.dbl)); // get environment - OldEnvironment.words.lo |= SET_INVALID; - __asm__ ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl )); - if ( target ) // pin result - return ( LONG_MAX ); - else - return ( LONG_MIN ); - } - - if ( target ) -/******************************************************************************* -* Is sign of x is "+"? * -*******************************************************************************/ - { - if ( x < 2147483647.5 ) -/******************************************************************************* -* x is in the range of a long. * -*******************************************************************************/ - { - y = ( x + doubleToLong ) - doubleToLong; // round at binary point - if ( y != x ) - { // inexact case - __asm__ ("mffs %0" : "=f" (OldEnvironment.dbl)); // save environment - __asm__ ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( kTZ.dbl )); // truncate rounding - z = x + 0.5; // truncate x + 0.5 - argument.dbl = z + doubleToLong; - __asm__ ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl )); - return ( ( long ) argument.words.lo ); - } - - argument.dbl = y + doubleToLong; // force result into argument.words.lo - return ( ( long ) argument.words.lo ); // return long result - } -/******************************************************************************* -* Rounded positive x is out of the range of a long. * -*******************************************************************************/ - __asm__ ("mffs %0" : "=f" (OldEnvironment.dbl)); - OldEnvironment.words.lo |= SET_INVALID; - __asm__ ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl )); - return ( LONG_MAX ); // return pinned result - } -/******************************************************************************* -* x < 0.0 and may or may not be out of the range of a long. * -*******************************************************************************/ - if ( x > -2147483648.5 ) -/******************************************************************************* -* x is in the range of a long. * -*******************************************************************************/ - { - y = ( x + doubleToLong ) - doubleToLong; // round at binary point - if ( y != x ) - { // inexact case - __asm__ ("mffs %0" : "=f" (OldEnvironment.dbl)); // save environment - __asm__ ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( kUP.dbl )); // round up - z = x - 0.5; // truncate x - 0.5 - argument.dbl = z + doubleToLong; - __asm__ ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl )); - return ( ( long ) argument.words.lo ); - } - - argument.dbl = y + doubleToLong; - return ( ( long ) argument.words.lo ); // return long result - } -/******************************************************************************* -* Rounded negative x is out of the range of a long. * -*******************************************************************************/ - __asm__ ("mffs %0" : "=f" (OldEnvironment.dbl)); - OldEnvironment.words.lo |= SET_INVALID; - __asm__ ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl )); - return ( LONG_MIN ); // return pinned result - } - -/******************************************************************************* -* The modf family of functions separate a floating-point number into its * -* fractional and integral parts, returning the fractional part and writing * -* the integral part in floating-point format to the object pointed to by a * -* pointer argument. If the input argument is integral or infinite in * -* value, the return value is a zero with the sign of the input argument. * -* The modf family of functions raises no floating-point exceptions. older * -* implemenation set the INVALID flag due to signaling NaN input. * -* * -*******************************************************************************/ - -/******************************************************************************* -* modf is the double implementation. * -*******************************************************************************/ - -libm_hidden_proto(modf) -double modf ( double x, double *iptr ) - { - register double OldEnvironment, xtrunc; - register unsigned long int xHead, signBit; - DblInHex argument; - - argument.dbl = x; - xHead = argument.words.hi & 0x7ffffffful; // |x| high bit pattern - signBit = ( argument.words.hi & 0x80000000ul ); // isolate sign bit - if (xHead == 0x7ff81fe0) - signBit = signBit | 0; - - if ( xHead < 0x43300000ul ) -/******************************************************************************* -* Is |x| < 2.0^53? * -*******************************************************************************/ - { - if ( xHead < 0x3ff00000ul ) -/******************************************************************************* -* Is |x| < 1.0? * -*******************************************************************************/ - { - argument.words.hi = signBit; // truncate to zero - argument.words.lo = 0ul; - *iptr = argument.dbl; - return ( x ); - } -/******************************************************************************* -* Is 1.0 < |x| < 2.0^52? * -*******************************************************************************/ - __asm__ ("mffs %0" : "=f" (OldEnvironment)); // save environment - // round toward zero - __asm__ ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( TOWARDZERO.dbl )); - if ( signBit == 0ul ) // truncate to integer - xtrunc = ( x + twoTo52 ) - twoTo52; - else - xtrunc = ( x - twoTo52 ) + twoTo52; - // restore caller's env - __asm__ ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment )); - *iptr = xtrunc; // store integral part - if ( x != xtrunc ) // nonzero fraction - return ( x - xtrunc ); - else - { // zero with x's sign - argument.words.hi = signBit; - argument.words.lo = 0ul; - return ( argument.dbl ); - } - } - - *iptr = x; // x is integral or NaN - if ( x != x ) // NaN is returned - return x; - else - { // zero with x's sign - argument.words.hi = signBit; - argument.words.lo = 0ul; - return ( argument.dbl ); - } - } -libm_hidden_def(modf) diff --git a/libm/powerpc/classic/s_nearbyint.c b/libm/powerpc/classic/s_nearbyint.c deleted file mode 100644 index d08430dc6..000000000 --- a/libm/powerpc/classic/s_nearbyint.c +++ /dev/null @@ -1,38 +0,0 @@ -#include <limits.h> -#include <math.h> - -/******************************************************************************* -* * -* The function nearbyint rounds its double argument to integral value * -* according to the current rounding direction and returns the result in * -* double format. This function does not signal inexact. * -* * -******************************************************************************** -* * -* This function calls fabs and copysign. * -* * -*******************************************************************************/ - -static const double twoTo52 = 4503599627370496.0; - -libm_hidden_proto(nearbyint) -double nearbyint ( double x ) - { - double y; - double OldEnvironment; - - y = twoTo52; - - __asm__ ("mffs %0" : "=f" (OldEnvironment)); /* get the environement */ - - if ( fabs ( x ) >= y ) /* huge case is exact */ - return x; - if ( x < 0 ) y = -y; /* negative case */ - y = ( x + y ) - y; /* force rounding */ - if ( y == 0.0 ) /* zero results mirror sign of x */ - y = copysign ( y, x ); -// restore old flags - __asm__ ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment )); - return ( y ); - } -libm_hidden_def(nearbyint) diff --git a/libm/powerpc/classic/s_rint.c b/libm/powerpc/classic/s_rint.c deleted file mode 100644 index dd2ae585e..000000000 --- a/libm/powerpc/classic/s_rint.c +++ /dev/null @@ -1,159 +0,0 @@ -/******************************************************************************* -** File: rndint.c -** -** Contains: C source code for implementations of floating-point -** functions which round to integral value or format, as -** defined in header <fp.h>. In particular, this file -** contains implementations of functions rint, nearbyint, -** rinttol, round, roundtol, trunc, modf and modfl. This file -** targets PowerPC or Power platforms. -** -** Written by: A. Sazegari, Apple AltiVec Group -** Created originally by Jon Okada, Apple Numerics Group -** -** Copyright: © 1992-2001 by Apple Computer, Inc., all rights reserved -** -** Change History (most recent first): -** -** 13 Jul 01 ram replaced --setflm calls with inline assembly -** 03 Mar 01 ali first port to os x using gcc, added the crucial __setflm -** definition. -** 1. removed double_t, put in double for now. -** 2. removed iclass from nearbyint. -** 3. removed wrong comments intrunc. -** 4. -** 13 May 97 ali made performance improvements in rint, rinttol, roundtol -** and trunc by folding some of the taligent ideas into this -** implementation. nearbyint is faster than the one in taligent, -** rint is more elegant, but slower by %30 than the taligent one. -** 09 Apr 97 ali deleted modfl and deferred to AuxiliaryDD.c -** 15 Sep 94 ali Major overhaul and performance improvements of all functions. -** 20 Jul 94 PAF New faster version -** 16 Jul 93 ali Added the modfl function. -** 18 Feb 93 ali Changed the return value of fenv functions -** feclearexcept and feraiseexcept to their new -** NCEG X3J11.1/93-001 definitions. -** 16 Dec 92 JPO Removed __itrunc implementation to a -** separate file. -** 15 Dec 92 JPO Added __itrunc implementation and modified -** rinttol to include conversion from double -** to long int format. Modified roundtol to -** call __itrunc. -** 10 Dec 92 JPO Added modf (double) implementation. -** 04 Dec 92 JPO First created. -** -*******************************************************************************/ - -#include <limits.h> -#include <math.h> -#include <endian.h> - -#define SET_INVALID 0x01000000UL - -typedef union - { - struct { -#if (__BYTE_ORDER == __BIG_ENDIAN) - unsigned long int hi; - unsigned long int lo; -#else - unsigned long int lo; - unsigned long int hi; -#endif - } words; - double dbl; - } DblInHex; - -static const unsigned long int signMask = 0x80000000ul; -static const double twoTo52 = 4503599627370496.0; -static const double doubleToLong = 4503603922337792.0; // 2^52 -static const DblInHex Huge = {{ 0x7FF00000, 0x00000000 }}; -static const DblInHex TOWARDZERO = {{ 0x00000000, 0x00000001 }}; - -libm_hidden_proto(rint) -/******************************************************************************* -* * -* The function rint rounds its double argument to integral value * -* according to the current rounding direction and returns the result in * -* double format. This function signals inexact if an ordered return * -* value is not equal to the operand. * -* * -******************************************************************************** -* * -* This function calls: fabs. * -* * -*******************************************************************************/ - -/******************************************************************************* -* First, an elegant implementation. * -******************************************************************************** -* -*double rint ( double x ) -* { -* double y; -* -* y = twoTo52.fval; -* -* if ( fabs ( x ) >= y ) // huge case is exact -* return x; -* if ( x < 0 ) y = -y; // negative case -* y = ( x + y ) - y; // force rounding -* if ( y == 0.0 ) // zero results mirror sign of x -* y = copysign ( y, x ); -* return ( y ); -* } -******************************************************************************** -* Now a bit twidling version that is about %30 faster. * -*******************************************************************************/ - -double rint ( double x ) - { - DblInHex argument; - register double y; - unsigned long int xHead; - register long int target; - - argument.dbl = x; - xHead = argument.words.hi & 0x7fffffffUL; // xHead <- high half of |x| - target = ( argument.words.hi < signMask ); // flags positive sign - - if ( xHead < 0x43300000ul ) -/******************************************************************************* -* Is |x| < 2.0^52? * -*******************************************************************************/ - { - if ( xHead < 0x3ff00000ul ) -/******************************************************************************* -* Is |x| < 1.0? * -*******************************************************************************/ - { - if ( target ) - y = ( x + twoTo52 ) - twoTo52; // round at binary point - else - y = ( x - twoTo52 ) + twoTo52; // round at binary point - if ( y == 0.0 ) - { // fix sign of zero result - if ( target ) - return ( 0.0 ); - else - return ( -0.0 ); - } - return y; - } - -/******************************************************************************* -* Is 1.0 < |x| < 2.0^52? * -*******************************************************************************/ - - if ( target ) - return ( ( x + twoTo52 ) - twoTo52 ); // round at binary pt. - else - return ( ( x - twoTo52 ) + twoTo52 ); - } - -/******************************************************************************* -* |x| >= 2.0^52 or x is a NaN. * -*******************************************************************************/ - return ( x ); - } -libm_hidden_def(rint) diff --git a/libm/powerpc/classic/s_round.c b/libm/powerpc/classic/s_round.c deleted file mode 100644 index 9fd73801d..000000000 --- a/libm/powerpc/classic/s_round.c +++ /dev/null @@ -1,115 +0,0 @@ -#include <limits.h> -#include <math.h> -#include <endian.h> - -typedef union - { - struct { -#if (__BYTE_ORDER == __BIG_ENDIAN) - unsigned long int hi; - unsigned long int lo; -#else - unsigned long int lo; - unsigned long int hi; -#endif - } words; - double dbl; - } DblInHex; - -static const unsigned long int signMask = 0x80000000ul; -static const double twoTo52 = 4503599627370496.0; - -/******************************************************************************* -* * -* The function round rounds its double argument to integral value * -* according to the "add half to the magnitude and truncate" rounding of * -* Pascal's Round function and FORTRAN's ANINT function and returns the * -* result in double format. This function signals inexact if an ordered * -* return value is not equal to the operand. * -* * -*******************************************************************************/ - -libm_hidden_proto(round) -double round ( double x ) - { - DblInHex argument, OldEnvironment; - register double y, z; - register unsigned long int xHead; - register long int target; - - argument.dbl = x; - xHead = argument.words.hi & 0x7fffffffUL; // xHead <- high half of |x| - target = ( argument.words.hi < signMask ); // flag positive sign - - if ( xHead < 0x43300000ul ) -/******************************************************************************* -* Is |x| < 2.0^52? * -*******************************************************************************/ - { - if ( xHead < 0x3ff00000ul ) -/******************************************************************************* -* Is |x| < 1.0? * -*******************************************************************************/ - { - __asm__ ("mffs %0" : "=f" (OldEnvironment.dbl)); // get environment - if ( xHead < 0x3fe00000ul ) -/******************************************************************************* -* Is |x| < 0.5? * -*******************************************************************************/ - { - if ( ( xHead | argument.words.lo ) != 0ul ) - OldEnvironment.words.lo |= 0x02000000ul; - __asm__ ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl )); - if ( target ) - return ( 0.0 ); - else - return ( -0.0 ); - } -/******************************************************************************* -* Is 0.5 ² |x| < 1.0? * -*******************************************************************************/ - OldEnvironment.words.lo |= 0x02000000ul; - __asm__ ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl )); - if ( target ) - return ( 1.0 ); - else - return ( -1.0 ); - } -/******************************************************************************* -* Is 1.0 < |x| < 2.0^52? * -*******************************************************************************/ - if ( target ) - { // positive x - y = ( x + twoTo52 ) - twoTo52; // round at binary point - if ( y == x ) // exact case - return ( x ); - z = x + 0.5; // inexact case - y = ( z + twoTo52 ) - twoTo52; // round at binary point - if ( y > z ) - return ( y - 1.0 ); - else - return ( y ); - } - -/******************************************************************************* -* Is x < 0? * -*******************************************************************************/ - else - { - y = ( x - twoTo52 ) + twoTo52; // round at binary point - if ( y == x ) - return ( x ); - z = x - 0.5; - y = ( z - twoTo52 ) + twoTo52; // round at binary point - if ( y < z ) - return ( y + 1.0 ); - else - return ( y ); - } - } -/******************************************************************************* -* |x| >= 2.0^52 or x is a NaN. * -*******************************************************************************/ - return ( x ); - } -libm_hidden_def(round) diff --git a/libm/powerpc/classic/s_trunc.c b/libm/powerpc/classic/s_trunc.c deleted file mode 100644 index e9c668127..000000000 --- a/libm/powerpc/classic/s_trunc.c +++ /dev/null @@ -1,89 +0,0 @@ -#include <limits.h> -#include <math.h> -#include <endian.h> - -typedef union - { - struct { -#if (__BYTE_ORDER == __BIG_ENDIAN) - unsigned long int hi; - unsigned long int lo; -#else - unsigned long int lo; - unsigned long int hi; -#endif - } words; - double dbl; - } DblInHex; - -static const unsigned long int signMask = 0x80000000ul; -static const double twoTo52 = 4503599627370496.0; - -/******************************************************************************* -* * -* The function trunc truncates its double argument to integral value * -* and returns the result in double format. This function signals * -* inexact if an ordered return value is not equal to the operand. * -* * -*******************************************************************************/ - -libm_hidden_proto(trunc) -double trunc ( double x ) - { - DblInHex argument,OldEnvironment; - register double y; - register unsigned long int xhi; - register long int target; - - argument.dbl = x; - xhi = argument.words.hi & 0x7fffffffUL; // xhi <- high half of |x| - target = ( argument.words.hi < signMask ); // flag positive sign - - if ( xhi < 0x43300000ul ) -/******************************************************************************* -* Is |x| < 2.0^53? * -*******************************************************************************/ - { - if ( xhi < 0x3ff00000ul ) -/******************************************************************************* -* Is |x| < 1.0? * -*******************************************************************************/ - { - if ( ( xhi | argument.words.lo ) != 0ul ) - { // raise deserved INEXACT - __asm__ ("mffs %0" : "=f" (OldEnvironment.dbl)); - OldEnvironment.words.lo |= 0x02000000ul; - __asm__ ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl )); - } - if ( target ) // return properly signed zero - return ( 0.0 ); - else - return ( -0.0 ); - } -/******************************************************************************* -* Is 1.0 < |x| < 2.0^52? * -*******************************************************************************/ - if ( target ) - { - y = ( x + twoTo52 ) - twoTo52; // round at binary point - if ( y > x ) - return ( y - 1.0 ); - else - return ( y ); - } - - else - { - y = ( x - twoTo52 ) + twoTo52; // round at binary point. - if ( y < x ) - return ( y + 1.0 ); - else - return ( y ); - } - } -/******************************************************************************* -* Is |x| >= 2.0^52 or x is a NaN. * -*******************************************************************************/ - return ( x ); - } -libm_hidden_def(trunc) diff --git a/libm/powerpc/classic/w_scalb.c b/libm/powerpc/classic/w_scalb.c deleted file mode 100644 index 408136001..000000000 --- a/libm/powerpc/classic/w_scalb.c +++ /dev/null @@ -1,94 +0,0 @@ -/*********************************************************************** -** File: scalb.c -** -** Contains: C source code for implementations of floating-point -** scalb functions defined in header <fp.h>. In -** particular, this file contains implementations of -** functions scalb and scalbl for double and long double -** formats on PowerPC platforms. -** -** Written by: Jon Okada, SANEitation Engineer, ext. 4-4838 -** -** Copyright: © 1992 by Apple Computer, Inc., all rights reserved -** -** Change History ( most recent first ): -** -** 28 May 97 ali made an speed improvement for large n, -** removed scalbl. -** 12 Dec 92 JPO First created. -** -***********************************************************************/ - -#include <math.h> -#include <endian.h> - -typedef union - { - struct { -#if (__BYTE_ORDER == __BIG_ENDIAN) - unsigned long int hi; - unsigned long int lo; -#else - unsigned long int lo; - unsigned long int hi; -#endif - } words; - double dbl; - } DblInHex; - -static const double twoTo1023 = 8.988465674311579539e307; // 0x1p1023 -static const double twoToM1022 = 2.225073858507201383e-308; // 0x1p-1022 - - -/*********************************************************************** - double scalb( double x, long int n ) returns its argument x scaled - by the factor 2^m. NaNs, signed zeros, and infinities are propagated - by this function regardless of the value of n. - - Exceptions: OVERFLOW/INEXACT or UNDERFLOW inexact may occur; - INVALID for signaling NaN inputs ( quiet NaN returned ). - - Calls: none. -***********************************************************************/ - -libm_hidden_proto(scalb) -#ifdef _SCALB_INT -double scalb ( double x, int n ) -#else -double scalb ( double x, double n ) -#endif - { - DblInHex xInHex; - - xInHex.words.lo = 0UL; // init. low half of xInHex - - if ( n > 1023 ) - { // large positive scaling - if ( n > 2097 ) // huge scaling - return ( ( x * twoTo1023 ) * twoTo1023 ) * twoTo1023; - while ( n > 1023 ) - { // scale reduction loop - x *= twoTo1023; // scale x by 2^1023 - n -= 1023; // reduce n by 1023 - } - } - - else if ( n < -1022 ) - { // large negative scaling - if ( n < -2098 ) // huge negative scaling - return ( ( x * twoToM1022 ) * twoToM1022 ) * twoToM1022; - while ( n < -1022 ) - { // scale reduction loop - x *= twoToM1022; // scale x by 2^( -1022 ) - n += 1022; // incr n by 1022 - } - } - -/******************************************************************************* -* -1022 <= n <= 1023; convert n to double scale factor. * -*******************************************************************************/ - - xInHex.words.hi = ( ( unsigned long ) ( n + 1023 ) ) << 20; - return ( x * xInHex.dbl ); - } -libm_hidden_def(scalb) diff --git a/libm/s_fpclassify.c b/libm/s_fpclassify.c new file mode 100644 index 000000000..6175f51ad --- /dev/null +++ b/libm/s_fpclassify.c @@ -0,0 +1,44 @@ +/* Return classification value corresponding to argument. + Copyright (C) 1997, 2002 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, write to the Free + Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA + 02111-1307 USA. */ + +#include <math.h> + +#include "math_private.h" + +libm_hidden_proto (__fpclassify) +int +__fpclassify (double x) +{ + u_int32_t hx, lx; + int retval = FP_NORMAL; + + EXTRACT_WORDS (hx, lx, x); + lx |= hx & 0xfffff; + hx &= 0x7ff00000; + if ((hx | lx) == 0) + retval = FP_ZERO; + else if (hx == 0) + retval = FP_SUBNORMAL; + else if (hx == 0x7ff00000) + retval = lx != 0 ? FP_NAN : FP_INFINITE; + + return retval; +} +libm_hidden_def (__fpclassify) diff --git a/libm/s_fpclassifyf.c b/libm/s_fpclassifyf.c new file mode 100644 index 000000000..3a8fd56df --- /dev/null +++ b/libm/s_fpclassifyf.c @@ -0,0 +1,44 @@ +/* Return classification value corresponding to argument. + Copyright (C) 1997, 2000, 2002 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, write to the Free + Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA + 02111-1307 USA. */ + +#include <math.h> + +#include "math_private.h" + + +libm_hidden_proto (__fpclassifyf) +int +__fpclassifyf (float x) +{ + u_int32_t wx; + int retval = FP_NORMAL; + + GET_FLOAT_WORD (wx, x); + wx &= 0x7fffffff; + if (wx == 0) + retval = FP_ZERO; + else if (wx < 0x800000) + retval = FP_SUBNORMAL; + else if (wx >= 0x7f800000) + retval = wx > 0x7f800000 ? FP_NAN : FP_INFINITE; + + return retval; +} +libm_hidden_def (__fpclassifyf) diff --git a/libm/s_isinf.c b/libm/s_isinf.c new file mode 100644 index 000000000..ce509f072 --- /dev/null +++ b/libm/s_isinf.c @@ -0,0 +1,34 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Changed to return -1 for -Inf by Ulrich Drepper <drepper@cygnus.com>. + * Public domain. + */ + +#if defined(LIBM_SCCS) && !defined(lint) +static char rcsid[] = "$NetBSD: s_isinf.c,v 1.3 1995/05/11 23:20:14 jtc Exp $"; +#endif + +/* + * isinf(x) returns 1 is x is inf, -1 if x is -inf, else 0; + * no branching! + */ + +#include "math.h" +#include "math_private.h" + +libm_hidden_proto(__isinf) +int +__isinf (double x) +{ + int32_t hx,lx; + EXTRACT_WORDS(hx,lx,x); + lx |= (hx & 0x7fffffff) ^ 0x7ff00000; + lx |= -lx; + return ~(lx >> 31) & (hx >> 30); +} +libm_hidden_def (__isinf) +weak_alias (__isinf, isinf) +#ifdef NO_LONG_DOUBLE +strong_alias (__isinf, __isinfl) +weak_alias (__isinf, isinfl) +#endif diff --git a/libm/s_isinff.c b/libm/s_isinff.c new file mode 100644 index 000000000..33e274947 --- /dev/null +++ b/libm/s_isinff.c @@ -0,0 +1,30 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + */ + +#if defined(LIBM_SCCS) && !defined(lint) +static char rcsid[] = "$NetBSD: s_isinff.c,v 1.3 1995/05/11 23:20:21 jtc Exp $"; +#endif + +/* + * isinff(x) returns 1 if x is inf, -1 if x is -inf, else 0; + * no branching! + */ + +#include "math.h" +#include "math_private.h" + +libm_hidden_proto(__isinff) +int +__isinff (float x) +{ + int32_t ix,t; + GET_FLOAT_WORD(ix,x); + t = ix & 0x7fffffff; + t ^= 0x7f800000; + t |= -t; + return ~(t >> 31) & (ix >> 30); +} +libm_hidden_def (__isinff) +weak_alias (__isinff, isinff) diff --git a/libm/s_isnan.c b/libm/s_isnan.c new file mode 100644 index 000000000..6263d1c3f --- /dev/null +++ b/libm/s_isnan.c @@ -0,0 +1,45 @@ +/* @(#)s_isnan.c 5.1 93/09/24 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +#if defined(LIBM_SCCS) && !defined(lint) +static char rcsid[] = "$NetBSD: s_isnan.c,v 1.8 1995/05/10 20:47:36 jtc Exp $"; +#endif + +/* + * isnan(x) returns 1 is x is nan, else 0; + * no branching! + */ + +#include "math.h" +#include "math_private.h" + +libm_hidden_proto (__isnan) +#ifdef __STDC__ + int __isnan(double x) +#else + int __isnan(x) + double x; +#endif +{ + int32_t hx,lx; + EXTRACT_WORDS(hx,lx,x); + hx &= 0x7fffffff; + hx |= (u_int32_t)(lx|(-lx))>>31; + hx = 0x7ff00000 - hx; + return (int)(((u_int32_t)hx)>>31); +} +libm_hidden_def (__isnan) +weak_alias (__isnan, isnan) +#ifdef NO_LONG_DOUBLE +strong_alias (__isnan, __isnanl) +weak_alias (__isnan, isnanl) +#endif diff --git a/libm/s_isnanf.c b/libm/s_isnanf.c new file mode 100644 index 000000000..fa2d2fe89 --- /dev/null +++ b/libm/s_isnanf.c @@ -0,0 +1,43 @@ +/* s_isnanf.c -- float version of s_isnan.c. + * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. + */ + +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +#if defined(LIBM_SCCS) && !defined(lint) +static char rcsid[] = "$NetBSD: s_isnanf.c,v 1.4 1995/05/10 20:47:38 jtc Exp $"; +#endif + +/* + * isnanf(x) returns 1 is x is nan, else 0; + * no branching! + */ + +#include "math.h" +#include "math_private.h" + +libm_hidden_proto (__isnanf) +#ifdef __STDC__ + int __isnanf(float x) +#else + int __isnanf(x) + float x; +#endif +{ + int32_t ix; + GET_FLOAT_WORD(ix,x); + ix &= 0x7fffffff; + ix = 0x7f800000 - ix; + return (int)(((u_int32_t)(ix))>>31); +} +libm_hidden_def (__isnanf) +weak_alias (__isnanf, isnanf) diff --git a/libm/s_signbit.c b/libm/s_signbit.c new file mode 100644 index 000000000..c8a72c539 --- /dev/null +++ b/libm/s_signbit.c @@ -0,0 +1,34 @@ +/* Return nonzero value if number is negative. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, write to the Free + Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA + 02111-1307 USA. */ + +#include <math.h> + +#include "math_private.h" + +libm_hidden_proto(__signbit) +int +__signbit (double x) +{ + int32_t hx; + + GET_HIGH_WORD (hx, x); + return hx & 0x80000000; +} +libm_hidden_def(__signbit) diff --git a/libm/s_signbitf.c b/libm/s_signbitf.c new file mode 100644 index 000000000..ef83d6b7f --- /dev/null +++ b/libm/s_signbitf.c @@ -0,0 +1,34 @@ +/* Return nonzero value if number is negative. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, write to the Free + Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA + 02111-1307 USA. */ + +#include <math.h> + +#include "math_private.h" + +libm_hidden_proto(__signbitf) +int +__signbitf (float x) +{ + int32_t hx; + + GET_FLOAT_WORD (hx, x); + return hx & 0x80000000; +} +libm_hidden_def(__signbitf) |