From d3937f3b373ee9b1af09176d8a7d86eb6fcf790f Mon Sep 17 00:00:00 2001 From: Andreas Sandberg Date: Mon, 30 Sep 2013 09:40:26 +0200 Subject: ext: Include libfputils This changeset includes libfputils from revision bbf0d61d75. This library can be used to convert to and from 80-bit floats and query the type of an 80-bit float, which is needed to support the x87 FPU. --- ext/fputils/fp80.c | 247 +++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 247 insertions(+) create mode 100644 ext/fputils/fp80.c (limited to 'ext/fputils/fp80.c') diff --git a/ext/fputils/fp80.c b/ext/fputils/fp80.c new file mode 100644 index 000000000..703b07ac9 --- /dev/null +++ b/ext/fputils/fp80.c @@ -0,0 +1,247 @@ +/* + * Copyright (c) 2013, Andreas Sandberg + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * + * 1. Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above + * copyright notice, this list of conditions and the following + * disclaimer in the documentation and/or other materials provided + * with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT + * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS + * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE + * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, + * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES + * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR + * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) + * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, + * STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED + * OF THE POSSIBILITY OF SUCH DAMAGE. + */ + +#include +#include "fpbits.h" + +#include +#include + +#include + +typedef union { + union { + uint64_t bits; + double value; + }; +} fp64_t; + + +const fp80_t fp80_pinf = BUILD_FP80(0, 0, FP80_EXP_SPECIAL); +const fp80_t fp80_ninf = BUILD_FP80(1, 0, FP80_EXP_SPECIAL); +const fp80_t fp80_qnan = BUILD_FP80(0, FP80_FRAC_QNAN, FP80_EXP_SPECIAL); +const fp80_t fp80_qnani = BUILD_FP80(1, FP80_FRAC_QNANI, FP80_EXP_SPECIAL); +const fp80_t fp80_snan = BUILD_FP80(0, FP80_FRAC_SNAN, FP80_EXP_SPECIAL); +const fp80_t fp80_nan = BUILD_FP80(0, FP80_FRAC_QNAN, FP80_EXP_SPECIAL); + +static const fp64_t fp64_pinf = BUILD_FP64(0, 0, FP64_EXP_SPECIAL); +static const fp64_t fp64_ninf = BUILD_FP64(1, 0, FP64_EXP_SPECIAL); +static const fp64_t fp64_qnan = BUILD_FP64(0, FP64_FRAC_QNAN, + FP64_EXP_SPECIAL); +static const fp64_t fp64_nqnan = BUILD_FP64(1, FP64_FRAC_QNAN, + FP64_EXP_SPECIAL); +static const fp64_t fp64_qnani = BUILD_FP64(1, FP64_FRAC_QNANI, + FP64_EXP_SPECIAL); +static const fp64_t fp64_snan = BUILD_FP64(0, FP64_FRAC_SNAN, + FP64_EXP_SPECIAL); +static const fp64_t fp64_nsnan = BUILD_FP64(1, FP64_FRAC_SNAN, + FP64_EXP_SPECIAL); + +static double +build_fp64(int sign, uint64_t frac, int exp) +{ + const fp64_t f = BUILD_FP64(sign, frac, exp); + + return f.value; +} + +int +fp80_sgn(fp80_t fp80) +{ + return (fp80.u.repr.se & FP80_SIGN_BIT) ? -1 : 1; +} + +int +fp80_isspecial(fp80_t fp80) +{ + const int exp = FP80_EXP(fp80); + + return exp == FP80_EXP_SPECIAL; +} + +int +fp80_isinf(fp80_t fp80) +{ + const uint64_t frac = FP80_FRAC(fp80); + + return fp80_isspecial(fp80) && frac == 0 ? fp80_sgn(fp80) : 0; +} + + +int +fp80_isqnan(fp80_t fp80) +{ + const uint64_t frac = FP80_FRAC(fp80); + + return fp80_isspecial(fp80) && (frac & FP80_QNAN_BIT); +} + +int +fp80_isqnani(fp80_t fp80) +{ + const uint64_t frac_low = fp80.u.repr.fi & (FP80_FRAC_MASK >> 1); + + return fp80_isqnan(fp80) && (fp80.u.repr.se & FP80_SIGN_BIT) && !frac_low; +} + +int +fp80_issnan(fp80_t fp80) +{ + const uint64_t frac = FP80_FRAC(fp80); + + return fp80_isspecial(fp80) && !(frac & FP80_QNAN_BIT) && frac; +} + +int +fp80_isfinite(fp80_t fp80) +{ + return !fp80_isnan(fp80) && !fp80_isinf(fp80); +} + +int +fp80_isnan(fp80_t fp80) +{ + return fp80_issnan(fp80) || fp80_isqnan(fp80) ? fp80_sgn(fp80) : 0; +} + +int +fp80_iszero(fp80_t fp80) +{ + return fp80.u.repr.fi == 0 && FP80_EXP(fp80) == 0 ? fp80_sgn(fp80) : 0; +} + +int +fp80_isnormal(fp80_t fp80) +{ + return FP80_EXP(fp80) != 0 && !fp80_isspecial(fp80) ? + fp80_sgn(fp80) : 0; +} + +int +fp80_issubnormal(fp80_t fp80) +{ + return FP80_FRAC(fp80) && FP80_EXP(fp80) == 0 ? fp80_sgn(fp80) : 0; +} + +int +fp80_classify(fp80_t fp80) +{ + if (fp80_issubnormal(fp80)) { + return FP_SUBNORMAL; + } else if (fp80_iszero(fp80)) { + return FP_ZERO; + } else if (fp80_isinf(fp80)) { + return FP_INFINITE; + } else if (fp80_isnan(fp80)) { + return FP_NAN; + } else { + assert(fp80_isfinite(fp80)); + return FP_NORMAL; + } +} + +double +fp80_cvtd(fp80_t fp80) +{ + const int sign = fp80.u.repr.se & FP80_SIGN_BIT; + + if (!fp80_isspecial(fp80)) { + const uint64_t frac = fp80.u.repr.fi; + const int unb_exp = FP80_EXP(fp80) - FP80_EXP_BIAS; + const int fp64_exp = unb_exp + FP64_EXP_BIAS; + const uint64_t fp64_frac = frac >> (FP80_FRAC_BITS - FP64_FRAC_BITS); + + if (fp64_exp > 0 && fp64_exp < FP64_EXP_SPECIAL) { + /* These numbers fall in the range of what we can express + * as normals */ + return build_fp64(sign, fp64_frac, fp64_exp); + } else if (fp64_exp <= 0) { + uint64_t fp64_denormal_frac = fp64_frac >> (-fp64_exp); + /* Generate a denormal or zero */ + return build_fp64(sign, fp64_denormal_frac, 0); + } else { + /* Infinity */ + return build_fp64(sign, 0, FP64_EXP_SPECIAL); + } + } else { + if (fp80_isinf(fp80)) { + return build_fp64(sign, 0, FP64_EXP_SPECIAL); + } else if (fp80_issnan(fp80)) { + return fp80_sgn(fp80) > 0 ? fp64_snan.value : fp64_nsnan.value; + } else if (fp80_isqnani(fp80)) { + return fp64_qnani.value; + } else { + assert(fp80_isqnan(fp80)); + return fp80_sgn(fp80) > 0 ? fp64_qnan.value : fp64_nqnan.value; + } + } +} + +fp80_t +fp80_cvfd(double value) +{ + const fp64_t fp64 = { .value = value }; + const uint64_t frac = FP64_FRAC(fp64); + const unsigned exp = FP64_EXP(fp64); + const int unb_exp = exp - FP64_EXP_BIAS; + const uint64_t fp80_frac = frac << (FP80_FRAC_BITS - FP64_FRAC_BITS); + + if (exp != 0) { + // Normal, inf, nan + const unsigned fp80_exp = exp == FP64_EXP_SPECIAL ? + FP80_EXP_SPECIAL : (unb_exp + FP80_EXP_BIAS); + const fp80_t fp80 = BUILD_FP80(fp64.bits & FP64_SIGN_BIT, + fp80_frac, fp80_exp); + return fp80; + } else if (exp == 0 && frac == 0) { + // Zero + const fp80_t fp80 = BUILD_FP80(fp64.bits & FP64_SIGN_BIT, 0, 0); + return fp80; + } else { + // Denormal + uint64_t fp80_fi = fp80_frac; + int shift_amt = 0; + while (!(fp80_fi & FP80_INT_BIT)) { + fp80_fi <<= 1; + ++shift_amt; + } + const unsigned fp80_exp = (unb_exp - shift_amt) + FP80_EXP_BIAS; + const fp80_t fp80 = BUILD_FP80(fp64.bits & FP64_SIGN_BIT, + fp80_fi, fp80_exp); + return fp80; + } +} + +void +fp80_debug_dump(FILE *fout, fp80_t fp80) +{ + fprintf(fout, "sgn: %i, int: %i, frac: 0x%llx, exp: 0x%x (%i)\n", + fp80_sgn(fp80), !!(fp80.u.repr.fi & FP80_INT_BIT), FP80_FRAC(fp80), + FP80_EXP(fp80), FP80_EXP(fp80) - FP80_EXP_BIAS); +} -- cgit v1.2.3