shithub: mc

ref: 1d6cb9bf584e35be0786467e1786e60ac909ae21
dir: /lib/math/ancillary/ulp.gp/

/*
  I always end up need this for debugging
 */

{ ulp32(a) = my(aa, q);
        aa = abs(a);
        if(aa < 2^(-150),return(2^(-126 - 23)),);
        q = floor(log(aa)/log(2));
        if(q < -126,q=-126,);
        return(2^(q-23));
}

{ ulp64(a) = my(aa, q);
        aa = abs(a);
        if(aa < 2^(-2000),return(2^(-1022 - 52)),);
        q = floor(log(aa)/log(2));
        if(q < -1022,q=-1022,);
        return(2^(q-52));
}

{ err32(x, y) =
        return(abs(x-y)/ulp32(x));
}

{ err64(x, y) =
        return(abs(x-y)/ulp64(x));
}

{ flt32fromuint32(a) = my(n, e, s);
        n = bitand(a, 0x80000000) != 0;
        e = bitand(a, 0x7f800000) >> 23;
        s = bitand(a, 0x007fffff);

        if(e != 0, s = bitor(s, 0x00800000), s = 2.0 * s);
        s = s * 2.0^(-23);
        e = e - 127;
        return((-1)^n * s * 2^(e));
}

{ flt64fromuint64(a) = my(n, e, s);
        n = bitand(a, 0x8000000000000000) != 0;
        e = bitand(a, 0x7ff0000000000000) >> 52;
        s = bitand(a, 0x000fffffffffffff);

        if(e != 0, s = bitor(s, 0x0010000000000000), s = 2.0 * s);
        s = s * 2.0^(-52);
        e = e - 1023;
        return((-1)^n * 2^(e) * s);
}