ref: c57f69776103121107ad70d86a7b0efbc2553daa
parent: bc7b9cf7c49595c758d68c30ee1179740e269121
author: S. Gilles <sgilles@math.umd.edu>
date: Thu Jun 28 06:37:51 EDT 2018
Add macros for relative error to lib/math/ancillary.
--- /dev/null
+++ b/lib/math/ancillary/ulp.gp
@@ -1,0 +1,27 @@
+/*
+ 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));
+}