6 #define ASUINT64(x) ((union {double f; uint64_t i;}){x}).i
7 #define ZEROINFNAN (0x7ff-0x3ff-52-1)
9 struct num
{ uint64_t m
; int e
; int sign
; };
11 static struct num
normalize(double x
)
13 uint64_t ix
= ASUINT64(x
);
18 ix
= ASUINT64(x
*0x1p
63);
26 return (struct num
){ix
,e
,sign
};
29 static void mul(uint64_t *hi
, uint64_t *lo
, uint64_t x
, uint64_t y
)
32 uint64_t xlo
= (uint32_t)x
, xhi
= x
>>32;
33 uint64_t ylo
= (uint32_t)y
, yhi
= y
>>32;
36 t2
= xlo
*yhi
+ xhi
*ylo
;
39 *hi
= t3
+ (t2
>>32) + (t1
> *lo
);
42 double fma(double x
, double y
, double z
)
44 #pragma STDC FENV_ACCESS ON
46 /* normalize so top 10bits and last bit are 0 */
47 struct num nx
, ny
, nz
;
52 if (nx
.e
>= ZEROINFNAN
|| ny
.e
>= ZEROINFNAN
)
54 if (nz
.e
>= ZEROINFNAN
) {
55 if (nz
.e
> ZEROINFNAN
) /* z==0 */
61 uint64_t rhi
, rlo
, zhi
, zlo
;
62 mul(&rhi
, &rlo
, nx
.m
, ny
.m
);
63 /* either top 20 or 21 bits of rhi and last 2 bits of rlo are 0 */
68 /* shift bits z<<=kz, r>>=kr, so kz+kr == d, set e = e+kr (== ez-kz) */
80 rlo
= rhi
<<64-d
| rlo
>>d
| !!(rlo
<<64-d
);
93 zlo
= nz
.m
>>d
| !!(nz
.m
<<64-d
);
100 int sign
= nx
.sign
^ny
.sign
;
101 int samesign
= !(sign
^nz
.sign
);
106 rhi
+= zhi
+ (rlo
< zlo
);
111 rhi
= rhi
- zhi
- (t
< rlo
);
120 /* set rhi to top 63bit of the result (last bit is sticky) */
125 rhi
= rhi
<<d
| rlo
>>64-d
| !!(rlo
<<d
);
129 rhi
= rlo
>>1 | (rlo
&1);
138 /* convert to double */
139 int64_t i
= rhi
; /* i is in [1<<62,(1<<63)-1] */
142 double r
= i
; /* |r| is in [0x1p62,0x1p63] */
145 /* result is subnormal before rounding */
151 /* min normal after rounding, underflow depends
152 on arch behaviour which can be imitated by
153 a double to float conversion */
154 float fltmin
= 0x0.ffffff8p
-63*FLT_MIN
* r
;
155 return DBL_MIN
/FLT_MIN
* fltmin
;
157 /* one bit is lost when scaled, add another top bit to
158 only round once at conversion if it is inexact */
160 i
= rhi
>>1 | (rhi
&1) | 1ull<<62;
164 r
= 2*r
- c
; /* remove top bit */
166 /* raise underflow portably, such that it
167 cannot be optimized away */
169 double_t tiny
= DBL_MIN
/FLT_MIN
* r
;
170 r
+= (double)(tiny
*tiny
) * (r
-r
);
174 /* only round once when scaled */
176 i
= ( rhi
>>d
| !!(rhi
<<64-d
) ) << d
;