111 lines
2.2 KiB
C
111 lines
2.2 KiB
C
|
#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
|
||
|
|
extern int64_t llabs(int64_t j);
|
||
|
|
#endif
|
||
|
|
|
||
|
|
int dbl_equals(double a, double b)
|
||
|
|
{
|
||
|
|
int64_t aint, bint;
|
||
|
|
|
||
|
|
if (a == b)
|
||
|
|
return 1;
|
||
|
|
aint = *(int64_t*)&a;
|
||
|
|
bint = *(int64_t*)&b;
|
||
|
|
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)
|
||
|
|
{
|
||
|
|
int32_t aint, bint;
|
||
|
|
|
||
|
|
if (a == b)
|
||
|
|
return 1;
|
||
|
|
aint = *(int32_t*)&a;
|
||
|
|
bint = *(int32_t*)&b;
|
||
|
|
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;
|
||
|
|
}
|
||
|
|
|
||
|
|
double randn()
|
||
|
|
{
|
||
|
|
double s, vre, vim, ure, uim;
|
||
|
|
static double next = -42;
|
||
|
|
|
||
|
|
if (next != -42) {
|
||
|
|
s = next;
|
||
|
|
next = -42;
|
||
|
|
return s;
|
||
|
|
}
|
||
|
|
do {
|
||
|
|
ure = rand_double();
|
||
|
|
uim = rand_double();
|
||
|
|
vre = 2*ure - 1;
|
||
|
|
vim = 2*uim - 1;
|
||
|
|
s = vre*vre + vim*vim;
|
||
|
|
} while (s >= 1);
|
||
|
|
s = sqrt(-2*log(s)/s);
|
||
|
|
next = s * vre;
|
||
|
|
return s * vim;
|
||
|
|
}
|