[BACK]Return to ldexp_ieee754.c CVS log [TXT][DIR] Up to [cvs.NetBSD.org] / src / lib / libc / gen

File: [cvs.NetBSD.org] / src / lib / libc / gen / ldexp_ieee754.c (download)

Revision 1.5, Mon Apr 28 20:22:59 2008 UTC (6 years, 3 months ago) by martin
Branch: MAIN
CVS Tags: yamt-pf42-base4, yamt-pf42-base3, yamt-pf42-base2, yamt-pagecache-tag8, yamt-pagecache-base9, yamt-pagecache-base8, yamt-pagecache-base7, yamt-pagecache-base6, yamt-pagecache-base5, yamt-pagecache-base4, yamt-pagecache-base3, yamt-pagecache-base2, yamt-pagecache-base, yamt-pagecache, wrstuden-revivesa-base-3, wrstuden-revivesa-base-2, wrstuden-revivesa-base-1, wrstuden-revivesa-base, wrstuden-revivesa, tls-maxphys-base, tls-maxphys, tls-earlyentropy-base, tls-earlyentropy, riastradh-xf86-video-intel-2-7-1-pre-2-21-15, riastradh-drm2-base3, riastradh-drm2-base2, riastradh-drm2-base1, riastradh-drm2-base, riastradh-drm2, netbsd-6-base, netbsd-6-1-RELEASE, netbsd-6-1-RC4, netbsd-6-1-RC3, netbsd-6-1-RC2, netbsd-6-1-RC1, netbsd-6-1-4-RELEASE, netbsd-6-1-3-RELEASE, netbsd-6-1-2-RELEASE, netbsd-6-1-1-RELEASE, netbsd-6-1, netbsd-6-0-RELEASE, netbsd-6-0-RC2, netbsd-6-0-RC1, netbsd-6-0-5-RELEASE, netbsd-6-0-4-RELEASE, netbsd-6-0-3-RELEASE, netbsd-6-0-2-RELEASE, netbsd-6-0-1-RELEASE, netbsd-6-0, netbsd-6, netbsd-5-base, netbsd-5-2-RELEASE, netbsd-5-2-RC1, netbsd-5-2-2-RELEASE, netbsd-5-2-1-RELEASE, netbsd-5-2, netbsd-5-1-RELEASE, netbsd-5-1-RC4, netbsd-5-1-RC3, netbsd-5-1-RC2, netbsd-5-1-RC1, netbsd-5-1-4-RELEASE, netbsd-5-1-3-RELEASE, netbsd-5-1-2-RELEASE, netbsd-5-1-1-RELEASE, netbsd-5-1, netbsd-5-0-RELEASE, netbsd-5-0-RC4, netbsd-5-0-RC3, netbsd-5-0-RC2, netbsd-5-0-RC1, netbsd-5-0-2-RELEASE, netbsd-5-0-1-RELEASE, netbsd-5-0, netbsd-5, matt-premerge-20091211, matt-nb6-plus-nbase, matt-nb6-plus-base, matt-nb6-plus, matt-nb5-pq3-base, matt-nb5-pq3, matt-nb5-mips64-u2-k2-k4-k7-k8-k9, matt-nb5-mips64-u1-k1-k5, matt-nb5-mips64-premerge-20101231, matt-nb5-mips64-premerge-20091211, matt-nb5-mips64-k15, matt-nb5-mips64, matt-nb4-mips64-k7-u2a-k9b, matt-mips64-premerge-20101231, matt-mips64-base2, jym-xensuspend-nbase, jym-xensuspend-base, jym-xensuspend, hpcarm-cleanup-nbase, christos-time_t-nbase, christos-time_t-base, cherry-xenmp-base, cherry-xenmp, bouyer-quota2-nbase, bouyer-quota2-base, bouyer-quota2, agc-symver-base, agc-symver, HEAD
Branch point for: christos-time_t
Changes since 1.4: +2 -9 lines

Remove clause 3 and 4 from TNF licenses

/*	$NetBSD: ldexp_ieee754.c,v 1.5 2008/04/28 20:22:59 martin Exp $	*/

/*-
 * 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.
 *
 * 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 defined(LIBC_SCCS) && !defined(lint)
__RCSID("$NetBSD: ldexp_ieee754.c,v 1.5 2008/04/28 20:22:59 martin Exp $");
#endif /* LIBC_SCCS and not lint */

#include <sys/types.h>
#include <machine/ieee.h>
#include <errno.h>
#include <math.h>

/*
 * Multiply the given value by 2^expon.
 */
double
ldexp(double val, int expon)
{
	int oldexp, newexp;
	union ieee_double_u u, mul;

	u.dblu_d = val;
	oldexp = u.dblu_dbl.dbl_exp;

	/*
	 * If input is zero, Inf or NaN, just return it.
	 */
	if (u.dblu_d == 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 (expon <= DBL_EXP_BIAS) {
			/*
			 * Optimization: if the scaling can be done in a single
			 * multiply, or underflows, just do it now.
			 */
			if (expon <= -DBL_FRACBITS) {
				errno = ERANGE;
				return (val < 0.0 ? -0.0 : 0.0);
			}
			mul.dblu_d = 0.0;
			mul.dblu_dbl.dbl_exp = expon + DBL_EXP_BIAS;
			u.dblu_d *= mul.dblu_d;
			if (u.dblu_d == 0.0) {
				errno = ERANGE;
				return (val < 0.0 ? -0.0 : 0.0);
			}
			return (u.dblu_d);
		} else {
			/*
			 * We know that expon 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.dblu_d = 0.0;
			mul.dblu_dbl.dbl_exp = DBL_FRACBITS + DBL_EXP_BIAS;
			u.dblu_d *= mul.dblu_d;
			expon -= DBL_FRACBITS;
			oldexp = u.dblu_dbl.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 + expon;

	if (newexp <= 0) {
		/*
		 * The output number is either denormal or underflows (see
		 * comments in machine/ieee.h).
		 */
		if (newexp <= -DBL_FRACBITS) {
			errno = ERANGE;
			return (val < 0.0 ? -0.0 : 0.0);
		}
		/*
		 * Denormalize the result.  We do this with a multiply.  If
		 * expon 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 (expon <= -DBL_EXP_BIAS) {
			u.dblu_dbl.dbl_exp = 1;
			expon += oldexp - 1;
		}
		mul.dblu_d = 0.0;
		mul.dblu_dbl.dbl_exp = expon + DBL_EXP_BIAS;
		u.dblu_d *= mul.dblu_d;
		return (u.dblu_d);
	} else if (newexp >= DBL_EXP_INFNAN) {
		/*
		 * The result overflowed; return +/-Inf.
		 */
		u.dblu_dbl.dbl_exp = DBL_EXP_INFNAN;
		u.dblu_dbl.dbl_frach = 0;
		u.dblu_dbl.dbl_fracl = 0;
		errno = ERANGE;
		return (u.dblu_d);
	} else {
		/*
		 * The result is normal; just replace the old exponent with the
		 * new one.
		 */
		u.dblu_dbl.dbl_exp = newexp;
		return (u.dblu_d);
	}
}