summaryrefslogtreecommitdiff
path: root/lib/math/ancillary/ulp.gp
blob: 9a4b15562438cf2380fdfdd76d4dec0c5b919af2 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
/*
  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 = 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 = s * 2.0^(-52);
        e = e - 1023;
        return((-1)^n * 2^(e) * s);
}