sysdeps/ieee754/ldbl-96/s_scalblnl.c

   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  */
   5
   6 /*
   7  * ====================================================
   8  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   9  *
  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  * ====================================================
  15  */
  16
  17 /*
  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.
  22  */
  23
  24 #include <math.h>
  25 #include <math_private.h>
  26
  27 static const long double
  28 two63   =  4.50359962737049600000e+15,
  29 twom63  =  1.08420217248550443400e-19,
  30 huge   = 1.0e+4900L,
  31 tiny   = 1.0e-4900L;
  32
  33 long double
  34 __scalblnl (long double x, long int n)
  35 {
  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 = (hx&0x7fff) - 63;
  44             }
  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 <= -63)
  56             return tiny*__copysignl(tiny,x);    /*underflow*/
  57         k += 63;                                /* subnormal result */
  58         SET_LDOUBLE_EXP(x,(es&0x8000)|k);
  59         return x*twom63;
  60 }
  61 weak_alias (__scalblnl, scalblnl)