From 178be753686094a0a998ab898e1f13d96626fbc9 Mon Sep 17 00:00:00 2001 From: "Peter S. Mazinger" Date: Thu, 22 Sep 2005 00:51:06 +0000 Subject: split out nearbyint, round, trunc from libm/powerpc/s_modf.c --- libm/powerpc/Makefile | 3 +- libm/powerpc/s_modf.c | 198 +-------------------------------------------- libm/powerpc/s_nearbyint.c | 36 +++++++++ libm/powerpc/s_round.c | 112 +++++++++++++++++++++++++ libm/powerpc/s_trunc.c | 86 ++++++++++++++++++++ 5 files changed, 238 insertions(+), 197 deletions(-) create mode 100644 libm/powerpc/s_nearbyint.c create mode 100644 libm/powerpc/s_round.c create mode 100644 libm/powerpc/s_trunc.c (limited to 'libm/powerpc') 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 . 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,42 +65,9 @@ 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 * @@ -195,99 +161,6 @@ long int rinttol ( double x ) return ( ( long ) argument.words.lo ); // return result & flags } -/******************************************************************************* -* * -* 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 * @@ -391,73 +264,6 @@ long int roundtol ( double x ) return ( LONG_MIN ); // return pinned result } -/******************************************************************************* -* * -* 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 * 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 +#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; + +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 +#include + +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 +#include + +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 ); + } -- cgit v1.2.3