summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
-rw-r--r--libm/powerpc/Makefile2
-rw-r--r--libm/powerpc/s_ceil.c109
-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.c47
-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.c157
-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