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 * frexpl/ldexpl implementation
26 */
27
28#include <ast.h>
29
30#include "FEATURE/float"
31
32#if _lib_frexpl && _lib_ldexpl
33
34NoN(frexpl)
35
36#else
37
38#ifndef LDBL_MAX_EXP
39#define LDBL_MAX_EXP	DBL_MAX_EXP
40#endif
41
42#if defined(_ast_fltmax_exp_index) && defined(_ast_fltmax_exp_shift)
43
44#define INIT()		_ast_fltmax_exp_t _pow_
45#define pow2(i)		(_pow_.f=1,_pow_.e[_ast_fltmax_exp_index]+=((i)<<_ast_fltmax_exp_shift),_pow_.f)
46
47#else
48
49static _ast_fltmax_t	pow2tab[LDBL_MAX_EXP + 1];
50
51static int
52init(void)
53{
54	register int		x;
55	_ast_fltmax_t		g;
56
57	g = 1;
58	for (x = 0; x < elementsof(pow2tab); x++)
59	{
60		pow2tab[x] = g;
61		g *= 2;
62	}
63	return 0;
64}
65
66#define INIT()		(pow2tab[0]?0:init())
67
68#define pow2(i)		(pow2tab[i])
69
70#endif
71
72#if !_lib_frexpl
73
74#undef	frexpl
75
76extern _ast_fltmax_t
77frexpl(_ast_fltmax_t f, int* p)
78{
79	register int		k;
80	register int		x;
81	_ast_fltmax_t		g;
82
83	INIT();
84
85	/*
86	 * normalize
87	 */
88
89	x = k = LDBL_MAX_EXP / 2;
90	if (f < 1)
91	{
92		g = 1.0L / f;
93		for (;;)
94		{
95			k = (k + 1) / 2;
96			if (g < pow2(x))
97				x -= k;
98			else if (k == 1 && g < pow2(x+1))
99				break;
100			else
101				x += k;
102		}
103		if (g == pow2(x))
104			x--;
105		x = -x;
106	}
107	else if (f > 1)
108	{
109		for (;;)
110		{
111			k = (k + 1) / 2;
112			if (f > pow2(x))
113				x += k;
114			else if (k == 1 && f > pow2(x-1))
115				break;
116			else
117				x -= k;
118		}
119		if (f == pow2(x))
120			x++;
121	}
122	else
123		x = 1;
124	*p = x;
125
126	/*
127	 * shift
128	 */
129
130	x = -x;
131	if (x < 0)
132		f /= pow2(-x);
133	else if (x < LDBL_MAX_EXP)
134		f *= pow2(x);
135	else
136		f = (f * pow2(LDBL_MAX_EXP - 1)) * pow2(x - (LDBL_MAX_EXP - 1));
137	return f;
138}
139
140#endif
141
142#if !_lib_ldexpl
143
144#undef	ldexpl
145
146extern _ast_fltmax_t
147ldexpl(_ast_fltmax_t f, register int x)
148{
149	INIT();
150	if (x < 0)
151		f /= pow2(-x);
152	else if (x < LDBL_MAX_EXP)
153		f *= pow2(x);
154	else
155		f = (f * pow2(LDBL_MAX_EXP - 1)) * pow2(x - (LDBL_MAX_EXP - 1));
156	return f;
157}
158
159#endif
160
161#endif
162