diff options
author | das <das@FreeBSD.org> | 2005-01-22 06:03:40 +0000 |
---|---|---|
committer | das <das@FreeBSD.org> | 2005-01-22 06:03:40 +0000 |
commit | 87888496cc9e2d5e2e57e2f77d2fa7f615a0e1ed (patch) | |
tree | 453d24603f8076b2c01c839cddf81fa32e3c6416 /lib | |
parent | cb779e8aab212ecfba1e0ec21a5362a2d5457c95 (diff) | |
download | FreeBSD-src-87888496cc9e2d5e2e57e2f77d2fa7f615a0e1ed.zip FreeBSD-src-87888496cc9e2d5e2e57e2f77d2fa7f615a0e1ed.tar.gz |
Replace the ldexp() implementation in libc with a renamed copy of the
scalbn() implementation from libm. (The two functions are defined to
be identical, but ldexp() lives in libc for backwards compatibility.)
The old ldexp() implementation...
- was more complicated than this one
- set errno instead of raising FP exceptions
- got some corner cases wrong
(e.g. ldexp(1.0, 2000) in round-to-zero mode)
The new implementation lives in libc/gen instead of
libc/$MACHINE_ARCH/gen, since we don't need N copies of a
machine-independent file. The amd64 and i386 platforms
retain their fast and correct MD implementations and
override this one.
Diffstat (limited to 'lib')
-rw-r--r-- | lib/libc/alpha/gen/ldexp.c | 137 | ||||
-rw-r--r-- | lib/libc/arm/gen/ldexp.c | 155 | ||||
-rw-r--r-- | lib/libc/gen/ldexp.c | 123 | ||||
-rw-r--r-- | lib/libc/ia64/gen/ldexp.c | 137 | ||||
-rw-r--r-- | lib/libc/powerpc/gen/ldexp.c | 155 | ||||
-rw-r--r-- | lib/libc/sparc64/gen/ldexp.c | 155 |
6 files changed, 123 insertions, 739 deletions
diff --git a/lib/libc/alpha/gen/ldexp.c b/lib/libc/alpha/gen/ldexp.c deleted file mode 100644 index 78decbb..0000000 --- a/lib/libc/alpha/gen/ldexp.c +++ /dev/null @@ -1,137 +0,0 @@ -/* - * Copyright (c) 1994, 1995 Carnegie-Mellon University. - * All rights reserved. - * - * Author: Chris G. Demetriou - * - * Permission to use, copy, modify and distribute this software and - * its documentation is hereby granted, provided that both the copyright - * notice and this permission notice appear in all copies of the - * software, derivative works or modified versions, and any portions - * thereof, and that both notices appear in supporting documentation. - * - * CARNEGIE MELLON ALLOWS FREE USE OF THIS SOFTWARE IN ITS "AS IS" - * CONDITION. CARNEGIE MELLON DISCLAIMS ANY LIABILITY OF ANY KIND - * FOR ANY DAMAGES WHATSOEVER RESULTING FROM THE USE OF THIS SOFTWARE. - * - * Carnegie Mellon requests users of this software to return to - * - * Software Distribution Coordinator or Software.Distribution@CS.CMU.EDU - * School of Computer Science - * Carnegie Mellon University - * Pittsburgh PA 15213-3890 - * - * any improvements or extensions that they make and grant Carnegie the - * rights to redistribute these changes. - * - * $NetBSD: ldexp.c,v 1.1 1995/02/10 17:50:24 cgd Exp $ - */ - -#include <sys/cdefs.h> -__FBSDID("$FreeBSD$"); - -#include <sys/types.h> -#include <machine/ieee.h> -#include <errno.h> -#include <math.h> - -/* - * double ldexp(double val, int exp) - * returns: val * (2**exp) - */ -double -ldexp(val, exp) - double val; - int exp; -{ - int oldexp, newexp, mulexp; - union doub { - double v; - struct ieee_double s; - } u, mul; - - /* - * If input is zero, or no change, just return input. - * Likewise, if input is Inf or NaN, just return it. - */ - u.v = val; - oldexp = u.s.dbl_exp; - if (val == 0 || exp == 0 || oldexp == DBL_EXP_INFNAN) - return (val); - - /* - * Compute new exponent and check for over/under flow. - * Underflow, unfortunately, could mean switching to denormal. - * If result out of range, set ERANGE and return 0 if too small - * or Inf if too big, with the same sign as the input value. - */ - newexp = oldexp + exp; - if (newexp >= DBL_EXP_INFNAN) { - /* u.s.dbl_sign = val < 0; -- already set */ - u.s.dbl_exp = DBL_EXP_INFNAN; - u.s.dbl_frach = u.s.dbl_fracl = 0; - errno = ERANGE; - return (u.v); /* Inf */ - } - if (newexp <= 0) { - /* - * The output number is either a denormal or underflows - * (see comments in machine/ieee.h). - */ - if (newexp <= -DBL_FRACBITS) { - /* u.s.dbl_sign = val < 0; -- already set */ - u.s.dbl_exp = 0; - u.s.dbl_frach = u.s.dbl_fracl = 0; - errno = ERANGE; - return (u.v); /* zero */ - } - /* - * We are going to produce a denorm. Our `exp' argument - * might be as small as -2097, and we cannot compute - * 2^-2097, so we may have to do this as many as three - * steps (not just two, as for positive `exp's below). - */ - mul.v = 0; - while (exp <= -DBL_EXP_BIAS) { - mul.s.dbl_exp = 1; - val *= mul.v; - exp += DBL_EXP_BIAS - 1; - } - mul.s.dbl_exp = exp + DBL_EXP_BIAS; - val *= mul.v; - return (val); - } - - /* - * Newexp is positive. - * - * If oldexp is zero, we are starting with a denorm, and simply - * adjusting the exponent will produce bogus answers. We need - * to fix that first. - */ - if (oldexp == 0) { - /* - * Multiply by 2^mulexp to make the number normalizable. - * We cannot multiply by more than 2^1023, but `exp' - * argument might be as large as 2046. A single - * adjustment, however, will normalize the number even - * for huge `exp's, and then we can use exponent - * arithmetic just as for normal `double's. - */ - mulexp = exp <= DBL_EXP_BIAS ? exp : DBL_EXP_BIAS; - mul.v = 0; - mul.s.dbl_exp = mulexp + DBL_EXP_BIAS; - val *= mul.v; - if (mulexp == exp) - return (val); - u.v = val; - newexp -= mulexp; - } - - /* - * Both oldexp and newexp are positive; just replace the - * old exponent with the new one. - */ - u.s.dbl_exp = newexp; - return (u.v); -} diff --git a/lib/libc/arm/gen/ldexp.c b/lib/libc/arm/gen/ldexp.c deleted file mode 100644 index 9aed74d..0000000 --- a/lib/libc/arm/gen/ldexp.c +++ /dev/null @@ -1,155 +0,0 @@ -/*- - * Copyright (c) 1999 The NetBSD Foundation, Inc. - * All rights reserved. - * - * This code is derived from software contributed to The NetBSD Foundation - * by Charles M. Hannum. - * - * 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 NetBSD - * Foundation, Inc. and its contributors. - * 4. Neither the name of The NetBSD Foundation 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 NETBSD FOUNDATION, INC. 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 FOUNDATION 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. - */ - -#include <sys/cdefs.h> -#if 0 -#if defined(LIBC_SCCS) && !defined(lint) -__RCSID("$NetBSD: ldexp.c,v 1.8 1999/08/30 18:28:26 mycroft Exp $"); -#endif /* LIBC_SCCS and not lint */ -#endif -__FBSDID("$FreeBSD$"); - -#include <sys/types.h> -#include <errno.h> -#include <math.h> -#include <machine/ieee.h> - -/* - * Multiply the given value by 2^exp. - */ -double -ldexp(val, exp) - double val; - int exp; -{ - int oldexp, newexp; - union { - double v; - struct ieee_double s; - } u, mul; - - u.v = val; - oldexp = u.s.dbl_exp; - - /* - * If input is zero, Inf or NaN, just return it. - */ - if (u.v == 0.0 || oldexp == DBL_EXP_INFNAN) - return (val); - - if (oldexp == 0) { - /* - * u.v is denormal. We must adjust it so that the exponent - * arithmetic below will work. - */ - if (exp <= DBL_EXP_BIAS) { - /* - * Optimization: if the scaling can be done in a single - * multiply, or underflows, just do it now. - */ - if (exp <= -DBL_FRACBITS) { - errno = ERANGE; - return (0.0); - } - mul.v = 0.0; - mul.s.dbl_exp = exp + DBL_EXP_BIAS; - u.v *= mul.v; - if (u.v == 0.0) { - errno = ERANGE; - return (0.0); - } - return (u.v); - } else { - /* - * We know that exp is very large, and therefore the - * result cannot be denormal (though it may be Inf). - * Shift u.v by just enough to make it normal. - */ - mul.v = 0.0; - mul.s.dbl_exp = DBL_FRACBITS + DBL_EXP_BIAS; - u.v *= mul.v; - exp -= DBL_FRACBITS; - oldexp = u.s.dbl_exp; - } - } - - /* - * u.v is now normalized and oldexp has been adjusted if necessary. - * Calculate the new exponent and check for underflow and overflow. - */ - newexp = oldexp + exp; - - if (newexp <= 0) { - /* - * The output number is either denormal or underflows (see - * comments in machine/ieee.h). - */ - if (newexp <= -DBL_FRACBITS) { - errno = ERANGE; - return (0.0); - } - /* - * Denormalize the result. We do this with a multiply. If exp - * is very large, it won't fit in a double, so we have to - * adjust the exponent first. This is safe because we know - * that u.v is normal at this point. - */ - if (exp <= -DBL_EXP_BIAS) { - u.s.dbl_exp = 1; - exp += oldexp - 1; - } - mul.v = 0.0; - mul.s.dbl_exp = exp + DBL_EXP_BIAS; - u.v *= mul.v; - return (u.v); - } else if (newexp >= DBL_EXP_INFNAN) { - /* - * The result overflowed; return +/-Inf. - */ - u.s.dbl_exp = DBL_EXP_INFNAN; - u.s.dbl_frach = 0; - u.s.dbl_fracl = 0; - errno = ERANGE; - return (u.v); - } else { - /* - * The result is normal; just replace the old exponent with the - * new one. - */ - u.s.dbl_exp = newexp; - return (u.v); - } -} diff --git a/lib/libc/gen/ldexp.c b/lib/libc/gen/ldexp.c new file mode 100644 index 0000000..887f673 --- /dev/null +++ b/lib/libc/gen/ldexp.c @@ -0,0 +1,123 @@ +/* @(#)s_scalbn.c 5.1 93/09/24 */ +/* @(#)fdlibm.h 5.1 93/09/24 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +#include <sys/cdefs.h> +__FBSDID("$FreeBSD$"); + +#include <sys/types.h> +#include <machine/endian.h> +#include <math.h> + +/* Bit fiddling routines copied from msun/src/math_private.h,v 1.15 */ + +#if BYTE_ORDER == BIG_ENDIAN + +typedef union +{ + double value; + struct + { + u_int32_t msw; + u_int32_t lsw; + } parts; +} ieee_double_shape_type; + +#endif + +#if BYTE_ORDER == LITTLE_ENDIAN + +typedef union +{ + double value; + struct + { + u_int32_t lsw; + u_int32_t msw; + } parts; +} ieee_double_shape_type; + +#endif + +/* Get two 32 bit ints from a double. */ + +#define EXTRACT_WORDS(ix0,ix1,d) \ +do { \ + ieee_double_shape_type ew_u; \ + ew_u.value = (d); \ + (ix0) = ew_u.parts.msw; \ + (ix1) = ew_u.parts.lsw; \ +} while (0) + +/* Get the more significant 32 bit int from a double. */ + +#define GET_HIGH_WORD(i,d) \ +do { \ + ieee_double_shape_type gh_u; \ + gh_u.value = (d); \ + (i) = gh_u.parts.msw; \ +} while (0) + +/* Set the more significant 32 bits of a double from an int. */ + +#define SET_HIGH_WORD(d,v) \ +do { \ + ieee_double_shape_type sh_u; \ + sh_u.value = (d); \ + sh_u.parts.msw = (v); \ + (d) = sh_u.value; \ +} while (0) + + +static const double +two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */ +twom54 = 5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */ +huge = 1.0e+300, +tiny = 1.0e-300; + +static double +_copysign(double x, double y) +{ + u_int32_t hx,hy; + GET_HIGH_WORD(hx,x); + GET_HIGH_WORD(hy,y); + SET_HIGH_WORD(x,(hx&0x7fffffff)|(hy&0x80000000)); + return x; +} + +double +ldexp(double x, int n) +{ + int32_t k,hx,lx; + EXTRACT_WORDS(hx,lx,x); + k = (hx&0x7ff00000)>>20; /* extract exponent */ + if (k==0) { /* 0 or subnormal x */ + if ((lx|(hx&0x7fffffff))==0) return x; /* +-0 */ + x *= two54; + GET_HIGH_WORD(hx,x); + k = ((hx&0x7ff00000)>>20) - 54; + if (n< -50000) return tiny*x; /*underflow*/ + } + if (k==0x7ff) return x+x; /* NaN or Inf */ + k = k+n; + if (k > 0x7fe) return huge*_copysign(huge,x); /* overflow */ + if (k > 0) /* normal result */ + {SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20)); return x;} + if (k <= -54) { + if (n > 50000) /* in case integer overflow in n+k */ + return huge*_copysign(huge,x); /*overflow*/ + else return tiny*_copysign(tiny,x); /*underflow*/ + } + k += 54; /* subnormal result */ + SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20)); + return x*twom54; +} diff --git a/lib/libc/ia64/gen/ldexp.c b/lib/libc/ia64/gen/ldexp.c deleted file mode 100644 index c431417..0000000 --- a/lib/libc/ia64/gen/ldexp.c +++ /dev/null @@ -1,137 +0,0 @@ -/* $NetBSD: ldexp.c,v 1.1 1995/02/10 17:50:24 cgd Exp $ */ - -/* - * Copyright (c) 1994, 1995 Carnegie-Mellon University. - * All rights reserved. - * - * Author: Chris G. Demetriou - * - * Permission to use, copy, modify and distribute this software and - * its documentation is hereby granted, provided that both the copyright - * notice and this permission notice appear in all copies of the - * software, derivative works or modified versions, and any portions - * thereof, and that both notices appear in supporting documentation. - * - * CARNEGIE MELLON ALLOWS FREE USE OF THIS SOFTWARE IN ITS "AS IS" - * CONDITION. CARNEGIE MELLON DISCLAIMS ANY LIABILITY OF ANY KIND - * FOR ANY DAMAGES WHATSOEVER RESULTING FROM THE USE OF THIS SOFTWARE. - * - * Carnegie Mellon requests users of this software to return to - * - * Software Distribution Coordinator or Software.Distribution@CS.CMU.EDU - * School of Computer Science - * Carnegie Mellon University - * Pittsburgh PA 15213-3890 - * - * any improvements or extensions that they make and grant Carnegie the - * rights to redistribute these changes. - */ - -#include <sys/cdefs.h> -__FBSDID("$FreeBSD$"); - -#include <sys/types.h> -#include <machine/ieee.h> -#include <errno.h> -#include <math.h> - -/* - * double ldexp(double val, int exp) - * returns: val * (2**exp) - */ -double -ldexp(val, exp) - double val; - int exp; -{ - int oldexp, newexp, mulexp; - union doub { - double v; - struct ieee_double s; - } u, mul; - - /* - * If input is zero, or no change, just return input. - * Likewise, if input is Inf or NaN, just return it. - */ - u.v = val; - oldexp = u.s.dbl_exp; - if (val == 0 || exp == 0 || oldexp == DBL_EXP_INFNAN) - return (val); - - /* - * Compute new exponent and check for over/under flow. - * Underflow, unfortunately, could mean switching to denormal. - * If result out of range, set ERANGE and return 0 if too small - * or Inf if too big, with the same sign as the input value. - */ - newexp = oldexp + exp; - if (newexp >= DBL_EXP_INFNAN) { - /* u.s.dbl_sign = val < 0; -- already set */ - u.s.dbl_exp = DBL_EXP_INFNAN; - u.s.dbl_frach = u.s.dbl_fracl = 0; - errno = ERANGE; - return (u.v); /* Inf */ - } - if (newexp <= 0) { - /* - * The output number is either a denormal or underflows - * (see comments in machine/ieee.h). - */ - if (newexp <= -DBL_FRACBITS) { - /* u.s.dbl_sign = val < 0; -- already set */ - u.s.dbl_exp = 0; - u.s.dbl_frach = u.s.dbl_fracl = 0; - errno = ERANGE; - return (u.v); /* zero */ - } - /* - * We are going to produce a denorm. Our `exp' argument - * might be as small as -2097, and we cannot compute - * 2^-2097, so we may have to do this as many as three - * steps (not just two, as for positive `exp's below). - */ - mul.v = 0; - while (exp <= -DBL_EXP_BIAS) { - mul.s.dbl_exp = 1; - val *= mul.v; - exp += DBL_EXP_BIAS - 1; - } - mul.s.dbl_exp = exp + DBL_EXP_BIAS; - val *= mul.v; - return (val); - } - - /* - * Newexp is positive. - * - * If oldexp is zero, we are starting with a denorm, and simply - * adjusting the exponent will produce bogus answers. We need - * to fix that first. - */ - if (oldexp == 0) { - /* - * Multiply by 2^mulexp to make the number normalizable. - * We cannot multiply by more than 2^1023, but `exp' - * argument might be as large as 2046. A single - * adjustment, however, will normalize the number even - * for huge `exp's, and then we can use exponent - * arithmetic just as for normal `double's. - */ - mulexp = exp <= DBL_EXP_BIAS ? exp : DBL_EXP_BIAS; - mul.v = 0; - mul.s.dbl_exp = mulexp + DBL_EXP_BIAS; - val *= mul.v; - if (mulexp == exp) - return (val); - u.v = val; - newexp -= mulexp; - } - - /* - * Both oldexp and newexp are positive; just replace the - * old exponent with the new one. - */ - u.s.dbl_exp = newexp; - return (u.v); -} diff --git a/lib/libc/powerpc/gen/ldexp.c b/lib/libc/powerpc/gen/ldexp.c deleted file mode 100644 index 3a85612..0000000 --- a/lib/libc/powerpc/gen/ldexp.c +++ /dev/null @@ -1,155 +0,0 @@ -/*- - * Copyright (c) 1999 The NetBSD Foundation, Inc. - * All rights reserved. - * - * This code is derived from software contributed to The NetBSD Foundation - * by Charles M. Hannum. - * - * 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 NetBSD - * Foundation, Inc. and its contributors. - * 4. Neither the name of The NetBSD Foundation 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 NETBSD FOUNDATION, INC. 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 FOUNDATION 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. - */ - -#include <sys/cdefs.h> -#if 0 -#if defined(LIBC_SCCS) && !defined(lint) -__RCSID("$NetBSD: ldexp.c,v 1.8 1999/08/30 18:28:26 mycroft Exp $"); -#endif /* LIBC_SCCS and not lint */ -#endif -__FBSDID("$FreeBSD$"); - -#include <sys/types.h> -#include <machine/ieee.h> -#include <errno.h> -#include <math.h> - -/* - * Multiply the given value by 2^exp. - */ -double -ldexp(val, exp) - double val; - int exp; -{ - int oldexp, newexp; - union { - double v; - struct ieee_double s; - } u, mul; - - u.v = val; - oldexp = u.s.dbl_exp; - - /* - * If input is zero, Inf or NaN, just return it. - */ - if (u.v == 0.0 || oldexp == DBL_EXP_INFNAN) - return (val); - - if (oldexp == 0) { - /* - * u.v is denormal. We must adjust it so that the exponent - * arithmetic below will work. - */ - if (exp <= DBL_EXP_BIAS) { - /* - * Optimization: if the scaling can be done in a single - * multiply, or underflows, just do it now. - */ - if (exp <= -DBL_FRACBITS) { - errno = ERANGE; - return (0.0); - } - mul.v = 0.0; - mul.s.dbl_exp = exp + DBL_EXP_BIAS; - u.v *= mul.v; - if (u.v == 0.0) { - errno = ERANGE; - return (0.0); - } - return (u.v); - } else { - /* - * We know that exp is very large, and therefore the - * result cannot be denormal (though it may be Inf). - * Shift u.v by just enough to make it normal. - */ - mul.v = 0.0; - mul.s.dbl_exp = DBL_FRACBITS + DBL_EXP_BIAS; - u.v *= mul.v; - exp -= DBL_FRACBITS; - oldexp = u.s.dbl_exp; - } - } - - /* - * u.v is now normalized and oldexp has been adjusted if necessary. - * Calculate the new exponent and check for underflow and overflow. - */ - newexp = oldexp + exp; - - if (newexp <= 0) { - /* - * The output number is either denormal or underflows (see - * comments in machine/ieee.h). - */ - if (newexp <= -DBL_FRACBITS) { - errno = ERANGE; - return (0.0); - } - /* - * Denormalize the result. We do this with a multiply. If exp - * is very large, it won't fit in a double, so we have to - * adjust the exponent first. This is safe because we know - * that u.v is normal at this point. - */ - if (exp <= -DBL_EXP_BIAS) { - u.s.dbl_exp = 1; - exp += oldexp - 1; - } - mul.v = 0.0; - mul.s.dbl_exp = exp + DBL_EXP_BIAS; - u.v *= mul.v; - return (u.v); - } else if (newexp >= DBL_EXP_INFNAN) { - /* - * The result overflowed; return +/-Inf. - */ - u.s.dbl_exp = DBL_EXP_INFNAN; - u.s.dbl_frach = 0; - u.s.dbl_fracl = 0; - errno = ERANGE; - return (u.v); - } else { - /* - * The result is normal; just replace the old exponent with the - * new one. - */ - u.s.dbl_exp = newexp; - return (u.v); - } -} diff --git a/lib/libc/sparc64/gen/ldexp.c b/lib/libc/sparc64/gen/ldexp.c deleted file mode 100644 index 3a85612..0000000 --- a/lib/libc/sparc64/gen/ldexp.c +++ /dev/null @@ -1,155 +0,0 @@ -/*- - * Copyright (c) 1999 The NetBSD Foundation, Inc. - * All rights reserved. - * - * This code is derived from software contributed to The NetBSD Foundation - * by Charles M. Hannum. - * - * 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 NetBSD - * Foundation, Inc. and its contributors. - * 4. Neither the name of The NetBSD Foundation 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 NETBSD FOUNDATION, INC. 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 FOUNDATION 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. - */ - -#include <sys/cdefs.h> -#if 0 -#if defined(LIBC_SCCS) && !defined(lint) -__RCSID("$NetBSD: ldexp.c,v 1.8 1999/08/30 18:28:26 mycroft Exp $"); -#endif /* LIBC_SCCS and not lint */ -#endif -__FBSDID("$FreeBSD$"); - -#include <sys/types.h> -#include <machine/ieee.h> -#include <errno.h> -#include <math.h> - -/* - * Multiply the given value by 2^exp. - */ -double -ldexp(val, exp) - double val; - int exp; -{ - int oldexp, newexp; - union { - double v; - struct ieee_double s; - } u, mul; - - u.v = val; - oldexp = u.s.dbl_exp; - - /* - * If input is zero, Inf or NaN, just return it. - */ - if (u.v == 0.0 || oldexp == DBL_EXP_INFNAN) - return (val); - - if (oldexp == 0) { - /* - * u.v is denormal. We must adjust it so that the exponent - * arithmetic below will work. - */ - if (exp <= DBL_EXP_BIAS) { - /* - * Optimization: if the scaling can be done in a single - * multiply, or underflows, just do it now. - */ - if (exp <= -DBL_FRACBITS) { - errno = ERANGE; - return (0.0); - } - mul.v = 0.0; - mul.s.dbl_exp = exp + DBL_EXP_BIAS; - u.v *= mul.v; - if (u.v == 0.0) { - errno = ERANGE; - return (0.0); - } - return (u.v); - } else { - /* - * We know that exp is very large, and therefore the - * result cannot be denormal (though it may be Inf). - * Shift u.v by just enough to make it normal. - */ - mul.v = 0.0; - mul.s.dbl_exp = DBL_FRACBITS + DBL_EXP_BIAS; - u.v *= mul.v; - exp -= DBL_FRACBITS; - oldexp = u.s.dbl_exp; - } - } - - /* - * u.v is now normalized and oldexp has been adjusted if necessary. - * Calculate the new exponent and check for underflow and overflow. - */ - newexp = oldexp + exp; - - if (newexp <= 0) { - /* - * The output number is either denormal or underflows (see - * comments in machine/ieee.h). - */ - if (newexp <= -DBL_FRACBITS) { - errno = ERANGE; - return (0.0); - } - /* - * Denormalize the result. We do this with a multiply. If exp - * is very large, it won't fit in a double, so we have to - * adjust the exponent first. This is safe because we know - * that u.v is normal at this point. - */ - if (exp <= -DBL_EXP_BIAS) { - u.s.dbl_exp = 1; - exp += oldexp - 1; - } - mul.v = 0.0; - mul.s.dbl_exp = exp + DBL_EXP_BIAS; - u.v *= mul.v; - return (u.v); - } else if (newexp >= DBL_EXP_INFNAN) { - /* - * The result overflowed; return +/-Inf. - */ - u.s.dbl_exp = DBL_EXP_INFNAN; - u.s.dbl_frach = 0; - u.s.dbl_fracl = 0; - errno = ERANGE; - return (u.v); - } else { - /* - * The result is normal; just replace the old exponent with the - * new one. - */ - u.s.dbl_exp = newexp; - return (u.v); - } -} |