7c52865d1b
This changeset updates the external library to git revision 52b6190b4e. This update includes changes that fix compilation errors on old gcc versions and fixes to test a case that affect ICC.
246 lines
7.2 KiB
C
246 lines
7.2 KiB
C
/*
|
|
* 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 <fputils/fp80.h>
|
|
#include "fpbits.h"
|
|
|
|
#include <assert.h>
|
|
#include <stdint.h>
|
|
|
|
#include <stdio.h>
|
|
|
|
typedef 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.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.repr.fi & (FP80_FRAC_MASK >> 1);
|
|
|
|
return fp80_isqnan(fp80) && (fp80.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.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.repr.se & FP80_SIGN_BIT;
|
|
|
|
if (!fp80_isspecial(fp80)) {
|
|
const uint64_t frac = fp80.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.repr.fi & FP80_INT_BIT), FP80_FRAC(fp80),
|
|
FP80_EXP(fp80), FP80_EXP(fp80) - FP80_EXP_BIAS);
|
|
}
|