2 * Copyright (c) 2013, Andreas Sandberg
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions
9 * 1. Redistributions of source code must retain the above copyright
10 * notice, this list of conditions and the following disclaimer.
11 * 2. Redistributions in binary form must reproduce the above
12 * copyright notice, this list of conditions and the following
13 * disclaimer in the documentation and/or other materials provided
14 * with the distribution.
16 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
17 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
18 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
19 * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
20 * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
21 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
22 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
23 * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
24 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
25 * STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
26 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
27 * OF THE POSSIBILITY OF SUCH DAMAGE.
30 #include <fputils/fp80.h>
31 #include <fputils/fp64.h>
40 const fp80_t fp80_pinf
= BUILD_FP80(0, 0, FP80_EXP_SPECIAL
);
41 const fp80_t fp80_ninf
= BUILD_FP80(1, 0, FP80_EXP_SPECIAL
);
42 const fp80_t fp80_qnan
= BUILD_FP80(0, FP80_FRAC_QNAN
, FP80_EXP_SPECIAL
);
43 const fp80_t fp80_qnani
= BUILD_FP80(1, FP80_FRAC_QNANI
, FP80_EXP_SPECIAL
);
44 const fp80_t fp80_snan
= BUILD_FP80(0, FP80_FRAC_SNAN
, FP80_EXP_SPECIAL
);
45 const fp80_t fp80_nan
= BUILD_FP80(0, FP80_FRAC_QNAN
, FP80_EXP_SPECIAL
);
50 return (fp80
.repr
.se
& FP80_SIGN_BIT
) ? -1 : 1;
54 fp80_isspecial(fp80_t fp80
)
56 const int exp
= FP80_EXP(fp80
);
58 return exp
== FP80_EXP_SPECIAL
;
62 fp80_isinf(fp80_t fp80
)
64 const uint64_t frac
= FP80_FRAC(fp80
);
66 return fp80_isspecial(fp80
) && frac
== 0 ? fp80_sgn(fp80
) : 0;
71 fp80_isqnan(fp80_t fp80
)
73 const uint64_t frac
= FP80_FRAC(fp80
);
75 return fp80_isspecial(fp80
) && (frac
& FP80_QNAN_BIT
);
79 fp80_isqnani(fp80_t fp80
)
81 const uint64_t frac_low
= fp80
.repr
.fi
& (FP80_FRAC_MASK
>> 1);
83 return fp80_isqnan(fp80
) && (fp80
.repr
.se
& FP80_SIGN_BIT
) && !frac_low
;
87 fp80_issnan(fp80_t fp80
)
89 const uint64_t frac
= FP80_FRAC(fp80
);
91 return fp80_isspecial(fp80
) && !(frac
& FP80_QNAN_BIT
) && frac
;
95 fp80_isfinite(fp80_t fp80
)
97 return !fp80_isnan(fp80
) && !fp80_isinf(fp80
);
101 fp80_isnan(fp80_t fp80
)
103 return fp80_issnan(fp80
) || fp80_isqnan(fp80
) ? fp80_sgn(fp80
) : 0;
107 fp80_iszero(fp80_t fp80
)
109 return fp80
.repr
.fi
== 0 && FP80_EXP(fp80
) == 0 ? fp80_sgn(fp80
) : 0;
113 fp80_isnormal(fp80_t fp80
)
115 return FP80_EXP(fp80
) != 0 && !fp80_isspecial(fp80
) ?
120 fp80_issubnormal(fp80_t fp80
)
122 return FP80_FRAC(fp80
) && FP80_EXP(fp80
) == 0 ? fp80_sgn(fp80
) : 0;
126 fp80_classify(fp80_t fp80
)
128 if (fp80_issubnormal(fp80
)) {
130 } else if (fp80_iszero(fp80
)) {
132 } else if (fp80_isinf(fp80
)) {
134 } else if (fp80_isnan(fp80
)) {
137 assert(fp80_isfinite(fp80
));
143 fp80_cvtd(fp80_t fp80
)
145 return fp80_cvtfp64(fp80
).value
;
149 fp80_cvtfp64(fp80_t fp80
)
151 const int sign
= fp80
.repr
.se
& FP80_SIGN_BIT
;
153 if (!fp80_isspecial(fp80
)) {
154 const uint64_t frac
= fp80
.repr
.fi
;
155 const int unb_exp
= FP80_EXP(fp80
) - FP80_EXP_BIAS
;
156 const int fp64_exp
= unb_exp
+ FP64_EXP_BIAS
;
157 const uint64_t fp64_frac
= frac
>> (FP80_FRAC_BITS
- FP64_FRAC_BITS
);
159 if (fp64_exp
> 0 && fp64_exp
< FP64_EXP_SPECIAL
) {
160 /* These numbers fall in the range of what we can express
162 return build_fp64(sign
, fp64_frac
, fp64_exp
);
163 } else if (fp64_exp
<= 0) {
164 uint64_t fp64_denormal_frac
= fp64_frac
>> (-fp64_exp
);
165 /* Generate a denormal or zero */
166 return build_fp64(sign
, fp64_denormal_frac
, 0);
169 return build_fp64(sign
, 0, FP64_EXP_SPECIAL
);
172 if (fp80_isinf(fp80
)) {
173 return build_fp64(sign
, 0, FP64_EXP_SPECIAL
);
174 } else if (fp80_issnan(fp80
)) {
175 return fp80_sgn(fp80
) > 0 ? fp64_snan
: fp64_nsnan
;
176 } else if (fp80_isqnani(fp80
)) {
179 assert(fp80_isqnan(fp80
));
180 return fp80_sgn(fp80
) > 0 ? fp64_qnan
: fp64_nqnan
;
186 fp80_cvfd(double value
)
188 const fp64_t fp64
= { .value
= value
};
190 return fp80_cvffp64(fp64
);
194 fp80_cvffp64(fp64_t fp64
)
196 const uint64_t frac
= FP64_FRAC(fp64
);
197 const unsigned exp
= FP64_EXP(fp64
);
198 const int unb_exp
= exp
- FP64_EXP_BIAS
;
199 const uint64_t fp80_frac
= frac
<< (FP80_FRAC_BITS
- FP64_FRAC_BITS
);
203 const unsigned fp80_exp
= exp
== FP64_EXP_SPECIAL
?
204 FP80_EXP_SPECIAL
: (unb_exp
+ FP80_EXP_BIAS
);
205 const fp80_t fp80
= BUILD_FP80(fp64
.bits
& FP64_SIGN_BIT
,
206 fp80_frac
, fp80_exp
);
208 } else if (exp
== 0 && frac
== 0) {
210 const fp80_t fp80
= BUILD_FP80(fp64
.bits
& FP64_SIGN_BIT
, 0, 0);
214 uint64_t fp80_fi
= fp80_frac
;
216 while (!(fp80_fi
& FP80_INT_BIT
)) {
220 const unsigned fp80_exp
= (unb_exp
- shift_amt
) + FP80_EXP_BIAS
;
221 const fp80_t fp80
= BUILD_FP80(fp64
.bits
& FP64_SIGN_BIT
,
228 fp80_debug_dump(FILE *fout
, fp80_t fp80
)
230 fprintf(fout
, "sgn: %i, int: %i, frac: 0x%llx, exp: 0x%x (%i)\n",
231 fp80_sgn(fp80
), !!(fp80
.repr
.fi
& FP80_INT_BIT
), FP80_FRAC(fp80
),
232 FP80_EXP(fp80
), FP80_EXP(fp80
) - FP80_EXP_BIAS
);