diff options
Diffstat (limited to 'sysdeps/ieee754')
| -rw-r--r-- | sysdeps/ieee754/Dist | 2 | ||||
| -rw-r--r-- | sysdeps/ieee754/cabs.c | 228 | ||||
| -rw-r--r-- | sysdeps/ieee754/cbrt.c | 120 | ||||
| -rw-r--r-- | sysdeps/ieee754/copysign.c | 37 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl2mpn.c | 107 | ||||
| -rw-r--r-- | sysdeps/ieee754/drem.c | 107 | ||||
| -rw-r--r-- | sysdeps/ieee754/frexp.c | 42 | ||||
| -rw-r--r-- | sysdeps/ieee754/huge_val.h | 87 | ||||
| -rw-r--r-- | sysdeps/ieee754/ieee754.h | 153 | ||||
| -rw-r--r-- | sysdeps/ieee754/infnan.c | 50 | ||||
| -rw-r--r-- | sysdeps/ieee754/isinf.c | 40 | ||||
| -rw-r--r-- | sysdeps/ieee754/isinfl.c | 44 | ||||
| -rw-r--r-- | sysdeps/ieee754/isnan.c | 37 | ||||
| -rw-r--r-- | sysdeps/ieee754/isnanl.c | 40 | ||||
| -rw-r--r-- | sysdeps/ieee754/ldbl2mpn.c | 93 | ||||
| -rw-r--r-- | sysdeps/ieee754/ldexp.c | 145 | ||||
| -rw-r--r-- | sysdeps/ieee754/log10.c | 30 | ||||
| -rw-r--r-- | sysdeps/ieee754/logb.c | 48 | ||||
| -rw-r--r-- | sysdeps/ieee754/mpn2dbl.c | 46 | ||||
| -rw-r--r-- | sysdeps/ieee754/mpn2flt.c | 42 | ||||
| -rw-r--r-- | sysdeps/ieee754/mpn2ldbl.c | 46 | ||||
| -rw-r--r-- | sysdeps/ieee754/nan.h | 46 | ||||
| -rw-r--r-- | sysdeps/ieee754/sqrt.c | 120 | ||||
| -rw-r--r-- | sysdeps/ieee754/support.c | 524 |
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 |