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