diff options
Diffstat (limited to 'libm/powerpc')
| -rw-r--r-- | libm/powerpc/Makefile | 60 | ||||
| -rw-r--r-- | libm/powerpc/ceilfloor.c | 174 | ||||
| -rw-r--r-- | libm/powerpc/frexpldexp.c | 72 | ||||
| -rw-r--r-- | libm/powerpc/logb.c | 103 | ||||
| -rw-r--r-- | libm/powerpc/rndint.c | 620 | ||||
| -rw-r--r-- | libm/powerpc/scalb.c | 85 | ||||
| -rw-r--r-- | libm/powerpc/sign.c | 56 |
7 files changed, 1170 insertions, 0 deletions
diff --git a/libm/powerpc/Makefile b/libm/powerpc/Makefile new file mode 100644 index 000000000..7464ee0a8 --- /dev/null +++ b/libm/powerpc/Makefile @@ -0,0 +1,60 @@ +# Makefile for uClibc's math library +# Copyright (C) 2001 by Lineo, inc. +# +# This math library is derived primarily from the Cephes Math Library, +# copyright by Stephen L. Moshier <moshier@world.std.com> +# +# This program is free software; you can redistribute it and/or modify it under +# the terms of the GNU Library General Public License as published by the Free +# Software Foundation; either version 2 of the License, or (at your option) any +# later version. +# +# This program is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +# FOR A PARTICULAR PURPOSE. See the GNU Library General Public License for more +# details. +# +# You should have received a copy of the GNU Library General Public License +# along with this program; if not, write to the Free Software Foundation, Inc., +# 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +# + +TOPDIR=../../ +include $(TOPDIR)Rules.mak + +LIBM=../libm.a +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 +else +CSRC = +endif +COBJS=$(patsubst %.c,%.o, $(CSRC)) +OBJS=$(COBJS) + + +ifneq ($(strip $(HAS_FLOATING_POINT)),true) +all: clean +else +all: $(OBJS) $(LIBM) +endif + +$(LIBM): ar-target + +ar-target: $(OBJS) + $(AR) $(ARFLAGS) $(LIBM) $(OBJS) + +$(COBJS): %.o : %.c + $(TARGET_CC) $(TARGET_CFLAGS) -c $< -o $@ + $(STRIPTOOL) -x -R .note -R .comment $*.o + +$(OBJ): Makefile + +tags: + ctags -R + +clean: + rm -f *.[oa] *~ core $(LIBM_SHARED)* $(LIBM_SHARED_FULLNAME)* + diff --git a/libm/powerpc/ceilfloor.c b/libm/powerpc/ceilfloor.c new file mode 100644 index 000000000..d5947ab50 --- /dev/null +++ b/libm/powerpc/ceilfloor.c @@ -0,0 +1,174 @@ +/******************************************************************************* +* * +* 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 ); + } + +/******************************************************************************* +* Floor(x) returns the largest integer not greater than x. * +*******************************************************************************/ + +double floor ( double x ) + { + DblInHex xInHex,OldEnvironment; + 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 ) +/******************************************************************************* +* 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 ( 0.0 ); + else + return ( -1.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/frexpldexp.c b/libm/powerpc/frexpldexp.c new file mode 100644 index 000000000..6b759f91b --- /dev/null +++ b/libm/powerpc/frexpldexp.c @@ -0,0 +1,72 @@ +/******************************************************************************* +* * +* 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> + +static const double two54 = 1.80143985094819840000e+16; /* 0x43500000, 0x00000000 */ + +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 ); + } + +double frexp ( double value, int *eptr ) + { + DblInHex argument; + unsigned long int valueHead; + + argument.dbl = value; + valueHead = argument.words.hi & 0x7fffffffUL; // valueHead <- |x| + + *eptr = 0; + if ( valueHead >= 0x7ff00000 || ( valueHead | argument.words.lo ) == 0 ) + return value; // 0, inf, or NaN + + if ( valueHead < 0x00100000 ) + { // denorm + argument.dbl = two54 * value; + valueHead = argument.words.hi &0x7fffffff; + *eptr = -54; + } + *eptr += ( valueHead >> 20 ) - 1022; + argument.words.hi = ( argument.words.hi & 0x800fffff ) | 0x3fe00000; + return argument.dbl; + } + diff --git a/libm/powerpc/logb.c b/libm/powerpc/logb.c new file mode 100644 index 000000000..3a410ba5c --- /dev/null +++ b/libm/powerpc/logb.c @@ -0,0 +1,103 @@ +/******************************************************************************* +* * +* File logb.c, * +* Functions logb. * +* Implementation of logb for the PowerPC. * +* * +* Copyright © 1991 Apple Computer, Inc. All rights reserved. * +* * +* Written by Ali Sazegari, started on June 1991, * +* * +* August 26 1991: removed CFront Version 1.1d17 warnings. * +* August 27 1991: no errors reported by the test suite. * +* November 11 1991: changed CLASSEXTENDED to the macro CLASSIFY and * +* + or - infinity to constants. * +* November 18 1991: changed the macro CLASSIFY to CLASSEXTENDEDint to * +* improve performance. * +* February 07 1992: changed bit operations to macros ( object size is * +* unchanged ). * +* September24 1992: took the "#include support.h" out. * +* December 03 1992: first rs/6000 port. * +* August 30 1992: set the divide by zero for the zero argument case. * +* October 05 1993: corrected the environment. * +* October 17 1994: replaced all environmental functions with __setflm. * +* May 28 1997: made speed improvements. * +* April 30 2001: forst mac os x port using gcc. * +* * +******************************************************************************** +* The C math library offers a similar function called "frexp". It is * +* different in details from logb, but similar in spirit. This current * +* implementation of logb follows the recommendation in IEEE Standard 854 * +* which is different in its handling of denormalized numbers from the IEEE * +* Standard 754. * +*******************************************************************************/ + +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 double twoTo52 = 4.50359962737049600e15; // 0x1p52 +static const double klTod = 4503601774854144.0; // 0x1.000008p52 +static const unsigned long int signMask = 0x80000000ul; +static const DblInHex minusInf = {{ 0xFFF00000, 0x00000000 }}; + + +/******************************************************************************* +******************************************************************************** +* L O G B * +******************************************************************************** +*******************************************************************************/ + +double logb ( double x ) + { + DblInHex xInHex; + long int shiftedExp; + + xInHex.dbl = x; + shiftedExp = ( xInHex.words.hi & 0x7ff00000UL ) >> 20; + + 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 ) + { // zero + xInHex.words.hi = 0x0UL; // return -infinity + return ( minusInf.dbl ); + } + + 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 + { // nonzero result + xInHex.dbl = klTod; + xInHex.words.lo += shiftedExp; + return ( xInHex.dbl - klTod ); + } + } + diff --git a/libm/powerpc/rndint.c b/libm/powerpc/rndint.c new file mode 100644 index 000000000..1b7a9ec81 --- /dev/null +++ b/libm/powerpc/rndint.c @@ -0,0 +1,620 @@ +/******************************************************************************* +** 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 ); + } + +/******************************************************************************* +* * +* 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 * +* NaN or the rounded intermediate result is out of range of the * +* destination long int format, and it delivers an unspecified result in * +* this case. This function signals inexact if the rounded result is * +* within range of the long int format but unequal to the operand. * +* * +*******************************************************************************/ + +long int rinttol ( double x ) + { + register double y; + 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 ) +/******************************************************************************* +* Sign of x is positive. * +*******************************************************************************/ + { + 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 ) + { // 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 ) + { // out of range of long + OldEnvironment.words.lo |= SET_INVALID; + asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl )); + return ( LONG_MAX ); + } + argument.dbl = y + doubleToLong; // in range + return ( ( long ) argument.words.lo ); // return result & flags + } + +/******************************************************************************* +* Sign of x is negative. * +*******************************************************************************/ + 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 ) + { // 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 ) + { // out of range of long + OldEnvironment.words.lo |= SET_INVALID; + asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl )); + return ( LONG_MIN ); + } + argument.dbl = y + doubleToLong; // in range + 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 * +* according to the "add half to the magnitude and chop" rounding mode of * +* Pascal's Round function and FORTRAN's NINT function. This conversion * +* signals invalid if the argument is a NaN or the rounded intermediate * +* result is out of range of the destination long int format, and it * +* delivers an unspecified result in this case. This function signals * +* inexact if the rounded result is within range of the long int format but * +* unequal to the operand. * +* * +*******************************************************************************/ + +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 ) +/******************************************************************************* +* Is x is out of long range or NaN? * +*******************************************************************************/ + { + asm ("mffs %0" : "=f" (OldEnvironment.dbl)); // get environment + OldEnvironment.words.lo |= SET_INVALID; + asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl )); + if ( target ) // pin result + return ( LONG_MAX ); + else + return ( LONG_MIN ); + } + + if ( target ) +/******************************************************************************* +* Is sign of x is "+"? * +*******************************************************************************/ + { + if ( x < 2147483647.5 ) +/******************************************************************************* +* x is in the range of a long. * +*******************************************************************************/ + { + y = ( x + doubleToLong ) - doubleToLong; // round at binary point + if ( y != x ) + { // inexact case + asm ("mffs %0" : "=f" (OldEnvironment.dbl)); // save environment + asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( kTZ.dbl )); // truncate rounding + z = x + 0.5; // truncate x + 0.5 + argument.dbl = z + doubleToLong; + 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 + } +/******************************************************************************* +* Rounded positive x is out of the range of a long. * +*******************************************************************************/ + asm ("mffs %0" : "=f" (OldEnvironment.dbl)); + OldEnvironment.words.lo |= SET_INVALID; + asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl )); + return ( LONG_MAX ); // return pinned result + } +/******************************************************************************* +* x < 0.0 and may or may not be out of the range of a long. * +*******************************************************************************/ + if ( x > -2147483648.5 ) +/******************************************************************************* +* x is in the range of a long. * +*******************************************************************************/ + { + y = ( x + doubleToLong ) - doubleToLong; // round at binary point + if ( y != x ) + { // inexact case + asm ("mffs %0" : "=f" (OldEnvironment.dbl)); // save environment + asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( kUP.dbl )); // round up + z = x - 0.5; // truncate x - 0.5 + argument.dbl = z + doubleToLong; + 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 + } +/******************************************************************************* +* Rounded negative x is out of the range of a long. * +*******************************************************************************/ + asm ("mffs %0" : "=f" (OldEnvironment.dbl)); + OldEnvironment.words.lo |= SET_INVALID; + asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl )); + 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 * +* the integral part in floating-point format to the object pointed to by a * +* pointer argument. If the input argument is integral or infinite in * +* value, the return value is a zero with the sign of the input argument. * +* The modf family of functions raises no floating-point exceptions. older * +* implemenation set the INVALID flag due to signaling NaN input. * +* * +*******************************************************************************/ + +/******************************************************************************* +* modf is the double implementation. * +*******************************************************************************/ + +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 ) +/******************************************************************************* +* Is |x| < 2.0^53? * +*******************************************************************************/ + { + if ( xHead < 0x3ff00000ul ) +/******************************************************************************* +* Is |x| < 1.0? * +*******************************************************************************/ + { + argument.words.hi = signBit; // truncate to zero + argument.words.lo = 0ul; + *iptr = argument.dbl; + return ( x ); + } +/******************************************************************************* +* Is 1.0 < |x| < 2.0^52? * +*******************************************************************************/ + asm ("mffs %0" : "=f" (OldEnvironment)); // save environment + // round toward zero + asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( TOWARDZERO.dbl )); + if ( signBit == 0ul ) // truncate to integer + xtrunc = ( x + twoTo52 ) - twoTo52; + else + xtrunc = ( x - twoTo52 ) + twoTo52; + // restore caller's env + asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment )); + *iptr = xtrunc; // store integral part + if ( x != xtrunc ) // nonzero fraction + return ( x - xtrunc ); + 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 + { // zero with x's sign + argument.words.hi = signBit; + argument.words.lo = 0ul; + return ( argument.dbl ); + } + } diff --git a/libm/powerpc/scalb.c b/libm/powerpc/scalb.c new file mode 100644 index 000000000..98a2b7d3f --- /dev/null +++ b/libm/powerpc/scalb.c @@ -0,0 +1,85 @@ +/*********************************************************************** +** 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 + { + 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 double twoTo1023 = 8.988465674311579539e307; // 0x1p1023 +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 ) + { // large positive scaling + if ( n > 2097 ) // huge scaling + return ( ( x * twoTo1023 ) * twoTo1023 ) * twoTo1023; + while ( n > 1023 ) + { // scale reduction loop + x *= twoTo1023; // scale x by 2^1023 + n -= 1023; // reduce n by 1023 + } + } + + else if ( n < -1022 ) + { // large negative scaling + if ( n < -2098 ) // huge negative scaling + return ( ( x * twoToM1022 ) * twoToM1022 ) * twoToM1022; + while ( n < -1022 ) + { // scale reduction loop + x *= twoToM1022; // scale x by 2^( -1022 ) + n += 1022; // incr n by 1022 + } + } + +/******************************************************************************* +* -1022 <= n <= 1023; convert n to double scale factor. * +*******************************************************************************/ + + xInHex.words.hi = ( ( unsigned long ) ( n + 1023 ) ) << 20; + return ( x * xInHex.dbl ); + } diff --git a/libm/powerpc/sign.c b/libm/powerpc/sign.c new file mode 100644 index 000000000..0bf2f3829 --- /dev/null +++ b/libm/powerpc/sign.c @@ -0,0 +1,56 @@ +/******************************************************************************* +* * +* File sign.c, * +* Functions copysign and __signbitd. * +* For PowerPC based machines. * +* * +* Copyright © 1991, 2001 Apple Computer, Inc. All rights reserved. * +* * +* Written by Ali Sazegari, started on June 1991. * +* * +* August 26 1991: no CFront Version 1.1d17 warnings. * +* September 06 1991: passes the test suite with invalid raised on * +* signaling nans. sane rom code behaves the same. * +* September 24 1992: took the ̉#include support.hÓ out. * +* Dcember 02 1992: PowerPC port. * +* July 20 1994: __fabs added * +* July 21 1994: deleted unnecessary functions: neg, COPYSIGNnew, * +* and SIGNNUMnew. * +* April 11 2001: first port to os x using gcc. * +* removed fabs and deffered to gcc for direct * +* instruction generation. * +* * +*******************************************************************************/ + +#include "fp_private.h" + +/******************************************************************************* +* * +* Function copysign. * +* Implementation of copysign for the PowerPC. * +* * +******************************************************************************** +* Note: The order of the operands in this function is reversed from that * +* suggested in the IEEE standard 754. * +*******************************************************************************/ + +double copysign ( double arg2, double arg1 ) + { + union + { + dHexParts hex; + double dbl; + } x, y; + +/******************************************************************************* +* 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; + } |