Import OpenBSD's libm (trunk, 4 July 2015) to a new vendor branch
[dragonfly.git] / contrib / openbsd_libm / src / s_scalbnf.c
1 /* s_scalbnf.c -- float version of s_scalbn.c.
2  * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
3  */
4
5 /*
6  * ====================================================
7  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
8  *
9  * Developed at SunPro, a Sun Microsystems, Inc. business.
10  * Permission to use, copy, modify, and distribute this
11  * software is freely granted, provided that this notice 
12  * is preserved.
13  * ====================================================
14  */
15
16 #include "math.h"
17 #include "math_private.h"
18
19 static const float
20 two25   =  3.355443200e+07,     /* 0x4c000000 */
21 twom25  =  2.9802322388e-08,    /* 0x33000000 */
22 huge   = 1.0e+30,
23 tiny   = 1.0e-30;
24
25 float
26 scalbnf(float x, int n)
27 {
28         int32_t k,ix;
29         GET_FLOAT_WORD(ix,x);
30         k = (ix&0x7f800000)>>23;                /* extract exponent */
31         if (k==0) {                             /* 0 or subnormal x */
32             if ((ix&0x7fffffff)==0) return x; /* +-0 */
33             x *= two25;
34             GET_FLOAT_WORD(ix,x);
35             k = ((ix&0x7f800000)>>23) - 25; 
36             if (n< -50000) return tiny*x;       /*underflow*/
37             }
38         if (k==0xff) return x+x;                /* NaN or Inf */
39         k = k+n; 
40         if (k >  0xfe) return huge*copysignf(huge,x); /* overflow  */
41         if (k > 0)                              /* normal result */
42             {SET_FLOAT_WORD(x,(ix&0x807fffff)|(k<<23)); return x;}
43         if (k <= -25)
44             if (n > 50000)      /* in case integer overflow in n+k */
45                 return huge*copysignf(huge,x);  /*overflow*/
46             else return tiny*copysignf(tiny,x); /*underflow*/
47         k += 25;                                /* subnormal result */
48         SET_FLOAT_WORD(x,(ix&0x807fffff)|(k<<23));
49         return x*twom25;
50 }
51
52 float
53 ldexpf(float x, int n)
54 {
55         return scalbnf(x, n);
56 }