initial import

24 files changed, 2234 insertions, 0 deletions

diff --git a/sysdeps/ieee754/Dist b/sysdeps/ieee754/Dist
new file mode 100644
index 0000000000..94cc5c9ac4
--- /dev/null
+++ b/sysdeps/ieee754/Dist
@@ -0,0 +1,2 @@
+support.c
+ieee754.h
diff --git a/sysdeps/ieee754/cabs.c b/sysdeps/ieee754/cabs.c
new file mode 100644
index 0000000000..6b0d4c4cde
--- /dev/null
+++ b/sysdeps/ieee754/cabs.c
@@ -0,0 +1,228 @@
+/*
+ * Copyright (c) 1985 Regents of the University of California.
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ * 3. All advertising materials mentioning features or use of this software
+ *    must display the following acknowledgement:
+ *	This product includes software developed by the University of
+ *	California, Berkeley and its contributors.
+ * 4. Neither the name of the University nor the names of its contributors
+ *    may be used to endorse or promote products derived from this software
+ *    without specific prior written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
+ */
+
+#ifndef lint
+static char sccsid[] = "@(#)cabs.c	5.6 (Berkeley) 10/9/90";
+#endif /* not lint */
+
+/* HYPOT(X,Y)
+ * RETURN THE SQUARE ROOT OF X^2 + Y^2  WHERE Z=X+iY
+ * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
+ * CODED IN C BY K.C. NG, 11/28/84; 
+ * REVISED BY K.C. NG, 7/12/85.
+ *
+ * Required system supported functions :
+ *	copysign(x,y)
+ *	finite(x)
+ *	scalb(x,N)
+ *	sqrt(x)
+ *
+ * Method :
+ *	1. replace x by |x| and y by |y|, and swap x and
+ *	   y if y > x (hence x is never smaller than y).
+ *	2. Hypot(x,y) is computed by:
+ *	   Case I, x/y > 2
+ *		
+ *				       y
+ *		hypot = x + -----------------------------
+ *			 		    2
+ *			    sqrt ( 1 + [x/y]  )  +  x/y
+ *
+ *	   Case II, x/y <= 2 
+ *				                   y
+ *		hypot = x + --------------------------------------------------
+ *				          		     2 
+ *				     			[x/y]   -  2
+ *			   (sqrt(2)+1) + (x-y)/y + -----------------------------
+ *			 		    			  2
+ *			    			  sqrt ( 1 + [x/y]  )  + sqrt(2)
+ *
+ *
+ *
+ * Special cases:
+ *	hypot(x,y) is INF if x or y is +INF or -INF; else
+ *	hypot(x,y) is NAN if x or y is NAN.
+ *
+ * Accuracy:
+ * 	hypot(x,y) returns the sqrt(x^2+y^2) with error less than 1 ulps (units
+ *	in the last place). See Kahan's "Interval Arithmetic Options in the
+ *	Proposed IEEE Floating Point Arithmetic Standard", Interval Mathematics
+ *      1980, Edited by Karl L.E. Nickel, pp 99-128. (A faster but less accurate
+ *	code follows in	comments.) In a test run with 500,000 random arguments
+ *	on a VAX, the maximum observed error was .959 ulps.
+ *
+ * Constants:
+ * The hexadecimal values are the intended ones for the following constants.
+ * The decimal values may be used, provided that the compiler will convert
+ * from decimal to binary accurately enough to produce the hexadecimal values
+ * shown.
+ */
+#include "mathimpl.h"
+
+vc(r2p1hi, 2.4142135623730950345E0   ,8279,411a,ef32,99fc,   2, .9A827999FCEF32)
+vc(r2p1lo, 1.4349369327986523769E-17 ,597d,2484,754b,89b3, -55, .84597D89B3754B)
+vc(sqrt2,  1.4142135623730950622E0   ,04f3,40b5,de65,33f9,   1, .B504F333F9DE65)
+
+ic(r2p1hi, 2.4142135623730949234E0   ,   1, 1.3504F333F9DE6)
+ic(r2p1lo, 1.2537167179050217666E-16 , -53, 1.21165F626CDD5)
+ic(sqrt2,  1.4142135623730951455E0   ,   0, 1.6A09E667F3BCD)
+
+#ifdef vccast
+#define	r2p1hi	vccast(r2p1hi)
+#define	r2p1lo	vccast(r2p1lo)
+#define	sqrt2	vccast(sqrt2)
+#endif
+
+double
+hypot(x,y)
+double x, y;
+{
+	static const double zero=0, one=1, 
+		      small=1.0E-18;	/* fl(1+small)==1 */
+	static const ibig=30;	/* fl(1+2**(2*ibig))==1 */
+	double t,r;
+	int exp;
+
+	if(finite(x))
+	    if(finite(y))
+	    {	
+		x=copysign(x,one);
+		y=copysign(y,one);
+		if(y > x) 
+		    { t=x; x=y; y=t; }
+		if(x == zero) return(zero);
+		if(y == zero) return(x);
+		exp= logb(x);
+		if(exp-(int)logb(y) > ibig ) 	
+			/* raise inexact flag and return |x| */
+		   { one+small; return(x); }
+
+	    /* start computing sqrt(x^2 + y^2) */
+		r=x-y;
+		if(r>y) { 	/* x/y > 2 */
+		    r=x/y;
+		    r=r+sqrt(one+r*r); }
+		else {		/* 1 <= x/y <= 2 */
+		    r/=y; t=r*(r+2.0);
+		    r+=t/(sqrt2+sqrt(2.0+t));
+		    r+=r2p1lo; r+=r2p1hi; }
+
+		r=y/r;
+		return(x+r);
+
+	    }
+
+	    else if(y==y)   	   /* y is +-INF */
+		     return(copysign(y,one));
+	    else 
+		     return(y);	   /* y is NaN and x is finite */
+
+	else if(x==x) 		   /* x is +-INF */
+	         return (copysign(x,one));
+	else if(finite(y))
+	         return(x);		   /* x is NaN, y is finite */
+#if !defined(vax)&&!defined(tahoe)
+	else if(y!=y) return(y);  /* x and y is NaN */
+#endif	/* !defined(vax)&&!defined(tahoe) */
+	else return(copysign(y,one));   /* y is INF */
+}
+
+/* CABS(Z)
+ * RETURN THE ABSOLUTE VALUE OF THE COMPLEX NUMBER  Z = X + iY
+ * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
+ * CODED IN C BY K.C. NG, 11/28/84.
+ * REVISED BY K.C. NG, 7/12/85.
+ *
+ * Required kernel function :
+ *	hypot(x,y)
+ *
+ * Method :
+ *	cabs(z) = hypot(x,y) .
+ */
+
+double
+cabs(z)
+struct __cabs_complex z;
+{
+	return hypot(z.__x,z.__y);
+}
+
+double
+z_abs(z)
+struct __cabs_complex *z;
+{
+	return hypot(z->__x,z->__y);
+}
+
+/* A faster but less accurate version of cabs(x,y) */
+#if 0
+double hypot(x,y)
+double x, y;
+{
+	static const double zero=0, one=1;
+		      small=1.0E-18;	/* fl(1+small)==1 */
+	static const ibig=30;	/* fl(1+2**(2*ibig))==1 */
+	double temp;
+	int exp;
+
+	if(finite(x))
+	    if(finite(y))
+	    {	
+		x=copysign(x,one);
+		y=copysign(y,one);
+		if(y > x) 
+		    { temp=x; x=y; y=temp; }
+		if(x == zero) return(zero);
+		if(y == zero) return(x);
+		exp= logb(x);
+		x=scalb(x,-exp);
+		if(exp-(int)logb(y) > ibig ) 
+			/* raise inexact flag and return |x| */
+		   { one+small; return(scalb(x,exp)); }
+		else y=scalb(y,-exp);
+		return(scalb(sqrt(x*x+y*y),exp));
+	    }
+
+	    else if(y==y)   	   /* y is +-INF */
+		     return(copysign(y,one));
+	    else 
+		     return(y);	   /* y is NaN and x is finite */
+
+	else if(x==x) 		   /* x is +-INF */
+	         return (copysign(x,one));
+	else if(finite(y))
+	         return(x);		   /* x is NaN, y is finite */
+	else if(y!=y) return(y);  	/* x and y is NaN */
+	else return(copysign(y,one));   /* y is INF */
+}
+#endif
diff --git a/sysdeps/ieee754/cbrt.c b/sysdeps/ieee754/cbrt.c
new file mode 100644
index 0000000000..fe5fb95511
--- /dev/null
+++ b/sysdeps/ieee754/cbrt.c
@@ -0,0 +1,120 @@
+/*
+ * Copyright (c) 1985, 1993
+ *	The Regents of the University of California.  All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ * 3. All advertising materials mentioning features or use of this software
+ *    must display the following acknowledgement:
+ *	This product includes software developed by the University of
+ *	California, Berkeley and its contributors.
+ * 4. Neither the name of the University nor the names of its contributors
+ *    may be used to endorse or promote products derived from this software
+ *    without specific prior written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
+ */
+
+#ifndef lint
+static char sccsid[] = "@(#)cbrt.c	8.1 (Berkeley) 6/4/93";
+#endif /* not lint */
+
+#include <sys/cdefs.h>
+
+/* kahan's cube root (53 bits IEEE double precision)
+ * for IEEE machines only
+ * coded in C by K.C. Ng, 4/30/85
+ *
+ * Accuracy:
+ *	better than 0.667 ulps according to an error analysis. Maximum
+ * error observed was 0.666 ulps in an 1,000,000 random arguments test.
+ *
+ * Warning: this code is semi machine dependent; the ordering of words in
+ * a floating point number must be known in advance. I assume that the
+ * long interger at the address of a floating point number will be the
+ * leading 32 bits of that floating point number (i.e., sign, exponent,
+ * and the 20 most significant bits).
+ * On a National machine, it has different ordering; therefore, this code 
+ * must be compiled with flag -DNATIONAL. 
+ */
+#if !defined(vax)&&!defined(tahoe)
+
+static const unsigned long
+		     B1 = 715094163, /* B1 = (682-0.03306235651)*2**20 */
+	             B2 = 696219795; /* B2 = (664-0.03306235651)*2**20 */
+static const double
+	    C= 19./35.,
+	    D= -864./1225.,
+	    E= 99./70.,
+	    F= 45./28.,
+	    G= 5./14.;
+
+double cbrt(x) 
+double x;
+{
+	double r,s,t=0.0,w;
+	unsigned long *px = (unsigned long *) &x,
+	              *pt = (unsigned long *) &t,
+		      mexp,sign;
+
+#ifdef national /* ordering of words in a floating points number */
+	const int n0=1,n1=0;
+#else	/* national */
+	const int n0=0,n1=1;
+#endif	/* national */
+
+	mexp=px[n0]&0x7ff00000;
+	if(mexp==0x7ff00000) return(x); /* cbrt(NaN,INF) is itself */
+	if(x==0.0) return(x);		/* cbrt(0) is itself */
+
+	sign=px[n0]&0x80000000; /* sign= sign(x) */
+	px[n0] ^= sign;		/* x=|x| */
+
+
+    /* rough cbrt to 5 bits */
+	if(mexp==0) 		/* subnormal number */
+	  {pt[n0]=0x43500000; 	/* set t= 2**54 */
+	   t*=x; pt[n0]=pt[n0]/3+B2;
+	  }
+	else
+	  pt[n0]=px[n0]/3+B1;	
+
+
+    /* new cbrt to 23 bits, may be implemented in single precision */
+	r=t*t/x;
+	s=C+r*t;
+	t*=G+F/(s+E+D/s);	
+
+    /* chopped to 20 bits and make it larger than cbrt(x) */ 
+	pt[n1]=0; pt[n0]+=0x00000001;
+
+
+    /* one step newton iteration to 53 bits with error less than 0.667 ulps */
+	s=t*t;		/* t*t is exact */
+	r=x/s;
+	w=t+t;
+	r=(r-t)/(w+r);	/* r-t is exact */
+	t=t+t*r;
+
+
+    /* retore the sign bit */
+	pt[n0] |= sign;
+	return(t);
+}
+#endif
diff --git a/sysdeps/ieee754/copysign.c b/sysdeps/ieee754/copysign.c
new file mode 100644
index 0000000000..f1f0591d5f
--- /dev/null
+++ b/sysdeps/ieee754/copysign.c
@@ -0,0 +1,37 @@
+/* Copyright (C) 1991, 1992, 1995 Free Software Foundation, Inc.
+This file is part of the GNU C Library.
+
+The GNU C Library 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.
+
+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
+Library General Public License for more details.
+
+You should have received a copy of the GNU Library General Public
+License along with the GNU C Library; see the file COPYING.LIB.  If
+not, write to the Free Software Foundation, Inc., 675 Mass Ave,
+Cambridge, MA 02139, USA.  */
+
+#include <ansidecl.h>
+#include <math.h>
+#include "ieee754.h"
+
+/* Return X with its signed changed to Y's.  */
+double
+DEFUN(__copysign, (x, y), double x AND double y)
+{
+  union ieee754_double ux, uy;
+
+  ux.d = x;
+  uy.d = y;
+
+  ux.ieee.negative = uy.ieee.negative;
+
+  return ux.d;
+}
+
+weak_alias (__copysign, copysign)
diff --git a/sysdeps/ieee754/dbl2mpn.c b/sysdeps/ieee754/dbl2mpn.c
new file mode 100644
index 0000000000..6b690f711e
--- /dev/null
+++ b/sysdeps/ieee754/dbl2mpn.c
@@ -0,0 +1,107 @@
+/* Copyright (C) 1993, 1994, 1995 Free Software Foundation, Inc.
+This file is part of the GNU C Library.
+
+The GNU C Library 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.
+
+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
+Library General Public License for more details.
+
+You should have received a copy of the GNU Library General Public
+License along with the GNU C Library; see the file COPYING.LIB.  If
+not, write to the Free Software Foundation, Inc., 675 Mass Ave,
+Cambridge, MA 02139, USA.  */
+
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "longlong.h"
+#include "ieee754.h"
+#include <float.h>
+#include <stdlib.h>
+
+/* Convert a `double' in IEEE754 standard double-precision format to a
+   multi-precision integer representing the significand scaled up by its
+   number of bits (52 for double) and an integral power of two (MPN frexp). */
+
+mp_size_t
+__mpn_extract_double (mp_ptr res_ptr, mp_size_t size,
+		      int *expt, int *is_neg,
+		      double value)
+{
+  union ieee754_double u;
+  u.d = value;
+
+  *is_neg = u.ieee.negative;
+  *expt = (int) u.ieee.exponent - IEEE754_DOUBLE_BIAS;
+
+#if BITS_PER_MP_LIMB == 32
+  res_ptr[0] = u.ieee.mantissa1; /* Low-order 32 bits of fraction.  */
+  res_ptr[1] = u.ieee.mantissa0; /* High-order 20 bits.  */
+  #define N 2
+#elif BITS_PER_MP_LIMB == 64
+  /* Hopefully the compiler will combine the two bitfield extracts
+     and this composition into just the original quadword extract.  */
+  res_ptr[0] = (u.ieee.mantissa0 << 32) | u.ieee.mantissa1;
+  #define N 1
+#else
+  #error "mp_limb size " BITS_PER_MP_LIMB "not accounted for"
+#endif
+/* The format does not fill the last limb.  There are some zeros.  */
+#define NUM_LEADING_ZEROS (BITS_PER_MP_LIMB \
+			   - (DBL_MANT_DIG - ((N - 1) * BITS_PER_MP_LIMB)))
+
+  if (u.ieee.exponent == 0)
+    {
+      /* A biased exponent of zero is a special case.
+	 Either it is a zero or it is a denormal number.  */
+      if (res_ptr[0] == 0 && res_ptr[N - 1] == 0) /* Assumes N<=2.  */
+	/* It's zero.  */
+	*expt = 0;
+      else
+	{
+          /* It is a denormal number, meaning it has no implicit leading
+	     one bit, and its exponent is in fact the format minimum.  */
+	  int cnt;
+
+	  if (res_ptr[N - 1] != 0)
+	    {
+	      count_leading_zeros (cnt, res_ptr[N - 1]);
+	      cnt -= NUM_LEADING_ZEROS;
+#if N == 2
+	      res_ptr[N - 1] = res_ptr[1] << cnt
+			       | (N - 1)
+			         * (res_ptr[0] >> (BITS_PER_MP_LIMB - cnt));
+	      res_ptr[0] <<= cnt;
+#else
+	      res_ptr[N - 1] <<= cnt;
+#endif
+	      *expt = DBL_MIN_EXP - 1 - cnt;
+	    }
+	  else
+	    {
+	      count_leading_zeros (cnt, res_ptr[0]);
+	      if (cnt >= NUM_LEADING_ZEROS)
+		{
+		  res_ptr[N - 1] = res_ptr[0] << (cnt - NUM_LEADING_ZEROS);
+		  res_ptr[0] = 0; 
+		}
+	      else
+		{
+		  res_ptr[N - 1] = res_ptr[0] >> (NUM_LEADING_ZEROS - cnt);
+		  res_ptr[0] <<= BITS_PER_MP_LIMB - (NUM_LEADING_ZEROS - cnt);
+		}
+	      *expt = DBL_MIN_EXP - 1
+		      - (BITS_PER_MP_LIMB - NUM_LEADING_ZEROS) - cnt;
+	    }
+	}
+    }
+  else
+    /* Add the implicit leading one bit for a normalized number.  */
+    res_ptr[N - 1] |= 1 << (DBL_MANT_DIG - 1 - ((N - 1) * BITS_PER_MP_LIMB));
+
+  return N;
+}
diff --git a/sysdeps/ieee754/drem.c b/sysdeps/ieee754/drem.c
new file mode 100644
index 0000000000..cab3a04535
--- /dev/null
+++ b/sysdeps/ieee754/drem.c
@@ -0,0 +1,107 @@
+/* Copyright (C) 1992, 1995 Free Software Foundation, Inc.
+This file is part of the GNU C Library.
+
+The GNU C Library 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.
+
+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
+Library General Public License for more details.
+
+You should have received a copy of the GNU Library General Public
+License along with the GNU C Library; see the file COPYING.LIB.  If
+not, write to the Free Software Foundation, Inc., 675 Mass Ave,
+Cambridge, MA 02139, USA.  */
+
+/*
+ * Copyright (c) 1985 Regents of the University of California.
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms are permitted provided
+ * that: (1) source distributions retain this entire copyright notice and
+ * comment, and (2) distributions including binaries display the following
+ * acknowledgement:  ``This product includes software developed by the
+ * University of California, Berkeley and its contributors'' in the
+ * documentation or other materials provided with the distribution and in
+ * all advertising materials mentioning features or use of this software.
+ * Neither the name of the University nor the names of its contributors may
+ * be used to endorse or promote products derived from this software without
+ * specific prior written permission.
+ * THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR IMPLIED
+ * WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED WARRANTIES OF
+ * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+ */
+
+#include <ansidecl.h>
+#include <math.h>
+#include <float.h>
+#include "ieee754.h"
+
+/* Return the remainder of X/Y.  */
+double
+DEFUN(__drem, (x, y),
+      double x AND double y)
+{
+  union ieee754_double ux, uy;
+
+  ux.d = x;
+  uy.d = y;
+#define x ux.d
+#define	y uy.d
+
+  uy.ieee.negative = 0;
+
+  if (!__finite (x) || y == 0.0)
+    return NAN;
+  else if (__isnan (y))
+    return y;
+  else if (__isinf (y))
+    return x;
+  else if (uy.ieee.exponent <= 1)
+    {
+      /* Subnormal (or almost subnormal) Y value.  */
+      double b = __scalb (1.0, 54);
+      y *= b;
+      x = __drem (x, y);
+      x *= b;
+      return __drem (x, y) / b;
+    }
+  else if (y >= 1.7e308 / 2)
+    {
+      y /= 2;
+      x /= 2;
+      return __drem (x, y) * 2;
+    }
+  else
+    {
+      union ieee754_double a;
+      double b;
+      unsigned int negative = ux.ieee.negative;
+      a.d = y + y;
+      b = y / 2;
+      ux.ieee.negative = 0;
+      while (x > a.d)
+	{
+	  unsigned short int k = ux.ieee.exponent - a.ieee.exponent;
+	  union ieee754_double tmp;
+	  tmp.d = a.d;
+	  tmp.ieee.exponent += k;
+	  if (x < tmp.d)
+	    --tmp.ieee.exponent;
+	  x -= tmp.d;
+	}
+      if (x > b)
+	{
+	  x -= y;
+	  if (x >= b)
+	    x -= y;
+	}
+      ux.ieee.negative ^= negative;
+      return x;
+    }
+}
+
+weak_alias (__drem, drem)
diff --git a/sysdeps/ieee754/frexp.c b/sysdeps/ieee754/frexp.c
new file mode 100644
index 0000000000..c56a17f6c1
--- /dev/null
+++ b/sysdeps/ieee754/frexp.c
@@ -0,0 +1,42 @@
+/* Copyright (C) 1993 Free Software Foundation, Inc.
+This file is part of the GNU C Library.
+
+The GNU C Library 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.
+
+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
+Library General Public License for more details.
+
+You should have received a copy of the GNU Library General Public
+License along with the GNU C Library; see the file COPYING.LIB.  If
+not, write to the Free Software Foundation, Inc., 675 Mass Ave,
+Cambridge, MA 02139, USA.  */
+
+#include <math.h>
+#include "ieee754.h"
+
+/* Break VALUE into a normalized fraction and an integral power of 2.  */
+
+double
+frexp (value, exp)
+     double value;
+     int *exp;
+{
+  if (value == 0)
+    {
+      *exp = 0;
+      return 0;
+    }
+  else
+    {
+      union ieee754_double u;
+      u.d = value;
+      *exp = u.ieee.exponent - 1022;
+      u.ieee.exponent = 1022;
+      return u.d;
+    }
+}
diff --git a/sysdeps/ieee754/huge_val.h b/sysdeps/ieee754/huge_val.h
new file mode 100644
index 0000000000..183f213b83
--- /dev/null
+++ b/sysdeps/ieee754/huge_val.h
@@ -0,0 +1,87 @@
+/* `HUGE_VAL' constant for IEEE 754 machines (where it is infinity).
+   Used by <stdlib.h> and <math.h> functions for overflow.
+
+Copyright (C) 1992, 1995 Free Software Foundation, Inc.
+This file is part of the GNU C Library.
+
+The GNU C Library 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.
+
+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
+Library General Public License for more details.
+
+You should have received a copy of the GNU Library General Public
+License along with the GNU C Library; see the file COPYING.LIB.  If
+not, write to the Free Software Foundation, Inc., 675 Mass Ave,
+Cambridge, MA 02139, USA.  */
+
+#ifndef	   _HUGE_VAL_H
+#define	   _HUGE_VAL_H	1
+
+#include <features.h>
+#include <sys/cdefs.h>
+#include <endian.h>
+
+/* IEEE positive infinity (-HUGE_VAL is negative infinity).  */
+
+#if __BYTE_ORDER == __BIG_ENDIAN
+#define	__HUGE_VAL_bytes	{ 0x7f, 0xf0, 0, 0, 0, 0, 0, 0 }
+#endif
+#if __BYTE_ORDER == __LITTLE_ENDIAN
+#define	__HUGE_VAL_bytes	{ 0, 0, 0, 0, 0, 0, 0xf0, 0x7f }
+#endif
+
+#define __huge_val_t	union { unsigned char __c[8]; double __d; }
+#ifdef	__GNUC__
+#define	HUGE_VAL	(__extension__ \
+			 ((__huge_val_t) { __c: __HUGE_VAL_bytes }).__d)
+#else	/* Not GCC.  */
+static __huge_val_t __huge_val = { __HUGE_VAL_bytes };
+#define	HUGE_VAL	(__huge_val.__d)
+#endif	/* GCC.  */
+
+
+/* GNU extensions: (float) HUGE_VALf and (long double) HUGE_VALl.  */
+
+#ifdef	__USE_GNU
+
+#if __BYTE_ORDER == __BIG_ENDIAN
+#define	__HUGE_VALf_bytes	{ 0x7f, 0x80, 0, 0 }
+#endif
+#if __BYTE_ORDER == __LITTLE_ENDIAN
+#define	__HUGE_VALf_bytes	{ 0, 0, 0x80, 0x7f }
+#endif
+
+#define __huge_valf_t	union { unsigned char __c[4]; float __f; }
+#ifdef	__GNUC__
+#define	HUGE_VALf	(__extension__ \
+			 ((__huge_valf_t) { __c: __HUGE_VALf_bytes }).__f)
+#else	/* Not GCC.  */
+static __huge_valf_t __huge_valf = { __HUGE_VALf_bytes };
+#define	HUGE_VALf	(__huge_valf.__f)
+#endif	/* GCC.  */
+
+#if __BYTE_ORDER == __BIG_ENDIAN
+#define	__HUGE_VALl_bytes	{ 0x7f, 0xff, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }
+#endif
+#if __BYTE_ORDER == __LITTLE_ENDIAN
+#define	__HUGE_VALl_bytes	{ 0, 0, 0, 0, 0, 0, 0, 0, 0xff, 0x7f, 0, 0 }
+#endif
+
+#define __huge_vall_t	union { unsigned char __c[12]; long double __ld; }
+#ifdef	__GNUC__
+#define	HUGE_VALl	(__extension__ \
+			 ((__huge_vall_t) { __c: __HUGE_VALl_bytes }).__ld)
+#else	/* Not GCC.  */
+static __huge_vall_t __huge_vall = { __HUGE_VALl_bytes };
+#define	HUGE_VALl	(__huge_vall.__ld)
+#endif	/* GCC.  */
+
+#endif	/* __USE_GNU.  */
+
+
+#endif	   /* huge_val.h */
diff --git a/sysdeps/ieee754/ieee754.h b/sysdeps/ieee754/ieee754.h
new file mode 100644
index 0000000000..9cc6f14707
--- /dev/null
+++ b/sysdeps/ieee754/ieee754.h
@@ -0,0 +1,153 @@
+/* Copyright (C) 1992, 1995 Free Software Foundation, Inc.
+This file is part of the GNU C Library.
+
+The GNU C Library 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.
+
+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
+Library General Public License for more details.
+
+You should have received a copy of the GNU Library General Public
+License along with the GNU C Library; see the file COPYING.LIB.  If
+not, write to the Free Software Foundation, Inc., 675 Mass Ave,
+Cambridge, MA 02139, USA.  */
+
+#include <endian.h>
+
+union ieee754_float
+  {
+    float f;
+    
+    /* This is the IEEE 754 single-precision format.  */
+    struct
+      {
+#if	__BYTE_ORDER == __BIG_ENDIAN
+	unsigned int negative:1;
+	unsigned int exponent:8;
+	unsigned int mantissa:23;
+#endif				/* Big endian.  */
+#if	__BYTE_ORDER == __LITTLE_ENDIAN
+	unsigned int mantissa:23;
+	unsigned int exponent:8;
+	unsigned int negative:1;
+#endif				/* Little endian.  */
+      } ieee;
+
+    /* This format makes it easier to see if a NaN is a signalling NaN.  */
+    struct
+      {
+#if	__BYTE_ORDER == __BIG_ENDIAN
+	unsigned int negative:1;
+	unsigned int exponent:8;
+	unsigned int quiet_nan:1;
+	unsigned int mantissa:22;
+#endif				/* Big endian.  */
+#if	__BYTE_ORDER == __LITTLE_ENDIAN
+	unsigned int mantissa:22;
+	unsigned int quiet_nan:1;
+	unsigned int exponent:8;
+	unsigned int negative:1;
+#endif				/* Little endian.  */
+      } ieee_nan;
+  };
+
+#define IEEE754_FLOAT_BIAS	0x7f /* Added to exponent.  */
+
+
+union ieee754_double
+  {
+    double d;
+    
+    /* This is the IEEE 754 double-precision format.  */
+    struct
+      {
+#if	__BYTE_ORDER == __BIG_ENDIAN
+	unsigned int negative:1;
+	unsigned int exponent:11;
+	/* Together these comprise the mantissa.  */
+	unsigned int mantissa0:20;
+	unsigned int mantissa1:32;
+#endif				/* Big endian.  */
+#if	__BYTE_ORDER == __LITTLE_ENDIAN
+	/* Together these comprise the mantissa.  */
+	unsigned int mantissa1:32;
+	unsigned int mantissa0:20;
+	unsigned int exponent:11;
+	unsigned int negative:1;
+#endif				/* Little endian.  */
+      } ieee;
+
+    /* This format makes it easier to see if a NaN is a signalling NaN.  */
+    struct
+      {
+#if	__BYTE_ORDER == __BIG_ENDIAN
+	unsigned int negative:1;
+	unsigned int exponent:11;
+	unsigned int quiet_nan:1;
+	/* Together these conprise the mantissa.  */
+	unsigned int mantissa0:19;
+	unsigned int mantissa1:32;
+#else
+	/* Together these conprise the mantissa.  */
+	unsigned int mantissa1:32;
+	unsigned int mantissa0:19;
+	unsigned int quiet_nan:1;
+	unsigned int exponent:11;
+	unsigned int negative:1;
+#endif
+      } ieee_nan;
+  };
+
+#define IEEE754_DOUBLE_BIAS	0x3ff /* Added to exponent.  */
+
+
+union ieee854_long_double
+  {
+    long double d;
+    
+    /* This is the IEEE 854 double-extended-precision format.  */
+    struct
+      {
+#if	__BYTE_ORDER == __BIG_ENDIAN
+	unsigned int negative:1;
+	unsigned int exponent:15;
+	unsigned int empty:16;
+	unsigned int mantissa0:32;
+	unsigned int mantissa1:32;
+#endif
+#if	__BYTE_ORDER == __LITTLE_ENDIAN
+	unsigned int mantissa1:32;
+	unsigned int mantissa0:32;
+	unsigned int exponent:15;
+	unsigned int negative:1;
+	unsigned int empty:16;
+#endif
+      } ieee;
+    
+    /* This is for NaNs in the IEEE 854 double-extended-precision format.  */
+    struct
+      {
+#if	__BYTE_ORDER == __BIG_ENDIAN
+	unsigned int negative:1;
+	unsigned int exponent:15;
+	unsigned int empty:16;
+	unsigned int quiet_nan:1;
+	unsigned int mantissa0:31;
+	unsigned int mantissa1:32;
+#endif
+#if	__BYTE_ORDER == __LITTLE_ENDIAN
+	unsigned int mantissa1:32;
+	unsigned int mantissa0:31;
+	unsigned int quiet_nan:1;
+	unsigned int exponent:15;
+	unsigned int negative:1;
+	unsigned int empty:16;
+#endif
+      } ieee_nan;
+  };
+
+#define IEEE854_LONG_DOUBLE_BIAS 0x3fff
diff --git a/sysdeps/ieee754/infnan.c b/sysdeps/ieee754/infnan.c
new file mode 100644
index 0000000000..89ab5b6e08
--- /dev/null
+++ b/sysdeps/ieee754/infnan.c
@@ -0,0 +1,50 @@
+/* Copyright (C) 1991, 1992, 1995 Free Software Foundation, Inc.
+This file is part of the GNU C Library.
+
+The GNU C Library 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.
+
+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
+Library General Public License for more details.
+
+You should have received a copy of the GNU Library General Public
+License along with the GNU C Library; see the file COPYING.LIB.  If
+not, write to the Free Software Foundation, Inc., 675 Mass Ave,
+Cambridge, MA 02139, USA.  */
+
+#include <ansidecl.h>
+#include <errno.h>
+#include <math.h>
+#include <float.h>
+#include "ieee754.h"
+
+/* Deal with an infinite or NaN result.
+   If ERROR is ERANGE, result is +Inf;
+   if ERROR is - ERANGE, result is -Inf;
+   otherwise result is NaN.
+   This will set `errno' to either ERANGE or EDOM,
+   and may return an infinity or NaN, or may do something else.  */
+double
+DEFUN(__infnan, (error), int error)
+{
+  switch (error)
+    {
+    case ERANGE:
+      errno = ERANGE;
+      return HUGE_VAL;
+
+    case - ERANGE:
+      errno = ERANGE;
+      return - HUGE_VAL;
+
+    default:
+      errno = EDOM;
+      return NAN;
+    }
+}
+
+weak_alias (__infnan, infnan)
diff --git a/sysdeps/ieee754/isinf.c b/sysdeps/ieee754/isinf.c
new file mode 100644
index 0000000000..79be916432
--- /dev/null
+++ b/sysdeps/ieee754/isinf.c
@@ -0,0 +1,40 @@
+/* Copyright (C) 1991, 1992, 1995 Free Software Foundation, Inc.
+This file is part of the GNU C Library.
+
+The GNU C Library 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.
+
+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
+Library General Public License for more details.
+
+You should have received a copy of the GNU Library General Public
+License along with the GNU C Library; see the file COPYING.LIB.  If
+not, write to the Free Software Foundation, Inc., 675 Mass Ave,
+Cambridge, MA 02139, USA.  */
+
+#include <ansidecl.h>
+#include <math.h>
+#include "ieee754.h"
+
+/* Return 0 if VALUE is finite or NaN, +1 if it
+   is +Infinity, -1 if it is -Infinity.  */
+int
+DEFUN(__isinf, (value), double value)
+{
+  union ieee754_double u;
+
+  u.d = value;
+  /* An IEEE 754 infinity has an exponent with the
+     maximum possible value and a zero mantissa.  */
+  if ((u.ieee.exponent & 0x7ff) == 0x7ff &&
+      u.ieee.mantissa0 == 0 && u.ieee.mantissa1 == 0)
+    return u.ieee.negative ? -1 : 1;
+
+  return 0;
+}
+
+weak_alias (__isinf, isinf)
diff --git a/sysdeps/ieee754/isinfl.c b/sysdeps/ieee754/isinfl.c
new file mode 100644
index 0000000000..ad03728ad8
--- /dev/null
+++ b/sysdeps/ieee754/isinfl.c
@@ -0,0 +1,44 @@
+/* Copyright (C) 1991, 1992, 1995 Free Software Foundation, Inc.
+This file is part of the GNU C Library.
+
+The GNU C Library 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.
+
+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
+Library General Public License for more details.
+
+You should have received a copy of the GNU Library General Public
+License along with the GNU C Library; see the file COPYING.LIB.  If
+not, write to the Free Software Foundation, Inc., 675 Mass Ave,
+Cambridge, MA 02139, USA.  */
+
+#include <math.h>
+#include "ieee754.h"
+
+#undef __isinfl
+#undef isinfl
+
+
+/* Return 0 if VALUE is finite or NaN, +1 if it
+   is +Infinity, -1 if it is -Infinity.  */
+int
+__isinfl (long double value)
+{
+  union ieee854_long_double u;
+
+  u.d = value;
+
+  /* An IEEE 854 infinity has an exponent with the
+     maximum possible value and a zero mantissa.  */
+  if ((u.ieee.exponent & 0x7fff) == 0x7fff &&
+      u.ieee.mantissa0 == 0 && u.ieee.mantissa1 == 0)
+    return u.ieee.negative ? -1 : 1;
+
+  return 0;
+}
+
+weak_alias (__isinfl, isinfl);
diff --git a/sysdeps/ieee754/isnan.c b/sysdeps/ieee754/isnan.c
new file mode 100644
index 0000000000..8a537cc0dd
--- /dev/null
+++ b/sysdeps/ieee754/isnan.c
@@ -0,0 +1,37 @@
+/* Copyright (C) 1991, 1992, 1995 Free Software Foundation, Inc.
+This file is part of the GNU C Library.
+
+The GNU C Library 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.
+
+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
+Library General Public License for more details.
+
+You should have received a copy of the GNU Library General Public
+License along with the GNU C Library; see the file COPYING.LIB.  If
+not, write to the Free Software Foundation, Inc., 675 Mass Ave,
+Cambridge, MA 02139, USA.  */
+
+#include <ansidecl.h>
+#include <math.h>
+#include "ieee754.h"
+
+/* Return nonzero if VALUE is not a number.  */
+int
+DEFUN(__isnan, (value), double value)
+{
+  union ieee754_double u;
+
+  u.d = value;
+
+  /* IEEE 754 NaN's have the maximum possible
+     exponent and a nonzero mantissa.  */
+  return ((u.ieee.exponent & 0x7ff) == 0x7ff &&
+	  (u.ieee.mantissa0 != 0 || u.ieee.mantissa1 != 0));
+}
+
+weak_alias (__isnan, isnan)
diff --git a/sysdeps/ieee754/isnanl.c b/sysdeps/ieee754/isnanl.c
new file mode 100644
index 0000000000..b06b31a3d2
--- /dev/null
+++ b/sysdeps/ieee754/isnanl.c
@@ -0,0 +1,40 @@
+/* Copyright (C) 1991, 1992, 1995 Free Software Foundation, Inc.
+This file is part of the GNU C Library.
+
+The GNU C Library 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.
+
+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
+Library General Public License for more details.
+
+You should have received a copy of the GNU Library General Public
+License along with the GNU C Library; see the file COPYING.LIB.  If
+not, write to the Free Software Foundation, Inc., 675 Mass Ave,
+Cambridge, MA 02139, USA.  */
+
+#include <math.h>
+#include "ieee754.h"
+
+#undef __isnanl
+#undef isnanl
+
+
+/* Return nonzero if VALUE is not a number.  */
+int
+__isnanl (long double value)
+{
+  union ieee854_long_double u;
+
+  u.d = value;
+
+  /* IEEE 854 NaN's have the maximum possible
+     exponent and a nonzero mantissa.  */
+  return ((u.ieee.exponent & 0x7fff) == 0x7fff &&
+	  (u.ieee.mantissa0 != 0 || u.ieee.mantissa1 != 0));
+}
+
+weak_alias (__isnanl, isnanl);
diff --git a/sysdeps/ieee754/ldbl2mpn.c b/sysdeps/ieee754/ldbl2mpn.c
new file mode 100644
index 0000000000..29e8a3b2ce
--- /dev/null
+++ b/sysdeps/ieee754/ldbl2mpn.c
@@ -0,0 +1,93 @@
+/* Copyright (C) 1995 Free Software Foundation, Inc.
+This file is part of the GNU C Library.
+
+The GNU C Library 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.
+
+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
+Library General Public License for more details.
+
+You should have received a copy of the GNU Library General Public
+License along with the GNU C Library; see the file COPYING.LIB.  If
+not, write to the Free Software Foundation, Inc., 675 Mass Ave,
+Cambridge, MA 02139, USA.  */
+
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "longlong.h"
+#include "ieee754.h"
+#include <float.h>
+#include <stdlib.h>
+
+/* Convert a `long double' in IEEE854 standard double-precision format to a
+   multi-precision integer representing the significand scaled up by its
+   number of bits (64 for long double) and an integral power of two
+   (MPN frexpl). */
+
+mp_size_t
+__mpn_extract_long_double (mp_ptr res_ptr, mp_size_t size,
+			   int *expt, int *is_neg,
+			   long double value)
+{
+  union ieee854_long_double u;
+  u.d = value;
+
+  *is_neg = u.ieee.negative;
+  *expt = (int) u.ieee.exponent - IEEE854_LONG_DOUBLE_BIAS;
+
+#if BITS_PER_MP_LIMB == 32
+  res_ptr[0] = u.ieee.mantissa1; /* Low-order 32 bits of fraction.  */
+  res_ptr[1] = u.ieee.mantissa0; /* High-order 32 bits.  */
+  #define N 2
+#elif BITS_PER_MP_LIMB == 64
+  /* Hopefully the compiler will combine the two bitfield extracts
+     and this composition into just the original quadword extract.  */
+  res_ptr[0] = (u.ieee.mantissa0 << 32) | u.ieee.mantissa1;
+  #define N 1
+#else
+  #error "mp_limb size " BITS_PER_MP_LIMB "not accounted for"
+#endif
+
+  if (u.ieee.exponent == 0)
+    {
+      /* A biased exponent of zero is a special case.
+	 Either it is a zero or it is a denormal number.  */
+      if (res_ptr[0] == 0 && res_ptr[N - 1] == 0) /* Assumes N<=2.  */
+	/* It's zero.  */
+	*expt = 0;
+      else
+	{
+          /* It is a denormal number, meaning it has no implicit leading
+  	     one bit, and its exponent is in fact the format minimum.  */
+	  int cnt;
+	  if (res_ptr[N - 1] != 0)
+	    {
+	      count_leading_zeros (cnt, res_ptr[N - 1]);
+	      if (cnt != 0)
+		{
+#if N == 2
+	          res_ptr[N - 1] = res_ptr[N - 1] << cnt
+			           | (res_ptr[0] >> (BITS_PER_MP_LIMB - cnt));
+	          res_ptr[0] <<= cnt;
+#else
+	          res_ptr[N - 1] <<= cnt;
+#endif
+		}
+	      *expt = LDBL_MIN_EXP - 2 - cnt;
+	    }
+	  else
+	    {
+	      count_leading_zeros (cnt, res_ptr[0]);
+	      res_ptr[N - 1] = res_ptr[0] << cnt;
+	      res_ptr[0] = 0;
+	      *expt = LDBL_MIN_EXP - 2 - BITS_PER_MP_LIMB - cnt;
+	    }
+	}
+    }
+
+  return N;
+}
diff --git a/sysdeps/ieee754/ldexp.c b/sysdeps/ieee754/ldexp.c
new file mode 100644
index 0000000000..7e1c74715b
--- /dev/null
+++ b/sysdeps/ieee754/ldexp.c
@@ -0,0 +1,145 @@
+/* Copyright (C) 1992, 1995 Free Software Foundation, Inc.
+This file is part of the GNU C Library.
+
+The GNU C Library 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.
+
+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
+Library General Public License for more details.
+
+You should have received a copy of the GNU Library General Public
+License along with the GNU C Library; see the file COPYING.LIB.  If
+not, write to the Free Software Foundation, Inc., 675 Mass Ave,
+Cambridge, MA 02139, USA.  */
+
+/*
+ * Copyright (c) 1985 Regents of the University of California.
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms are permitted provided
+ * that: (1) source distributions retain this entire copyright notice and
+ * comment, and (2) distributions including binaries display the following
+ * acknowledgement:  ``This product includes software developed by the
+ * University of California, Berkeley and its contributors'' in the
+ * documentation or other materials provided with the distribution and in
+ * all advertising materials mentioning features or use of this software.
+ * Neither the name of the University nor the names of its contributors may
+ * be used to endorse or promote products derived from this software without
+ * specific prior written permission.
+ * THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR IMPLIED
+ * WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED WARRANTIES OF
+ * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+ */
+
+#include <ansidecl.h>
+#include <math.h>
+#include <float.h>
+#include "ieee754.h"
+
+double
+DEFUN(ldexp, (x, exp),
+      double x AND int exp)
+{
+  union ieee754_double u;
+  unsigned int exponent;
+
+  u.d = x;
+#define	x u.d
+
+  exponent = u.ieee.exponent;
+
+  /* The order of the tests is carefully chosen to handle
+     the usual case first, with no branches taken.  */
+
+  if (exponent != 0)
+    {
+      /* X is nonzero and not denormalized.  */
+
+      if (exponent <= DBL_MAX_EXP - DBL_MIN_EXP + 1)
+  	{
+	  /* X is finite.  When EXP < 0, overflow is actually underflow.  */
+
+	  exponent += exp;
+
+	  if (exponent != 0)
+	    {
+	      if (exponent <= DBL_MAX_EXP - DBL_MIN_EXP + 1)
+		{
+		  /* In range.  */
+		  u.ieee.exponent = exponent;
+		  return x;
+		}
+
+	      if (exp >= 0)
+	      overflow:
+		{
+		  CONST int negative = u.ieee.negative;
+		  u.d = HUGE_VAL;
+		  u.ieee.negative = negative;
+		  errno = ERANGE;
+		  return u.d;
+		}
+
+	      if (exponent <= - (unsigned int) (DBL_MANT_DIG + 1))
+		{
+		  /* Underflow.  */
+		  CONST int negative = u.ieee.negative;
+		  u.d = 0.0;
+		  u.ieee.negative = negative;
+		  errno = ERANGE;
+		  return u.d;
+		}
+	    }
+
+	  /* Gradual underflow.  */
+	  u.ieee.exponent = 1;
+	  u.d *= ldexp (1.0, (int) exponent - 1);
+	  if (u.ieee.mantissa0 == 0 && u.ieee.mantissa1 == 0)
+	    /* Underflow.  */
+	    errno = ERANGE;
+	  return u.d;
+  	}
+
+      /* X is +-infinity or NaN.  */
+      if (u.ieee.mantissa0 == 0 && u.ieee.mantissa1 == 0)
+  	{
+	  /* X is +-infinity.  */
+	  if (exp >= 0)
+	    goto overflow;
+	  else
+	    {
+	      /* (infinity * number < 1).  With infinite precision,
+		 (infinity / finite) would be infinity, but otherwise it's
+		 safest to regard (infinity / 2) as indeterminate.  The
+		 infinity might be (2 * finite).  */
+	      CONST int negative = u.ieee.negative;
+	      u.d = NAN;
+	      u.ieee.negative = negative;
+	      errno = EDOM;
+	      return u.d;
+	    }
+	}
+
+      /* X is NaN.  */
+      errno = EDOM;
+      return u.d;
+    }
+
+  /* X is zero or denormalized.  */
+  if (u.ieee.mantissa0 == 0 && u.ieee.mantissa1 == 0)
+    /* X is +-0.0. */
+    return x;
+
+  /* X is denormalized.
+     Multiplying by 2 ** DBL_MANT_DIG normalizes it;
+     we then subtract the DBL_MANT_DIG we added to the exponent.  */
+  return ldexp (x * ldexp (1.0, DBL_MANT_DIG), exp - DBL_MANT_DIG);
+}
+
+/* Compatibility names for the same function.  */
+weak_alias (ldexp, __scalb)
+weak_alias (ldexp, scalb)
diff --git a/sysdeps/ieee754/log10.c b/sysdeps/ieee754/log10.c
new file mode 100644
index 0000000000..da2f5b4d2b
--- /dev/null
+++ b/sysdeps/ieee754/log10.c
@@ -0,0 +1,30 @@
+/* Copyright (C) 1991, 1995 Free Software Foundation, Inc.
+This file is part of the GNU C Library.
+
+The GNU C Library 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.
+
+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
+Library General Public License for more details.
+
+You should have received a copy of the GNU Library General Public
+License along with the GNU C Library; see the file COPYING.LIB.  If
+not, write to the Free Software Foundation, Inc., 675 Mass Ave,
+Cambridge, MA 02139, USA.  */
+
+#include <ansidecl.h>
+#include <errno.h>
+#include <math.h>
+
+/* Return the base-ten logarithm of X.  */
+double
+DEFUN(log10, (x), double x)
+{
+  CONST double inverse_ln10 = 4.3429448190325181667e-1;	/* 1 / log(10) */
+
+  return inverse_ln10 * log(x);
+}
diff --git a/sysdeps/ieee754/logb.c b/sysdeps/ieee754/logb.c
new file mode 100644
index 0000000000..918de9ebad
--- /dev/null
+++ b/sysdeps/ieee754/logb.c
@@ -0,0 +1,48 @@
+/* Copyright (C) 1992, 1995 Free Software Foundation, Inc.
+This file is part of the GNU C Library.
+
+The GNU C Library 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.
+
+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
+Library General Public License for more details.
+
+You should have received a copy of the GNU Library General Public
+License along with the GNU C Library; see the file COPYING.LIB.  If
+not, write to the Free Software Foundation, Inc., 675 Mass Ave,
+Cambridge, MA 02139, USA.  */
+
+#include <ansidecl.h>
+#include <math.h>
+#include <float.h>
+#include "ieee754.h"
+
+/* Return the base 2 signed integral exponent of X.  */
+double
+DEFUN(__logb, (x), double x)
+{
+  union ieee754_double u;
+
+  if (__isnan (x))
+    return x;
+  else if (__isinf (x))
+    return HUGE_VAL;
+  else if (x == 0.0)
+    return - HUGE_VAL;
+
+  u.d = x;
+
+  if (u.ieee.exponent == 0)
+    /* A denormalized number.
+       Multiplying by 2 ** DBL_MANT_DIG normalizes it;
+       we then subtract the DBL_MANT_DIG we added to the exponent.  */
+    return (__logb (x * ldexp (1.0, DBL_MANT_DIG)) - DBL_MANT_DIG);
+
+  return (int) u.ieee.exponent - (DBL_MAX_EXP - 1);
+}
+
+weak_alias (__logb, logb)
diff --git a/sysdeps/ieee754/mpn2dbl.c b/sysdeps/ieee754/mpn2dbl.c
new file mode 100644
index 0000000000..01e0019b8f
--- /dev/null
+++ b/sysdeps/ieee754/mpn2dbl.c
@@ -0,0 +1,46 @@
+/* Copyright (C) 1995 Free Software Foundation, Inc.
+This file is part of the GNU C Library.
+
+The GNU C Library 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.
+
+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
+Library General Public License for more details.
+
+You should have received a copy of the GNU Library General Public
+License along with the GNU C Library; see the file COPYING.LIB.  If
+not, write to the Free Software Foundation, Inc., 675 Mass Ave,
+Cambridge, MA 02139, USA.  */
+
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "ieee754.h"
+#include <float.h>
+
+/* Convert a multi-precision integer of the needed number of bits (53 for
+   double) and an integral power of two to a `double' in IEEE754 double-
+   precision format.  */ 
+
+double
+__mpn_construct_double (mp_srcptr frac_ptr, int expt, int negative)
+{
+  union ieee754_double u;
+  
+  u.ieee.negative = negative;
+  u.ieee.exponent = expt + IEEE754_DOUBLE_BIAS;
+#if BITS_PER_MP_LIMB == 32
+  u.ieee.mantissa1 = frac_ptr[0];
+  u.ieee.mantissa0 = frac_ptr[1] & ((1 << (DBL_MANT_DIG - 32)) - 1);
+#elif BITS_PER_MP_LIMB == 64
+  u.ieee.mantissa1 = frac_ptr[0] & ((1 << 32) - 1);
+  i.ieee.mantissa0 = (frac_ptr[0] >> 32) & ((1 << (DBL_MANT_DIG - 32)) - 1);
+#else
+  #error "mp_limb size " BITS_PER_MP_LIMB "not accounted for"
+#endif
+
+  return u.d;
+}
diff --git a/sysdeps/ieee754/mpn2flt.c b/sysdeps/ieee754/mpn2flt.c
new file mode 100644
index 0000000000..e9bbc490ad
--- /dev/null
+++ b/sysdeps/ieee754/mpn2flt.c
@@ -0,0 +1,42 @@
+/* Copyright (C) 1995 Free Software Foundation, Inc.
+This file is part of the GNU C Library.
+
+The GNU C Library 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.
+
+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
+Library General Public License for more details.
+
+You should have received a copy of the GNU Library General Public
+License along with the GNU C Library; see the file COPYING.LIB.  If
+not, write to the Free Software Foundation, Inc., 675 Mass Ave,
+Cambridge, MA 02139, USA.  */
+
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "ieee754.h"
+#include <float.h>
+
+/* Convert a multi-precision integer of the needed number of bits (24 for
+   float) and an integral power of two to a `float' in IEEE754 single-
+   precision format.  */ 
+
+float
+__mpn_construct_float (mp_srcptr frac_ptr, int expt, int sign)
+{
+  union ieee754_float u;
+  
+  u.ieee.negative = sign;
+  u.ieee.exponent = expt + IEEE754_FLOAT_BIAS;
+#if BITS_PER_MP_LIMB > FLT_MANT_DIG
+  u.ieee.mantissa = frac_ptr[0] & ((1 << FLT_MANT_DIG) - 1);
+#else
+  #error "mp_limb size " BITS_PER_MP_LIMB "not accounted for"
+#endif
+
+  return u.f;
+}
diff --git a/sysdeps/ieee754/mpn2ldbl.c b/sysdeps/ieee754/mpn2ldbl.c
new file mode 100644
index 0000000000..831f2e5f7b
--- /dev/null
+++ b/sysdeps/ieee754/mpn2ldbl.c
@@ -0,0 +1,46 @@
+/* Copyright (C) 1995 Free Software Foundation, Inc.
+This file is part of the GNU C Library.
+
+The GNU C Library 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.
+
+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
+Library General Public License for more details.
+
+You should have received a copy of the GNU Library General Public
+License along with the GNU C Library; see the file COPYING.LIB.  If
+not, write to the Free Software Foundation, Inc., 675 Mass Ave,
+Cambridge, MA 02139, USA.  */
+
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "ieee754.h"
+#include <float.h>
+
+/* Convert a multi-precision integer of the needed number of bits (64 for
+   long double) and an integral power of two to a `long double' in IEEE854
+   extended-precision format.  */ 
+
+long double
+__mpn_construct_long_double (mp_srcptr frac_ptr, int expt, int sign)
+{
+  union ieee854_long_double u;
+  
+  u.ieee.negative = sign;
+  u.ieee.exponent = expt + IEEE854_LONG_DOUBLE_BIAS;
+#if BITS_PER_MP_LIMB == 32
+  u.ieee.mantissa1 = frac_ptr[0];
+  u.ieee.mantissa0 = frac_ptr[1];
+#elif BITS_PER_MP_LIMB == 64
+  u.ieee.mantissa1 = frac_ptr[0] & ((1 << 32) - 1);
+  i.ieee.mantissa0 = frac_ptr[0] >> 32;
+#else
+  #error "mp_limb size " BITS_PER_MP_LIMB "not accounted for"
+#endif
+
+  return u.d;
+}
diff --git a/sysdeps/ieee754/nan.h b/sysdeps/ieee754/nan.h
new file mode 100644
index 0000000000..e34c172dcb
--- /dev/null
+++ b/sysdeps/ieee754/nan.h
@@ -0,0 +1,46 @@
+/* `NAN' constant for IEEE 754 machines.
+
+Copyright (C) 1992 Free Software Foundation, Inc.
+This file is part of the GNU C Library.
+
+The GNU C Library 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.
+
+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
+Library General Public License for more details.
+
+You should have received a copy of the GNU Library General Public
+License along with the GNU C Library; see the file COPYING.LIB.  If
+not, write to the Free Software Foundation, Inc., 675 Mass Ave,
+Cambridge, MA 02139, USA.  */
+
+#ifndef	_NAN_H
+
+#define	_NAN_H	1
+
+/* IEEE Not A Number.  */
+
+#include <endian.h>
+
+#if __BYTE_ORDER == __BIG_ENDIAN
+#define	__nan_bytes		{ 0x7f, 0xf8, 0, 0, 0, 0, 0, 0 }
+#endif
+#if __BYTE_ORDER == __LITTLE_ENDIAN
+#define	__nan_bytes		{ 0, 0, 0, 0, 0, 0, 0xf8, 0x7f }
+#endif
+
+#ifdef	__GNUC__
+#define	NAN \
+  (__extension__ ((union { unsigned char __c[8];			      \
+			   double __d; })				      \
+		  { __nan_bytes }).__d)
+#else	/* Not GCC.  */
+static CONST char __nan[8] = __nan_bytes;
+#define	NAN	(*(CONST double *) __nan)
+#endif	/* GCC.  */
+
+#endif	/* nan.h */
diff --git a/sysdeps/ieee754/sqrt.c b/sysdeps/ieee754/sqrt.c
new file mode 100644
index 0000000000..7e350e0d91
--- /dev/null
+++ b/sysdeps/ieee754/sqrt.c
@@ -0,0 +1,120 @@
+/*
+ * Copyright (c) 1985 Regents of the University of California.
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms are permitted provided
+ * that: (1) source distributions retain this entire copyright notice and
+ * comment, and (2) distributions including binaries display the following
+ * acknowledgement:  ``This product includes software developed by the
+ * University of California, Berkeley and its contributors'' in the
+ * documentation or other materials provided with the distribution and in
+ * all advertising materials mentioning features or use of this software.
+ * Neither the name of the University nor the names of its contributors may
+ * be used to endorse or promote products derived from this software without
+ * specific prior written permission.
+ * THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR IMPLIED
+ * WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED WARRANTIES OF
+ * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+ */
+
+#include <ansidecl.h>
+#include <errno.h>
+#include <math.h>
+
+/* Return the square root of X.  */
+double
+DEFUN (sqrt, (x), double x)
+{
+  double q, s, b, r, t;
+  CONST double zero = 0.0;
+  int m, n, i;
+
+  /* sqrt (NaN) is NaN; sqrt (+-0) is +-0.  */
+  if (__isnan (x) || x == zero)
+    return x;
+
+  if (x < zero)
+    return zero / zero;
+
+  /* sqrt (Inf) is Inf.  */
+  if (__isinf (x))
+    return x;
+
+  /* Scale X to [1,4).  */
+  n = __logb (x);
+  x = __scalb (x, -n);
+  m = __logb (x);
+  if (m != 0)
+    /* Subnormal number.  */
+    x = __scalb (x, -m);
+
+  m += n;
+  n = m / 2;
+
+  if ((n + n) != m)
+    {
+      x *= 2;
+      --m;
+      n = m / 2;
+    }
+
+  /* Generate sqrt (X) bit by bit (accumulating in Q).  */
+  q = 1.0;
+  s = 4.0;
+  x -= 1.0;
+  r = 1;
+  for (i = 1; i <= 51; i++)
+    {
+      t = s + 1;
+      x *= 4;
+      r /= 2;
+      if (t <= x)
+	{
+	  s = t + t + 2, x -= t;
+	  q += r;
+	}
+      else
+	s *= 2;
+    }
+
+  /* Generate the last bit and determine the final rounding.  */
+  r /= 2;
+  x *= 4;
+  if (x == zero)
+    goto end;
+  (void) (100 + r);		/* Trigger inexact flag.  */
+  if (s < x)
+    {
+      q += r;
+      x -= s;
+      s += 2;
+      s *= 2;
+      x *= 4;
+      t = (x - s) - 5;
+      b = 1.0 + 3 * r / 4;
+      if (b == 1.0)
+	goto end;		/* B == 1: Round to zero.  */
+      b = 1.0 + r / 4;
+      if (b > 1.0)
+	t = 1;			/* B > 1: Round to +Inf.  */
+      if (t >= 0)
+	q += r;
+    }				/* Else round to nearest.  */
+  else
+    {
+      s *= 2;
+      x *= 4;
+      t = (x - s) - 1;
+      b = 1.0 + 3 * r / 4;
+      if (b == 1.0)
+	goto end;
+      b = 1.0 + r / 4;
+      if (b > 1.0)
+	t = 1;
+      if (t >= 0)
+	q += r;
+    }
+
+end:
+  return __scalb (q, n);
+}
diff --git a/sysdeps/ieee754/support.c b/sysdeps/ieee754/support.c
new file mode 100644
index 0000000000..e976839421
--- /dev/null
+++ b/sysdeps/ieee754/support.c
@@ -0,0 +1,524 @@
+/*
+ * Copyright (c) 1985, 1993
+ *	The Regents of the University of California.  All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ * 3. All advertising materials mentioning features or use of this software
+ *    must display the following acknowledgement:
+ *	This product includes software developed by the University of
+ *	California, Berkeley and its contributors.
+ * 4. Neither the name of the University nor the names of its contributors
+ *    may be used to endorse or promote products derived from this software
+ *    without specific prior written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
+ */
+
+#ifndef lint
+static char sccsid[] = "@(#)support.c	8.1 (Berkeley) 6/4/93";
+#endif /* not lint */
+
+/* 
+ * Some IEEE standard 754 recommended functions and remainder and sqrt for 
+ * supporting the C elementary functions.
+ ******************************************************************************
+ * WARNING:
+ *      These codes are developed (in double) to support the C elementary
+ * functions temporarily. They are not universal, and some of them are very
+ * slow (in particular, drem and sqrt is extremely inefficient). Each 
+ * computer system should have its implementation of these functions using 
+ * its own assembler.
+ ******************************************************************************
+ *
+ * IEEE 754 required operations:
+ *     drem(x,p) 
+ *              returns  x REM y  =  x - [x/y]*y , where [x/y] is the integer
+ *              nearest x/y; in half way case, choose the even one.
+ *     sqrt(x) 
+ *              returns the square root of x correctly rounded according to 
+ *		the rounding mod.
+ *
+ * IEEE 754 recommended functions:
+ * (a) copysign(x,y) 
+ *              returns x with the sign of y. 
+ * (b) scalb(x,N) 
+ *              returns  x * (2**N), for integer values N.
+ * (c) logb(x) 
+ *              returns the unbiased exponent of x, a signed integer in 
+ *              double precision, except that logb(0) is -INF, logb(INF) 
+ *              is +INF, and logb(NAN) is that NAN.
+ * (d) finite(x) 
+ *              returns the value TRUE if -INF < x < +INF and returns 
+ *              FALSE otherwise.
+ *
+ *
+ * CODED IN C BY K.C. NG, 11/25/84;
+ * REVISED BY K.C. NG on 1/22/85, 2/13/85, 3/24/85.
+ */
+
+#include "mathimpl.h"
+
+#if defined(vax)||defined(tahoe)      /* VAX D format */
+#include <errno.h>
+    static const unsigned short msign=0x7fff , mexp =0x7f80 ;
+    static const short  prep1=57, gap=7, bias=129           ;   
+    static const double novf=1.7E38, nunf=3.0E-39, zero=0.0 ;
+#else	/* defined(vax)||defined(tahoe) */
+    static const unsigned short msign=0x7fff, mexp =0x7ff0  ;
+    static const short prep1=54, gap=4, bias=1023           ;
+    static const double novf=1.7E308, nunf=3.0E-308,zero=0.0;
+#endif	/* defined(vax)||defined(tahoe) */
+
+double scalb(x,N)
+double x; int N;
+{
+        int k;
+
+#ifdef national
+        unsigned short *px=(unsigned short *) &x + 3;
+#else	/* national */
+        unsigned short *px=(unsigned short *) &x;
+#endif	/* national */
+
+        if( x == zero )  return(x); 
+
+#if defined(vax)||defined(tahoe)
+        if( (k= *px & mexp ) != ~msign ) {
+            if (N < -260)
+		return(nunf*nunf);
+	    else if (N > 260) {
+		return(copysign(infnan(ERANGE),x));
+	    }
+#else	/* defined(vax)||defined(tahoe) */
+        if( (k= *px & mexp ) != mexp ) {
+            if( N<-2100) return(nunf*nunf); else if(N>2100) return(novf+novf);
+            if( k == 0 ) {
+                 x *= scalb(1.0,(int)prep1);  N -= prep1; return(scalb(x,N));}
+#endif	/* defined(vax)||defined(tahoe) */
+
+            if((k = (k>>gap)+ N) > 0 )
+                if( k < (mexp>>gap) ) *px = (*px&~mexp) | (k<<gap);
+                else x=novf+novf;               /* overflow */
+            else
+                if( k > -prep1 ) 
+                                        /* gradual underflow */
+                    {*px=(*px&~mexp)|(short)(1<<gap); x *= scalb(1.0,k-1);}
+                else
+                return(nunf*nunf);
+            }
+        return(x);
+}
+
+
+double copysign(x,y)
+double x,y;
+{
+#ifdef national
+        unsigned short  *px=(unsigned short *) &x+3,
+                        *py=(unsigned short *) &y+3;
+#else	/* national */
+        unsigned short  *px=(unsigned short *) &x,
+                        *py=(unsigned short *) &y;
+#endif	/* national */
+
+#if defined(vax)||defined(tahoe)
+        if ( (*px & mexp) == 0 ) return(x);
+#endif	/* defined(vax)||defined(tahoe) */
+
+        *px = ( *px & msign ) | ( *py & ~msign );
+        return(x);
+}
+
+double logb(x)
+double x; 
+{
+
+#ifdef national
+        short *px=(short *) &x+3, k;
+#else	/* national */
+        short *px=(short *) &x, k;
+#endif	/* national */
+
+#if defined(vax)||defined(tahoe)
+        return (int)(((*px&mexp)>>gap)-bias);
+#else	/* defined(vax)||defined(tahoe) */
+        if( (k= *px & mexp ) != mexp )
+            if ( k != 0 )
+                return ( (k>>gap) - bias );
+            else if( x != zero)
+                return ( -1022.0 );
+            else        
+                return(-(1.0/zero));    
+        else if(x != x)
+            return(x);
+        else
+            {*px &= msign; return(x);}
+#endif	/* defined(vax)||defined(tahoe) */
+}
+
+finite(x)
+double x;    
+{
+#if defined(vax)||defined(tahoe)
+        return(1);
+#else	/* defined(vax)||defined(tahoe) */
+#ifdef national
+        return( (*((short *) &x+3 ) & mexp ) != mexp );
+#else	/* national */
+        return( (*((short *) &x ) & mexp ) != mexp );
+#endif	/* national */
+#endif	/* defined(vax)||defined(tahoe) */
+}
+
+double drem(x,p)
+double x,p;
+{
+        short sign;
+        double hp,dp,tmp;
+        unsigned short  k; 
+#ifdef national
+        unsigned short
+              *px=(unsigned short *) &x  +3, 
+              *pp=(unsigned short *) &p  +3,
+              *pd=(unsigned short *) &dp +3,
+              *pt=(unsigned short *) &tmp+3;
+#else	/* national */
+        unsigned short
+              *px=(unsigned short *) &x  , 
+              *pp=(unsigned short *) &p  ,
+              *pd=(unsigned short *) &dp ,
+              *pt=(unsigned short *) &tmp;
+#endif	/* national */
+
+        *pp &= msign ;
+
+#if defined(vax)||defined(tahoe)
+        if( ( *px & mexp ) == ~msign )	/* is x a reserved operand? */
+#else	/* defined(vax)||defined(tahoe) */
+        if( ( *px & mexp ) == mexp )
+#endif	/* defined(vax)||defined(tahoe) */
+		return  (x-p)-(x-p);	/* create nan if x is inf */
+	if (p == zero) {
+#if defined(vax)||defined(tahoe)
+		return(infnan(EDOM));
+#else	/* defined(vax)||defined(tahoe) */
+		return zero/zero;
+#endif	/* defined(vax)||defined(tahoe) */
+	}
+
+#if defined(vax)||defined(tahoe)
+        if( ( *pp & mexp ) == ~msign )	/* is p a reserved operand? */
+#else	/* defined(vax)||defined(tahoe) */
+        if( ( *pp & mexp ) == mexp )
+#endif	/* defined(vax)||defined(tahoe) */
+		{ if (p != p) return p; else return x;}
+
+        else  if ( ((*pp & mexp)>>gap) <= 1 ) 
+                /* subnormal p, or almost subnormal p */
+            { double b; b=scalb(1.0,(int)prep1);
+              p *= b; x = drem(x,p); x *= b; return(drem(x,p)/b);}
+        else  if ( p >= novf/2)
+            { p /= 2 ; x /= 2; return(drem(x,p)*2);}
+        else 
+            {
+                dp=p+p; hp=p/2;
+                sign= *px & ~msign ;
+                *px &= msign       ;
+                while ( x > dp )
+                    {
+                        k=(*px & mexp) - (*pd & mexp) ;
+                        tmp = dp ;
+                        *pt += k ;
+
+#if defined(vax)||defined(tahoe)
+                        if( x < tmp ) *pt -= 128 ;
+#else	/* defined(vax)||defined(tahoe) */
+                        if( x < tmp ) *pt -= 16 ;
+#endif	/* defined(vax)||defined(tahoe) */
+
+                        x -= tmp ;
+                    }
+                if ( x > hp )
+                    { x -= p ;  if ( x >= hp ) x -= p ; }
+
+#if defined(vax)||defined(tahoe)
+		if (x)
+#endif	/* defined(vax)||defined(tahoe) */
+			*px ^= sign;
+                return( x);
+
+            }
+}
+
+
+double sqrt(x)
+double x;
+{
+        double q,s,b,r;
+        double t;
+	double const zero=0.0;
+        int m,n,i;
+#if defined(vax)||defined(tahoe)
+        int k=54;
+#else	/* defined(vax)||defined(tahoe) */
+        int k=51;
+#endif	/* defined(vax)||defined(tahoe) */
+
+    /* sqrt(NaN) is NaN, sqrt(+-0) = +-0 */
+        if(x!=x||x==zero) return(x);
+
+    /* sqrt(negative) is invalid */
+        if(x<zero) {
+#if defined(vax)||defined(tahoe)
+		return (infnan(EDOM));	/* NaN */
+#else	/* defined(vax)||defined(tahoe) */
+		return(zero/zero);
+#endif	/* defined(vax)||defined(tahoe) */
+	}
+
+    /* sqrt(INF) is INF */
+        if(!finite(x)) return(x);               
+
+    /* scale x to [1,4) */
+        n=logb(x);
+        x=scalb(x,-n);
+        if((m=logb(x))!=0) x=scalb(x,-m);       /* subnormal number */
+        m += n; 
+        n = m/2;
+        if((n+n)!=m) {x *= 2; m -=1; n=m/2;}
+
+    /* generate sqrt(x) bit by bit (accumulating in q) */
+            q=1.0; s=4.0; x -= 1.0; r=1;
+            for(i=1;i<=k;i++) {
+                t=s+1; x *= 4; r /= 2;
+                if(t<=x) {
+                    s=t+t+2, x -= t; q += r;}
+                else
+                    s *= 2;
+                }
+            
+    /* generate the last bit and determine the final rounding */
+            r/=2; x *= 4; 
+            if(x==zero) goto end; 100+r; /* trigger inexact flag */
+            if(s<x) {
+                q+=r; x -=s; s += 2; s *= 2; x *= 4;
+                t = (x-s)-5; 
+                b=1.0+3*r/4; if(b==1.0) goto end; /* b==1 : Round-to-zero */
+                b=1.0+r/4;   if(b>1.0) t=1;	/* b>1 : Round-to-(+INF) */
+                if(t>=0) q+=r; }	      /* else: Round-to-nearest */
+            else { 
+                s *= 2; x *= 4; 
+                t = (x-s)-1; 
+                b=1.0+3*r/4; if(b==1.0) goto end;
+                b=1.0+r/4;   if(b>1.0) t=1;
+                if(t>=0) q+=r; }
+            
+end:        return(scalb(q,n));
+}
+
+#if 0
+/* DREM(X,Y)
+ * RETURN X REM Y =X-N*Y, N=[X/Y] ROUNDED (ROUNDED TO EVEN IN THE HALF WAY CASE)
+ * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
+ * INTENDED FOR ASSEMBLY LANGUAGE
+ * CODED IN C BY K.C. NG, 3/23/85, 4/8/85.
+ *
+ * Warning: this code should not get compiled in unless ALL of
+ * the following machine-dependent routines are supplied.
+ * 
+ * Required machine dependent functions (not on a VAX):
+ *     swapINX(i): save inexact flag and reset it to "i"
+ *     swapENI(e): save inexact enable and reset it to "e"
+ */
+
+double drem(x,y)	
+double x,y;
+{
+
+#ifdef national		/* order of words in floating point number */
+	static const n0=3,n1=2,n2=1,n3=0;
+#else /* VAX, SUN, ZILOG, TAHOE */
+	static const n0=0,n1=1,n2=2,n3=3;
+#endif
+
+    	static const unsigned short mexp =0x7ff0, m25 =0x0190, m57 =0x0390;
+	static const double zero=0.0;
+	double hy,y1,t,t1;
+	short k;
+	long n;
+	int i,e; 
+	unsigned short xexp,yexp, *px  =(unsigned short *) &x  , 
+	      		nx,nf,	  *py  =(unsigned short *) &y  ,
+	      		sign,	  *pt  =(unsigned short *) &t  ,
+	      			  *pt1 =(unsigned short *) &t1 ;
+
+	xexp = px[n0] & mexp ;	/* exponent of x */
+	yexp = py[n0] & mexp ;	/* exponent of y */
+	sign = px[n0] &0x8000;	/* sign of x     */
+
+/* return NaN if x is NaN, or y is NaN, or x is INF, or y is zero */
+	if(x!=x) return(x); if(y!=y) return(y);	     /* x or y is NaN */
+	if( xexp == mexp )   return(zero/zero);      /* x is INF */
+	if(y==zero) return(y/y);
+
+/* save the inexact flag and inexact enable in i and e respectively
+ * and reset them to zero
+ */
+	i=swapINX(0);	e=swapENI(0);	
+
+/* subnormal number */
+	nx=0;
+	if(yexp==0) {t=1.0,pt[n0]+=m57; y*=t; nx=m57;}
+
+/* if y is tiny (biased exponent <= 57), scale up y to y*2**57 */
+	if( yexp <= m57 ) {py[n0]+=m57; nx+=m57; yexp+=m57;}
+
+	nf=nx;
+	py[n0] &= 0x7fff;	
+	px[n0] &= 0x7fff;
+
+/* mask off the least significant 27 bits of y */
+	t=y; pt[n3]=0; pt[n2]&=0xf800; y1=t;
+
+/* LOOP: argument reduction on x whenever x > y */
+loop:
+	while ( x > y )
+	{
+	    t=y;
+	    t1=y1;
+	    xexp=px[n0]&mexp;	  /* exponent of x */
+	    k=xexp-yexp-m25;
+	    if(k>0) 	/* if x/y >= 2**26, scale up y so that x/y < 2**26 */
+		{pt[n0]+=k;pt1[n0]+=k;}
+	    n=x/t; x=(x-n*t1)-n*(t-t1);
+	}	
+    /* end while (x > y) */
+
+	if(nx!=0) {t=1.0; pt[n0]+=nx; x*=t; nx=0; goto loop;}
+
+/* final adjustment */
+
+	hy=y/2.0;
+	if(x>hy||((x==hy)&&n%2==1)) x-=y; 
+	px[n0] ^= sign;
+	if(nf!=0) { t=1.0; pt[n0]-=nf; x*=t;}
+
+/* restore inexact flag and inexact enable */
+	swapINX(i); swapENI(e);	
+
+	return(x);	
+}
+#endif
+
+#if 0
+/* SQRT
+ * RETURN CORRECTLY ROUNDED (ACCORDING TO THE ROUNDING MODE) SQRT
+ * FOR IEEE DOUBLE PRECISION ONLY, INTENDED FOR ASSEMBLY LANGUAGE
+ * CODED IN C BY K.C. NG, 3/22/85.
+ *
+ * Warning: this code should not get compiled in unless ALL of
+ * the following machine-dependent routines are supplied.
+ * 
+ * Required machine dependent functions:
+ *     swapINX(i)  ...return the status of INEXACT flag and reset it to "i"
+ *     swapRM(r)   ...return the current Rounding Mode and reset it to "r"
+ *     swapENI(e)  ...return the status of inexact enable and reset it to "e"
+ *     addc(t)     ...perform t=t+1 regarding t as a 64 bit unsigned integer
+ *     subc(t)     ...perform t=t-1 regarding t as a 64 bit unsigned integer
+ */
+
+static const unsigned long table[] = {
+0, 1204, 3062, 5746, 9193, 13348, 18162, 23592, 29598, 36145, 43202, 50740,
+58733, 67158, 75992, 85215, 83599, 71378, 60428, 50647, 41945, 34246, 27478,
+21581, 16499, 12183, 8588, 5674, 3403, 1742, 661, 130, };
+
+double newsqrt(x)
+double x;
+{
+        double y,z,t,addc(),subc()
+	double const b54=134217728.*134217728.; /* b54=2**54 */
+        long mx,scalx;
+	long const mexp=0x7ff00000;
+        int i,j,r,e,swapINX(),swapRM(),swapENI();       
+        unsigned long *py=(unsigned long *) &y   ,
+                      *pt=(unsigned long *) &t   ,
+                      *px=(unsigned long *) &x   ;
+#ifdef national         /* ordering of word in a floating point number */
+        const int n0=1, n1=0; 
+#else
+        const int n0=0, n1=1; 
+#endif
+/* Rounding Mode:  RN ...round-to-nearest 
+ *                 RZ ...round-towards 0
+ *                 RP ...round-towards +INF
+ *		   RM ...round-towards -INF
+ */
+        const int RN=0,RZ=1,RP=2,RM=3;
+				/* machine dependent: work on a Zilog Z8070
+                                 * and a National 32081 & 16081
+                                 */
+
+/* exceptions */
+	if(x!=x||x==0.0) return(x);  /* sqrt(NaN) is NaN, sqrt(+-0) = +-0 */
+	if(x<0) return((x-x)/(x-x)); /* sqrt(negative) is invalid */
+        if((mx=px[n0]&mexp)==mexp) return(x);  /* sqrt(+INF) is +INF */
+
+/* save, reset, initialize */
+        e=swapENI(0);   /* ...save and reset the inexact enable */
+        i=swapINX(0);   /* ...save INEXACT flag */
+        r=swapRM(RN);   /* ...save and reset the Rounding Mode to RN */
+        scalx=0;
+
+/* subnormal number, scale up x to x*2**54 */
+        if(mx==0) {x *= b54 ; scalx-=0x01b00000;}
+
+/* scale x to avoid intermediate over/underflow:
+ * if (x > 2**512) x=x/2**512; if (x < 2**-512) x=x*2**512 */
+        if(mx>0x5ff00000) {px[n0] -= 0x20000000; scalx+= 0x10000000;}
+        if(mx<0x1ff00000) {px[n0] += 0x20000000; scalx-= 0x10000000;}
+
+/* magic initial approximation to almost 8 sig. bits */
+        py[n0]=(px[n0]>>1)+0x1ff80000;
+        py[n0]=py[n0]-table[(py[n0]>>15)&31];
+
+/* Heron's rule once with correction to improve y to almost 18 sig. bits */
+        t=x/y; y=y+t; py[n0]=py[n0]-0x00100006; py[n1]=0;
+
+/* triple to almost 56 sig. bits; now y approx. sqrt(x) to within 1 ulp */
+        t=y*y; z=t;  pt[n0]+=0x00100000; t+=z; z=(x-z)*y; 
+        t=z/(t+x) ;  pt[n0]+=0x00100000; y+=t;
+
+/* twiddle last bit to force y correctly rounded */ 
+        swapRM(RZ);     /* ...set Rounding Mode to round-toward-zero */
+        swapINX(0);     /* ...clear INEXACT flag */
+        swapENI(e);     /* ...restore inexact enable status */
+        t=x/y;          /* ...chopped quotient, possibly inexact */
+        j=swapINX(i);   /* ...read and restore inexact flag */
+        if(j==0) { if(t==y) goto end; else t=subc(t); }  /* ...t=t-ulp */
+        b54+0.1;        /* ..trigger inexact flag, sqrt(x) is inexact */
+        if(r==RN) t=addc(t);            /* ...t=t+ulp */
+        else if(r==RP) { t=addc(t);y=addc(y);}/* ...t=t+ulp;y=y+ulp; */
+        y=y+t;                          /* ...chopped sum */
+        py[n0]=py[n0]-0x00100000;       /* ...correctly rounded sqrt(x) */
+end:    py[n0]=py[n0]+scalx;            /* ...scale back y */
+        swapRM(r);                      /* ...restore Rounding Mode */
+        return(y);
+}
+#endif