blob: 7d119bb25a872870d25da5757d705f78550c79ca (
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
|
#include <u.h>
#include <base.h>
/*
* this is big/little endian non-portable
* it gets the endian from the FPdbleword
* union in u.h.
*/
#define MASK 0x7ffL
#define SHIFT 20
#define BIAS 1022L
#define SIG 52
double
math·frexp(double d, int *ep)
{
sys·DoubleWord x, a;
*ep = 0;
/* order matters: only isNaN can operate on NaN */
if(math·isNaN(d) || math·isInf(d, 0) || d == 0)
return d;
x.x = d;
a.x = math·fabs(d);
if((a.hi >> SHIFT) == 0){ /* normalize subnormal numbers */
x.x = (double)(1ULL<<SIG) * d;
*ep = -SIG;
}
*ep += ((x.hi >> SHIFT) & MASK) - BIAS;
x.hi &= ~(MASK << SHIFT);
x.hi |= BIAS << SHIFT;
return x.x;
}
double
math·ldexp(double d, int deltae)
{
int e, bits;
sys·DoubleWord x;
ulong z;
if(d == 0)
return 0;
x.x = d;
e = (x.hi >> SHIFT) & MASK;
if(deltae >= 0 || e+deltae >= 1){ /* no underflow */
e += deltae;
if(e >= MASK){ /* overflow */
if(d < 0)
return math·Inf(-1);
return math·Inf(1);
}
}else{ /* underflow gracefully */
deltae = -deltae;
/* need to shift d right deltae */
if(e > 1){ /* shift e-1 by exponent manipulation */
deltae -= e-1;
e = 1;
}
if(deltae > 0 && e==1){ /* shift 1 by switch from 1.fff to 0.1ff */
deltae--;
e = 0;
x.lo >>= 1;
x.lo |= (x.hi&1)<<31;
z = x.hi & ((1<<SHIFT)-1);
x.hi &= ~((1<<SHIFT)-1);
x.hi |= (1<<(SHIFT-1)) | (z>>1);
}
while(deltae > 0){ /* shift bits down */
bits = deltae;
if(bits > SHIFT)
bits = SHIFT;
x.lo >>= bits;
x.lo |= (x.hi&((1<<bits)-1)) << (32-bits);
z = x.hi & ((1<<SHIFT)-1);
x.hi &= ~((1<<SHIFT)-1);
x.hi |= z>>bits;
deltae -= bits;
}
}
x.hi &= ~(MASK << SHIFT);
x.hi |= (long)e << SHIFT;
return x.x;
}
double
math·modf(double d, double *ip)
{
sys·DoubleWord x;
int e;
x.x = d;
e = (x.hi >> SHIFT) & MASK;
if(e == MASK){
*ip = d;
if(x.lo != 0 || (x.hi & 0xfffffL) != 0) /* NaN */
return d;
/* ±Inf */
x.hi &= 0x80000000L;
return x.x;
}
if(d < 1) {
if(d < 0) {
x.x = math·modf(-d, ip);
*ip = -*ip;
return -x.x;
}
*ip = 0;
return d;
}
e -= BIAS;
if(e <= SHIFT+1) {
x.hi &= ~(0x1fffffL >> e);
x.lo = 0;
} else
if(e <= SHIFT+33)
x.lo &= ~(0x7fffffffL >> (e-SHIFT-2));
*ip = x.x;
return d - x.x;
}
|
