aboutsummaryrefslogtreecommitdiff
path: root/sysdeps/ieee754
diff options
context:
space:
mode:
Diffstat (limited to 'sysdeps/ieee754')
-rw-r--r--sysdeps/ieee754/Dist2
-rw-r--r--sysdeps/ieee754/cabs.c228
-rw-r--r--sysdeps/ieee754/cbrt.c120
-rw-r--r--sysdeps/ieee754/copysign.c37
-rw-r--r--sysdeps/ieee754/dbl2mpn.c107
-rw-r--r--sysdeps/ieee754/drem.c107
-rw-r--r--sysdeps/ieee754/frexp.c42
-rw-r--r--sysdeps/ieee754/huge_val.h87
-rw-r--r--sysdeps/ieee754/ieee754.h153
-rw-r--r--sysdeps/ieee754/infnan.c50
-rw-r--r--sysdeps/ieee754/isinf.c40
-rw-r--r--sysdeps/ieee754/isinfl.c44
-rw-r--r--sysdeps/ieee754/isnan.c37
-rw-r--r--sysdeps/ieee754/isnanl.c40
-rw-r--r--sysdeps/ieee754/ldbl2mpn.c93
-rw-r--r--sysdeps/ieee754/ldexp.c145
-rw-r--r--sysdeps/ieee754/log10.c30
-rw-r--r--sysdeps/ieee754/logb.c48
-rw-r--r--sysdeps/ieee754/mpn2dbl.c46
-rw-r--r--sysdeps/ieee754/mpn2flt.c42
-rw-r--r--sysdeps/ieee754/mpn2ldbl.c46
-rw-r--r--sysdeps/ieee754/nan.h46
-rw-r--r--sysdeps/ieee754/sqrt.c120
-rw-r--r--sysdeps/ieee754/support.c524
24 files changed, 2234 insertions, 0 deletions
diff --git a/sysdeps/ieee754/Dist b/sysdeps/ieee754/Dist
new file mode 100644
index 0000000..94cc5c9
--- /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 0000000..6b0d4c4
--- /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 0000000..fe5fb95
--- /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 0000000..f1f0591
--- /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 0000000..6b690f7
--- /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 0000000..cab3a04
--- /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 0000000..c56a17f
--- /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 0000000..183f213
--- /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 0000000..9cc6f14
--- /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 0000000..89ab5b6
--- /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 0000000..79be916
--- /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 0000000..ad03728
--- /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 0000000..8a537cc
--- /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 0000000..b06b31a
--- /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 0000000..29e8a3b
--- /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 0000000..7e1c747
--- /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 0000000..da2f5b4
--- /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 0000000..918de9e
--- /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 0000000..01e0019
--- /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 0000000..e9bbc49
--- /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 0000000..831f2e5
--- /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 0000000..e34c172
--- /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 0000000..7e350e0
--- /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 0000000..e976839
--- /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