From e7bcf43b6440ac9fc61a0eef5591393810daafb5 Mon Sep 17 00:00:00 2001 From: Joakim Tjernlund Date: Sat, 31 Mar 2007 13:28:15 +0000 Subject: From Steve Papacharalambous: Add math support for PowerPC e500. --- libm/Makefile.in | 21 +++ libm/powerpc/Makefile.arch | 18 -- libm/powerpc/classic/Makefile.arch | 18 ++ libm/powerpc/classic/s_ceil.c | 113 ++++++++++++ libm/powerpc/classic/s_copysign.c | 59 ++++++ libm/powerpc/classic/s_floor.c | 113 ++++++++++++ libm/powerpc/classic/s_frexp.c | 65 +++++++ libm/powerpc/classic/s_ldexp.c | 49 +++++ libm/powerpc/classic/s_logb.c | 107 +++++++++++ libm/powerpc/classic/s_modf.c | 344 +++++++++++++++++++++++++++++++++++ libm/powerpc/classic/s_nearbyint.c | 38 ++++ libm/powerpc/classic/s_rint.c | 159 ++++++++++++++++ libm/powerpc/classic/s_round.c | 115 ++++++++++++ libm/powerpc/classic/s_trunc.c | 89 +++++++++ libm/powerpc/classic/w_scalb.c | 94 ++++++++++ libm/powerpc/e500/Makefile.arch | 9 + libm/powerpc/e500/README.txt | 8 + libm/powerpc/e500/fpu/Makefile.arch | 19 ++ libm/powerpc/e500/fpu/fclrexcpt.c | 39 ++++ libm/powerpc/e500/fpu/fe_nomask.c | 32 ++++ libm/powerpc/e500/fpu/fedisblxcpt.c | 60 ++++++ libm/powerpc/e500/fpu/feenablxcpt.c | 60 ++++++ libm/powerpc/e500/fpu/fegetenv.c | 38 ++++ libm/powerpc/e500/fpu/fegetexcept.c | 31 ++++ libm/powerpc/e500/fpu/fegetround.c | 31 ++++ libm/powerpc/e500/fpu/feholdexcpt.c | 45 +++++ libm/powerpc/e500/fpu/fenv_const.c | 27 +++ libm/powerpc/e500/fpu/fenv_libc.h | 77 ++++++++ libm/powerpc/e500/fpu/fesetenv.c | 38 ++++ libm/powerpc/e500/fpu/fesetround.c | 37 ++++ libm/powerpc/e500/fpu/feupdateenv.c | 49 +++++ libm/powerpc/e500/fpu/fgetexcptflg.c | 39 ++++ libm/powerpc/e500/fpu/fraiseexcpt.c | 29 +++ libm/powerpc/e500/fpu/fsetexcptflg.c | 46 +++++ libm/powerpc/e500/fpu/ftestexcept.c | 32 ++++ libm/powerpc/e500/spe-raise.c | 67 +++++++ libm/powerpc/s_ceil.c | 113 ------------ libm/powerpc/s_copysign.c | 59 ------ libm/powerpc/s_floor.c | 113 ------------ libm/powerpc/s_frexp.c | 65 ------- libm/powerpc/s_ldexp.c | 49 ----- libm/powerpc/s_logb.c | 107 ----------- libm/powerpc/s_modf.c | 344 ----------------------------------- libm/powerpc/s_nearbyint.c | 38 ---- libm/powerpc/s_rint.c | 159 ---------------- libm/powerpc/s_round.c | 115 ------------ libm/powerpc/s_trunc.c | 89 --------- libm/powerpc/w_scalb.c | 94 ---------- 48 files changed, 2197 insertions(+), 1363 deletions(-) delete mode 100644 libm/powerpc/Makefile.arch create mode 100644 libm/powerpc/classic/Makefile.arch create mode 100644 libm/powerpc/classic/s_ceil.c create mode 100644 libm/powerpc/classic/s_copysign.c create mode 100644 libm/powerpc/classic/s_floor.c create mode 100644 libm/powerpc/classic/s_frexp.c create mode 100644 libm/powerpc/classic/s_ldexp.c create mode 100644 libm/powerpc/classic/s_logb.c create mode 100644 libm/powerpc/classic/s_modf.c create mode 100644 libm/powerpc/classic/s_nearbyint.c create mode 100644 libm/powerpc/classic/s_rint.c create mode 100644 libm/powerpc/classic/s_round.c create mode 100644 libm/powerpc/classic/s_trunc.c create mode 100644 libm/powerpc/classic/w_scalb.c create mode 100644 libm/powerpc/e500/Makefile.arch create mode 100644 libm/powerpc/e500/README.txt create mode 100644 libm/powerpc/e500/fpu/Makefile.arch create mode 100644 libm/powerpc/e500/fpu/fclrexcpt.c create mode 100644 libm/powerpc/e500/fpu/fe_nomask.c create mode 100644 libm/powerpc/e500/fpu/fedisblxcpt.c create mode 100644 libm/powerpc/e500/fpu/feenablxcpt.c create mode 100644 libm/powerpc/e500/fpu/fegetenv.c create mode 100644 libm/powerpc/e500/fpu/fegetexcept.c create mode 100644 libm/powerpc/e500/fpu/fegetround.c create mode 100644 libm/powerpc/e500/fpu/feholdexcpt.c create mode 100644 libm/powerpc/e500/fpu/fenv_const.c create mode 100644 libm/powerpc/e500/fpu/fenv_libc.h create mode 100644 libm/powerpc/e500/fpu/fesetenv.c create mode 100644 libm/powerpc/e500/fpu/fesetround.c create mode 100644 libm/powerpc/e500/fpu/feupdateenv.c create mode 100644 libm/powerpc/e500/fpu/fgetexcptflg.c create mode 100644 libm/powerpc/e500/fpu/fraiseexcpt.c create mode 100644 libm/powerpc/e500/fpu/fsetexcptflg.c create mode 100644 libm/powerpc/e500/fpu/ftestexcept.c create mode 100644 libm/powerpc/e500/spe-raise.c delete mode 100644 libm/powerpc/s_ceil.c delete mode 100644 libm/powerpc/s_copysign.c delete mode 100644 libm/powerpc/s_floor.c delete mode 100644 libm/powerpc/s_frexp.c delete mode 100644 libm/powerpc/s_ldexp.c delete mode 100644 libm/powerpc/s_logb.c delete mode 100644 libm/powerpc/s_modf.c delete mode 100644 libm/powerpc/s_nearbyint.c delete mode 100644 libm/powerpc/s_rint.c delete mode 100644 libm/powerpc/s_round.c delete mode 100644 libm/powerpc/s_trunc.c delete mode 100644 libm/powerpc/w_scalb.c (limited to 'libm') diff --git a/libm/Makefile.in b/libm/Makefile.in index 5b2fbff36..a972530cf 100644 --- a/libm/Makefile.in +++ b/libm/Makefile.in @@ -29,8 +29,21 @@ LIBS-libm.so := $(LIBS) libm_FULL_NAME := libm-$(VERSION).so + +# Fix builds for powerpc as there are different cores in this +# section now.` +ifeq ($(TARGET_ARCH),powerpc) +ifeq ($(CONFIG_E500),y) +libm_ARCH_DIR:=$(top_srcdir)libm/$(TARGET_ARCH)/e500 +libm_ARCH_OUT:=$(top_builddir)libm/$(TARGET_ARCH)/e500 +else +libm_ARCH_DIR:=$(top_srcdir)libm/$(TARGET_ARCH)/classic +libm_ARCH_OUT:=$(top_builddir)libm/$(TARGET_ARCH)/classic +endif +else libm_ARCH_DIR:=$(top_srcdir)libm/$(TARGET_ARCH) libm_ARCH_OUT:=$(top_builddir)libm/$(TARGET_ARCH) +endif libm_ARCH_fpu_DIR:=$(libm_ARCH_DIR)/fpu libm_ARCH_fpu_OUT:=$(libm_ARCH_OUT)/fpu @@ -91,7 +104,15 @@ libm_OUT := $(top_builddir)libm ifeq ($(UCLIBC_HAS_FPU),y) ifeq ($(DO_C99_MATH),y) ifneq ($(strip $(libm_ARCH_OBJS)),) +ifeq ($(TARGET_ARCH),powerpc) +ifeq ($(CONFIG_E500),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 # remove generic sources, if arch specific version is present ifneq ($(strip $(libm_ARCH_SRC)),) diff --git a/libm/powerpc/Makefile.arch b/libm/powerpc/Makefile.arch deleted file mode 100644 index 7c7600f80..000000000 --- a/libm/powerpc/Makefile.arch +++ /dev/null @@ -1,18 +0,0 @@ -# Makefile for uClibc -# -# Copyright (C) 2000-2006 Erik Andersen -# -# Licensed under the LGPL v2.1, see the file COPYING.LIB in this tarball. -# - -libm_ARCH_SRC:=$(wildcard $(libm_ARCH_DIR)/*.c) -libm_ARCH_OBJ:=$(patsubst $(libm_ARCH_DIR)/%.c,$(libm_ARCH_OUT)/%.o,$(libm_ARCH_SRC)) - -libm_ARCH_OBJS:=$(libm_ARCH_OBJ) - -ifeq ($(DOPIC),y) -libm-a-y+=$(libm_ARCH_OBJS:.o=.os) -else -libm-a-y+=$(libm_ARCH_OBJS) -endif -libm-so-y+=$(libm_ARCH_OBJS:.o=.os) diff --git a/libm/powerpc/classic/Makefile.arch b/libm/powerpc/classic/Makefile.arch new file mode 100644 index 000000000..7c7600f80 --- /dev/null +++ b/libm/powerpc/classic/Makefile.arch @@ -0,0 +1,18 @@ +# Makefile for uClibc +# +# Copyright (C) 2000-2006 Erik Andersen +# +# Licensed under the LGPL v2.1, see the file COPYING.LIB in this tarball. +# + +libm_ARCH_SRC:=$(wildcard $(libm_ARCH_DIR)/*.c) +libm_ARCH_OBJ:=$(patsubst $(libm_ARCH_DIR)/%.c,$(libm_ARCH_OUT)/%.o,$(libm_ARCH_SRC)) + +libm_ARCH_OBJS:=$(libm_ARCH_OBJ) + +ifeq ($(DOPIC),y) +libm-a-y+=$(libm_ARCH_OBJS:.o=.os) +else +libm-a-y+=$(libm_ARCH_OBJS) +endif +libm-so-y+=$(libm_ARCH_OBJS:.o=.os) diff --git a/libm/powerpc/classic/s_ceil.c b/libm/powerpc/classic/s_ceil.c new file mode 100644 index 000000000..8db5ce537 --- /dev/null +++ b/libm/powerpc/classic/s_ceil.c @@ -0,0 +1,113 @@ +/******************************************************************************* +* * +* 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 +#include + +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 new file mode 100644 index 000000000..1a988ac87 --- /dev/null +++ b/libm/powerpc/classic/s_copysign.c @@ -0,0 +1,59 @@ +/******************************************************************************* +* * +* 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 +#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 new file mode 100644 index 000000000..ff3436707 --- /dev/null +++ b/libm/powerpc/classic/s_floor.c @@ -0,0 +1,113 @@ +/******************************************************************************* +* * +* 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 +#include + +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 new file mode 100644 index 000000000..001aaf708 --- /dev/null +++ b/libm/powerpc/classic/s_frexp.c @@ -0,0 +1,65 @@ +/******************************************************************************* +* * +* 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 +#include +#include + +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 new file mode 100644 index 000000000..10100d7c2 --- /dev/null +++ b/libm/powerpc/classic/s_ldexp.c @@ -0,0 +1,49 @@ +/******************************************************************************* +* * +* 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 +#include +#include + +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 new file mode 100644 index 000000000..81daa412e --- /dev/null +++ b/libm/powerpc/classic/s_logb.c @@ -0,0 +1,107 @@ +/******************************************************************************* +* * +* 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 +#include + +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 new file mode 100644 index 000000000..b9d69445d --- /dev/null +++ b/libm/powerpc/classic/s_modf.c @@ -0,0 +1,344 @@ +/******************************************************************************* +** File: rndint.c +** +** Contains: C source code for implementations of floating-point +** functions which round to integral value or format, as +** defined in header . 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 +#include +#include + +#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 new file mode 100644 index 000000000..068e21378 --- /dev/null +++ b/libm/powerpc/classic/s_nearbyint.c @@ -0,0 +1,38 @@ +#include +#include + +/******************************************************************************* +* * +* 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 new file mode 100644 index 000000000..dd2ae585e --- /dev/null +++ b/libm/powerpc/classic/s_rint.c @@ -0,0 +1,159 @@ +/******************************************************************************* +** File: rndint.c +** +** Contains: C source code for implementations of floating-point +** functions which round to integral value or format, as +** defined in header . 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 +#include +#include + +#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 new file mode 100644 index 000000000..62d5936d9 --- /dev/null +++ b/libm/powerpc/classic/s_round.c @@ -0,0 +1,115 @@ +#include +#include +#include + +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 new file mode 100644 index 000000000..f793992a7 --- /dev/null +++ b/libm/powerpc/classic/s_trunc.c @@ -0,0 +1,89 @@ +#include +#include +#include + +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 new file mode 100644 index 000000000..408136001 --- /dev/null +++ b/libm/powerpc/classic/w_scalb.c @@ -0,0 +1,94 @@ +/*********************************************************************** +** File: scalb.c +** +** Contains: C source code for implementations of floating-point +** scalb functions defined in header . 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 +#include + +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/powerpc/e500/Makefile.arch b/libm/powerpc/e500/Makefile.arch new file mode 100644 index 000000000..bec21caef --- /dev/null +++ b/libm/powerpc/e500/Makefile.arch @@ -0,0 +1,9 @@ +# Makefile for uClibc +# +# Copyright (C) 2000-2006 Erik Andersen +# +# Licensed under the LGPL v2.1, see the file COPYING.LIB in this tarball. +# + +-include $(libm_ARCH_fpu_DIR)/Makefile.arch + diff --git a/libm/powerpc/e500/README.txt b/libm/powerpc/e500/README.txt new file mode 100644 index 000000000..354d7921a --- /dev/null +++ b/libm/powerpc/e500/README.txt @@ -0,0 +1,8 @@ +The routines this math library have been derived from the e500 math +library in eglibc. The original e500 port to glibc was done by +Aldy Hernandez, , and the port to eglibc was done by +Joseph S. Myers, + +It has been ported to uClibc by Steve Papacharalambous + 19 December, 2006 + diff --git a/libm/powerpc/e500/fpu/Makefile.arch b/libm/powerpc/e500/fpu/Makefile.arch new file mode 100644 index 000000000..8f00e0916 --- /dev/null +++ b/libm/powerpc/e500/fpu/Makefile.arch @@ -0,0 +1,19 @@ +# Makefile for uClibc +# +# Copyright (C) 2000-2006 Erik Andersen +# +# Licensed under the LGPL v2.1, see the file COPYING.LIB in this tarball. +# + + +libm_ARCH_SRC:=$(wildcard $(libm_ARCH_fpu_DIR)/*.c) +libm_ARCH_OBJ:=$(patsubst $(libm_ARCH_fpu_DIR)/%.c,$(libm_ARCH_fpu_OUT)/%.o,$(libm_ARCH_SRC)) + +libm_ARCH_OBJS:=$(libm_ARCH_OBJ) + +ifeq ($(DOPIC),y) +libm-a-y+=$(libm_ARCH_OBJS:.o=.os) +else +libm-a-y+=$(libm_ARCH_OBJS) +endif +libm-so-y+=$(libm_ARCH_OBJS:.o=.os) diff --git a/libm/powerpc/e500/fpu/fclrexcpt.c b/libm/powerpc/e500/fpu/fclrexcpt.c new file mode 100644 index 000000000..42e8b6bab --- /dev/null +++ b/libm/powerpc/e500/fpu/fclrexcpt.c @@ -0,0 +1,39 @@ +/* Clear given exceptions in current floating-point environment. + Copyright (C) 2004 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + 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 "fenv_libc.h" +#undef feclearexcept +int +__feclearexcept (int excepts) +{ + unsigned int fpescr; + + /* Get the current state. */ + fpescr = fegetenv_register (); + + /* Clear the relevant bits. */ + fpescr &= ~(excepts & FE_ALL_EXCEPT); + + /* Put the new state in effect. */ + fesetenv_register (fpescr); + + /* Success. */ + return 0; +} + diff --git a/libm/powerpc/e500/fpu/fe_nomask.c b/libm/powerpc/e500/fpu/fe_nomask.c new file mode 100644 index 000000000..ad7b324c0 --- /dev/null +++ b/libm/powerpc/e500/fpu/fe_nomask.c @@ -0,0 +1,32 @@ +/* Procedure definition for FE_NOMASK_ENV. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + 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 +#include + +/* This is presently a stub, until it's decided how the kernels should + support this. */ + +const fenv_t * +__fe_nomask_env(void) +{ + __set_errno (ENOSYS); + return FE_ENABLED_ENV; +} + diff --git a/libm/powerpc/e500/fpu/fedisblxcpt.c b/libm/powerpc/e500/fpu/fedisblxcpt.c new file mode 100644 index 000000000..d2c5caba7 --- /dev/null +++ b/libm/powerpc/e500/fpu/fedisblxcpt.c @@ -0,0 +1,60 @@ +/* Disable floating-point exceptions. + Copyright (C) 2004 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Aldy Hernandez , 2004. + + 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 "fenv_libc.h" +#include +#include + +int +fedisableexcept (int excepts) +{ + unsigned int result = 0, pflags, r; + INTERNAL_SYSCALL_DECL (err); + + INTERNAL_SYSCALL (prctl, err, 2, PR_GET_FPEXC, &pflags); + + /* Save old enable bits. */ + if (pflags & PR_FP_EXC_OVF) + result |= FE_OVERFLOW; + if (pflags & PR_FP_EXC_UND) + result |= FE_UNDERFLOW; + if (pflags & PR_FP_EXC_INV) + result |= FE_INVALID; + if (pflags & PR_FP_EXC_DIV) + result |= FE_DIVBYZERO; + if (pflags & PR_FP_EXC_RES) + result |= FE_INEXACT; + + if (excepts & FE_INEXACT) + pflags &= ~PR_FP_EXC_RES; + if (excepts & FE_DIVBYZERO) + pflags &= ~PR_FP_EXC_DIV; + if (excepts & FE_UNDERFLOW) + pflags &= ~PR_FP_EXC_UND; + if (excepts & FE_OVERFLOW) + pflags &= ~PR_FP_EXC_OVF; + if (excepts & FE_INVALID) + pflags &= ~PR_FP_EXC_INV; + r = INTERNAL_SYSCALL (prctl, err, 2, PR_SET_FPEXC, pflags); + if (INTERNAL_SYSCALL_ERROR_P (r, err)) + return -1; + + return result; +} diff --git a/libm/powerpc/e500/fpu/feenablxcpt.c b/libm/powerpc/e500/fpu/feenablxcpt.c new file mode 100644 index 000000000..4f00a918d --- /dev/null +++ b/libm/powerpc/e500/fpu/feenablxcpt.c @@ -0,0 +1,60 @@ +/* Enable floating-point exceptions. + Copyright (C) 2004 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Aldy Hernandez , 2004. + + 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 "fenv_libc.h" +#include +#include + +int +feenableexcept (int excepts) +{ + unsigned int result = 0, pflags, r; + INTERNAL_SYSCALL_DECL (err); + + INTERNAL_SYSCALL (prctl, err, 2, PR_GET_FPEXC, &pflags); + + /* Save old enable bits. */ + if (pflags & PR_FP_EXC_OVF) + result |= FE_OVERFLOW; + if (pflags & PR_FP_EXC_UND) + result |= FE_UNDERFLOW; + if (pflags & PR_FP_EXC_INV) + result |= FE_INVALID; + if (pflags & PR_FP_EXC_DIV) + result |= FE_DIVBYZERO; + if (pflags & PR_FP_EXC_RES) + result |= FE_INEXACT; + + if (excepts & FE_INEXACT) + pflags |= PR_FP_EXC_RES; + if (excepts & FE_DIVBYZERO) + pflags |= PR_FP_EXC_DIV; + if (excepts & FE_UNDERFLOW) + pflags |= PR_FP_EXC_UND; + if (excepts & FE_OVERFLOW) + pflags |= PR_FP_EXC_OVF; + if (excepts & FE_INVALID) + pflags |= PR_FP_EXC_INV; + r = INTERNAL_SYSCALL (prctl, err, 2, PR_SET_FPEXC, pflags); + if (INTERNAL_SYSCALL_ERROR_P (r, err)) + return -1; + + return result; +} diff --git a/libm/powerpc/e500/fpu/fegetenv.c b/libm/powerpc/e500/fpu/fegetenv.c new file mode 100644 index 000000000..fbafdc208 --- /dev/null +++ b/libm/powerpc/e500/fpu/fegetenv.c @@ -0,0 +1,38 @@ +/* Store current floating-point environment. + Copyright (C) 2004 Free Software Foundation, Inc. + Contributed by Aldy Hernandez , 2004 + This file is part of the GNU C Library. + + 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 "fenv_libc.h" +#include +#include + +int +__fegetenv (fenv_t *envp) +{ + fenv_union_t u; + INTERNAL_SYSCALL_DECL (err); + + INTERNAL_SYSCALL (prctl, err, 2, PR_GET_FPEXC, &u.l[0]); + u.l[1] = fegetenv_register (); + *envp = u.fenv; + + /* Success. */ + return 0; +} + diff --git a/libm/powerpc/e500/fpu/fegetexcept.c b/libm/powerpc/e500/fpu/fegetexcept.c new file mode 100644 index 000000000..c0398461c --- /dev/null +++ b/libm/powerpc/e500/fpu/fegetexcept.c @@ -0,0 +1,31 @@ +/* Get floating-point exceptions. + Copyright (C) 2004 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Aldy Hernandez , 2004. + + 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 "fenv_libc.h" + +int +fegetexcept (void) +{ + unsigned long fpescr; + + fpescr = fegetenv_register (); + + return fpescr & FE_ALL_EXCEPT; +} diff --git a/libm/powerpc/e500/fpu/fegetround.c b/libm/powerpc/e500/fpu/fegetround.c new file mode 100644 index 000000000..6e568fa23 --- /dev/null +++ b/libm/powerpc/e500/fpu/fegetround.c @@ -0,0 +1,31 @@ +/* Return current rounding direction. + Copyright (C) 2004 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Aldy Hernandez , 2004. + + 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 "fenv_libc.h" + +#undef fegetround +int +fegetround (void) +{ + unsigned long fpescr; + + fpescr = fegetenv_register (); + return fpescr & 3; +} diff --git a/libm/powerpc/e500/fpu/feholdexcpt.c b/libm/powerpc/e500/fpu/feholdexcpt.c new file mode 100644 index 000000000..eccde4217 --- /dev/null +++ b/libm/powerpc/e500/fpu/feholdexcpt.c @@ -0,0 +1,45 @@ +/* Store current floating-point environment and clear exceptions. + Copyright (C) 2004, 2006 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Aldy Hernandez , 2004. + + 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 "fenv_libc.h" +#include +#include + +int +feholdexcept (fenv_t *envp) +{ + fenv_union_t u; + INTERNAL_SYSCALL_DECL (err); + + + /* Get the current state. */ + INTERNAL_SYSCALL (prctl, err, 2, PR_GET_FPEXC, &u.l[0]); + u.l[1] = fegetenv_register (); + *envp = u.fenv; + + /* Clear everything except for the rounding mode. */ + u.l[1] &= 3; + + /* Put the new state in effect. */ + fesetenv_register (u.l[1]); + + return 0; +} +libm_hidden_def (feholdexcept) diff --git a/libm/powerpc/e500/fpu/fenv_const.c b/libm/powerpc/e500/fpu/fenv_const.c new file mode 100644 index 000000000..073fc9277 --- /dev/null +++ b/libm/powerpc/e500/fpu/fenv_const.c @@ -0,0 +1,27 @@ +/* Constants for fenv_bits.h. + Copyright (C) 2004 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Aldy Hernandez , 2004. + + 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. */ + +/* If the default argument is used we use this value. */ +const unsigned long long __fe_dfl_env __attribute__ ((aligned (8))) = +0x0ULL; + +/* Floating-point environment where none of the exceptions are masked. */ +const unsigned long long __fe_enabled_env __attribute__ ((aligned (8))) = +0xF00000000ULL; diff --git a/libm/powerpc/e500/fpu/fenv_libc.h b/libm/powerpc/e500/fpu/fenv_libc.h new file mode 100644 index 000000000..cd003eab7 --- /dev/null +++ b/libm/powerpc/e500/fpu/fenv_libc.h @@ -0,0 +1,77 @@ +/* Internal libc stuff for floating point environment routines. + Copyright (C) 2004 Free Software Foundation, Inc. + Contributed by Aldy Hernandez , 2004 + This file is part of the GNU C Library. + + 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. */ + +#ifndef _FENV_LIBC_H +#define _FENV_LIBC_H 1 + +#include + +extern int __feraiseexcept_internal (int __excepts); + +/* Equivalent to fegetenv, but returns a fenv_t instead of taking a + pointer. */ +#define fegetenv_register() \ + ({ unsigned fscr; asm volatile ("mfspefscr %0" : "=r" (fscr)); fscr; }) + +/* Equivalent to fesetenv, but takes a fenv_t instead of a pointer. */ +#define fesetenv_register(fscr) \ + ({ asm volatile ("mtspefscr %0" : : "r" (fscr)); }) + +typedef union +{ + fenv_t fenv; + unsigned int l[2]; +} fenv_union_t; + +/* Definitions of all the SPEFSCR bit numbers. */ +enum { + SPEFSCR_SOVH = 0x80000000, + SPEFSCR_OVH = 0x40000000, + SPEFSCR_FGH = 0x20000000, + SPEFSCR_FXH = 0x10000000, + SPEFSCR_FINVH = 0x08000000, + SPEFSCR_FDBZH = 0x04000000, + SPEFSCR_FUNFH = 0x02000000, + SPEFSCR_FOVFH = 0x01000000, + /* 2 unused bits. */ + SPEFSCR_FINXS = 0x00200000, + SPEFSCR_FINVS = 0x00100000, + SPEFSCR_FDBZS = 0x00080000, + SPEFSCR_FUNFS = 0x00040000, + SPEFSCR_FOVFS = 0x00020000, + SPEFSCR_MODE = 0x00010000, + SPEFSCR_SOV = 0x00008000, + SPEFSCR_OV = 0x00004000, + SPEFSCR_FG = 0x00002000, + SPEFSCR_FX = 0x00001000, + SPEFSCR_FINV = 0x00000800, + SPEFSCR_FDBZ = 0x00000400, + SPEFSCR_FUNF = 0x00000200, + SPEFSCR_FOVF = 0x00000100, + /* 1 unused bit. */ + SPEFSCR_FINXE = 0x00000040, + SPEFSCR_FINVE = 0x00000020, + SPEFSCR_FDBZE = 0x00000010, + SPEFSCR_FUNFE = 0x00000008, + SPEFSCR_FOVFE = 0x00000004, + SPEFSCR_FRMC = 0x00000003 +}; + +#endif /* fenv_libc.h */ diff --git a/libm/powerpc/e500/fpu/fesetenv.c b/libm/powerpc/e500/fpu/fesetenv.c new file mode 100644 index 000000000..cc0cc18be --- /dev/null +++ b/libm/powerpc/e500/fpu/fesetenv.c @@ -0,0 +1,38 @@ +/* Install given floating-point environment. + Copyright (C) 1997,99,2000,01,02 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + 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 "fenv_libc.h" +#include +#include + +int +__fesetenv (const fenv_t *envp) +{ + fenv_union_t u; + INTERNAL_SYSCALL_DECL (err); + + u.fenv = *envp; + INTERNAL_SYSCALL (prctl, err, 2, PR_SET_FPEXC, &u.l[0]); + fesetenv_register (u.l[1]); + + /* Success. */ + return 0; +} + +libm_hidden_ver (__fesetenv, fesetenv) diff --git a/libm/powerpc/e500/fpu/fesetround.c b/libm/powerpc/e500/fpu/fesetround.c new file mode 100644 index 000000000..0f29a91a6 --- /dev/null +++ b/libm/powerpc/e500/fpu/fesetround.c @@ -0,0 +1,37 @@ +/* Set current rounding direction. + Copyright (C) 2004, 2006 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Aldy Hernandez , 2004. + + 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 "fenv_libc.h" + +int +fesetround (int round) +{ + unsigned long fpescr; + + if ((unsigned int) round > 3) + return 1; + + fpescr = fegetenv_register (); + fpescr = (fpescr & ~SPEFSCR_FRMC) | (round & 3); + fesetenv_register (fpescr); + + return 0; +} +libm_hidden_def (fesetround) diff --git a/libm/powerpc/e500/fpu/feupdateenv.c b/libm/powerpc/e500/fpu/feupdateenv.c new file mode 100644 index 000000000..03f3af8d9 --- /dev/null +++ b/libm/powerpc/e500/fpu/feupdateenv.c @@ -0,0 +1,49 @@ +/* Install given floating-point environment and raise exceptions. + Copyright (C) 2004 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Aldy Hernandez , 2004. + + 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 "fenv_libc.h" +#include +#include + +int +__feupdateenv (const fenv_t *envp) +{ + unsigned long fpescr, old, new, pflags; + fenv_union_t u; + INTERNAL_SYSCALL_DECL (err); + + /* Save the currently set exceptions. */ + u.fenv = *envp; + new = u.l[1]; + old = fegetenv_register (); + new |= (old & FE_ALL_EXCEPT); + + INTERNAL_SYSCALL (prctl, err, 2, PR_GET_FPEXC, &pflags); + pflags |= u.l[0]; + INTERNAL_SYSCALL (prctl, err, 2, PR_SET_FPEXC, pflags); + + /* Enable and raise (if appropriate) exceptions set in `new'. */ + fesetenv_register (new); + feraiseexcept (new & FE_ALL_EXCEPT); + + /* Success. */ + return 0; +} + diff --git a/libm/powerpc/e500/fpu/fgetexcptflg.c b/libm/powerpc/e500/fpu/fgetexcptflg.c new file mode 100644 index 000000000..62d7f2c32 --- /dev/null +++ b/libm/powerpc/e500/fpu/fgetexcptflg.c @@ -0,0 +1,39 @@ +/* Store current representation for exceptions. + Copyright (C) 2004 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Aldy Hernandez , 2004. + + 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 "fenv_libc.h" + +int +__fegetexceptflag (fexcept_t *flagp, int excepts) +{ + unsigned long fpescr; + + /* Get the current state. */ + fpescr = fegetenv_register (); + + /* ?? Classic PPC doesn't do anything with `excepts', so we'll do + the same here. (We should really be ignoring exceptions in + excepts) ?? */ + *flagp = fpescr & FE_ALL_EXCEPT; + + /* Success. */ + return 0; +} + diff --git a/libm/powerpc/e500/fpu/fraiseexcpt.c b/libm/powerpc/e500/fpu/fraiseexcpt.c new file mode 100644 index 000000000..39caf79d3 --- /dev/null +++ b/libm/powerpc/e500/fpu/fraiseexcpt.c @@ -0,0 +1,29 @@ +/* Raise given exceptions. + Copyright (C) 2004, 2006 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Aldy Hernandez , 2004. + + 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 "fenv_libc.h" + +#undef feraiseexcept + +#define __FERAISEEXCEPT_INTERNAL feraiseexcept + +#include "../spe-raise.c" + +libm_hidden_def (feraiseexcept) diff --git a/libm/powerpc/e500/fpu/fsetexcptflg.c b/libm/powerpc/e500/fpu/fsetexcptflg.c new file mode 100644 index 000000000..8e4bb895d --- /dev/null +++ b/libm/powerpc/e500/fpu/fsetexcptflg.c @@ -0,0 +1,46 @@ +/* Set floating-point environment exception handling. + Copyright (C) 1997,99,2000,01,04 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + 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 "fenv_libc.h" + +int +__fesetexceptflag (const fexcept_t *flagp, int excepts) +{ + unsigned long spefscr; + fexcept_t flag; + + /* Get the current state. */ + spefscr = fegetenv_register (); + + /* Ignore exceptions not listed in 'excepts'. */ + flag = *flagp & excepts; + + /* Replace the exception status */ + spefscr = (spefscr & ~FE_ALL_EXCEPT) | flag; + + /* Store the new status word (along with the rest of the environment). + This may cause floating-point exceptions if the restored state + requests it. */ + fesetenv_register (spefscr); + feraiseexcept (spefscr & FE_ALL_EXCEPT); + + /* Success. */ + return 0; +} + diff --git a/libm/powerpc/e500/fpu/ftestexcept.c b/libm/powerpc/e500/fpu/ftestexcept.c new file mode 100644 index 000000000..971c51965 --- /dev/null +++ b/libm/powerpc/e500/fpu/ftestexcept.c @@ -0,0 +1,32 @@ +/* Test exception in current environment. + Copyright (C) 2004 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Aldy Hernandez , 2004. + + 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 "fenv_libc.h" + +int +fetestexcept (int excepts) +{ + unsigned long f; + + /* Get the current state. */ + f = fegetenv_register (); + + return f & excepts; +} diff --git a/libm/powerpc/e500/spe-raise.c b/libm/powerpc/e500/spe-raise.c new file mode 100644 index 000000000..fb53dcec7 --- /dev/null +++ b/libm/powerpc/e500/spe-raise.c @@ -0,0 +1,67 @@ +/* Raise given exceptions. + Copyright (C) 1997,99,2000,01,02,04 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + 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 "fpu/fenv_libc.h" + +int +__FERAISEEXCEPT_INTERNAL (int excepts) +{ + unsigned long f; + + f = fegetenv_register (); + f |= (excepts & FE_ALL_EXCEPT); + fesetenv_register (f); + + /* Force the operations that cause the exceptions. */ + if ((FE_INVALID & excepts) != 0) + { + /* ?? Does not set sticky bit ?? */ + /* 0 / 0 */ + asm volatile ("efsdiv %0,%0,%1" : : "r" (0), "r" (0)); + } + + if ((FE_DIVBYZERO & excepts) != 0) + { + /* 1.0 / 0.0 */ + asm volatile ("efsdiv %0,%0,%1" : : "r" (1.0F), "r" (0)); + } + + if ((FE_OVERFLOW & excepts) != 0) + { + /* ?? Does not set sticky bit ?? */ + /* Largest normalized number plus itself. */ + asm volatile ("efsadd %0,%0,%1" : : "r" (0x7f7fffff), "r" (0x7f7fffff)); + } + + if ((FE_UNDERFLOW & excepts) != 0) + { + /* ?? Does not set sticky bit ?? */ + /* Smallest normalized number times itself. */ + asm volatile ("efsmul %0,%0,%1" : : "r" (0x800000), "r" (0x800000)); + } + + if ((FE_INEXACT & excepts) != 0) + { + /* Smallest normalized minus 1.0 raises the inexact flag. */ + asm volatile ("efssub %0,%0,%1" : : "r" (0x00800000), "r" (1.0F)); + } + + /* Success. */ + return 0; +} diff --git a/libm/powerpc/s_ceil.c b/libm/powerpc/s_ceil.c deleted file mode 100644 index 8db5ce537..000000000 --- a/libm/powerpc/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 -#include - -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/s_copysign.c b/libm/powerpc/s_copysign.c deleted file mode 100644 index 1a988ac87..000000000 --- a/libm/powerpc/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 -#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/s_floor.c b/libm/powerpc/s_floor.c deleted file mode 100644 index ff3436707..000000000 --- a/libm/powerpc/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 -#include - -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/s_frexp.c b/libm/powerpc/s_frexp.c deleted file mode 100644 index 001aaf708..000000000 --- a/libm/powerpc/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 -#include -#include - -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/s_ldexp.c b/libm/powerpc/s_ldexp.c deleted file mode 100644 index 10100d7c2..000000000 --- a/libm/powerpc/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 -#include -#include - -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/s_logb.c b/libm/powerpc/s_logb.c deleted file mode 100644 index 81daa412e..000000000 --- a/libm/powerpc/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 -#include - -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/s_modf.c b/libm/powerpc/s_modf.c deleted file mode 100644 index b9d69445d..000000000 --- a/libm/powerpc/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 . 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 -#include -#include - -#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/s_nearbyint.c b/libm/powerpc/s_nearbyint.c deleted file mode 100644 index 068e21378..000000000 --- a/libm/powerpc/s_nearbyint.c +++ /dev/null @@ -1,38 +0,0 @@ -#include -#include - -/******************************************************************************* -* * -* 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/s_rint.c b/libm/powerpc/s_rint.c deleted file mode 100644 index dd2ae585e..000000000 --- a/libm/powerpc/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 . 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 -#include -#include - -#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/s_round.c b/libm/powerpc/s_round.c deleted file mode 100644 index 62d5936d9..000000000 --- a/libm/powerpc/s_round.c +++ /dev/null @@ -1,115 +0,0 @@ -#include -#include -#include - -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/s_trunc.c b/libm/powerpc/s_trunc.c deleted file mode 100644 index f793992a7..000000000 --- a/libm/powerpc/s_trunc.c +++ /dev/null @@ -1,89 +0,0 @@ -#include -#include -#include - -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/w_scalb.c b/libm/powerpc/w_scalb.c deleted file mode 100644 index 408136001..000000000 --- a/libm/powerpc/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 . 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 -#include - -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) -- cgit v1.2.3