Use round-to-nearest internally in jn, test with ALL_RM_TEST (bug 18602).
[glibc.git] / sysdeps / ieee754 / ldbl-96 / s_scalblnl.c
blob457e999c6c9c36737e4159d909d2d8b003cd3c0a
1 /* s_scalbnl.c -- long double version of s_scalbn.c.
2 * Conversion to long double by Ulrich Drepper,
3 * Cygnus Support, drepper@cygnus.com.
4 */
6 /*
7 * ====================================================
8 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
10 * Developed at SunPro, a Sun Microsystems, Inc. business.
11 * Permission to use, copy, modify, and distribute this
12 * software is freely granted, provided that this notice
13 * is preserved.
14 * ====================================================
18 * scalbnl (long double x, int n)
19 * scalbnl(x,n) returns x* 2**n computed by exponent
20 * manipulation rather than by actually performing an
21 * exponentiation or a multiplication.
24 #include <math.h>
25 #include <math_private.h>
27 static const long double
28 two63 = 0x1p63L,
29 twom64 = 0x1p-64L,
30 huge = 1.0e+4900L,
31 tiny = 1.0e-4900L;
33 long double
34 __scalblnl (long double x, long int n)
36 int32_t k,es,hx,lx;
37 GET_LDOUBLE_WORDS(es,hx,lx,x);
38 k = es&0x7fff; /* extract exponent */
39 if (__builtin_expect(k==0, 0)) { /* 0 or subnormal x */
40 if ((lx|(hx&0x7fffffff))==0) return x; /* +-0 */
41 x *= two63;
42 GET_LDOUBLE_EXP(es,x);
43 k = (es&0x7fff) - 63;
45 if (__builtin_expect(k==0x7fff, 0)) return x+x; /* NaN or Inf */
46 if (__builtin_expect(n< -50000, 0))
47 return tiny*__copysignl(tiny,x);
48 if (__builtin_expect(n> 50000 || k+n > 0x7ffe, 0))
49 return huge*__copysignl(huge,x); /* overflow */
50 /* Now k and n are bounded we know that k = k+n does not
51 overflow. */
52 k = k+n;
53 if (__builtin_expect(k > 0, 1)) /* normal result */
54 {SET_LDOUBLE_EXP(x,(es&0x8000)|k); return x;}
55 if (k <= -64)
56 return tiny*__copysignl(tiny,x); /*underflow*/
57 k += 64; /* subnormal result */
58 SET_LDOUBLE_EXP(x,(es&0x8000)|k);
59 return x*twom64;