Home | History | Annotate | Line # | Download | only in gen 
      1 /* $NetBSD: compat_ldexp_ieee754.c,v 1.9 2024/06/17 18:16:58 riastradh 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  *
     19  * THIS SOFTWARE IS PROVIDED BY THE NETBSD FOUNDATION, INC. AND CONTRIBUTORS
     20  * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
     21  * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
     22  * PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR CONTRIBUTORS
     23  * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
     24  * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
     25  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
     26  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
     27  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
     28  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
     29  * POSSIBILITY OF SUCH DAMAGE.
     30  */
     31 
     32 #include <sys/cdefs.h>
     33 #if defined(LIBC_SCCS) && !defined(lint)
     34 __RCSID("$NetBSD: compat_ldexp_ieee754.c,v 1.9 2024/06/17 18:16:58 riastradh Exp $");
     35 #endif /* LIBC_SCCS and not lint */
     36 
     37 #include <sys/types.h>
     38 
     39 #include <machine/ieee.h>
     40 
     41 #include <errno.h>
     42 #include <float.h>
     43 #include <math.h>
     44 
     45 static volatile const double tiny = DBL_MIN, huge = DBL_MAX;
     46 
     47 static double
     48 underflow(double val)
     49 {
     50 
     51 	errno = ERANGE;
     52 	return (val < 0 ? -tiny*tiny : tiny*tiny);
     53 }
     54 
     55 static double
     56 overflow(double val)
     57 {
     58 
     59 	errno = ERANGE;
     60 	return (val < 0 ? -huge*huge : huge*huge);
     61 }
     62 
     63 /*
     64  * Multiply the given value by 2^expon.
     65  */
     66 double
     67 ldexp(double val, int expon)
     68 {
     69 	int oldexp, newexp;
     70 	union ieee_double_u u, mul;
     71 
     72 	u.dblu_d = val;
     73 	oldexp = u.dblu_dbl.dbl_exp;
     74 
     75 	/*
     76 	 * If input is zero, Inf or NaN, just return it, but raise
     77 	 * invalid exception if it is a signalling NaN: adding any of
     78 	 * these inputs to itself gives itself as output; arithmetic on
     79 	 * a signalling NaN additionally raises invalid-operation.
     80 	 */
     81 	if (u.dblu_d == 0.0 || oldexp == DBL_EXP_INFNAN)
     82 		return (val + val);
     83 
     84 	if (oldexp == 0) {
     85 		/*
     86 		 * u.v is denormal.  We must adjust it so that the exponent
     87 		 * arithmetic below will work.
     88 		 */
     89 		if (expon <= DBL_EXP_BIAS) {
     90 			/*
     91 			 * Optimization: if the scaling can be done in a single
     92 			 * multiply, or underflows, just do it now.
     93 			 */
     94 			if (expon <= -DBL_FRACBITS)
     95 				return underflow(val);
     96 			mul.dblu_d = 0.0;
     97 			mul.dblu_dbl.dbl_exp = expon + DBL_EXP_BIAS;
     98 			u.dblu_d *= mul.dblu_d;
     99 			if (u.dblu_d == 0.0)
    100 				return underflow(val);
    101 			return (u.dblu_d);
    102 		} else {
    103 			/*
    104 			 * We know that expon is very large, and therefore the
    105 			 * result cannot be denormal (though it may be Inf).
    106 			 * Shift u.v by just enough to make it normal.
    107 			 */
    108 			mul.dblu_d = 0.0;
    109 			mul.dblu_dbl.dbl_exp = DBL_FRACBITS + DBL_EXP_BIAS;
    110 			u.dblu_d *= mul.dblu_d;
    111 			expon -= DBL_FRACBITS;
    112 			oldexp = u.dblu_dbl.dbl_exp;
    113 		}
    114 	}
    115 
    116 	/*
    117 	 * u.v is now normalized and oldexp has been adjusted if necessary.
    118 	 * We have
    119 	 *
    120 	 *	0 <= oldexp <= DBL_EXP_INFNAN,
    121 	 *
    122 	 * but
    123 	 *
    124 	 *	INT_MIN <= expon <= INT_MAX.
    125 	 *
    126 	 * Check for underflow and overflow, and if none, calculate the
    127 	 * new exponent.
    128 	 */
    129 	if (expon >= DBL_EXP_INFNAN - oldexp) {
    130 		/*
    131 		 * The result overflowed; return +/-Inf.
    132 		 */
    133 		return overflow(val);
    134 	}
    135 
    136 	/*
    137 	 * We now have INT_MIN <= oldexp + expon <= DBL_EXP_INFNAN <= INT_MAX,
    138 	 * so the arithmetic is safe.
    139 	 */
    140 	newexp = oldexp + expon;
    141 
    142 	if (newexp <= 0) {
    143 		/*
    144 		 * The output number is either denormal or underflows (see
    145 		 * comments in machine/ieee.h).
    146 		 */
    147 		if (newexp <= -DBL_FRACBITS)
    148 			return underflow(val);
    149 		/*
    150 		 * Denormalize the result.  We do this with a multiply.  If
    151 		 * expon is very large, it won't fit in a double, so we have
    152 		 * to adjust the exponent first.  This is safe because we know
    153 		 * that u.v is normal at this point.
    154 		 */
    155 		if (expon <= -DBL_EXP_BIAS) {
    156 			u.dblu_dbl.dbl_exp = 1;
    157 			expon += oldexp - 1;
    158 		}
    159 		mul.dblu_d = 0.0;
    160 		mul.dblu_dbl.dbl_exp = expon + DBL_EXP_BIAS;
    161 		u.dblu_d *= mul.dblu_d;
    162 		return (u.dblu_d);
    163 	} else {
    164 		/*
    165 		 * The result is normal; just replace the old exponent with the
    166 		 * new one.
    167 		 */
    168 		u.dblu_dbl.dbl_exp = newexp;
    169 		return (u.dblu_d);
    170 	}
    171 }
    172