Subversion Repositories planix.SVN

Rev

Rev 2 | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
2 - 1 
#include <math.h>
2 
#include <errno.h>
3 
#define _RESEARCH_SOURCE
4 
#include <float.h>
5 
 
6 
#define	MASK	0x7ffL
7 
#define	SHIFT	20
8 
#define	BIAS	1022L
9 
#define	SIG	52
10 
 
11 
typedef	union
12 
{
13 
	double	d;
14 
	struct
15 
	{
16 
#ifdef IEEE_8087
17 
		long	ls;
18 
		long	ms;
19 
#else
20 
		long	ms;
21 
		long	ls;
22 
#endif
23 
	};
24 
} Cheat;
25 
 
26 
double
27 
frexp(double d, int *ep)
28 
{
29 
	Cheat x, a;
30 
 
31 
	*ep = 0;
32 
	/* order matters: only isNaN can operate on NaN */
33 
	if(isNaN(d) || isInf(d, 0) || d == 0)
34 
		return d;
35 
	x.d = d;
36 
	a.d = fabs(d);
37 
	if((a.ms >> SHIFT) == 0){	/* normalize subnormal numbers */
38 
		x.d = (double)(1ULL<<SIG) * d;
39 
		*ep = -SIG;
40 
	}
41 
	*ep = ((x.ms >> SHIFT) & MASK) - BIAS;
42 
	x.ms &= ~(MASK << SHIFT);
43 
	x.ms |= BIAS << SHIFT;
44 
	return x.d;
45 
}
46 
 
47 
double
48 
ldexp(double d, int e)
49 
{
50 
	Cheat x;
51 
 
52 
	if(d == 0)
53 
		return 0;
54 
	x.d = d;
55 
	e += (x.ms >> SHIFT) & MASK;
56 
	if(e <= 0)
57 
		return 0;
58 
	if(e >= MASK){
59 
		errno = ERANGE;
60 
		if(d < 0)
61 
			return -HUGE_VAL;
62 
		return HUGE_VAL;
63 
	}
64 
	x.ms &= ~(MASK << SHIFT);
65 
	x.ms |= (long)e << SHIFT;
66 
	return x.d;
67 
}
68 
 
69 
double
70 
modf(double d, double *ip)
71 
{
72 
	double f;
73 
	Cheat x;
74 
	int e;
75 
 
76 
	x.d = d;
77 
	e = (x.ms >> SHIFT) & MASK;
78 
	if(e == MASK){
79 
		*ip = d;
80 
		if(x.ls != 0 || (x.ms & 0xfffffL) != 0)	/* NaN */
81 
			return d;
82 
		/* ±Inf */
83 
		x.ms &= 0x80000000L;
84 
		return x.d;
85 
	}
86 
	if(d < 1) {
87 
		if(d < 0) {
88 
			f = modf(-d, ip);
89 
			*ip = -*ip;
90 
			return -f;
91 
		}
92 
		*ip = 0;
93 
		return d;
94 
	}
95 
	e -= BIAS;
96 
	if(e <= SHIFT+1) {
97 
		x.ms &= ~(0x1fffffL >> e);
98 
		x.ls = 0;
99 
	} else
100 
	if(e <= SHIFT+33)
101 
		x.ls &= ~(0x7fffffffL >> (e-SHIFT-2));
102 
	*ip = x.d;
103 
	return d - x.d;
104 
}