femtolisp/llt/fp.c

92 lines
1.9 KiB
C
Raw Normal View History

#include <math.h>
#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include "ieee754.h"
#include "dtypes.h"
#include "hashing.h"
static uint64_t max_ulps;
static uint32_t flt_max_ulps;
static uint64_t nexti64pow2(uint64_t i)
{
if (i==0) return 1;
if ((i&(i-1))==0) return i;
if (i&BIT63) return BIT63;
// repeatedly clear bottom bit
while (i&(i-1))
i = i&(i-1);
return i<<1;
}
static uint32_t nexti32pow2(uint32_t i)
{
if (i==0) return 1;
if ((i&(i-1))==0) return i;
if (i&BIT31) return BIT31;
// repeatedly clear bottom bit
while (i&(i-1))
i = i&(i-1);
return i<<1;
}
void dbl_tolerance(double tol)
{
max_ulps = nexti64pow2((uint64_t)(tol/DBL_EPSILON));
}
void flt_tolerance(float tol)
{
flt_max_ulps = nexti32pow2((uint32_t)(tol/FLT_EPSILON));
}
#ifdef __INTEL_COMPILER
static inline int64_t llabs(int64_t j)
{
return NBABS(j, 64);
}
#else
#endif
int dbl_equals(double a, double b)
{
union { double d; int64_t i; } ua;
union { double d; int64_t i; } ub;
int64_t aint, bint;
if (a == b)
return 1;
ua.d = a; aint = ua.i;
ub.d = b; bint = ub.i;
if (aint < 0)
aint = BIT63 - aint;
if (bint < 0)
bint = BIT63 - bint;
/* you'd think it makes no difference whether the result of llabs is
signed or unsigned, but if it's signed then the case of
0x8000000000000000 blows up, making 4 == -1 :) */
if ((uint64_t)llabs(aint-bint) <= max_ulps)
return 1;
return 0;
}
int flt_equals(float a, float b)
{
union { float f; int32_t i; } ua;
union { float f; int32_t i; } ub;
int32_t aint, bint;
if (a == b)
return 1;
ua.f = a; aint = ua.i;
ub.f = b; bint = ub.i;
if (aint < 0)
aint = BIT31 - aint;
if (bint < 0)
bint = BIT31 - bint;
if ((uint32_t)abs(aint-bint) <= flt_max_ulps)
return 1;
return 0;
}