- remove files that are not either LGPL or Public Domain.

- pull replacement funcs for fpmacros.c from glibc

  This removes the powerpc/classic implementation which did not state any
  license but read:
  Copyright © 1991 Apple Computer, Inc.  All rights reserved.
  and thus was dubious (and not needed).

libm/Makefile.in

@@ -34,15 +34,10 @@ libm_OUT:=$(top_builddir)libm
 
 # Fix builds for powerpc as there are different cores in this 
 # section now.`
-ifeq ($(TARGET_ARCH),powerpc)
-ifeq ($(CONFIG_E500),y)
+ifeq ($(TARGET_ARCH)-$(CONFIG_E500),powerpc-y)
 libm_ARCH_DIR:=$(libm_DIR)/$(TARGET_ARCH)/e500
 libm_ARCH_OUT:=$(libm_OUT)/$(TARGET_ARCH)/e500
 else
-libm_ARCH_DIR:=$(libm_DIR)/$(TARGET_ARCH)/classic
-libm_ARCH_OUT:=$(libm_OUT)/$(TARGET_ARCH)/classic
-endif
-else
 libm_ARCH_DIR:=$(libm_DIR)/$(TARGET_ARCH)
 libm_ARCH_OUT:=$(libm_OUT)/$(TARGET_ARCH)
 endif
@@ -74,7 +69,9 @@ libm_CSRC := \
 	w_cabs.c w_cosh.c w_drem.c w_exp.c w_fmod.c w_gamma.c w_gamma_r.c \
 	w_hypot.c w_j0.c w_j1.c w_jn.c w_lgamma.c w_lgamma_r.c \
 	w_log.c w_log10.c w_pow.c w_remainder.c w_scalb.c w_sinh.c \
-	w_sqrt.c fpmacros.c nan.c carg.c s_llrint.c
+	w_sqrt.c nan.c carg.c s_llrint.c \
+	s_fpclassify.c s_fpclassifyf.c s_signbit.c s_signbitf.c \
+	s_isnan.c s_isnanf.c s_isnf.c s_isinff.c
 FL_MOBJ := \
 	acosf.o acoshf.o asinf.o asinhf.o atan2f.o atanf.o atanhf.o cbrtf.o \
 	ceilf.o copysignf.o cosf.o coshf.o erfcf.o erff.o exp2f.o expf.o \
@@ -103,13 +100,9 @@ endif
 ifeq ($(UCLIBC_HAS_FPU),y)
 ifeq ($(DO_C99_MATH),y)
 ifneq ($(strip $(libm_ARCH_OBJS)),)
-ifeq ($(TARGET_ARCH),powerpc)
-ifeq ($(CONFIG_E500),y)
+ifeq ($(TARGET_ARCH)-$(CONFIG_E500),powerpc-y)
 CFLAGS-libm/$(TARGET_ARCH)/e500/ := $(CFLAGS-libm)
 else
-CFLAGS-libm/$(TARGET_ARCH)/classic/ := $(CFLAGS-libm)
-endif
-else
 CFLAGS-libm/$(TARGET_ARCH)/ := $(CFLAGS-libm)
 endif
 

libm/fp_private.h

@@ -1,90 +0,0 @@
-/*******************************************************************************
-*                                                                              *
-*      File fp_private.h,                                                      *
-*      All pack 4 dependencies for the MathLib elems plus some defines used    *
-*      throughout MathLib.                                                     *
-*                                                                              *
-*      Copyright © 1991 Apple Computer, Inc.  All rights reserved.             *
-*                                                                              *
-*      Written by Ali Sazegari, started on October 1991,                       *
-*                                                                              *
-*      W A R N I N G:  This routine expects a 64 bit double model.             *
-*                                                                              *
-*******************************************************************************/
-
-#define      NoException            0
-
-/*******************************************************************************
-*                              Values of constants.                            *
-*******************************************************************************/
-
-//#define    SgnMask            0x8000
-#define      dSgnMask           0x80000000
-#define      sSgnMask           0x7FFFFFFF
-
-//#define    ExpMask            0x7FFF
-#define      dExpMask           0x7FF00000
-#define      sExpMask           0xFF000000
-
-                                          /* according to rounding BIG & SMALL are:  */
-#define      BIG               1.1e+300   /* used to deliver ±° or largest number,   */
-#define      SMALL             1.1e-300   /* used to deliver ±0 or smallest number.  */
-#define      InfExp            0x7FF
-#define      dMaxExp           0x7FF00000
-
-#define      MaxExpP1          1024
-#define      MaxExp            1023
-
-#define      DenormLimit       -52
-
-//#define    ManMask           0x80000000
-#define      dManMask          0x00080000
-
-//#define    IsItDenorm         0x80000000
-#define      dIsItDenorm        0x00080000
-
-//#define    xIsItSNaN          0x40000000
-#define      dIsItSNaN          0x00080000
-
-#define      dHighMan           0x000FFFFF
-#define      dFirstBitSet       0x00080000
-#define      BIAS               0x3FF
-
-//#define    GetSign            0x8000
-#define      dGetSign           0x80000000
-#define      sGetSign           0x80000000
-
-//#define    Infinity(x)       ( x.hex.exponent & ExpMask ) == ExpMask
-#define      dInfinity(x)      ( x.hex.high & dExpMask ) == dExpMask
-#define      sInfinity(x)      ( ( x.hexsgl << 1 ) & sExpMask ) == sExpMask
-
-//#define    Exponent(x)       x.hex.exponent & ExpMask
-#define      dExponent(x)      x.hex.high & dExpMask
-#define      sExponent(x)      ( ( x.hexsgl << 1 ) & sExpMask )
-
-#define      sZero(x)          ( x.hexsgl & sSgnMask ) == 0
-//#define    Sign(x)           ( x.hex.exponent & SgnMask ) == SgnMask
-
-/*******************************************************************************
-*                        Types used in the auxiliary functions.                *
-*******************************************************************************/
-
-#include <stdint.h>
-#include <endian.h>
-
-typedef struct                   /*      Hex representation of a double.      */
-      {
-#if (__BYTE_ORDER == __BIG_ENDIAN)
-      uint32_t high;
-      uint32_t low;
-#else
-      uint32_t low;
-      uint32_t high;
-#endif
-      } dHexParts;
-
-typedef union
-      {
-      unsigned char byties[8];
-      double dbl;
-      } DblInHex;

libm/fpmacros.c

@@ -1,303 +0,0 @@
-/***********************************************************************
-**  File:  fpmacros.c
-**
-**  Contains:  C source code for implementations of floating-point
-**             functions which involve float format numbers, as
-**             defined in header <fp.h>.  In particular, this file
-**             contains implementations of functions
-**              __fpclassify(d,f), __isnormal(d,f), __isfinite(d,f),
-**             __isnan(d,f), and __signbit(d,f).  This file targets
-**             PowerPC platforms.
-**
-**  Written by:   Robert A. Murley, Ali Sazegari
-**
-**  Copyright:   c 2001 by Apple Computer, Inc., all rights reserved
-**
-**  Change History (most recent first):
-**
-**     07 Jul 01   ram      First created from fpfloatfunc.c, fp.c,
-**				classify.c and sign.c in MathLib v3 Mac OS9.
-**
-***********************************************************************/
-
-#include     <features.h>
-#include     <sys/types.h>
-#include     <math.h>
-#include     "fp_private.h"
-
-#define SIGN_MASK 0x80000000
-#define NSIGN_MASK 0x7fffffff
-#define FEXP_MASK 0x7f800000
-#define FFRAC_MASK 0x007fffff
-
-/***********************************************************************
-   int __fpclassifyf(float x) returns the classification code of the
-   argument x, as defined in <fp.h>.
-
-   Exceptions:  INVALID signaled if x is a signaling NaN; in this case,
-                the FP_QNAN code is returned.
-
-   Calls:  none
-***********************************************************************/
-
-libm_hidden_proto(__fpclassifyf)
-int __fpclassifyf ( float x )
-{
-   unsigned int iexp;
-
-   union {
-      u_int32_t lval;
-      float fval;
-   } z;
-
-   z.fval = x;
-   iexp = z.lval & FEXP_MASK;                 /* isolate float exponent */
-
-   if (iexp == FEXP_MASK) {                   /* NaN or INF case */
-      if ((z.lval & 0x007fffff) == 0)
-         return FP_INFINITE;
-	return FP_NAN;
-   }
-
-   if (iexp != 0)                             /* normal float */
-      return FP_NORMAL;
-
-   if (x == 0.0)
-      return FP_ZERO;             /* zero */
-   else
-      return FP_SUBNORMAL;        /* must be subnormal */
-}
-libm_hidden_def(__fpclassifyf)
-
-
-/***********************************************************************
-      Function __fpclassify,
-      Implementation of classify of a double number for the PowerPC.
-
-   Exceptions:  INVALID signaled if x is a signaling NaN; in this case,
-                the FP_QNAN code is returned.
-
-   Calls:  none
-***********************************************************************/
-
-libm_hidden_proto(__fpclassify)
-int __fpclassify ( double arg )
-{
-	register unsigned int exponent;
-      union
-            {
-            dHexParts hex;
-            double dbl;
-            } x;
-
-	x.dbl = arg;
-
-	exponent = x.hex.high & dExpMask;
-	if ( exponent == dExpMask )
-		{
-		if ( ( ( x.hex.high & dHighMan ) | x.hex.low ) == 0 )
-			return FP_INFINITE;
-		else
-            	return FP_NAN;
-		}
-	else if ( exponent != 0)
-		return FP_NORMAL;
-	else {
-		if ( arg == 0.0 )
-			return FP_ZERO;
-		else
-			return FP_SUBNORMAL;
-		}
-}
-libm_hidden_def(__fpclassify)
-
-
-/***********************************************************************
-   int __isnormalf(float x) returns nonzero if and only if x is a
-   normalized float number and zero otherwise.
-
-   Exceptions:  INVALID is raised if x is a signaling NaN; in this case,
-                zero is returned.
-
-   Calls:  none
-***********************************************************************/
-
-int __isnormalf ( float x );
-int __isnormalf ( float x )
-{
-   unsigned int iexp;
-   union {
-      u_int32_t lval;
-      float fval;
-   } z;
-
-   z.fval = x;
-   iexp = z.lval & FEXP_MASK;                 /* isolate float exponent */
-   return ((iexp != FEXP_MASK) && (iexp != 0));
-}
-
-
-int __isnormal ( double x );
-int __isnormal ( double x )
-{
-	return ( __fpclassify ( x ) == FP_NORMAL );
-}
-
-
-/***********************************************************************
-   int __isfinitef(float x) returns nonzero if and only if x is a
-   finite (normal, subnormal, or zero) float number and zero otherwise.
-
-   Exceptions:  INVALID is raised if x is a signaling NaN; in this case,
-                zero is returned.
-
-   Calls:  none
-***********************************************************************/
-
-int __finitef ( float x )
-{
-   union {
-      u_int32_t lval;
-      float fval;
-   } z;
-
-   z.fval = x;
-   return ((z.lval & FEXP_MASK) != FEXP_MASK);
-}
-strong_alias(__finitef,finitef)
-
-#if 0 /* use __finite in s_finite.c */
-int __finite ( double x )
-{
-	return ( __fpclassify ( x ) >= FP_ZERO );
-}
-strong_alias(__finite,finite)
-#endif
-
-
-/***********************************************************************
-   int __signbitf(float x) returns nonzero if and only if the sign
-   bit of x is set and zero otherwise.
-
-   Exceptions:  INVALID is raised if x is a signaling NaN.
-
-   Calls:  none
-***********************************************************************/
-
-libm_hidden_proto(__signbitf)
-int __signbitf ( float x )
-{
-   union {
-      u_int32_t lval;
-      float fval;
-   } z;
-
-   z.fval = x;
-   return ((z.lval & SIGN_MASK) != 0);
-}
-libm_hidden_def(__signbitf)
-
-
-/***********************************************************************
-      Function sign of a double.
-      Implementation of sign bit for the PowerPC.
-
-   Calls:  none
-***********************************************************************/
-
-libm_hidden_proto(__signbit)
-int __signbit ( double arg )
-{
-      union
-            {
-            dHexParts hex;
-            double dbl;
-            } x;
-      int sign;
-
-      x.dbl = arg;
-      sign = ( ( x.hex.high & dSgnMask ) == dSgnMask ) ? 1 : 0;
-      return sign;
-}
-libm_hidden_def(__signbit)
-
-
-/***********************************************************************
-* int __isinff(float x) returns -1 if value represents  negative
-*	infinity,  1  if value represents positive infinity,
-*	and 0 otherwise.
-*
-* Calls:  __signbit
-* +***********************************************************************/
-int __isinff ( float x )
-{
-    int class = __fpclassifyf(x);
-    if ( class == FP_INFINITE ) {
-	return ( (__signbitf(x)) ? -1 : 1);
-    }
-    return 0;
-}
-strong_alias(__isinff,isinff)
-
-int __isinf ( double x )
-{
-    int class = __fpclassify(x);
-    if ( class == FP_INFINITE ) {
-	return ( (__signbit(x)) ? -1 : 1);
-    }
-    return 0;
-}
-strong_alias(__isinf,isinf)
-
-#if 0
-int __isinfl ( long double x )
-{
-    int class = __fpclassify(x);
-    if ( class == FP_INFINITE ) {
-	return ( (__signbit(x)) ? -1 : 1);
-    }
-    return 0;
-}
-strong_alias(__isinfl,isinfl)
-#endif
-
-/***********************************************************************
-   int __isnanf(float x) returns nonzero if and only if x is a
-   NaN and zero otherwise.
-
-   Exceptions:  INVALID is raised if x is a signaling NaN; in this case,
-                nonzero is returned.
-
-   Calls:  none
-***********************************************************************/
-
-int __isnanf ( float x )
-{
-   union {
-      u_int32_t lval;
-      float fval;
-   } z;
-
-   z.fval = x;
-   return (((z.lval&FEXP_MASK) == FEXP_MASK) && ((z.lval&FFRAC_MASK) != 0));
-}
-strong_alias(__isnanf,isnanf)
-
-libm_hidden_proto(__isnan)
-int __isnan ( double x )
-{
-	int class = __fpclassify(x);
-	return ( class == FP_NAN );
-}
-libm_hidden_def(__isnan)
-strong_alias(__isnan,isnan)
-
-#if 0
-int __isnanl ( long double x )
-{
-	int class = __fpclassify(x);
-	return ( class == FP_NAN );
-}
-strong_alias(__isnanl,isnanl)
-#endif
-

libm/powerpc/classic/s_ceil.c

@@ -1,113 +0,0 @@
-/*******************************************************************************
-*                                                                              *
-*      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)

libm/powerpc/classic/s_copysign.c

@@ -1,59 +0,0 @@
-/*******************************************************************************
-*                                                                              *
-*      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)

libm/powerpc/classic/s_floor.c

@@ -1,113 +0,0 @@
-/*******************************************************************************
-*                                                                              *
-*      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)

libm/powerpc/classic/s_frexp.c

@@ -1,65 +0,0 @@
-/*******************************************************************************
-*                                                                              *
-*      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)

libm/powerpc/classic/s_ldexp.c

@@ -1,49 +0,0 @@
-/*******************************************************************************
-*                                                                              *
-*      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)

libm/powerpc/classic/s_logb.c

@@ -1,107 +0,0 @@
-/*******************************************************************************
-*                                                                              *
-*      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)

libm/powerpc/classic/s_modf.c

@@ -1,344 +0,0 @@
-/*******************************************************************************
-**      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)

libm/powerpc/classic/s_nearbyint.c

@@ -1,38 +0,0 @@
-#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)

libm/powerpc/classic/s_rint.c

@@ -1,159 +0,0 @@
-/*******************************************************************************
-**      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)

libm/powerpc/classic/s_round.c

@@ -1,115 +0,0 @@
-#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)

libm/powerpc/classic/s_trunc.c

@@ -1,89 +0,0 @@
-#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)

libm/powerpc/classic/w_scalb.c

@@ -1,94 +0,0 @@
-/***********************************************************************
-**      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)

libm/s_fpclassify.c

@@ -0,0 +1,44 @@
+/* Return classification value corresponding to argument.
+   Copyright (C) 1997, 2002 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+   Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C Library 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
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, write to the Free
+   Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
+   02111-1307 USA.  */
+
+#include <math.h>
+
+#include "math_private.h"
+
+libm_hidden_proto (__fpclassify)
+int
+__fpclassify (double x)
+{
+  u_int32_t hx, lx;
+  int retval = FP_NORMAL;
+
+  EXTRACT_WORDS (hx, lx, x);
+  lx |= hx & 0xfffff;
+  hx &= 0x7ff00000;
+  if ((hx | lx) == 0)
+    retval = FP_ZERO;
+  else if (hx == 0)
+    retval = FP_SUBNORMAL;
+  else if (hx == 0x7ff00000)
+    retval = lx != 0 ? FP_NAN : FP_INFINITE;
+
+  return retval;
+}
+libm_hidden_def (__fpclassify)

libm/s_fpclassifyf.c

@@ -0,0 +1,44 @@
+/* Return classification value corresponding to argument.
+   Copyright (C) 1997, 2000, 2002 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+   Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C Library 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
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, write to the Free
+   Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
+   02111-1307 USA.  */
+
+#include <math.h>
+
+#include "math_private.h"
+
+
+libm_hidden_proto (__fpclassifyf)
+int
+__fpclassifyf (float x)
+{
+  u_int32_t wx;
+  int retval = FP_NORMAL;
+
+  GET_FLOAT_WORD (wx, x);
+  wx &= 0x7fffffff;
+  if (wx == 0)
+    retval = FP_ZERO;
+  else if (wx < 0x800000)
+    retval = FP_SUBNORMAL;
+  else if (wx >= 0x7f800000)
+    retval = wx > 0x7f800000 ? FP_NAN : FP_INFINITE;
+
+  return retval;
+}
+libm_hidden_def (__fpclassifyf)

libm/s_isinf.c

@@ -0,0 +1,34 @@
+/*
+ * Written by J.T. Conklin <jtc@netbsd.org>.
+ * Changed to return -1 for -Inf by Ulrich Drepper <drepper@cygnus.com>.
+ * Public domain.
+ */
+
+#if defined(LIBM_SCCS) && !defined(lint)
+static char rcsid[] = "$NetBSD: s_isinf.c,v 1.3 1995/05/11 23:20:14 jtc Exp $";
+#endif
+
+/*
+ * isinf(x) returns 1 is x is inf, -1 if x is -inf, else 0;
+ * no branching!
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+libm_hidden_proto(__isinf)
+int
+__isinf (double x)
+{
+	int32_t hx,lx;
+	EXTRACT_WORDS(hx,lx,x);
+	lx |= (hx & 0x7fffffff) ^ 0x7ff00000;
+	lx |= -lx;
+	return ~(lx >> 31) & (hx >> 30);
+}
+libm_hidden_def (__isinf)
+weak_alias (__isinf, isinf)
+#ifdef NO_LONG_DOUBLE
+strong_alias (__isinf, __isinfl)
+weak_alias (__isinf, isinfl)
+#endif

libm/s_isinff.c

@@ -0,0 +1,30 @@
+/*
+ * Written by J.T. Conklin <jtc@netbsd.org>.
+ * Public domain.
+ */
+
+#if defined(LIBM_SCCS) && !defined(lint)
+static char rcsid[] = "$NetBSD: s_isinff.c,v 1.3 1995/05/11 23:20:21 jtc Exp $";
+#endif
+
+/*
+ * isinff(x) returns 1 if x is inf, -1 if x is -inf, else 0;
+ * no branching!
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+libm_hidden_proto(__isinff)
+int
+__isinff (float x)
+{
+	int32_t ix,t;
+	GET_FLOAT_WORD(ix,x);
+	t = ix & 0x7fffffff;
+	t ^= 0x7f800000;
+	t |= -t;
+	return ~(t >> 31) & (ix >> 30);
+}
+libm_hidden_def (__isinff)
+weak_alias (__isinff, isinff)

libm/s_isnan.c

@@ -0,0 +1,45 @@
+/* @(#)s_isnan.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+
+#if defined(LIBM_SCCS) && !defined(lint)
+static char rcsid[] = "$NetBSD: s_isnan.c,v 1.8 1995/05/10 20:47:36 jtc Exp $";
+#endif
+
+/*
+ * isnan(x) returns 1 is x is nan, else 0;
+ * no branching!
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+libm_hidden_proto (__isnan)
+#ifdef __STDC__
+	int __isnan(double x)
+#else
+	int __isnan(x)
+	double x;
+#endif
+{
+	int32_t hx,lx;
+	EXTRACT_WORDS(hx,lx,x);
+	hx &= 0x7fffffff;
+	hx |= (u_int32_t)(lx|(-lx))>>31;
+	hx = 0x7ff00000 - hx;
+	return (int)(((u_int32_t)hx)>>31);
+}
+libm_hidden_def (__isnan)
+weak_alias (__isnan, isnan)
+#ifdef NO_LONG_DOUBLE
+strong_alias (__isnan, __isnanl)
+weak_alias (__isnan, isnanl)
+#endif

libm/s_isnanf.c

@@ -0,0 +1,43 @@
+/* s_isnanf.c -- float version of s_isnan.c.
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ */
+
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+
+#if defined(LIBM_SCCS) && !defined(lint)
+static char rcsid[] = "$NetBSD: s_isnanf.c,v 1.4 1995/05/10 20:47:38 jtc Exp $";
+#endif
+
+/*
+ * isnanf(x) returns 1 is x is nan, else 0;
+ * no branching!
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+libm_hidden_proto (__isnanf)
+#ifdef __STDC__
+	int __isnanf(float x)
+#else
+	int __isnanf(x)
+	float x;
+#endif
+{
+	int32_t ix;
+	GET_FLOAT_WORD(ix,x);
+	ix &= 0x7fffffff;
+	ix = 0x7f800000 - ix;
+	return (int)(((u_int32_t)(ix))>>31);
+}
+libm_hidden_def (__isnanf)
+weak_alias (__isnanf, isnanf)

libm/s_signbit.c

@@ -0,0 +1,34 @@
+/* Return nonzero value if number is negative.
+   Copyright (C) 1997 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+   Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C Library 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
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, write to the Free
+   Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
+   02111-1307 USA.  */
+
+#include <math.h>
+
+#include "math_private.h"
+
+libm_hidden_proto(__signbit)
+int
+__signbit (double x)
+{
+  int32_t hx;
+
+  GET_HIGH_WORD (hx, x);
+  return hx & 0x80000000;
+}
+libm_hidden_def(__signbit)

libm/s_signbitf.c

@@ -0,0 +1,34 @@
+/* Return nonzero value if number is negative.
+   Copyright (C) 1997 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+   Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C Library 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
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, write to the Free
+   Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
+   02111-1307 USA.  */
+
+#include <math.h>
+
+#include "math_private.h"
+
+libm_hidden_proto(__signbitf)
+int
+__signbitf (float x)
+{
+  int32_t hx;
+
+  GET_FLOAT_WORD (hx, x);
+  return hx & 0x80000000;
+}
+libm_hidden_def(__signbitf)