1 /***********************************************************************
2 *                                                                      *
3 *               This software is part of the ast package               *
4 *          Copyright (c) 1985-2011 AT&T Intellectual Property          *
5 *                      and is licensed under the                       *
6 *                 Eclipse Public License, Version 1.0                  *
7 *                    by AT&T Intellectual Property                     *
8 *                                                                      *
9 *                A copy of the License is available at                 *
10 *          http://www.eclipse.org/org/documents/epl-v10.html           *
11 *         (with md5 checksum b35adb5213ca9657e911e9befb180842)         *
12 *                                                                      *
13 *              Information and Software Systems Research               *
14 *                            AT&T Research                             *
15 *                           Florham Park NJ                            *
16 *                                                                      *
17 *                 Glenn Fowler <gsf@research.att.com>                  *
18 *                  David Korn <dgk@research.att.com>                   *
19 *                   Phong Vo <kpv@research.att.com>                    *
20 *                                                                      *
21 ***********************************************************************/
22 #pragma prototyped
23 
24 /*
25  * frexp/ldexp implementation
26  */
27 
28 #include <ast.h>
29 
30 #include "FEATURE/float"
31 
32 #if _lib_frexp && _lib_ldexp
33 
34 NoN(frexp)
35 
36 #else
37 
38 #if defined(_ast_dbl_exp_index) && defined(_ast_dbl_exp_shift)
39 
40 #define INIT()		_ast_dbl_exp_t _pow_
41 #define pow2(i)		(_pow_.f=1,_pow_.e[_ast_dbl_exp_index]+=((i)<<_ast_dbl_exp_shift),_pow_.f)
42 
43 #else
44 
45 static double		pow2tab[DBL_MAX_EXP + 1];
46 
47 static int
48 init(void)
49 {
50 	register int	x;
51 	double		g;
52 
53 	g = 1;
54 	for (x = 0; x < elementsof(pow2tab); x++)
55 	{
56 		pow2tab[x] = g;
57 		g *= 2;
58 	}
59 	return 0;
60 }
61 
62 #define INIT()		(pow2tab[0]?0:init())
63 
64 #define pow2(i)		pow2tab[i]
65 
66 #endif
67 
68 #if !_lib_frexp
69 
70 extern double
71 frexp(double f, int* p)
72 {
73 	register int	k;
74 	register int	x;
75 	double		g;
76 
77 	INIT();
78 
79 	/*
80 	 * normalize
81 	 */
82 
83 	x = k = DBL_MAX_EXP / 2;
84 	if (f < 1)
85 	{
86 		g = 1.0L / f;
87 		for (;;)
88 		{
89 			k = (k + 1) / 2;
90 			if (g < pow2(x))
91 				x -= k;
92 			else if (k == 1 && g < pow2(x+1))
93 				break;
94 			else
95 				x += k;
96 		}
97 		if (g == pow2(x))
98 			x--;
99 		x = -x;
100 	}
101 	else if (f > 1)
102 	{
103 		for (;;)
104 		{
105 			k = (k + 1) / 2;
106 			if (f > pow2(x))
107 				x += k;
108 			else if (k == 1 && f > pow2(x-1))
109 				break;
110 			else
111 				x -= k;
112 		}
113 		if (f == pow2(x))
114 			x++;
115 	}
116 	else
117 		x = 1;
118 	*p = x;
119 
120 	/*
121 	 * shift
122 	 */
123 
124 	x = -x;
125 	if (x < 0)
126 		f /= pow2(-x);
127 	else if (x < DBL_MAX_EXP)
128 		f *= pow2(x);
129 	else
130 		f = (f * pow2(DBL_MAX_EXP - 1)) * pow2(x - (DBL_MAX_EXP - 1));
131 	return f;
132 }
133 
134 #endif
135 
136 #if !_lib_ldexp
137 
138 extern double
139 ldexp(double f, register int x)
140 {
141 	INIT();
142 	if (x < 0)
143 		f /= pow2(-x);
144 	else if (x < DBL_MAX_EXP)
145 		f *= pow2(x);
146 	else
147 		f = (f * pow2(DBL_MAX_EXP - 1)) * pow2(x - (DBL_MAX_EXP - 1));
148 	return f;
149 }
150 
151 #endif
152 
153 #endif
154