Deleted Added
sdiff udiff text old ( 9888:68d6b600d51f ) new ( 9899:0392ef94d766 )
full compact
1/*
2 * Copyright (c) 2013, Andreas Sandberg
3 * All rights reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions
7 * are met:
8 *
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.
15 *
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.
28 */
29
30#include <fputils/fp80.h>
31#include "fpbits.h"
32
33#include <assert.h>
34#include <stdint.h>
35
36#include <stdio.h>
37
38typedef union {
39 union {
40 uint64_t bits;
41 double value;
42 };
43} fp64_t;
44
45
46const fp80_t fp80_pinf = BUILD_FP80(0, 0, FP80_EXP_SPECIAL);
47const fp80_t fp80_ninf = BUILD_FP80(1, 0, FP80_EXP_SPECIAL);
48const fp80_t fp80_qnan = BUILD_FP80(0, FP80_FRAC_QNAN, FP80_EXP_SPECIAL);
49const fp80_t fp80_qnani = BUILD_FP80(1, FP80_FRAC_QNANI, FP80_EXP_SPECIAL);
50const fp80_t fp80_snan = BUILD_FP80(0, FP80_FRAC_SNAN, FP80_EXP_SPECIAL);
51const fp80_t fp80_nan = BUILD_FP80(0, FP80_FRAC_QNAN, FP80_EXP_SPECIAL);
52
53static const fp64_t fp64_pinf = BUILD_FP64(0, 0, FP64_EXP_SPECIAL);
54static const fp64_t fp64_ninf = BUILD_FP64(1, 0, FP64_EXP_SPECIAL);
55static const fp64_t fp64_qnan = BUILD_FP64(0, FP64_FRAC_QNAN,
56 FP64_EXP_SPECIAL);
57static const fp64_t fp64_nqnan = BUILD_FP64(1, FP64_FRAC_QNAN,
58 FP64_EXP_SPECIAL);
59static const fp64_t fp64_qnani = BUILD_FP64(1, FP64_FRAC_QNANI,
60 FP64_EXP_SPECIAL);
61static const fp64_t fp64_snan = BUILD_FP64(0, FP64_FRAC_SNAN,
62 FP64_EXP_SPECIAL);
63static const fp64_t fp64_nsnan = BUILD_FP64(1, FP64_FRAC_SNAN,
64 FP64_EXP_SPECIAL);
65
66static double
67build_fp64(int sign, uint64_t frac, int exp)
68{
69 const fp64_t f = BUILD_FP64(sign, frac, exp);
70
71 return f.value;
72}
73
74int
75fp80_sgn(fp80_t fp80)
76{
77 return (fp80.u.repr.se & FP80_SIGN_BIT) ? -1 : 1;
78}
79
80int
81fp80_isspecial(fp80_t fp80)
82{
83 const int exp = FP80_EXP(fp80);
84
85 return exp == FP80_EXP_SPECIAL;
86}
87
88int
89fp80_isinf(fp80_t fp80)
90{
91 const uint64_t frac = FP80_FRAC(fp80);
92
93 return fp80_isspecial(fp80) && frac == 0 ? fp80_sgn(fp80) : 0;
94}
95
96
97int
98fp80_isqnan(fp80_t fp80)
99{
100 const uint64_t frac = FP80_FRAC(fp80);
101
102 return fp80_isspecial(fp80) && (frac & FP80_QNAN_BIT);
103}
104
105int
106fp80_isqnani(fp80_t fp80)
107{
108 const uint64_t frac_low = fp80.u.repr.fi & (FP80_FRAC_MASK >> 1);
109
110 return fp80_isqnan(fp80) && (fp80.u.repr.se & FP80_SIGN_BIT) && !frac_low;
111}
112
113int
114fp80_issnan(fp80_t fp80)
115{
116 const uint64_t frac = FP80_FRAC(fp80);
117
118 return fp80_isspecial(fp80) && !(frac & FP80_QNAN_BIT) && frac;
119}
120
121int
122fp80_isfinite(fp80_t fp80)
123{
124 return !fp80_isnan(fp80) && !fp80_isinf(fp80);
125}
126
127int
128fp80_isnan(fp80_t fp80)
129{
130 return fp80_issnan(fp80) || fp80_isqnan(fp80) ? fp80_sgn(fp80) : 0;
131}
132
133int
134fp80_iszero(fp80_t fp80)
135{
136 return fp80.u.repr.fi == 0 && FP80_EXP(fp80) == 0 ? fp80_sgn(fp80) : 0;
137}
138
139int
140fp80_isnormal(fp80_t fp80)
141{
142 return FP80_EXP(fp80) != 0 && !fp80_isspecial(fp80) ?
143 fp80_sgn(fp80) : 0;
144}
145
146int
147fp80_issubnormal(fp80_t fp80)
148{
149 return FP80_FRAC(fp80) && FP80_EXP(fp80) == 0 ? fp80_sgn(fp80) : 0;
150}
151
152int
153fp80_classify(fp80_t fp80)
154{
155 if (fp80_issubnormal(fp80)) {
156 return FP_SUBNORMAL;
157 } else if (fp80_iszero(fp80)) {
158 return FP_ZERO;
159 } else if (fp80_isinf(fp80)) {
160 return FP_INFINITE;
161 } else if (fp80_isnan(fp80)) {
162 return FP_NAN;
163 } else {
164 assert(fp80_isfinite(fp80));
165 return FP_NORMAL;
166 }
167}
168
169double
170fp80_cvtd(fp80_t fp80)
171{
172 const int sign = fp80.u.repr.se & FP80_SIGN_BIT;
173
174 if (!fp80_isspecial(fp80)) {
175 const uint64_t frac = fp80.u.repr.fi;
176 const int unb_exp = FP80_EXP(fp80) - FP80_EXP_BIAS;
177 const int fp64_exp = unb_exp + FP64_EXP_BIAS;
178 const uint64_t fp64_frac = frac >> (FP80_FRAC_BITS - FP64_FRAC_BITS);
179
180 if (fp64_exp > 0 && fp64_exp < FP64_EXP_SPECIAL) {
181 /* These numbers fall in the range of what we can express
182 * as normals */
183 return build_fp64(sign, fp64_frac, fp64_exp);
184 } else if (fp64_exp <= 0) {
185 uint64_t fp64_denormal_frac = fp64_frac >> (-fp64_exp);
186 /* Generate a denormal or zero */
187 return build_fp64(sign, fp64_denormal_frac, 0);
188 } else {
189 /* Infinity */
190 return build_fp64(sign, 0, FP64_EXP_SPECIAL);
191 }
192 } else {
193 if (fp80_isinf(fp80)) {
194 return build_fp64(sign, 0, FP64_EXP_SPECIAL);
195 } else if (fp80_issnan(fp80)) {
196 return fp80_sgn(fp80) > 0 ? fp64_snan.value : fp64_nsnan.value;
197 } else if (fp80_isqnani(fp80)) {
198 return fp64_qnani.value;
199 } else {
200 assert(fp80_isqnan(fp80));
201 return fp80_sgn(fp80) > 0 ? fp64_qnan.value : fp64_nqnan.value;
202 }
203 }
204}
205
206fp80_t
207fp80_cvfd(double value)
208{
209 const fp64_t fp64 = { .value = value };
210 const uint64_t frac = FP64_FRAC(fp64);
211 const unsigned exp = FP64_EXP(fp64);
212 const int unb_exp = exp - FP64_EXP_BIAS;
213 const uint64_t fp80_frac = frac << (FP80_FRAC_BITS - FP64_FRAC_BITS);
214
215 if (exp != 0) {
216 // Normal, inf, nan
217 const unsigned fp80_exp = exp == FP64_EXP_SPECIAL ?
218 FP80_EXP_SPECIAL : (unb_exp + FP80_EXP_BIAS);
219 const fp80_t fp80 = BUILD_FP80(fp64.bits & FP64_SIGN_BIT,
220 fp80_frac, fp80_exp);
221 return fp80;
222 } else if (exp == 0 && frac == 0) {
223 // Zero
224 const fp80_t fp80 = BUILD_FP80(fp64.bits & FP64_SIGN_BIT, 0, 0);
225 return fp80;
226 } else {
227 // Denormal
228 uint64_t fp80_fi = fp80_frac;
229 int shift_amt = 0;
230 while (!(fp80_fi & FP80_INT_BIT)) {
231 fp80_fi <<= 1;
232 ++shift_amt;
233 }
234 const unsigned fp80_exp = (unb_exp - shift_amt) + FP80_EXP_BIAS;
235 const fp80_t fp80 = BUILD_FP80(fp64.bits & FP64_SIGN_BIT,
236 fp80_fi, fp80_exp);
237 return fp80;
238 }
239}
240
241void
242fp80_debug_dump(FILE *fout, fp80_t fp80)
243{
244 fprintf(fout, "sgn: %i, int: %i, frac: 0x%llx, exp: 0x%x (%i)\n",
245 fp80_sgn(fp80), !!(fp80.u.repr.fi & FP80_INT_BIT), FP80_FRAC(fp80),
246 FP80_EXP(fp80), FP80_EXP(fp80) - FP80_EXP_BIAS);
247}