src/libm/s_scalbn.c
changeset 11683 48bcba563d9c
parent 7678 286c42d7c5ed
child 11684 eccdf37f8996
     1.1 --- a/src/libm/s_scalbn.c	Sat Nov 04 15:34:14 2017 -0700
     1.2 +++ b/src/libm/s_scalbn.c	Sat Nov 04 15:53:19 2017 -0700
     1.3 @@ -1,4 +1,3 @@
     1.4 -/* @(#)s_scalbn.c 5.1 93/09/24 */
     1.5  /*
     1.6   * ====================================================
     1.7   * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
     1.8 @@ -10,70 +9,69 @@
     1.9   * ====================================================
    1.10   */
    1.11  
    1.12 -#if defined(LIBM_SCCS) && !defined(lint)
    1.13 -static const char rcsid[] =
    1.14 -    "$NetBSD: s_scalbn.c,v 1.8 1995/05/10 20:48:08 jtc Exp $";
    1.15 -#endif
    1.16 -
    1.17  /*
    1.18 - * scalbn (double x, int n)
    1.19 - * scalbn(x,n) returns x* 2**n  computed by  exponent
    1.20 + * scalbln(double x, long n)
    1.21 + * scalbln(x,n) returns x * 2**n computed by exponent
    1.22   * manipulation rather than by actually performing an
    1.23   * exponentiation or a multiplication.
    1.24   */
    1.25  
    1.26  #include "math_libm.h"
    1.27  #include "math_private.h"
    1.28 +#include <limits.h>
    1.29  
    1.30 -libm_hidden_proto(copysign)
    1.31 -#ifdef __STDC__
    1.32 -     static const double
    1.33 +static const double
    1.34 +two54  = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
    1.35 +twom54 = 5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */
    1.36 +huge   = 1.0e+300,
    1.37 +tiny   = 1.0e-300;
    1.38 +
    1.39 +double scalbln(double x, long n)
    1.40 +{
    1.41 +	int32_t k, hx, lx;
    1.42 +
    1.43 +	EXTRACT_WORDS(hx, lx, x);
    1.44 +	k = (hx & 0x7ff00000) >> 20; /* extract exponent */
    1.45 +	if (k == 0) { /* 0 or subnormal x */
    1.46 +		if ((lx | (hx & 0x7fffffff)) == 0)
    1.47 +			return x; /* +-0 */
    1.48 +		x *= two54;
    1.49 +		GET_HIGH_WORD(hx, x);
    1.50 +		k = ((hx & 0x7ff00000) >> 20) - 54;
    1.51 +	}
    1.52 +	if (k == 0x7ff)
    1.53 +		return x + x; /* NaN or Inf */
    1.54 +	k = k + n;
    1.55 +	if (k > 0x7fe)
    1.56 +		return huge * copysign(huge, x); /* overflow */
    1.57 +	if (n < -50000)
    1.58 +		return tiny * copysign(tiny, x); /* underflow */
    1.59 +	if (k > 0) { /* normal result */
    1.60 +		SET_HIGH_WORD(x, (hx & 0x800fffff) | (k << 20));
    1.61 +		return x;
    1.62 +	}
    1.63 +	if (k <= -54) {
    1.64 +		if (n > 50000) /* in case integer overflow in n+k */
    1.65 +			return huge * copysign(huge, x); /* overflow */
    1.66 +		return tiny * copysign(tiny, x); /* underflow */
    1.67 +	}
    1.68 +	k += 54; /* subnormal result */
    1.69 +	SET_HIGH_WORD(x, (hx & 0x800fffff) | (k << 20));
    1.70 +	return x * twom54;
    1.71 +}
    1.72 +libm_hidden_def(scalbln)
    1.73 +
    1.74 +#if LONG_MAX == INT_MAX
    1.75 +/* strong_alias(scalbln, scalbn) - "error: conflicting types for 'scalbn'"
    1.76 + * because it tries to declare "typeof(scalbln) scalbn;"
    1.77 + * which tries to give "long" parameter to scalbn.
    1.78 + * Doing it by hand:
    1.79 + */
    1.80 +__typeof(scalbn) scalbn __attribute__((alias("scalbln")));
    1.81  #else
    1.82 -     static double
    1.83 +double scalbn(double x, int n)
    1.84 +{
    1.85 +	return scalbln(x, n);
    1.86 +}
    1.87  #endif
    1.88 -       two54 = 1.80143985094819840000e+16,      /* 0x43500000, 0x00000000 */
    1.89 -         twom54 = 5.55111512312578270212e-17,   /* 0x3C900000, 0x00000000 */
    1.90 -         huge_val = 1.0e+300, tiny = 1.0e-300;
    1.91 -
    1.92 -libm_hidden_proto(scalbn)
    1.93 -#ifdef __STDC__
    1.94 -     double scalbn(double x, int n)
    1.95 -#else
    1.96 -     double scalbn(x, n)
    1.97 -     double x;
    1.98 -     int n;
    1.99 -#endif
   1.100 -{
   1.101 -    int32_t k, hx, lx;
   1.102 -    EXTRACT_WORDS(hx, lx, x);
   1.103 -    k = (hx & 0x7ff00000) >> 20;        /* extract exponent */
   1.104 -    if (k == 0) {               /* 0 or subnormal x */
   1.105 -        if ((lx | (hx & 0x7fffffff)) == 0)
   1.106 -            return x;           /* +-0 */
   1.107 -        x *= two54;
   1.108 -        GET_HIGH_WORD(hx, x);
   1.109 -        k = ((hx & 0x7ff00000) >> 20) - 54;
   1.110 -        if (n < -50000)
   1.111 -            return tiny * x;    /* underflow */
   1.112 -    }
   1.113 -    if (k == 0x7ff)
   1.114 -        return x + x;           /* NaN or Inf */
   1.115 -    k = k + n;
   1.116 -    if (k > 0x7fe)
   1.117 -        return huge_val * copysign(huge_val, x);        /* overflow  */
   1.118 -    if (k > 0) {                /* normal result */
   1.119 -        SET_HIGH_WORD(x, (hx & 0x800fffff) | (k << 20));
   1.120 -        return x;
   1.121 -    }
   1.122 -    if (k <= -54) {
   1.123 -        if (n > 50000)          /* in case integer overflow in n+k */
   1.124 -            return huge_val * copysign(huge_val, x);    /* overflow */
   1.125 -        else
   1.126 -            return tiny * copysign(tiny, x);    /* underflow */
   1.127 -    }
   1.128 -    k += 54;                    /* subnormal result */
   1.129 -    SET_HIGH_WORD(x, (hx & 0x800fffff) | (k << 20));
   1.130 -    return x * twom54;
   1.131 -}
   1.132 -
   1.133  libm_hidden_def(scalbn)