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