1 #include <u.h>
2 #include <libc.h>
3
4 /*
5 * this is big/little endian non-portable
6 * it gets the endian from the FPdbleword
7 * union in u.h.
8 */
9 #define MASK 0x7ffL
10 #define SHIFT 20
11 #define BIAS 1022L
12 #define SIG 52
13
14 double
frexp(double d,int * ep)15 frexp(double d, int *ep)
16 {
17 FPdbleword x, a;
18
19 *ep = 0;
20 /* order matters: only isNaN can operate on NaN */
21 if(isNaN(d) || isInf(d, 0) || d == 0)
22 return d;
23 x.x = d;
24 a.x = fabs(d);
25 if((a.hi >> SHIFT) == 0){ /* normalize subnormal numbers */
26 x.x = (double)(1ULL<<SIG) * d;
27 *ep = -SIG;
28 }
29 *ep += ((x.hi >> SHIFT) & MASK) - BIAS;
30 x.hi &= ~(MASK << SHIFT);
31 x.hi |= BIAS << SHIFT;
32 return x.x;
33 }
34
35 double
ldexp(double d,int deltae)36 ldexp(double d, int deltae)
37 {
38 int e, bits;
39 FPdbleword x;
40 ulong z;
41
42 if(d == 0)
43 return 0;
44 x.x = d;
45 e = (x.hi >> SHIFT) & MASK;
46 if(deltae >= 0 || e+deltae >= 1){ /* no underflow */
47 e += deltae;
48 if(e >= MASK){ /* overflow */
49 if(d < 0)
50 return Inf(-1);
51 return Inf(1);
52 }
53 }else{ /* underflow gracefully */
54 deltae = -deltae;
55 /* need to shift d right deltae */
56 if(e > 1){ /* shift e-1 by exponent manipulation */
57 deltae -= e-1;
58 e = 1;
59 }
60 if(deltae > 0 && e==1){ /* shift 1 by switch from 1.fff to 0.1ff */
61 deltae--;
62 e = 0;
63 x.lo >>= 1;
64 x.lo |= (x.hi&1)<<31;
65 z = x.hi & ((1<<SHIFT)-1);
66 x.hi &= ~((1<<SHIFT)-1);
67 x.hi |= (1<<(SHIFT-1)) | (z>>1);
68 }
69 while(deltae > 0){ /* shift bits down */
70 bits = deltae;
71 if(bits > SHIFT)
72 bits = SHIFT;
73 x.lo >>= bits;
74 x.lo |= (x.hi&((1<<bits)-1)) << (32-bits);
75 z = x.hi & ((1<<SHIFT)-1);
76 x.hi &= ~((1<<SHIFT)-1);
77 x.hi |= z>>bits;
78 deltae -= bits;
79 }
80 }
81 x.hi &= ~(MASK << SHIFT);
82 x.hi |= (long)e << SHIFT;
83 return x.x;
84 }
85
86 double
modf(double d,double * ip)87 modf(double d, double *ip)
88 {
89 FPdbleword x;
90 int e;
91
92 x.x = d;
93 e = (x.hi >> SHIFT) & MASK;
94 if(e == MASK){
95 *ip = d;
96 if(x.lo != 0 || (x.hi & 0xfffffL) != 0) /* NaN */
97 return d;
98 /* ±Inf */
99 x.hi &= 0x80000000L;
100 return x.x;
101 }
102 if(d < 1) {
103 if(d < 0) {
104 x.x = modf(-d, ip);
105 *ip = -*ip;
106 return -x.x;
107 }
108 *ip = 0;
109 return d;
110 }
111 e -= BIAS;
112 if(e <= SHIFT+1) {
113 x.hi &= ~(0x1fffffL >> e);
114 x.lo = 0;
115 } else
116 if(e <= SHIFT+33)
117 x.lo &= ~(0x7fffffffL >> (e-SHIFT-2));
118 *ip = x.x;
119 return d - x.x;
120 }
121