1*6ff43c94SPeter Avalos /* @(#)s_scalbn.c 5.1 93/09/24 */
2*6ff43c94SPeter Avalos /* @(#)fdlibm.h 5.1 93/09/24 */
3*6ff43c94SPeter Avalos /* $FreeBSD: head/lib/libc/gen/ldexp.c 140607 2005-01-22 06:03:40Z das $ */
4*6ff43c94SPeter Avalos /*
5*6ff43c94SPeter Avalos * ====================================================
6*6ff43c94SPeter Avalos * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
7*6ff43c94SPeter Avalos *
8*6ff43c94SPeter Avalos * Developed at SunPro, a Sun Microsystems, Inc. business.
9*6ff43c94SPeter Avalos * Permission to use, copy, modify, and distribute this
10*6ff43c94SPeter Avalos * software is freely granted, provided that this notice
11*6ff43c94SPeter Avalos * is preserved.
12*6ff43c94SPeter Avalos * ====================================================
13*6ff43c94SPeter Avalos */
14*6ff43c94SPeter Avalos
15*6ff43c94SPeter Avalos #include <sys/types.h>
16*6ff43c94SPeter Avalos #include <machine/endian.h>
17*6ff43c94SPeter Avalos #include <math.h>
18*6ff43c94SPeter Avalos
19*6ff43c94SPeter Avalos /* Bit fiddling routines copied from msun/src/math_private.h,v 1.15 */
20*6ff43c94SPeter Avalos
21*6ff43c94SPeter Avalos #if BYTE_ORDER == BIG_ENDIAN
22*6ff43c94SPeter Avalos
23*6ff43c94SPeter Avalos typedef union
24*6ff43c94SPeter Avalos {
25*6ff43c94SPeter Avalos double value;
26*6ff43c94SPeter Avalos struct
27*6ff43c94SPeter Avalos {
28*6ff43c94SPeter Avalos u_int32_t msw;
29*6ff43c94SPeter Avalos u_int32_t lsw;
30*6ff43c94SPeter Avalos } parts;
31*6ff43c94SPeter Avalos } ieee_double_shape_type;
32*6ff43c94SPeter Avalos
33*6ff43c94SPeter Avalos #endif
34*6ff43c94SPeter Avalos
35*6ff43c94SPeter Avalos #if BYTE_ORDER == LITTLE_ENDIAN
36*6ff43c94SPeter Avalos
37*6ff43c94SPeter Avalos typedef union
38*6ff43c94SPeter Avalos {
39*6ff43c94SPeter Avalos double value;
40*6ff43c94SPeter Avalos struct
41*6ff43c94SPeter Avalos {
42*6ff43c94SPeter Avalos u_int32_t lsw;
43*6ff43c94SPeter Avalos u_int32_t msw;
44*6ff43c94SPeter Avalos } parts;
45*6ff43c94SPeter Avalos } ieee_double_shape_type;
46*6ff43c94SPeter Avalos
47*6ff43c94SPeter Avalos #endif
48*6ff43c94SPeter Avalos
49*6ff43c94SPeter Avalos /* Get two 32 bit ints from a double. */
50*6ff43c94SPeter Avalos
51*6ff43c94SPeter Avalos #define EXTRACT_WORDS(ix0,ix1,d) \
52*6ff43c94SPeter Avalos do { \
53*6ff43c94SPeter Avalos ieee_double_shape_type ew_u; \
54*6ff43c94SPeter Avalos ew_u.value = (d); \
55*6ff43c94SPeter Avalos (ix0) = ew_u.parts.msw; \
56*6ff43c94SPeter Avalos (ix1) = ew_u.parts.lsw; \
57*6ff43c94SPeter Avalos } while (0)
58*6ff43c94SPeter Avalos
59*6ff43c94SPeter Avalos /* Get the more significant 32 bit int from a double. */
60*6ff43c94SPeter Avalos
61*6ff43c94SPeter Avalos #define GET_HIGH_WORD(i,d) \
62*6ff43c94SPeter Avalos do { \
63*6ff43c94SPeter Avalos ieee_double_shape_type gh_u; \
64*6ff43c94SPeter Avalos gh_u.value = (d); \
65*6ff43c94SPeter Avalos (i) = gh_u.parts.msw; \
66*6ff43c94SPeter Avalos } while (0)
67*6ff43c94SPeter Avalos
68*6ff43c94SPeter Avalos /* Set the more significant 32 bits of a double from an int. */
69*6ff43c94SPeter Avalos
70*6ff43c94SPeter Avalos #define SET_HIGH_WORD(d,v) \
71*6ff43c94SPeter Avalos do { \
72*6ff43c94SPeter Avalos ieee_double_shape_type sh_u; \
73*6ff43c94SPeter Avalos sh_u.value = (d); \
74*6ff43c94SPeter Avalos sh_u.parts.msw = (v); \
75*6ff43c94SPeter Avalos (d) = sh_u.value; \
76*6ff43c94SPeter Avalos } while (0)
77*6ff43c94SPeter Avalos
78*6ff43c94SPeter Avalos
79*6ff43c94SPeter Avalos static const double
80*6ff43c94SPeter Avalos two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
81*6ff43c94SPeter Avalos twom54 = 5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */
82*6ff43c94SPeter Avalos huge = 1.0e+300,
83*6ff43c94SPeter Avalos tiny = 1.0e-300;
84*6ff43c94SPeter Avalos
85*6ff43c94SPeter Avalos static double
_copysign(double x,double y)86*6ff43c94SPeter Avalos _copysign(double x, double y)
87*6ff43c94SPeter Avalos {
88*6ff43c94SPeter Avalos u_int32_t hx,hy;
89*6ff43c94SPeter Avalos GET_HIGH_WORD(hx,x);
90*6ff43c94SPeter Avalos GET_HIGH_WORD(hy,y);
91*6ff43c94SPeter Avalos SET_HIGH_WORD(x,(hx&0x7fffffff)|(hy&0x80000000));
92*6ff43c94SPeter Avalos return x;
93*6ff43c94SPeter Avalos }
94*6ff43c94SPeter Avalos
95*6ff43c94SPeter Avalos double
ldexp(double x,int n)96*6ff43c94SPeter Avalos ldexp(double x, int n)
97*6ff43c94SPeter Avalos {
98*6ff43c94SPeter Avalos int32_t k,hx,lx;
99*6ff43c94SPeter Avalos EXTRACT_WORDS(hx,lx,x);
100*6ff43c94SPeter Avalos k = (hx&0x7ff00000)>>20; /* extract exponent */
101*6ff43c94SPeter Avalos if (k==0) { /* 0 or subnormal x */
102*6ff43c94SPeter Avalos if ((lx|(hx&0x7fffffff))==0) return x; /* +-0 */
103*6ff43c94SPeter Avalos x *= two54;
104*6ff43c94SPeter Avalos GET_HIGH_WORD(hx,x);
105*6ff43c94SPeter Avalos k = ((hx&0x7ff00000)>>20) - 54;
106*6ff43c94SPeter Avalos if (n< -50000) return tiny*x; /*underflow*/
107*6ff43c94SPeter Avalos }
108*6ff43c94SPeter Avalos if (k==0x7ff) return x+x; /* NaN or Inf */
109*6ff43c94SPeter Avalos k = k+n;
110*6ff43c94SPeter Avalos if (k > 0x7fe) return huge*_copysign(huge,x); /* overflow */
111*6ff43c94SPeter Avalos if (k > 0) /* normal result */
112*6ff43c94SPeter Avalos {SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20)); return x;}
113*6ff43c94SPeter Avalos if (k <= -54) {
114*6ff43c94SPeter Avalos if (n > 50000) /* in case integer overflow in n+k */
115*6ff43c94SPeter Avalos return huge*_copysign(huge,x); /*overflow*/
116*6ff43c94SPeter Avalos else return tiny*_copysign(tiny,x); /*underflow*/
117*6ff43c94SPeter Avalos }
118*6ff43c94SPeter Avalos k += 54; /* subnormal result */
119*6ff43c94SPeter Avalos SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20));
120*6ff43c94SPeter Avalos return x*twom54;
121*6ff43c94SPeter Avalos }
122