xref: /plan9-contrib/sys/src/ape/lib/ap/plan9/frexp.c (revision eaba85aa6b158bdf68fdb77f770e3ba0899a8b5e)
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 	if(d < 1) {
77 		if(d < 0) {
78 			f = modf(-d, ip);
79 			*ip = -*ip;
80 			return -f;
81 		}
82 		*ip = 0;
83 		return d;
84 	}
85 	x.d = d;
86 	e = ((x.ms >> SHIFT) & MASK) - BIAS;
87 	if(e <= SHIFT+1) {
88 		x.ms &= ~(0x1fffffL >> e);
89 		x.ls = 0;
90 	} else
91 	if(e <= SHIFT+33)
92 		x.ls &= ~(0x7fffffffL >> (e-SHIFT-2));
93 	*ip = x.d;
94 	return d - x.d;
95 }
96