shithub: femtolisp

ref: 209b77a534b953d8197b6d0dc2ddb4db16f8fb80
dir: /llt/fp.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;
}