lib/libc/compat/gen/compat_ldexp_ieee754.c

   1 /* $NetBSD: compat_ldexp_ieee754.c,v 1.3 2008/09/28 15:19:09 christos 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 #ifdef USE_LIBM
  33 #include "s_finite.c"
  34 #include "s_scalbn.c"
  35 #include "s_ldexp.c"
  36 #else
  37 #include <sys/cdefs.h>
  38 #if defined(LIBC_SCCS) && !defined(lint)
  39 __RCSID("$NetBSD: compat_ldexp_ieee754.c,v 1.3 2008/09/28 15:19:09 christos Exp $");
  40 #endif /* LIBC_SCCS and not lint */
  41
  42 #include <sys/types.h>
  43 #include <machine/ieee.h>
  44 #include <errno.h>
  45 #include <math.h>
  46
  47 /*
  48  * Multiply the given value by 2^expon.
  49  */
  50 double
  51 ldexp(double val, int expon)
  52 {
  53         int oldexp, newexp;
  54         union ieee_double_u u, mul;
  55
  56         u.dblu_d = val;
  57         oldexp = u.dblu_dbl.dbl_exp;
  58
  59         /*
  60          * If input is zero, Inf or NaN, just return it.
  61          */
  62         if (u.dblu_d == 0.0 || oldexp == DBL_EXP_INFNAN)
  63                 return (val);
  64
  65         if (oldexp == 0) {
  66                 /*
  67                  * u.v is denormal.  We must adjust it so that the exponent
  68                  * arithmetic below will work.
  69                  */
  70                 if (expon <= DBL_EXP_BIAS) {
  71                         /*
  72                          * Optimization: if the scaling can be done in a single
  73                          * multiply, or underflows, just do it now.
  74                          */
  75                         if (expon <= -DBL_FRACBITS) {
  76                                 errno = ERANGE;
  77                                 return (val < 0.0 ? -0.0 : 0.0);
  78                         }
  79                         mul.dblu_d = 0.0;
  80                         mul.dblu_dbl.dbl_exp = expon + DBL_EXP_BIAS;
  81                         u.dblu_d *= mul.dblu_d;
  82                         if (u.dblu_d == 0.0) {
  83                                 errno = ERANGE;
  84                                 return (val < 0.0 ? -0.0 : 0.0);
  85                         }
  86                         return (u.dblu_d);
  87                 } else {
  88                         /*
  89                          * We know that expon is very large, and therefore the
  90                          * result cannot be denormal (though it may be Inf).
  91                          * Shift u.v by just enough to make it normal.
  92                          */
  93                         mul.dblu_d = 0.0;
  94                         mul.dblu_dbl.dbl_exp = DBL_FRACBITS + DBL_EXP_BIAS;
  95                         u.dblu_d *= mul.dblu_d;
  96                         expon -= DBL_FRACBITS;
  97                         oldexp = u.dblu_dbl.dbl_exp;
  98                 }
  99         }
 100
 101         /*
 102          * u.v is now normalized and oldexp has been adjusted if necessary.
 103          * Calculate the new exponent and check for underflow and overflow.
 104          */
 105         newexp = oldexp + expon;
 106
 107         if (newexp <= 0) {
 108                 /*
 109                  * The output number is either denormal or underflows (see
 110                  * comments in machine/ieee.h).
 111                  */
 112                 if (newexp <= -DBL_FRACBITS) {
 113                         errno = ERANGE;
 114                         return (val < 0.0 ? -0.0 : 0.0);
 115                 }
 116                 /*
 117                  * Denormalize the result.  We do this with a multiply.  If
 118                  * expon is very large, it won't fit in a double, so we have
 119                  * to adjust the exponent first.  This is safe because we know
 120                  * that u.v is normal at this point.
 121                  */
 122                 if (expon <= -DBL_EXP_BIAS) {
 123                         u.dblu_dbl.dbl_exp = 1;
 124                         expon += oldexp - 1;
 125                 }
 126                 mul.dblu_d = 0.0;
 127                 mul.dblu_dbl.dbl_exp = expon + DBL_EXP_BIAS;
 128                 u.dblu_d *= mul.dblu_d;
 129                 return (u.dblu_d);
 130         } else if (newexp >= DBL_EXP_INFNAN) {
 131                 /*
 132                  * The result overflowed; return +/-Inf.
 133                  */
 134                 u.dblu_dbl.dbl_exp = DBL_EXP_INFNAN;
 135                 u.dblu_dbl.dbl_frach = 0;
 136                 u.dblu_dbl.dbl_fracl = 0;
 137                 errno = ERANGE;
 138                 return (u.dblu_d);
 139         } else {
 140                 /*
 141                  * The result is normal; just replace the old exponent with the
 142                  * new one.
 143                  */
 144                 u.dblu_dbl.dbl_exp = newexp;
 145                 return (u.dblu_d);
 146         }
 147 }
 148 #endif