diff options
author | Peter S. Mazinger <ps.m@gmx.net> | 2005-09-22 00:51:06 +0000 |
---|---|---|
committer | Peter S. Mazinger <ps.m@gmx.net> | 2005-09-22 00:51:06 +0000 |
commit | 178be753686094a0a998ab898e1f13d96626fbc9 (patch) | |
tree | 94782c457501e0d48317753fc4e4e5196a673ab9 /libm/powerpc/s_modf.c | |
parent | f8d5244380826053b8c75b3c302d39bdd1f9a121 (diff) |
split out nearbyint, round, trunc from libm/powerpc/s_modf.c
Diffstat (limited to 'libm/powerpc/s_modf.c')
-rw-r--r-- | libm/powerpc/s_modf.c | 198 |
1 files changed, 2 insertions, 196 deletions
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 * |