diff options
Diffstat (limited to 'libm/powerpc')
| -rw-r--r-- | libm/powerpc/Makefile | 8 | ||||
| -rw-r--r-- | libm/powerpc/s_ceil.c | 16 | ||||
| -rw-r--r-- | libm/powerpc/s_copysign.c | 6 | ||||
| -rw-r--r-- | libm/powerpc/s_floor.c | 16 | ||||
| -rw-r--r-- | libm/powerpc/s_frexp.c | 8 | ||||
| -rw-r--r-- | libm/powerpc/s_ldexp.c | 12 | ||||
| -rw-r--r-- | libm/powerpc/s_logb.c | 26 | ||||
| -rw-r--r-- | libm/powerpc/s_modf.c | 138 | ||||
| -rw-r--r-- | libm/powerpc/s_rint.c | 44 | ||||
| -rw-r--r-- | libm/powerpc/w_scalb.c | 32 |
10 files changed, 153 insertions, 153 deletions
diff --git a/libm/powerpc/Makefile b/libm/powerpc/Makefile index 84e0b94c9..9e4a309da 100644 --- a/libm/powerpc/Makefile +++ b/libm/powerpc/Makefile @@ -6,14 +6,14 @@ # math library for Apple's MacOS X/Darwin math library, which was # itself swiped from FreeBSD. The original copyright information # is as follows: -# +# # Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. -# +# # Developed at SunPro, a Sun Microsystems, Inc. business. # Permission to use, copy, modify, and distribute this # software is freely granted, provided that this notice # is preserved. -# +# # It has been ported to work with uClibc and generally behave # by Erik Andersen <andersen@codepoet.org> # @@ -66,6 +66,6 @@ $(OBJ): Makefile tags: ctags -R -clean: +clean: $(RM) *.[oa] *~ core diff --git a/libm/powerpc/s_ceil.c b/libm/powerpc/s_ceil.c index 790dd7ba0..fd073de7b 100644 --- a/libm/powerpc/s_ceil.c +++ b/libm/powerpc/s_ceil.c @@ -52,24 +52,24 @@ double ceil ( double x ) 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 ) + + if ( xhi < 0x43300000ul ) /******************************************************************************* * Is |x| < 2.0^52? * *******************************************************************************/ { - if ( xhi < 0x3ff00000ul ) + if ( xhi < 0x3ff00000ul ) /******************************************************************************* * Is |x| < 1.0? * *******************************************************************************/ { if ( ( xhi | xInHex.words.lo ) == 0ul ) // zero x is exact case return ( x ); - else + else { // inexact case asm ("mffs %0" : "=f" (OldEnvironment.dbl)); OldEnvironment.words.lo |= 0x02000000ul; @@ -83,7 +83,7 @@ double ceil ( double x ) /******************************************************************************* * Is 1.0 < |x| < 2.0^52? * *******************************************************************************/ - if ( target ) + if ( target ) { y = ( x + twoTo52 ) - twoTo52; // round at binary pt. if ( y < x ) @@ -91,8 +91,8 @@ double ceil ( double x ) else return ( y ); } - - else + + else { y = ( x - twoTo52 ) + twoTo52; // round at binary pt. if ( y < x ) diff --git a/libm/powerpc/s_copysign.c b/libm/powerpc/s_copysign.c index 18e91998a..d3fb5f40b 100644 --- a/libm/powerpc/s_copysign.c +++ b/libm/powerpc/s_copysign.c @@ -45,12 +45,12 @@ double copysign ( double arg2, double arg1 ) /******************************************************************************* * No need to flush NaNs out. * *******************************************************************************/ - + x.dbl = arg1; y.dbl = arg2; - + y.hex.high = y.hex.high & 0x7FFFFFFF; y.hex.high = ( y.hex.high | ( x.hex.high & dSgnMask ) ); - + return y.dbl; } diff --git a/libm/powerpc/s_floor.c b/libm/powerpc/s_floor.c index f7b4534bb..94677b4d2 100644 --- a/libm/powerpc/s_floor.c +++ b/libm/powerpc/s_floor.c @@ -52,24 +52,24 @@ double floor ( double x ) register double y; register unsigned long int xhi; register long 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 ) + + if ( xhi < 0x43300000ul ) /******************************************************************************* * Is |x| < 2.0^52? * *******************************************************************************/ { - if ( xhi < 0x3ff00000ul ) + if ( xhi < 0x3ff00000ul ) /******************************************************************************* * Is |x| < 1.0? * *******************************************************************************/ { if ( ( xhi | xInHex.words.lo ) == 0ul ) // zero x is exact case return ( x ); - else + else { // inexact case asm ("mffs %0" : "=f" (OldEnvironment.dbl)); OldEnvironment.words.lo |= 0x02000000ul; @@ -83,7 +83,7 @@ double floor ( double x ) /******************************************************************************* * Is 1.0 < |x| < 2.0^52? * *******************************************************************************/ - if ( target ) + if ( target ) { y = ( x + twoTo52 ) - twoTo52; // round at binary pt. if ( y > x ) @@ -91,8 +91,8 @@ double floor ( double x ) else return ( y ); } - - else + + else { y = ( x - twoTo52 ) + twoTo52; // round at binary pt. if ( y > x ) diff --git a/libm/powerpc/s_frexp.c b/libm/powerpc/s_frexp.c index 0ca7c1514..9909f2ce7 100644 --- a/libm/powerpc/s_frexp.c +++ b/libm/powerpc/s_frexp.c @@ -31,12 +31,12 @@ typedef union unsigned long int hi; unsigned long int lo; #else - unsigned long int lo; - unsigned long int hi; + unsigned long int lo; + unsigned long int hi; #endif } words; double dbl; - } DblInHex; + } DblInHex; double frexp ( double value, int *eptr ) { @@ -49,7 +49,7 @@ double frexp ( double value, int *eptr ) *eptr = 0; if ( valueHead >= 0x7ff00000 || ( valueHead | argument.words.lo ) == 0 ) return value; // 0, inf, or NaN - + if ( valueHead < 0x00100000 ) { // denorm argument.dbl = two54 * value; diff --git a/libm/powerpc/s_ldexp.c b/libm/powerpc/s_ldexp.c index 0e9a66611..ce9ec8b1b 100644 --- a/libm/powerpc/s_ldexp.c +++ b/libm/powerpc/s_ldexp.c @@ -29,18 +29,18 @@ typedef union unsigned long int hi; unsigned long int lo; #else - unsigned long int lo; - unsigned long int hi; + unsigned long int lo; + unsigned long int hi; #endif } words; double dbl; - } DblInHex; + } DblInHex; -double ldexp ( double value, int exp ) +double ldexp ( double value, int exp ) { - if ( exp > SHRT_MAX ) + if ( exp > SHRT_MAX ) exp = SHRT_MAX; - else if ( exp < -SHRT_MAX ) + else if ( exp < -SHRT_MAX ) exp = -SHRT_MAX; return scalb ( value, exp ); } diff --git a/libm/powerpc/s_logb.c b/libm/powerpc/s_logb.c index 3a410ba5c..23c7270f9 100644 --- a/libm/powerpc/s_logb.c +++ b/libm/powerpc/s_logb.c @@ -32,8 +32,8 @@ * Standard 754. * *******************************************************************************/ -typedef union - { +typedef union + { struct { #if defined(__BIG_ENDIAN__) unsigned long int hi; @@ -62,38 +62,38 @@ double logb ( double x ) { DblInHex xInHex; long int shiftedExp; - + xInHex.dbl = x; shiftedExp = ( xInHex.words.hi & 0x7ff00000UL ) >> 20; - - if ( shiftedExp == 2047 ) + + if ( shiftedExp == 2047 ) { // NaN or INF if ( ( ( xInHex.words.hi & signMask ) == 0 ) || ( x != x ) ) return x; // NaN or +INF return x else return -x; // -INF returns +INF } - + if ( shiftedExp != 0 ) // normal number shiftedExp -= 1023; // unbias exponent - - else if ( x == 0.0 ) + + else if ( x == 0.0 ) { // zero xInHex.words.hi = 0x0UL; // return -infinity return ( minusInf.dbl ); } - - else + + else { // subnormal number xInHex.dbl *= twoTo52; // scale up shiftedExp = ( xInHex.words.hi & 0x7ff00000UL ) >> 20; shiftedExp -= 1075; // unbias exponent } - + if ( shiftedExp == 0 ) // zero result return ( 0.0 ); - - else + + else { // nonzero result xInHex.dbl = klTod; xInHex.words.lo += shiftedExp; diff --git a/libm/powerpc/s_modf.c b/libm/powerpc/s_modf.c index 66c0655f5..403c54b79 100644 --- a/libm/powerpc/s_modf.c +++ b/libm/powerpc/s_modf.c @@ -1,18 +1,18 @@ /******************************************************************************* ** 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 @@ -21,7 +21,7 @@ ** 1. removed double_t, put in double for now. ** 2. removed iclass from nearbyint. ** 3. removed wrong comments intrunc. -** 4. +** 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, @@ -33,7 +33,7 @@ ** 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 +** 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 @@ -41,7 +41,7 @@ ** call __itrunc. ** 10 Dec 92 JPO Added modf (double) implementation. ** 04 Dec 92 JPO First created. -** +** *******************************************************************************/ #include <limits.h> @@ -81,14 +81,14 @@ static const DblInHex TOWARDZERO = {{ 0x00000000, 0x00000001 }}; * 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 */ @@ -99,9 +99,9 @@ double nearbyint ( double x ) y = copysign ( y, x ); // restore old flags asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment )); - return ( y ); + return ( y ); } - + /******************************************************************************* * * * The function rinttol converts its double argument to integral value * @@ -120,38 +120,38 @@ long int rinttol ( double x ) DblInHex argument, OldEnvironment; unsigned long int xHead; register long int target; - + argument.dbl = x; target = ( argument.words.hi < signMask ); // flag positive sign xHead = argument.words.hi & 0x7ffffffful; // high 32 bits of x - - if ( target ) + + if ( target ) /******************************************************************************* * Sign of x is positive. * *******************************************************************************/ { - if ( xHead < 0x41dffffful ) + if ( xHead < 0x41dffffful ) { // x is safely in long range y = ( x + twoTo52 ) - twoTo52; // round at binary point argument.dbl = y + doubleToLong; // force result into argument.words.lo return ( ( long ) argument.words.lo ); } - + asm ("mffs %0" : "=f" (OldEnvironment.dbl)); // get environment - - if ( xHead > 0x41dffffful ) + + if ( xHead > 0x41dffffful ) { // x is safely out of long range OldEnvironment.words.lo |= SET_INVALID; asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl )); return ( LONG_MAX ); } - + /******************************************************************************* * x > 0.0 and may or may not be out of range of long. * *******************************************************************************/ y = ( x + twoTo52 ) - twoTo52; // do rounding - if ( y > ( double ) LONG_MAX ) + if ( y > ( double ) LONG_MAX ) { // out of range of long OldEnvironment.words.lo |= SET_INVALID; asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl )); @@ -160,32 +160,32 @@ long int rinttol ( double x ) argument.dbl = y + doubleToLong; // in range return ( ( long ) argument.words.lo ); // return result & flags } - + /******************************************************************************* * Sign of x is negative. * *******************************************************************************/ - if ( xHead < 0x41e00000ul ) + if ( xHead < 0x41e00000ul ) { // x is safely in long range y = ( x - twoTo52 ) + twoTo52; argument.dbl = y + doubleToLong; return ( ( long ) argument.words.lo ); } - + asm ("mffs %0" : "=f" (OldEnvironment.dbl)); // get environment - - if ( xHead > 0x41e00000ul ) + + if ( xHead > 0x41e00000ul ) { // x is safely out of long range OldEnvironment.words.lo |= SET_INVALID; asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl )); return ( LONG_MIN ); } - + /******************************************************************************* * x < 0.0 and may or may not be out of range of long. * *******************************************************************************/ y = ( x - twoTo52 ) + twoTo52; // do rounding - if ( y < ( double ) LONG_MIN ) + if ( y < ( double ) LONG_MIN ) { // out of range of long OldEnvironment.words.lo |= SET_INVALID; asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl )); @@ -204,30 +204,30 @@ long int rinttol ( double x ) * 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 ) + + if ( xHead < 0x43300000ul ) /******************************************************************************* * Is |x| < 2.0^52? * *******************************************************************************/ { - if ( xHead < 0x3ff00000ul ) + if ( xHead < 0x3ff00000ul ) /******************************************************************************* * Is |x| < 1.0? * *******************************************************************************/ { asm ("mffs %0" : "=f" (OldEnvironment.dbl)); // get environment - if ( xHead < 0x3fe00000ul ) + if ( xHead < 0x3fe00000ul ) /******************************************************************************* * Is |x| < 0.5? * *******************************************************************************/ @@ -235,7 +235,7 @@ double round ( double x ) if ( ( xHead | argument.words.lo ) != 0ul ) OldEnvironment.words.lo |= 0x02000000ul; asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl )); - if ( target ) + if ( target ) return ( 0.0 ); else return ( -0.0 ); @@ -253,7 +253,7 @@ double round ( double x ) /******************************************************************************* * Is 1.0 < |x| < 2.0^52? * *******************************************************************************/ - if ( target ) + if ( target ) { // positive x y = ( x + twoTo52 ) - twoTo52; // round at binary point if ( y == x ) // exact case @@ -265,11 +265,11 @@ double round ( double x ) else return ( y ); } - + /******************************************************************************* * Is x < 0? * *******************************************************************************/ - else + else { y = ( x - twoTo52 ) + twoTo52; // round at binary point if ( y == x ) @@ -302,19 +302,19 @@ double round ( double x ) *******************************************************************************/ long int roundtol ( double x ) - { + { register double y, z; DblInHex argument, OldEnvironment; register unsigned long int xhi; register long int target; const DblInHex kTZ = {{ 0x0, 0x1 }}; const DblInHex kUP = {{ 0x0, 0x2 }}; - + argument.dbl = x; xhi = argument.words.hi & 0x7ffffffful; // high 32 bits of x target = ( argument.words.hi < signMask ); // flag positive sign - - if ( xhi > 0x41e00000ul ) + + if ( xhi > 0x41e00000ul ) /******************************************************************************* * Is x is out of long range or NaN? * *******************************************************************************/ @@ -327,19 +327,19 @@ long int roundtol ( double x ) else return ( LONG_MIN ); } - - if ( target ) + + if ( target ) /******************************************************************************* * Is sign of x is "+"? * *******************************************************************************/ { - if ( x < 2147483647.5 ) + if ( x < 2147483647.5 ) /******************************************************************************* * x is in the range of a long. * *******************************************************************************/ { y = ( x + doubleToLong ) - doubleToLong; // round at binary point - if ( y != x ) + if ( y != x ) { // inexact case asm ("mffs %0" : "=f" (OldEnvironment.dbl)); // save environment asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( kTZ.dbl )); // truncate rounding @@ -348,7 +348,7 @@ long int roundtol ( double x ) asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl )); return ( ( long ) argument.words.lo ); } - + argument.dbl = y + doubleToLong; // force result into argument.words.lo return ( ( long ) argument.words.lo ); // return long result } @@ -363,13 +363,13 @@ long int roundtol ( double x ) /******************************************************************************* * x < 0.0 and may or may not be out of the range of a long. * *******************************************************************************/ - if ( x > -2147483648.5 ) + if ( x > -2147483648.5 ) /******************************************************************************* * x is in the range of a long. * *******************************************************************************/ { y = ( x + doubleToLong ) - doubleToLong; // round at binary point - if ( y != x ) + if ( y != x ) { // inexact case asm ("mffs %0" : "=f" (OldEnvironment.dbl)); // save environment asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( kUP.dbl )); // round up @@ -378,7 +378,7 @@ long int roundtol ( double x ) asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl )); return ( ( long ) argument.words.lo ); } - + argument.dbl = y + doubleToLong; return ( ( long ) argument.words.lo ); // return long result } @@ -398,29 +398,29 @@ long int roundtol ( double x ) * 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 ) + + if ( xhi < 0x43300000ul ) /******************************************************************************* * Is |x| < 2.0^53? * *******************************************************************************/ { - if ( xhi < 0x3ff00000ul ) + if ( xhi < 0x3ff00000ul ) /******************************************************************************* * Is |x| < 1.0? * *******************************************************************************/ { - if ( ( xhi | argument.words.lo ) != 0ul ) + if ( ( xhi | argument.words.lo ) != 0ul ) { // raise deserved INEXACT asm ("mffs %0" : "=f" (OldEnvironment.dbl)); OldEnvironment.words.lo |= 0x02000000ul; @@ -434,7 +434,7 @@ double trunc ( double x ) /******************************************************************************* * Is 1.0 < |x| < 2.0^52? * *******************************************************************************/ - if ( target ) + if ( target ) { y = ( x + twoTo52 ) - twoTo52; // round at binary point if ( y > x ) @@ -442,8 +442,8 @@ double trunc ( double x ) else return ( y ); } - - else + + else { y = ( x - twoTo52 ) + twoTo52; // round at binary point. if ( y < x ) @@ -470,7 +470,7 @@ double trunc ( double x ) *******************************************************************************/ /******************************************************************************* -* modf is the double implementation. * +* modf is the double implementation. * *******************************************************************************/ double modf ( double x, double *iptr ) @@ -478,19 +478,19 @@ double modf ( double x, double *iptr ) register double OldEnvironment, xtrunc; register unsigned long int xHead, signBit; DblInHex argument; - + argument.dbl = x; xHead = argument.words.hi & 0x7ffffffful; // |x| high bit pattern signBit = ( argument.words.hi & 0x80000000ul ); // isolate sign bit if (xHead == 0x7ff81fe0) signBit = signBit | 0; - - if ( xHead < 0x43300000ul ) + + if ( xHead < 0x43300000ul ) /******************************************************************************* * Is |x| < 2.0^53? * *******************************************************************************/ { - if ( xHead < 0x3ff00000ul ) + if ( xHead < 0x3ff00000ul ) /******************************************************************************* * Is |x| < 1.0? * *******************************************************************************/ @@ -515,18 +515,18 @@ double modf ( double x, double *iptr ) *iptr = xtrunc; // store integral part if ( x != xtrunc ) // nonzero fraction return ( x - xtrunc ); - else + else { // zero with x's sign argument.words.hi = signBit; argument.words.lo = 0ul; return ( argument.dbl ); } } - + *iptr = x; // x is integral or NaN if ( x != x ) // NaN is returned return x; - else + else { // zero with x's sign argument.words.hi = signBit; argument.words.lo = 0ul; diff --git a/libm/powerpc/s_rint.c b/libm/powerpc/s_rint.c index 47a79ae5e..72c4834d0 100644 --- a/libm/powerpc/s_rint.c +++ b/libm/powerpc/s_rint.c @@ -1,18 +1,18 @@ /******************************************************************************* ** 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 @@ -21,7 +21,7 @@ ** 1. removed double_t, put in double for now. ** 2. removed iclass from nearbyint. ** 3. removed wrong comments intrunc. -** 4. +** 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, @@ -33,7 +33,7 @@ ** 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 +** 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 @@ -41,7 +41,7 @@ ** call __itrunc. ** 10 Dec 92 JPO Added modf (double) implementation. ** 04 Dec 92 JPO First created. -** +** *******************************************************************************/ #include <limits.h> @@ -73,7 +73,7 @@ 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 * +* double format. This function signals inexact if an ordered return * * value is not equal to the operand. * * * ******************************************************************************** @@ -89,16 +89,16 @@ static const DblInHex TOWARDZERO = {{ 0x00000000, 0x00000001 }}; *double rint ( double x ) * { * double y; -* +* * y = twoTo52.fval; -* -* if ( fabs ( x ) >= y ) // huge case is exact +* +* 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 +* 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 ); +* return ( y ); * } ******************************************************************************** * Now a bit twidling version that is about %30 faster. * @@ -110,17 +110,17 @@ double rint ( double x ) 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 ) + + if ( xHead < 0x43300000ul ) /******************************************************************************* * Is |x| < 2.0^52? * *******************************************************************************/ { - if ( xHead < 0x3ff00000ul ) + if ( xHead < 0x3ff00000ul ) /******************************************************************************* * Is |x| < 1.0? * *******************************************************************************/ @@ -129,7 +129,7 @@ double rint ( double x ) y = ( x + twoTo52 ) - twoTo52; // round at binary point else y = ( x - twoTo52 ) + twoTo52; // round at binary point - if ( y == 0.0 ) + if ( y == 0.0 ) { // fix sign of zero result if ( target ) return ( 0.0 ); @@ -138,7 +138,7 @@ double rint ( double x ) } return y; } - + /******************************************************************************* * Is 1.0 < |x| < 2.0^52? * *******************************************************************************/ @@ -148,7 +148,7 @@ double rint ( double x ) else return ( ( x - twoTo52 ) + twoTo52 ); } - + /******************************************************************************* * |x| >= 2.0^52 or x is a NaN. * *******************************************************************************/ diff --git a/libm/powerpc/w_scalb.c b/libm/powerpc/w_scalb.c index 98a2b7d3f..c93c74b68 100644 --- a/libm/powerpc/w_scalb.c +++ b/libm/powerpc/w_scalb.c @@ -1,26 +1,26 @@ /*********************************************************************** ** File: scalb.c -** +** ** Contains: C source code for implementations of floating-point ** scalb functions defined in header <fp.h>. In ** particular, this file contains implementations of ** functions scalb and scalbl for double and long double ** formats on PowerPC platforms. -** +** ** Written by: Jon Okada, SANEitation Engineer, ext. 4-4838 -** +** ** Copyright: © 1992 by Apple Computer, Inc., all rights reserved -** +** ** Change History ( most recent first ): ** ** 28 May 97 ali made an speed improvement for large n, ** removed scalbl. ** 12 Dec 92 JPO First created. -** +** ***********************************************************************/ -typedef union - { +typedef union + { struct { #if defined(__BIG_ENDIAN__) unsigned long int hi; @@ -41,35 +41,35 @@ static const double twoToM1022 = 2.225073858507201383e-308; // 0x1p-1022 double scalb( double x, long int n ) returns its argument x scaled by the factor 2^m. NaNs, signed zeros, and infinities are propagated by this function regardless of the value of n. - + Exceptions: OVERFLOW/INEXACT or UNDERFLOW inexact may occur; INVALID for signaling NaN inputs ( quiet NaN returned ). - + Calls: none. ***********************************************************************/ double scalb ( double x, int n ) { DblInHex xInHex; - + xInHex.words.lo = 0UL; // init. low half of xInHex - - if ( n > 1023 ) + + if ( n > 1023 ) { // large positive scaling if ( n > 2097 ) // huge scaling return ( ( x * twoTo1023 ) * twoTo1023 ) * twoTo1023; - while ( n > 1023 ) + while ( n > 1023 ) { // scale reduction loop x *= twoTo1023; // scale x by 2^1023 n -= 1023; // reduce n by 1023 } } - - else if ( n < -1022 ) + + else if ( n < -1022 ) { // large negative scaling if ( n < -2098 ) // huge negative scaling return ( ( x * twoToM1022 ) * twoToM1022 ) * twoToM1022; - while ( n < -1022 ) + while ( n < -1022 ) { // scale reduction loop x *= twoToM1022; // scale x by 2^( -1022 ) n += 1022; // incr n by 1022 |