ref: 156c048f0239b99b870cf1694c6886e39b85ae09
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); }