Move powerpc specific optimizations (courtesy of apple) to powerpc

subdir.  Put together a theoretical framework for adding arch specific
optimizations.  Havn't tried this yet but it looks correct...
 -Erik

7 files changed, 1170 insertions, 0 deletions

diff --git a/libm/powerpc/Makefile b/libm/powerpc/Makefile
new file mode 100644
index 000000000..7464ee0a8
--- /dev/null
+++ b/libm/powerpc/Makefile
@@ -0,0 +1,60 @@
+# Makefile for uClibc's math library
+# Copyright (C) 2001 by Lineo, inc.
+#
+# This math library is derived primarily from the Cephes Math Library,
+# copyright by Stephen L. Moshier <moshier@world.std.com>
+#
+# This program is free software; you can redistribute it and/or modify it under
+# the terms of the GNU Library General Public License as published by the Free
+# Software Foundation; either version 2 of the License, or (at your option) any
+# later version.
+#
+# This program is distributed in the hope that it will be useful, but WITHOUT
+# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+# FOR A PARTICULAR PURPOSE. See the GNU Library General Public License for more
+# details.
+#
+# You should have received a copy of the GNU Library General Public License
+# along with this program; if not, write to the Free Software Foundation, Inc.,
+# 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+
+TOPDIR=../../
+include $(TOPDIR)Rules.mak
+
+LIBM=../libm.a
+TARGET_CC= $(TOPDIR)extra/gcc-uClibc/$(TARGET_ARCH)-uclibc-gcc
+TARGET_CFLAGS+=-D_IEEE_LIBM -D_ISOC99_SOURCE -D_SVID_SOURCE
+
+ifeq ($(strip $(DO_C99_MATH)),true)
+CSRC =   ceilfloor.c  frexpldexp.c  logb.c  rndint.c  scalb.c  sign.c
+else
+CSRC =
+endif
+COBJS=$(patsubst %.c,%.o, $(CSRC))
+OBJS=$(COBJS)
+
+
+ifneq ($(strip $(HAS_FLOATING_POINT)),true)
+all:  clean
+else
+all: $(OBJS) $(LIBM)
+endif
+
+$(LIBM): ar-target
+
+ar-target: $(OBJS)
+	$(AR) $(ARFLAGS) $(LIBM) $(OBJS)
+
+$(COBJS): %.o : %.c
+	$(TARGET_CC) $(TARGET_CFLAGS) -c $< -o $@
+	$(STRIPTOOL) -x -R .note -R .comment $*.o
+
+$(OBJ): Makefile
+
+tags:
+	ctags -R
+
+clean: 
+	rm -f *.[oa] *~ core $(LIBM_SHARED)* $(LIBM_SHARED_FULLNAME)*
+
diff --git a/libm/powerpc/ceilfloor.c b/libm/powerpc/ceilfloor.c
new file mode 100644
index 000000000..d5947ab50
--- /dev/null
+++ b/libm/powerpc/ceilfloor.c
@@ -0,0 +1,174 @@
+/*******************************************************************************
+*                                                                              *
+*      File ceilfloor.c,                                                       *
+*      Function ceil(x) and floor(x),                                          *
+*      Implementation of ceil and floor for the PowerPC.                       *
+*                                                                              *
+*      Copyright © 1991 Apple Computer, Inc.  All rights reserved.             *
+*                                                                              *
+*      Written by Ali Sazegari, started on November 1991,                      *
+*                                                                              *
+*      based on math.h, library code for Macintoshes with a 68881/68882        *
+*      by Jim Thomas.                                                          *
+*                                                                              *
+*      W A R N I N G:  This routine expects a 64 bit double model.             *
+*                                                                              *
+*      December 03 1992: first rs6000 port.                                    *
+*      July     14 1993: comment changes and addition of #pragma fenv_access.  *
+*	 May      06 1997: port of the ibm/taligent ceil and floor routines.     *
+*	 April    11 2001: first port to os x using gcc.				 *
+*	 June	    13 2001: replaced __setflm with in-line assembly			 *
+*                                                                              *
+*******************************************************************************/
+
+static const double        twoTo52  = 4503599627370496.0;
+static const unsigned long signMask = 0x80000000ul;
+
+typedef union
+      {
+      struct {
+#if defined(__BIG_ENDIAN__)
+	unsigned long int hi;
+	unsigned long int lo;
+#else
+	unsigned long int lo;
+	unsigned long int hi;
+#endif
+      } words;
+      double dbl;
+      } DblInHex;
+
+/*******************************************************************************
+*            Functions needed for the computation.                             *
+*******************************************************************************/
+
+/*******************************************************************************
+*      Ceil(x) returns the smallest integer not less than x.                   *
+*******************************************************************************/
+
+double ceil ( double x )
+	{
+	DblInHex xInHex,OldEnvironment;
+	register double y;
+	register unsigned long int xhi;
+	register int target;
+	
+	xInHex.dbl = x;
+	xhi = xInHex.words.hi & 0x7fffffffUL;	  // xhi is the high half of |x|
+	target = ( xInHex.words.hi < signMask );
+	
+	if ( xhi < 0x43300000ul ) 
+/*******************************************************************************
+*      Is |x| < 2.0^52?                                                        *
+*******************************************************************************/
+		{
+		if ( xhi < 0x3ff00000ul ) 
+/*******************************************************************************
+*      Is |x| < 1.0?                                                           *
+*******************************************************************************/
+			{
+			if ( ( xhi | xInHex.words.lo ) == 0ul )  // zero x is exact case
+				return ( x );
+			else 
+				{			                // inexact case
+				asm ("mffs %0" : "=f" (OldEnvironment.dbl));
+				OldEnvironment.words.lo |= 0x02000000ul;
+				asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
+				if ( target )
+					return ( 1.0 );
+				else
+					return ( -0.0 );
+				}
+			}
+/*******************************************************************************
+*      Is 1.0 < |x| < 2.0^52?                                                  *
+*******************************************************************************/
+		if ( target ) 
+			{
+			y = ( x + twoTo52 ) - twoTo52;          // round at binary pt.
+			if ( y < x )
+				return ( y + 1.0 );
+			else
+				return ( y );
+			}
+		
+		else 
+			{
+			y = ( x - twoTo52 ) + twoTo52;          // round at binary pt.
+			if ( y < x )
+				return ( y + 1.0 );
+			else
+				return ( y );
+			}
+		}
+/*******************************************************************************
+*      |x| >= 2.0^52 or x is a NaN.                                            *
+*******************************************************************************/
+	return ( x );
+	}
+
+/*******************************************************************************
+*      Floor(x) returns the largest integer not greater than x.                *
+*******************************************************************************/
+
+double floor ( double x )
+	{
+	DblInHex xInHex,OldEnvironment;
+	register double y;
+	register unsigned long int xhi;
+	register long int target;
+	
+	xInHex.dbl = x;
+	xhi = xInHex.words.hi & 0x7fffffffUL;	  // xhi is the high half of |x|
+	target = ( xInHex.words.hi < signMask );
+	
+	if ( xhi < 0x43300000ul ) 
+/*******************************************************************************
+*      Is |x| < 2.0^52?                                                        *
+*******************************************************************************/
+		{
+		if ( xhi < 0x3ff00000ul ) 
+/*******************************************************************************
+*      Is |x| < 1.0?                                                           *
+*******************************************************************************/
+			{
+			if ( ( xhi | xInHex.words.lo ) == 0ul )  // zero x is exact case
+				return ( x );
+			else 
+				{			                // inexact case
+				asm ("mffs %0" : "=f" (OldEnvironment.dbl));
+				OldEnvironment.words.lo |= 0x02000000ul;
+				asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
+				if ( target )
+					return ( 0.0 );
+				else
+					return ( -1.0 );
+				}
+			}
+/*******************************************************************************
+*      Is 1.0 < |x| < 2.0^52?                                                  *
+*******************************************************************************/
+		if ( target ) 
+			{
+			y = ( x + twoTo52 ) - twoTo52;          // round at binary pt.
+			if ( y > x )
+				return ( y - 1.0 );
+			else
+				return ( y );
+			}
+		
+		else 
+			{
+			y = ( x - twoTo52 ) + twoTo52;          // round at binary pt.
+			if ( y > x )
+				return ( y - 1.0 );
+			else
+				return ( y );
+			}
+		}
+/*******************************************************************************
+*      |x| >= 2.0^52 or x is a NaN.                                            *
+*******************************************************************************/
+	return ( x );
+	}
+
diff --git a/libm/powerpc/frexpldexp.c b/libm/powerpc/frexpldexp.c
new file mode 100644
index 000000000..6b759f91b
--- /dev/null
+++ b/libm/powerpc/frexpldexp.c
@@ -0,0 +1,72 @@
+/*******************************************************************************
+*                                                                              *
+*      File frexpldexp.c,                                                      *
+*      Functions frexp(x) and ldexp(x),                                        *
+*      Implementation of frexp and ldexp functions for the PowerPC.            *
+*                                                                              *
+*      Copyright © 1991 Apple Computer, Inc.  All rights reserved.             *
+*                                                                              *
+*      Written by Ali Sazegari, started on January 1991,                       *
+*                                                                              *
+*      W A R N I N G:  This routine expects a 64 bit double model.             *
+*                                                                              *
+*      December03 1992: first rs6000 implementation.                           *
+*      October 05 1993: added special cases for NaN and ° in frexp.            *
+*      May     27 1997: improved the performance of frexp by eliminating the   *
+*                       switch statement.                                      *
+*	 June	   13 2001: (ram) rewrote frexp to eliminate calls to scalb and    *
+*				logb.									 *
+*                                                                              *
+*******************************************************************************/
+
+#include <limits.h>
+#include <math.h>
+
+static const double two54 =  1.80143985094819840000e+16; /* 0x43500000, 0x00000000 */
+
+typedef union
+      {
+      struct {
+#if defined(__BIG_ENDIAN__)
+        unsigned long int hi;
+        unsigned long int lo;
+#else
+        unsigned long int lo; 
+        unsigned long int hi; 
+#endif
+      } words;
+      double dbl;
+      } DblInHex;       
+
+double ldexp ( double value, int exp ) 
+      {
+      if ( exp > SHRT_MAX ) 
+            exp = SHRT_MAX;
+      else if ( exp < -SHRT_MAX ) 
+            exp = -SHRT_MAX;
+      return scalb ( value, exp  );
+      }
+
+double frexp ( double value, int *eptr )
+      {
+      DblInHex argument;
+      unsigned long int valueHead;
+
+      argument.dbl = value;
+      valueHead = argument.words.hi & 0x7fffffffUL; // valueHead <- |x|
+
+      *eptr = 0;
+	if ( valueHead >= 0x7ff00000 || ( valueHead | argument.words.lo ) == 0 )
+		return value;		// 0, inf, or NaN
+	
+	if ( valueHead < 0x00100000 )
+		{	// denorm
+		argument.dbl = two54 * value;
+		valueHead = argument.words.hi &0x7fffffff;
+		*eptr = -54;
+		}
+	*eptr += ( valueHead >> 20 ) - 1022;
+	argument.words.hi = ( argument.words.hi & 0x800fffff ) | 0x3fe00000;
+	return argument.dbl;
+	}
+
diff --git a/libm/powerpc/logb.c b/libm/powerpc/logb.c
new file mode 100644
index 000000000..3a410ba5c
--- /dev/null
+++ b/libm/powerpc/logb.c
@@ -0,0 +1,103 @@
+/*******************************************************************************
+*                                                                              *
+*      File logb.c,                                                            *
+*      Functions logb.                                                         *
+*      Implementation of logb for the PowerPC.                                 *
+*                                                                              *
+*      Copyright © 1991 Apple Computer, Inc.  All rights reserved.             *
+*                                                                              *
+*      Written by Ali Sazegari, started on June 1991,                          *
+*                                                                              *
+*      August   26 1991: removed CFront Version 1.1d17 warnings.               *
+*      August   27 1991: no errors reported by the test suite.                 *
+*      November 11 1991: changed CLASSEXTENDED to the macro CLASSIFY and       *
+*                        + or - infinity to constants.                         *
+*      November 18 1991: changed the macro CLASSIFY to CLASSEXTENDEDint to     *
+*                        improve performance.                                  *
+*      February 07 1992: changed bit operations to macros (  object size is    *
+*                        unchanged  ).                                         *
+*      September24 1992: took the "#include support.h" out.                    *
+*      December 03 1992: first rs/6000 port.                                   *
+*      August   30 1992: set the divide by zero for the zero argument case.    *
+*      October  05 1993: corrected the environment.                            *
+*      October  17 1994: replaced all environmental functions with __setflm.   *
+*      May      28 1997: made speed improvements.                              *
+*      April    30 2001: forst mac os x port using gcc.                        *
+*                                                                              *
+********************************************************************************
+*     The C math library offers a similar function called "frexp".  It is      *
+*     different in details from logb, but similar in spirit.  This current     *
+*     implementation of logb follows the recommendation in IEEE Standard 854   *
+*     which is different in its handling of denormalized numbers from the IEEE *
+*     Standard 754.                                                            *
+*******************************************************************************/
+
+typedef union           
+      { 
+      struct {
+#if defined(__BIG_ENDIAN__)
+        unsigned long int hi;
+        unsigned long int lo;
+#else
+        unsigned long int lo;
+        unsigned long int hi;
+#endif
+      } words;
+      double dbl;
+      } DblInHex;
+
+static const double twoTo52 = 4.50359962737049600e15;              // 0x1p52
+static const double klTod = 4503601774854144.0;                    // 0x1.000008p52
+static const unsigned long int signMask = 0x80000000ul;
+static const DblInHex minusInf  = {{ 0xFFF00000, 0x00000000 }};
+
+
+/*******************************************************************************
+********************************************************************************
+*                                    L  O  G  B                                *
+********************************************************************************
+*******************************************************************************/
+
+double logb (  double x  )
+      {
+      DblInHex xInHex;
+      long int shiftedExp;
+      
+      xInHex.dbl = x;
+      shiftedExp = ( xInHex.words.hi & 0x7ff00000UL ) >> 20;
+      
+      if ( shiftedExp == 2047 ) 
+            {                                            // NaN or INF
+            if ( ( ( xInHex.words.hi & signMask ) == 0 ) || ( x != x ) )
+                  return x;                              // NaN or +INF return x
+            else
+                  return -x;                             // -INF returns +INF
+            }
+      
+      if ( shiftedExp != 0 )                             // normal number
+            shiftedExp -= 1023;                          // unbias exponent
+      
+      else if ( x == 0.0 ) 
+            {                                            // zero
+            xInHex.words.hi = 0x0UL;                      // return -infinity
+            return (  minusInf.dbl  );
+            }
+      
+      else 
+            {                                            // subnormal number
+            xInHex.dbl *= twoTo52;                       // scale up
+            shiftedExp = ( xInHex.words.hi & 0x7ff00000UL ) >> 20;
+            shiftedExp -= 1075;                          // unbias exponent
+            }
+      
+      if ( shiftedExp == 0 )                             // zero result
+            return ( 0.0 );
+      
+      else 
+            {                                            // nonzero result
+            xInHex.dbl = klTod;
+            xInHex.words.lo += shiftedExp;
+            return ( xInHex.dbl - klTod );
+            }
+      }
+
diff --git a/libm/powerpc/rndint.c b/libm/powerpc/rndint.c
new file mode 100644
index 000000000..1b7a9ec81
--- /dev/null
+++ b/libm/powerpc/rndint.c
@@ -0,0 +1,620 @@
+/*******************************************************************************
+**      File:   rndint.c
+**      
+**      Contains: C source code for implementations of floating-point
+**                functions which round to integral value or format, as
+**                defined in header <fp.h>.  In particular, this file
+**                contains implementations of functions rint, nearbyint,
+**                rinttol, round, roundtol, trunc, modf and modfl.  This file
+**                targets PowerPC or Power platforms.
+**                        
+**      Written by: A. Sazegari, Apple AltiVec Group
+**	    Created originally by Jon Okada, Apple Numerics Group
+**      
+**      Copyright: © 1992-2001 by Apple Computer, Inc., all rights reserved
+**      
+**      Change History (most recent first):
+**
+**      13 Jul 01  ram  replaced --setflm calls with inline assembly
+**      03 Mar 01  ali	first port to os x using gcc, added the crucial __setflm
+**                      definition.
+**				1. removed double_t, put in double for now.
+**				2. removed iclass from nearbyint.
+**				3. removed wrong comments intrunc.
+**				4. 
+**      13 May 97  ali  made performance improvements in rint, rinttol, roundtol
+**                      and trunc by folding some of the taligent ideas into this
+**                      implementation.  nearbyint is faster than the one in taligent,
+**                      rint is more elegant, but slower by %30 than the taligent one.
+**      09 Apr 97  ali  deleted modfl and deferred to AuxiliaryDD.c
+**      15 Sep 94  ali  Major overhaul and performance improvements of all functions.
+**      20 Jul 94  PAF  New faster version
+**      16 Jul 93  ali  Added the modfl function.
+**      18 Feb 93  ali  Changed the return value of fenv functions
+**                      feclearexcept and feraiseexcept to their new
+**                      NCEG X3J11.1/93-001 definitions.
+**      16 Dec 92  JPO  Removed __itrunc implementation to a 
+**                      separate file.
+**      15 Dec 92  JPO  Added __itrunc implementation and modified
+**                      rinttol to include conversion from double
+**                      to long int format.  Modified roundtol to
+**                      call __itrunc.
+**      10 Dec 92  JPO  Added modf (double) implementation.
+**      04 Dec 92  JPO  First created.
+**                        
+*******************************************************************************/
+
+#include <limits.h>
+#include <math.h>
+
+#define      SET_INVALID      0x01000000UL
+
+typedef union
+      {
+      struct {
+#if defined(__BIG_ENDIAN__)
+        unsigned long int hi;
+        unsigned long int lo;
+#else
+        unsigned long int lo;
+        unsigned long int hi;
+#endif
+      } words;
+      double dbl;
+      } DblInHex;
+
+static const unsigned long int signMask = 0x80000000ul;
+static const double twoTo52      = 4503599627370496.0;
+static const double doubleToLong = 4503603922337792.0;	            // 2^52
+static const DblInHex Huge       = {{ 0x7FF00000, 0x00000000 }};
+static const DblInHex TOWARDZERO = {{ 0x00000000, 0x00000001 }};
+
+/*******************************************************************************
+*                                                                              *
+*     The function rint rounds its double argument to integral value           *
+*     according to the current rounding direction and returns the result in    *
+*     double format.  This function signals inexact if an ordered return       * 
+*     value is not equal to the operand.                                       *
+*                                                                              *
+********************************************************************************
+*                                                                              *
+*     This function calls:  fabs. 	                                           *
+*                                                                              *
+*******************************************************************************/
+
+/*******************************************************************************
+*     First, an elegant implementation.                                        *
+********************************************************************************
+*
+*double rint ( double x )
+*      {
+*      double y;
+*      
+*      y = twoTo52.fval;
+*      
+*      if ( fabs ( x ) >= y )                          // huge case is exact 
+*            return x;
+*      if ( x < 0 ) y = -y;                            // negative case 
+*      y = ( x + y ) - y;                              // force rounding 
+*      if ( y == 0.0 )                                 // zero results mirror sign of x 
+*            y = copysign ( y, x );
+*      return ( y );      
+*      }
+********************************************************************************
+*     Now a bit twidling version that is about %30 faster.                     *
+*******************************************************************************/
+
+double rint ( double x )
+      {
+      DblInHex argument;
+      register double y;
+      unsigned long int xHead;
+      register long int target;
+      
+      argument.dbl = x;
+      xHead = argument.words.hi & 0x7fffffffUL;          // xHead <- high half of |x|
+      target = ( argument.words.hi < signMask );         // flags positive sign
+      
+      if ( xHead < 0x43300000ul ) 
+/*******************************************************************************
+*     Is |x| < 2.0^52?                                                         *
+*******************************************************************************/
+            {
+            if ( xHead < 0x3ff00000ul ) 
+/*******************************************************************************
+*     Is |x| < 1.0?                                                            *
+*******************************************************************************/
+                  {
+                  if ( target )
+                        y = ( x + twoTo52 ) - twoTo52;  // round at binary point
+                  else
+                        y = ( x - twoTo52 ) + twoTo52;  // round at binary point
+                  if ( y == 0.0 ) 
+                        {                               // fix sign of zero result
+                        if ( target )
+                              return ( 0.0 );
+                        else
+                              return ( -0.0 );
+                        }
+                  return y;
+                  }
+            
+/*******************************************************************************
+*     Is 1.0 < |x| < 2.0^52?                                                   *
+*******************************************************************************/
+
+            if ( target )
+                  return ( ( x + twoTo52 ) - twoTo52 ); //   round at binary pt.
+            else
+                  return ( ( x - twoTo52 ) + twoTo52 );
+            }
+      
+/*******************************************************************************
+*     |x| >= 2.0^52 or x is a NaN.                                             *
+*******************************************************************************/
+      return ( x );
+      }
+
+/*******************************************************************************
+*                                                                              *
+*     The function nearbyint rounds its double argument to integral value      *
+*     according to the current rounding direction and returns the result in    *
+*     double format.  This function does not signal inexact.                   *
+*                                                                              *
+********************************************************************************
+*                                                                              *
+*     This function calls fabs and copysign.		                         *
+*                                                                              *
+*******************************************************************************/
+   
+double nearbyint ( double x )
+      {
+	double y;
+	double OldEnvironment;
+      
+	y = twoTo52;
+	
+	asm ("mffs %0" : "=f" (OldEnvironment));	/* get the environement */
+
+      if ( fabs ( x ) >= y )                          /* huge case is exact */
+            return x;
+      if ( x < 0 ) y = -y;                                   /* negative case */
+      y = ( x + y ) - y;                                    /* force rounding */
+      if ( y == 0.0 )                        /* zero results mirror sign of x */
+            y = copysign ( y, x );
+//	restore old flags
+	asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment ));
+      return ( y );      
+	}
+      
+/*******************************************************************************
+*                                                                              *
+*     The function rinttol converts its double argument to integral value      *
+*     according to the current rounding direction and returns the result in    *
+*     long int format.  This conversion signals invalid if the argument is a   *
+*     NaN or the rounded intermediate result is out of range of the            *
+*     destination long int format, and it delivers an unspecified result in    *
+*     this case.  This function signals inexact if the rounded result is       *
+*     within range of the long int format but unequal to the operand.          *
+*                                                                              *
+*******************************************************************************/
+
+long int rinttol ( double x )
+      {
+      register double y;
+      DblInHex argument, OldEnvironment;
+      unsigned long int xHead;
+      register long int target;
+      
+      argument.dbl = x;
+      target = ( argument.words.hi < signMask );        // flag positive sign
+      xHead = argument.words.hi & 0x7ffffffful;         // high 32 bits of x
+      
+      if ( target ) 
+/*******************************************************************************
+*    Sign of x is positive.                                                    *
+*******************************************************************************/
+            {
+            if ( xHead < 0x41dffffful ) 
+                  {                                    // x is safely in long range
+                  y = ( x + twoTo52 ) - twoTo52;       // round at binary point
+                  argument.dbl = y + doubleToLong;     // force result into argument.words.lo
+                  return ( ( long ) argument.words.lo );
+                  }
+            
+		asm ("mffs %0" : "=f" (OldEnvironment.dbl));	// get environment
+            
+            if ( xHead > 0x41dffffful ) 
+                  {                                    // x is safely out of long range
+                  OldEnvironment.words.lo |= SET_INVALID;
+			asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
+                  return ( LONG_MAX );
+                  }
+            
+/*******************************************************************************
+*     x > 0.0 and may or may not be out of range of long.                      *
+*******************************************************************************/
+
+            y = ( x + twoTo52 ) - twoTo52;             // do rounding
+            if ( y > ( double ) LONG_MAX ) 
+                  {                                    // out of range of long
+                  OldEnvironment.words.lo |= SET_INVALID;
+			asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
+                  return ( LONG_MAX );
+                  }
+            argument.dbl = y + doubleToLong;           // in range
+            return ( ( long ) argument.words.lo );      // return result & flags
+            }
+      
+/*******************************************************************************
+*    Sign of x is negative.                                                    *
+*******************************************************************************/
+      if ( xHead < 0x41e00000ul ) 
+            {                                          // x is safely in long range
+            y = ( x - twoTo52 ) + twoTo52;
+            argument.dbl = y + doubleToLong;
+            return ( ( long ) argument.words.lo );
+            }
+      
+	asm ("mffs %0" : "=f" (OldEnvironment.dbl));	// get environment
+      
+      if ( xHead > 0x41e00000ul ) 
+            {                                          // x is safely out of long range
+            OldEnvironment.words.lo |= SET_INVALID;
+		asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
+            return ( LONG_MIN );
+            }
+      
+/*******************************************************************************
+*    x < 0.0 and may or may not be out of range of long.                       *
+*******************************************************************************/
+
+      y = ( x - twoTo52 ) + twoTo52;                   // do rounding
+      if ( y < ( double ) LONG_MIN ) 
+            {                                          // out of range of long
+            OldEnvironment.words.lo |= SET_INVALID;
+		asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
+            return ( LONG_MIN );
+            }
+      argument.dbl = y + doubleToLong;                       // in range
+      return ( ( long ) argument.words.lo );           // return result & flags
+      }
+
+/*******************************************************************************
+*                                                                              *
+*     The function round rounds its double argument to integral value          *
+*     according to the "add half to the magnitude and truncate" rounding of    *
+*     Pascal's Round function and FORTRAN's ANINT function and returns the     *
+*     result in double format.  This function signals inexact if an ordered    *
+*     return value is not equal to the operand.                                *
+*                                                                              *
+*******************************************************************************/
+   
+double round ( double x )
+      {      
+      DblInHex argument, OldEnvironment;
+      register double y, z;
+      register unsigned long int xHead;
+      register long int target;
+      
+      argument.dbl = x;
+      xHead = argument.words.hi & 0x7fffffffUL;      // xHead <- high half of |x|
+      target = ( argument.words.hi < signMask );     // flag positive sign
+      
+      if ( xHead < 0x43300000ul ) 
+/*******************************************************************************
+*     Is |x| < 2.0^52?                                                        *
+*******************************************************************************/
+            {
+            if ( xHead < 0x3ff00000ul ) 
+/*******************************************************************************
+*     Is |x| < 1.0?                                                           *
+*******************************************************************************/
+                  {
+			asm ("mffs %0" : "=f" (OldEnvironment.dbl));	// get environment
+                  if ( xHead < 0x3fe00000ul ) 
+/*******************************************************************************
+*     Is |x| < 0.5?                                                           *
+*******************************************************************************/
+                        {
+                        if ( ( xHead | argument.words.lo ) != 0ul )
+                              OldEnvironment.words.lo |= 0x02000000ul;
+				asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
+                        if ( target ) 
+                              return ( 0.0 );
+                        else
+                              return ( -0.0 );
+                        }
+/*******************************************************************************
+*     Is 0.5 ² |x| < 1.0?                                                      *
+*******************************************************************************/
+                  OldEnvironment.words.lo |= 0x02000000ul;
+			asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
+                  if ( target )
+                        return ( 1.0 );
+                  else
+                        return ( -1.0 );
+                  }
+/*******************************************************************************
+*     Is 1.0 < |x| < 2.0^52?                                                   *
+*******************************************************************************/
+            if ( target ) 
+                  {                                     // positive x
+                  y = ( x + twoTo52 ) - twoTo52;        // round at binary point
+                  if ( y == x )                         // exact case
+                        return ( x );
+                  z = x + 0.5;                          // inexact case
+                  y = ( z + twoTo52 ) - twoTo52;        // round at binary point
+                  if ( y > z )
+                        return ( y - 1.0 );
+                  else
+                        return ( y );
+                  }
+            
+/*******************************************************************************
+*     Is x < 0?                                                                *
+*******************************************************************************/
+            else 
+                  {
+                  y = ( x - twoTo52 ) + twoTo52;        // round at binary point
+                  if ( y == x )
+                        return ( x );
+                  z = x - 0.5;
+                  y = ( z - twoTo52 ) + twoTo52;        // round at binary point
+                  if ( y < z )
+                        return ( y + 1.0 );
+                  else
+                  return ( y );
+                  }
+            }
+/*******************************************************************************
+*      |x| >= 2.0^52 or x is a NaN.                                            *
+*******************************************************************************/
+      return ( x );
+      }
+
+/*******************************************************************************
+*                                                                              *
+*     The function roundtol converts its double argument to integral format    *
+*     according to the "add half to the magnitude and chop" rounding mode of   *
+*     Pascal's Round function and FORTRAN's NINT function.  This conversion    *
+*     signals invalid if the argument is a NaN or the rounded intermediate     *
+*     result is out of range of the destination long int format, and it        *
+*     delivers an unspecified result in this case.  This function signals      *
+*     inexact if the rounded result is within range of the long int format but *
+*     unequal to the operand.                                                  *
+*                                                                              *
+*******************************************************************************/
+
+long int roundtol ( double x )
+	{	
+	register double y, z;
+	DblInHex argument, OldEnvironmen