Home | History | Annotate | Line # | Download | only in gen
compat_ldexp_ieee754.c revision 1.1 
      1 /* $NetBSD: compat_ldexp_ieee754.c,v 1.1 2006/06/27 18:16:47 drochner Exp $ */
      2 
      3 /*-
      4  * Copyright (c) 1999 The NetBSD Foundation, Inc.
      5  * All rights reserved.
      6  *
      7  * This code is derived from software contributed to The NetBSD Foundation
      8  * by Charles M. Hannum.
      9  *
     10  * Redistribution and use in source and binary forms, with or without
     11  * modification, are permitted provided that the following conditions
     12  * are met:
     13  * 1. Redistributions of source code must retain the above copyright
     14  *    notice, this list of conditions and the following disclaimer.
     15  * 2. Redistributions in binary form must reproduce the above copyright
     16  *    notice, this list of conditions and the following disclaimer in the
     17  *    documentation and/or other materials provided with the distribution.
     18  * 3. All advertising materials mentioning features or use of this software
     19  *    must display the following acknowledgement:
     20  *        This product includes software developed by the NetBSD
     21  *        Foundation, Inc. and its contributors.
     22  * 4. Neither the name of The NetBSD Foundation nor the names of its
     23  *    contributors may be used to endorse or promote products derived
     24  *    from this software without specific prior written permission.
     25  *
     26  * THIS SOFTWARE IS PROVIDED BY THE NETBSD FOUNDATION, INC. AND CONTRIBUTORS
     27  * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
     28  * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
     29  * PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR CONTRIBUTORS
     30  * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
     31  * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
     32  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
     33  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
     34  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
     35  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
     36  * POSSIBILITY OF SUCH DAMAGE.
     37  */
     38 
     39 #include <sys/cdefs.h>
     40 #if defined(LIBC_SCCS) && !defined(lint)
     41 __RCSID("$NetBSD: compat_ldexp_ieee754.c,v 1.1 2006/06/27 18:16:47 drochner Exp $");
     42 #endif /* LIBC_SCCS and not lint */
     43 
     44 #include <sys/types.h>
     45 #include <machine/ieee.h>
     46 #include <errno.h>
     47 #include <math.h>
     48 
     49 /*
     50  * Multiply the given value by 2^expon.
     51  */
     52 double
     53 ldexp(double val, int expon)
     54 {
     55 	int oldexp, newexp;
     56 	union ieee_double_u u, mul;
     57 
     58 	u.dblu_d = val;
     59 	oldexp = u.dblu_dbl.dbl_exp;
     60 
     61 	/*
     62 	 * If input is zero, Inf or NaN, just return it.
     63 	 */
     64 	if (u.dblu_d == 0.0 || oldexp == DBL_EXP_INFNAN)
     65 		return (val);
     66 
     67 	if (oldexp == 0) {
     68 		/*
     69 		 * u.v is denormal.  We must adjust it so that the exponent
     70 		 * arithmetic below will work.
     71 		 */
     72 		if (expon <= DBL_EXP_BIAS) {
     73 			/*
     74 			 * Optimization: if the scaling can be done in a single
     75 			 * multiply, or underflows, just do it now.
     76 			 */
     77 			if (expon <= -DBL_FRACBITS) {
     78 				errno = ERANGE;
     79 				return (val < 0.0 ? -0.0 : 0.0);
     80 			}
     81 			mul.dblu_d = 0.0;
     82 			mul.dblu_dbl.dbl_exp = expon + DBL_EXP_BIAS;
     83 			u.dblu_d *= mul.dblu_d;
     84 			if (u.dblu_d == 0.0) {
     85 				errno = ERANGE;
     86 				return (val < 0.0 ? -0.0 : 0.0);
     87 			}
     88 			return (u.dblu_d);
     89 		} else {
     90 			/*
     91 			 * We know that expon is very large, and therefore the
     92 			 * result cannot be denormal (though it may be Inf).
     93 			 * Shift u.v by just enough to make it normal.
     94 			 */
     95 			mul.dblu_d = 0.0;
     96 			mul.dblu_dbl.dbl_exp = DBL_FRACBITS + DBL_EXP_BIAS;
     97 			u.dblu_d *= mul.dblu_d;
     98 			expon -= DBL_FRACBITS;
     99 			oldexp = u.dblu_dbl.dbl_exp;
    100 		}
    101 	}
    102 
    103 	/*
    104 	 * u.v is now normalized and oldexp has been adjusted if necessary.
    105 	 * Calculate the new exponent and check for underflow and overflow.
    106 	 */
    107 	newexp = oldexp + expon;
    108 
    109 	if (newexp <= 0) {
    110 		/*
    111 		 * The output number is either denormal or underflows (see
    112 		 * comments in machine/ieee.h).
    113 		 */
    114 		if (newexp <= -DBL_FRACBITS) {
    115 			errno = ERANGE;
    116 			return (val < 0.0 ? -0.0 : 0.0);
    117 		}
    118 		/*
    119 		 * Denormalize the result.  We do this with a multiply.  If
    120 		 * expon is very large, it won't fit in a double, so we have
    121 		 * to adjust the exponent first.  This is safe because we know
    122 		 * that u.v is normal at this point.
    123 		 */
    124 		if (expon <= -DBL_EXP_BIAS) {
    125 			u.dblu_dbl.dbl_exp = 1;
    126 			expon += oldexp - 1;
    127 		}
    128 		mul.dblu_d = 0.0;
    129 		mul.dblu_dbl.dbl_exp = expon + DBL_EXP_BIAS;
    130 		u.dblu_d *= mul.dblu_d;
    131 		return (u.dblu_d);
    132 	} else if (newexp >= DBL_EXP_INFNAN) {
    133 		/*
    134 		 * The result overflowed; return +/-Inf.
    135 		 */
    136 		u.dblu_dbl.dbl_exp = DBL_EXP_INFNAN;
    137 		u.dblu_dbl.dbl_frach = 0;
    138 		u.dblu_dbl.dbl_fracl = 0;
    139 		errno = ERANGE;
    140 		return (u.dblu_d);
    141 	} else {
    142 		/*
    143 		 * The result is normal; just replace the old exponent with the
    144 		 * new one.
    145 		 */
    146 		u.dblu_dbl.dbl_exp = newexp;
    147 		return (u.dblu_d);
    148 	}
    149 }
    150