diff options
-rw-r--r-- | libm/powerpc/Makefile | 2 | ||||
-rw-r--r-- | libm/powerpc/s_ceil.c | 109 | ||||
-rw-r--r-- | libm/powerpc/s_copysign.c (renamed from libm/powerpc/sign.c) | 0 | ||||
-rw-r--r-- | libm/powerpc/s_floor.c (renamed from libm/powerpc/ceilfloor.c) | 65 | ||||
-rw-r--r-- | libm/powerpc/s_frexp.c (renamed from libm/powerpc/frexpldexp.c) | 9 | ||||
-rw-r--r-- | libm/powerpc/s_ldexp.c | 47 | ||||
-rw-r--r-- | libm/powerpc/s_logb.c (renamed from libm/powerpc/logb.c) | 0 | ||||
-rw-r--r-- | libm/powerpc/s_modf.c (renamed from libm/powerpc/rndint.c) | 85 | ||||
-rw-r--r-- | libm/powerpc/s_rint.c | 157 | ||||
-rw-r--r-- | libm/powerpc/w_scalb.c (renamed from libm/powerpc/scalb.c) | 0 |
10 files changed, 314 insertions, 160 deletions
diff --git a/libm/powerpc/Makefile b/libm/powerpc/Makefile index 9872d6077..dec4c49f3 100644 --- a/libm/powerpc/Makefile +++ b/libm/powerpc/Makefile @@ -27,7 +27,7 @@ TARGET_CC= $(TOPDIR)extra/gcc-uClibc/$(TARGET_ARCH)-uclibc-gcc TARGET_CFLAGS+=-D_IEEE_LIBM -D_ISOC99_SOURCE -D_SVID_SOURCE ifeq ($(strip $(DO_C99_MATH)),true) -CSRC = ceilfloor.c frexpldexp.c logb.c rndint.c scalb.c sign.c +CSRC = s_ceil.c s_floor.c s_ldexp.c s_frexp.c s_logb.c s_modf.c w_scalb.c s_copysign.c s_rint.c else CSRC = endif diff --git a/libm/powerpc/s_ceil.c b/libm/powerpc/s_ceil.c new file mode 100644 index 000000000..790dd7ba0 --- /dev/null +++ b/libm/powerpc/s_ceil.c @@ -0,0 +1,109 @@ +/******************************************************************************* +* * +* 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 * +* * +*******************************************************************************/ + +static const double twoTo52 = 4503599627370496.0; +static const unsigned long signMask = 0x80000000ul; + +typedef union + { + struct { +#if defined(__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. * +*******************************************************************************/ + +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 ); + } + diff --git a/libm/powerpc/sign.c b/libm/powerpc/s_copysign.c index 18e91998a..18e91998a 100644 --- a/libm/powerpc/sign.c +++ b/libm/powerpc/s_copysign.c diff --git a/libm/powerpc/ceilfloor.c b/libm/powerpc/s_floor.c index d5947ab50..f7b4534bb 100644 --- a/libm/powerpc/ceilfloor.c +++ b/libm/powerpc/s_floor.c @@ -43,71 +43,6 @@ typedef union *******************************************************************************/ /******************************************************************************* -* Ceil(x) returns the smallest integer not less than x. * -*******************************************************************************/ - -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 ); - } - -/******************************************************************************* * Floor(x) returns the largest integer not greater than x. * *******************************************************************************/ diff --git a/libm/powerpc/frexpldexp.c b/libm/powerpc/s_frexp.c index 6b759f91b..0ca7c1514 100644 --- a/libm/powerpc/frexpldexp.c +++ b/libm/powerpc/s_frexp.c @@ -38,15 +38,6 @@ typedef union double dbl; } DblInHex; -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 ); - } - double frexp ( double value, int *eptr ) { DblInHex argument; diff --git a/libm/powerpc/s_ldexp.c b/libm/powerpc/s_ldexp.c new file mode 100644 index 000000000..0e9a66611 --- /dev/null +++ b/libm/powerpc/s_ldexp.c @@ -0,0 +1,47 @@ +/******************************************************************************* +* * +* File frexpldexp.c, * +* Functions frexp(x) and ldexp(x), * +* Implementation of frexp and ldexp functions for the PowerPC. * +* * +* Copyright © 1991 Apple Computer, Inc. All rights reserved. * +* * +* Written by Ali Sazegari, started on January 1991, * +* * +* W A R N I N G: This routine expects a 64 bit double model. * +* * +* December03 1992: first rs6000 implementation. * +* October 05 1993: added special cases for NaN and ° in frexp. * +* May 27 1997: improved the performance of frexp by eliminating the * +* switch statement. * +* June 13 2001: (ram) rewrote frexp to eliminate calls to scalb and * +* logb. * +* * +*******************************************************************************/ + +#include <limits.h> +#include <math.h> + +typedef union + { + struct { +#if defined(__BIG_ENDIAN__) + unsigned long int hi; + unsigned long int lo; +#else + unsigned long int lo; + unsigned long int hi; +#endif + } words; + double dbl; + } DblInHex; + +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 ); + } + diff --git a/libm/powerpc/logb.c b/libm/powerpc/s_logb.c index 3a410ba5c..3a410ba5c 100644 --- a/libm/powerpc/logb.c +++ b/libm/powerpc/s_logb.c diff --git a/libm/powerpc/rndint.c b/libm/powerpc/s_modf.c index 1b7a9ec81..66c0655f5 100644 --- a/libm/powerpc/rndint.c +++ b/libm/powerpc/s_modf.c @@ -69,91 +69,6 @@ static const double doubleToLong = 4503603922337792.0; // 2^52 static const DblInHex Huge = {{ 0x7FF00000, 0x00000000 }}; static const DblInHex TOWARDZERO = {{ 0x00000000, 0x00000001 }}; -/******************************************************************************* -* * -* 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 ); - } /******************************************************************************* * * diff --git a/libm/powerpc/s_rint.c b/libm/powerpc/s_rint.c new file mode 100644 index 000000000..47a79ae5e --- /dev/null +++ b/libm/powerpc/s_rint.c @@ -0,0 +1,157 @@ +/******************************************************************************* +** File: rndint.c +** +** Contains: C source code for implementations of floating-point +** functions which round to integral value or format, as +** defined in header <fp.h>. In particular, this file +** contains implementations of functions rint, nearbyint, +** rinttol, round, roundtol, trunc, modf and modfl. This file +** targets PowerPC or Power platforms. +** +** Written by: A. Sazegari, Apple AltiVec Group +** Created originally by Jon Okada, Apple Numerics Group +** +** Copyright: © 1992-2001 by Apple Computer, Inc., all rights reserved +** +** Change History (most recent first): +** +** 13 Jul 01 ram replaced --setflm calls with inline assembly +** 03 Mar 01 ali first port to os x using gcc, added the crucial __setflm +** definition. +** 1. removed double_t, put in double for now. +** 2. removed iclass from nearbyint. +** 3. removed wrong comments intrunc. +** 4. +** 13 May 97 ali made performance improvements in rint, rinttol, roundtol +** and trunc by folding some of the taligent ideas into this +** implementation. nearbyint is faster than the one in taligent, +** rint is more elegant, but slower by %30 than the taligent one. +** 09 Apr 97 ali deleted modfl and deferred to AuxiliaryDD.c +** 15 Sep 94 ali Major overhaul and performance improvements of all functions. +** 20 Jul 94 PAF New faster version +** 16 Jul 93 ali Added the modfl function. +** 18 Feb 93 ali Changed the return value of fenv functions +** feclearexcept and feraiseexcept to their new +** NCEG X3J11.1/93-001 definitions. +** 16 Dec 92 JPO Removed __itrunc implementation to a +** separate file. +** 15 Dec 92 JPO Added __itrunc implementation and modified +** rinttol to include conversion from double +** to long int format. Modified roundtol to +** call __itrunc. +** 10 Dec 92 JPO Added modf (double) implementation. +** 04 Dec 92 JPO First created. +** +*******************************************************************************/ + +#include <limits.h> +#include <math.h> + +#define SET_INVALID 0x01000000UL + +typedef union + { + struct { +#if defined(__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 }}; + +/******************************************************************************* +* * +* 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 ); + } + diff --git a/libm/powerpc/scalb.c b/libm/powerpc/w_scalb.c index 98a2b7d3f..98a2b7d3f 100644 --- a/libm/powerpc/scalb.c +++ b/libm/powerpc/w_scalb.c |