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*                  Common Public License, Version 1.0                  *
7*                    by AT&T Intellectual Property                     *
8*                                                                      *
9*                A copy of the License is available at                 *
10*            http://www.opensource.org/licenses/cpl1.0.txt             *
11*         (with md5 checksum 059e8cd6165cb4c31e351f2b69388fd9)         *
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
34NoN(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
45static double		pow2tab[DBL_MAX_EXP + 1];
46
47static int
48init(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
70extern double
71frexp(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
138extern double
139ldexp(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