summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--libm/powerpc/Makefile3
-rw-r--r--libm/powerpc/s_modf.c198
-rw-r--r--libm/powerpc/s_nearbyint.c36
-rw-r--r--libm/powerpc/s_round.c112
-rw-r--r--libm/powerpc/s_trunc.c86
5 files changed, 238 insertions, 197 deletions
diff --git a/libm/powerpc/Makefile b/libm/powerpc/Makefile
index 9e4a309da..342c182ba 100644
--- a/libm/powerpc/Makefile
+++ b/libm/powerpc/Makefile
@@ -38,7 +38,8 @@ LIBM=../libm.a
CFLAGS+=-D_IEEE_LIBM -D_ISOC99_SOURCE -D_SVID_SOURCE
ifeq ($(strip $(DO_C99_MATH)),y)
-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
+CSRC = s_ceil.c s_copysign.c s_floor.c s_frexp.c s_ldexp.c s_logb.c s_modf.c \
+ s_nearbyint.c s_rint.c s_round.c s_trunc.c w_scalb.c
else
CSRC =
endif
diff --git a/libm/powerpc/s_modf.c b/libm/powerpc/s_modf.c
index 403c54b79..f4344bda8 100644
--- a/libm/powerpc/s_modf.c
+++ b/libm/powerpc/s_modf.c
@@ -4,9 +4,8 @@
** 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.
+** 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
@@ -66,44 +65,11 @@ typedef union
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 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. *
-* *
-*******************************************************************************/
-
-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 );
- }
-
-/*******************************************************************************
-* *
* 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 *
@@ -197,99 +163,6 @@ long int rinttol ( double x )
/*******************************************************************************
* *
-* 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. *
-* *
-*******************************************************************************/
-
-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 );
- }
-
-/*******************************************************************************
-* *
* 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 *
@@ -392,73 +265,6 @@ long int roundtol ( double x )
}
/*******************************************************************************
-* *
-* 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. *
-* *
-*******************************************************************************/
-
-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 );
- }
-
-/*******************************************************************************
* 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 *
diff --git a/libm/powerpc/s_nearbyint.c b/libm/powerpc/s_nearbyint.c
new file mode 100644
index 000000000..f2d7ded35
--- /dev/null
+++ b/libm/powerpc/s_nearbyint.c
@@ -0,0 +1,36 @@
+#include <limits.h>
+#include <math.h>
+
+/*******************************************************************************
+* *
+* The function nearbyint rounds its double argument to integral value *
+* according to the current rounding direction and returns the result in *
+* double format. This function does not signal inexact. *
+* *
+********************************************************************************
+* *
+* This function calls fabs and copysign. *
+* *
+*******************************************************************************/
+
+static const double twoTo52 = 4503599627370496.0;
+
+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 );
+ }
diff --git a/libm/powerpc/s_round.c b/libm/powerpc/s_round.c
new file mode 100644
index 000000000..81f4d0fec
--- /dev/null
+++ b/libm/powerpc/s_round.c
@@ -0,0 +1,112 @@
+#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;
+
+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. *
+* *
+*******************************************************************************/
+
+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 );
+ }
diff --git a/libm/powerpc/s_trunc.c b/libm/powerpc/s_trunc.c
new file mode 100644
index 000000000..4b61355ea
--- /dev/null
+++ b/libm/powerpc/s_trunc.c
@@ -0,0 +1,86 @@
+#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;
+
+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. *
+* *
+*******************************************************************************/
+
+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 );
+ }
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-01 06:07:07 +0000