summaryrefslogtreecommitdiff
path: root/libm/powerpc/classic
diff options
context:
space:
mode:
Diffstat (limited to 'libm/powerpc/classic')
-rw-r--r--libm/powerpc/classic/Makefile.arch18
-rw-r--r--libm/powerpc/classic/s_ceil.c113
-rw-r--r--libm/powerpc/classic/s_copysign.c59
-rw-r--r--libm/powerpc/classic/s_floor.c113
-rw-r--r--libm/powerpc/classic/s_frexp.c65
-rw-r--r--libm/powerpc/classic/s_ldexp.c49
-rw-r--r--libm/powerpc/classic/s_logb.c107
-rw-r--r--libm/powerpc/classic/s_modf.c344
-rw-r--r--libm/powerpc/classic/s_nearbyint.c38
-rw-r--r--libm/powerpc/classic/s_rint.c159
-rw-r--r--libm/powerpc/classic/s_round.c115
-rw-r--r--libm/powerpc/classic/s_trunc.c89
-rw-r--r--libm/powerpc/classic/w_scalb.c94
13 files changed, 1363 insertions, 0 deletions
diff --git a/libm/powerpc/classic/Makefile.arch b/libm/powerpc/classic/Makefile.arch
new file mode 100644
index 000000000..7c7600f80
--- /dev/null
+++ b/libm/powerpc/classic/Makefile.arch
@@ -0,0 +1,18 @@
+# Makefile for uClibc
+#
+# Copyright (C) 2000-2006 Erik Andersen <andersen@uclibc.org>
+#
+# Licensed under the LGPL v2.1, see the file COPYING.LIB in this tarball.
+#
+
+libm_ARCH_SRC:=$(wildcard $(libm_ARCH_DIR)/*.c)
+libm_ARCH_OBJ:=$(patsubst $(libm_ARCH_DIR)/%.c,$(libm_ARCH_OUT)/%.o,$(libm_ARCH_SRC))
+
+libm_ARCH_OBJS:=$(libm_ARCH_OBJ)
+
+ifeq ($(DOPIC),y)
+libm-a-y+=$(libm_ARCH_OBJS:.o=.os)
+else
+libm-a-y+=$(libm_ARCH_OBJS)
+endif
+libm-so-y+=$(libm_ARCH_OBJS:.o=.os)
diff --git a/libm/powerpc/classic/s_ceil.c b/libm/powerpc/classic/s_ceil.c
new file mode 100644
index 000000000..8db5ce537
--- /dev/null
+++ b/libm/powerpc/classic/s_ceil.c
@@ -0,0 +1,113 @@
+/*******************************************************************************
+* *
+* 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 *
+* *
+*******************************************************************************/
+
+#include <math.h>
+#include <endian.h>
+
+static const double twoTo52 = 4503599627370496.0;
+static const unsigned long signMask = 0x80000000ul;
+
+typedef union
+ {
+ struct {
+#if (__BYTE_ORDER == __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. *
+*******************************************************************************/
+
+libm_hidden_proto(ceil)
+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 );
+ }
+libm_hidden_def(ceil)
diff --git a/libm/powerpc/classic/s_copysign.c b/libm/powerpc/classic/s_copysign.c
new file mode 100644
index 000000000..1a988ac87
--- /dev/null
+++ b/libm/powerpc/classic/s_copysign.c
@@ -0,0 +1,59 @@
+/*******************************************************************************
+* *
+* 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 <math.h>
+#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. *
+*******************************************************************************/
+
+libm_hidden_proto(copysign)
+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;
+ }
+libm_hidden_def(copysign)
diff --git a/libm/powerpc/classic/s_floor.c b/libm/powerpc/classic/s_floor.c
new file mode 100644
index 000000000..ff3436707
--- /dev/null
+++ b/libm/powerpc/classic/s_floor.c
@@ -0,0 +1,113 @@
+/*******************************************************************************
+* *
+* 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 *
+* *
+*******************************************************************************/
+
+#include <math.h>
+#include <endian.h>
+
+static const double twoTo52 = 4503599627370496.0;
+static const unsigned long signMask = 0x80000000ul;
+
+typedef union
+ {
+ struct {
+#if (__BYTE_ORDER == __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. *
+*******************************************************************************/
+
+/*******************************************************************************
+* Floor(x) returns the largest integer not greater than x. *
+*******************************************************************************/
+
+libm_hidden_proto(floor)
+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 );
+ }
+libm_hidden_def(floor)
diff --git a/libm/powerpc/classic/s_frexp.c b/libm/powerpc/classic/s_frexp.c
new file mode 100644
index 000000000..001aaf708
--- /dev/null
+++ b/libm/powerpc/classic/s_frexp.c
@@ -0,0 +1,65 @@
+/*******************************************************************************
+* *
+* 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>
+#include <endian.h>
+
+static const double two54 = 1.80143985094819840000e+16; /* 0x43500000, 0x00000000 */
+
+typedef union
+ {
+ struct {
+#if (__BYTE_ORDER == __BIG_ENDIAN)
+ unsigned long int hi;
+ unsigned long int lo;
+#else
+ unsigned long int lo;
+ unsigned long int hi;
+#endif
+ } words;
+ double dbl;
+ } DblInHex;
+
+libm_hidden_proto(frexp)
+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;
+ }
+libm_hidden_def(frexp)
diff --git a/libm/powerpc/classic/s_ldexp.c b/libm/powerpc/classic/s_ldexp.c
new file mode 100644
index 000000000..10100d7c2
--- /dev/null
+++ b/libm/powerpc/classic/s_ldexp.c
@@ -0,0 +1,49 @@
+/*******************************************************************************
+* *
+* 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>
+#include <endian.h>
+
+typedef union
+ {
+ struct {
+#if (__BYTE_ORDER == __BIG_ENDIAN)
+ unsigned long int hi;
+ unsigned long int lo;
+#else
+ unsigned long int lo;
+ unsigned long int hi;
+#endif
+ } words;
+ double dbl;
+ } DblInHex;
+
+libm_hidden_proto(ldexp)
+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 );
+ }
+libm_hidden_def(ldexp)
diff --git a/libm/powerpc/classic/s_logb.c b/libm/powerpc/classic/s_logb.c
new file mode 100644
index 000000000..81daa412e
--- /dev/null
+++ b/libm/powerpc/classic/s_logb.c
@@ -0,0 +1,107 @@
+/*******************************************************************************
+* *
+* 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. *
+*******************************************************************************/
+
+#include <math.h>
+#include <endian.h>
+
+typedef union
+ {
+ struct {
+#if (__BYTE_ORDER == __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 *
+********************************************************************************
+*******************************************************************************/
+
+libm_hidden_proto(logb)
+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 );
+ }
+ }
+libm_hidden_def(logb)
diff --git a/libm/powerpc/classic/s_modf.c b/libm/powerpc/classic/s_modf.c
new file mode 100644
index 000000000..b9d69445d
--- /dev/null
+++ b/libm/powerpc/classic/s_modf.c
@@ -0,0 +1,344 @@
+/*******************************************************************************
+** 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 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
+**
+** 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>
+#include <endian.h>
+
+#define SET_INVALID 0x01000000UL
+
+typedef union
+ {
+ struct {
+#if (__BYTE_ORDER == __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 TOWARDZERO = {{ 0x00000000, 0x00000001 }};
+
+
+/*******************************************************************************
+* *
+* 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 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 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. *
+*******************************************************************************/
+
+libm_hidden_proto(modf)
+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 );
+ }
+ }
+libm_hidden_def(modf)
diff --git a/libm/powerpc/classic/s_nearbyint.c b/libm/powerpc/classic/s_nearbyint.c
new file mode 100644
index 000000000..068e21378
--- /dev/null
+++ b/libm/powerpc/classic/s_nearbyint.c
@@ -0,0 +1,38 @@
+#include <limits.h>
+#include <math.h>
+
+/*******************************************************************************
+* *
+* The function nearbyint rounds its double argument to integral value *
+* according to the current rounding direction and returns the result in *
+* double format. This function does not signal inexact. *
+* *
+********************************************************************************
+* *
+* This function calls fabs and copysign. *
+* *
+*******************************************************************************/
+
+static const double twoTo52 = 4503599627370496.0;
+
+libm_hidden_proto(nearbyint)
+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 );
+ }
+libm_hidden_def(nearbyint)
diff --git a/libm/powerpc/classic/s_rint.c b/libm/powerpc/classic/s_rint.c
new file mode 100644
index 000000000..dd2ae585e
--- /dev/null
+++ b/libm/powerpc/classic/s_rint.c
@@ -0,0 +1,159 @@
+/*******************************************************************************
+** 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>
+#include <endian.h>
+
+#define SET_INVALID 0x01000000UL
+
+typedef union
+ {
+ struct {
+#if (__BYTE_ORDER == __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 }};
+
+libm_hidden_proto(rint)
+/*******************************************************************************
+* *
+* 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 );
+ }
+libm_hidden_def(rint)
diff --git a/libm/powerpc/classic/s_round.c b/libm/powerpc/classic/s_round.c
new file mode 100644
index 000000000..62d5936d9
--- /dev/null
+++ b/libm/powerpc/classic/s_round.c
@@ -0,0 +1,115 @@
+#include <limits.h>
+#include <math.h>
+#include <endian.h>
+
+typedef union
+ {
+ struct {
+#if (__BYTE_ORDER == __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. *
+* *
+*******************************************************************************/
+
+libm_hidden_proto(round)
+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 );
+ }
+libm_hidden_def(round)
diff --git a/libm/powerpc/classic/s_trunc.c b/libm/powerpc/classic/s_trunc.c
new file mode 100644
index 000000000..f793992a7
--- /dev/null
+++ b/libm/powerpc/classic/s_trunc.c
@@ -0,0 +1,89 @@
+#include <limits.h>
+#include <math.h>
+#include <endian.h>
+
+typedef union
+ {
+ struct {
+#if (__BYTE_ORDER == __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. *
+* *
+*******************************************************************************/
+
+libm_hidden_proto(trunc)
+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 );
+ }
+libm_hidden_def(trunc)
diff --git a/libm/powerpc/classic/w_scalb.c b/libm/powerpc/classic/w_scalb.c
new file mode 100644
index 000000000..408136001
--- /dev/null
+++ b/libm/powerpc/classic/w_scalb.c
@@ -0,0 +1,94 @@
+/***********************************************************************
+** 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.
+**
+***********************************************************************/
+
+#include <math.h>
+#include <endian.h>
+
+typedef union
+ {
+ struct {
+#if (__BYTE_ORDER == __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.
+***********************************************************************/
+
+libm_hidden_proto(scalb)
+#ifdef _SCALB_INT
+double scalb ( double x, int n )
+#else
+double scalb ( double x, double n )
+#endif
+ {
+ 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 );
+ }
+libm_hidden_def(scalb)
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-03 11:27:41 +0000