src/libm/s_scalbn.c
author Sam Lantinga <slouken@libsdl.org>
Sat, 04 Nov 2017 15:53:19 -0700
changeset 11683 48bcba563d9c
parent 7678 286c42d7c5ed
child 11684 eccdf37f8996
permissions -rw-r--r--
Updated math code from the uClibc 0.9.33 release
slouken@2756
     1
/*
slouken@2756
     2
 * ====================================================
slouken@2756
     3
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
slouken@2756
     4
 *
slouken@2756
     5
 * Developed at SunPro, a Sun Microsystems, Inc. business.
slouken@2756
     6
 * Permission to use, copy, modify, and distribute this
slouken@2756
     7
 * software is freely granted, provided that this notice
slouken@2756
     8
 * is preserved.
slouken@2756
     9
 * ====================================================
slouken@2756
    10
 */
slouken@2756
    11
slouken@2756
    12
/*
slouken@11683
    13
 * scalbln(double x, long n)
slouken@11683
    14
 * scalbln(x,n) returns x * 2**n computed by exponent
slouken@2756
    15
 * manipulation rather than by actually performing an
slouken@2756
    16
 * exponentiation or a multiplication.
slouken@2756
    17
 */
slouken@2756
    18
slouken@6044
    19
#include "math_libm.h"
slouken@2756
    20
#include "math_private.h"
slouken@11683
    21
#include <limits.h>
slouken@2756
    22
slouken@11683
    23
static const double
slouken@11683
    24
two54  = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
slouken@11683
    25
twom54 = 5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */
slouken@11683
    26
huge   = 1.0e+300,
slouken@11683
    27
tiny   = 1.0e-300;
slouken@11683
    28
slouken@11683
    29
double scalbln(double x, long n)
slouken@11683
    30
{
slouken@11683
    31
	int32_t k, hx, lx;
slouken@11683
    32
slouken@11683
    33
	EXTRACT_WORDS(hx, lx, x);
slouken@11683
    34
	k = (hx & 0x7ff00000) >> 20; /* extract exponent */
slouken@11683
    35
	if (k == 0) { /* 0 or subnormal x */
slouken@11683
    36
		if ((lx | (hx & 0x7fffffff)) == 0)
slouken@11683
    37
			return x; /* +-0 */
slouken@11683
    38
		x *= two54;
slouken@11683
    39
		GET_HIGH_WORD(hx, x);
slouken@11683
    40
		k = ((hx & 0x7ff00000) >> 20) - 54;
slouken@11683
    41
	}
slouken@11683
    42
	if (k == 0x7ff)
slouken@11683
    43
		return x + x; /* NaN or Inf */
slouken@11683
    44
	k = k + n;
slouken@11683
    45
	if (k > 0x7fe)
slouken@11683
    46
		return huge * copysign(huge, x); /* overflow */
slouken@11683
    47
	if (n < -50000)
slouken@11683
    48
		return tiny * copysign(tiny, x); /* underflow */
slouken@11683
    49
	if (k > 0) { /* normal result */
slouken@11683
    50
		SET_HIGH_WORD(x, (hx & 0x800fffff) | (k << 20));
slouken@11683
    51
		return x;
slouken@11683
    52
	}
slouken@11683
    53
	if (k <= -54) {
slouken@11683
    54
		if (n > 50000) /* in case integer overflow in n+k */
slouken@11683
    55
			return huge * copysign(huge, x); /* overflow */
slouken@11683
    56
		return tiny * copysign(tiny, x); /* underflow */
slouken@11683
    57
	}
slouken@11683
    58
	k += 54; /* subnormal result */
slouken@11683
    59
	SET_HIGH_WORD(x, (hx & 0x800fffff) | (k << 20));
slouken@11683
    60
	return x * twom54;
slouken@11683
    61
}
slouken@11683
    62
libm_hidden_def(scalbln)
slouken@11683
    63
slouken@11683
    64
#if LONG_MAX == INT_MAX
slouken@11683
    65
/* strong_alias(scalbln, scalbn) - "error: conflicting types for 'scalbn'"
slouken@11683
    66
 * because it tries to declare "typeof(scalbln) scalbn;"
slouken@11683
    67
 * which tries to give "long" parameter to scalbn.
slouken@11683
    68
 * Doing it by hand:
slouken@11683
    69
 */
slouken@11683
    70
__typeof(scalbn) scalbn __attribute__((alias("scalbln")));
slouken@2756
    71
#else
slouken@11683
    72
double scalbn(double x, int n)
slouken@11683
    73
{
slouken@11683
    74
	return scalbln(x, n);
slouken@11683
    75
}
slouken@2756
    76
#endif
slouken@2756
    77
libm_hidden_def(scalbn)