On Saturday, 7 November, 2015 09:08, James K. Lowden <jklowden at
schemamania.org> wrote:
> On Fri, 06 Nov 2015 22:16:57 -0700
> "Keith Medcalf" <kmedcalf at dessus.com> wrote:
> > I wrote a function called "ulps" which can be used as an extension to
> > SQLite3
> Bravo, Keith!
> One suggestion, if I may. If you expect "ulps" to be used to test for
> equality, perhaps it would be more convenient to ignore the sign.
> Something like
> fequal(x, y, delta) === abs(ulps(x -y)) < delta
> might express the idea more directly?
> --jkl
Simple wrappers can be created that emulate all the equality tests.
feq(x, y, delta) --> abs(ulps(x, y)) <= delta : x == y
fne(x, y, delta) --> abs(ulps(x, y)) > delta : x != y
fgt(x, y, delta) --> ulps(x, y) > delta : x > y
fge(x, y, delta) --> ulps(x, y) => -delta : x >= y
flt(x, y, delta) --> ulps(x, y) < -delta : x < y
fle(x, y, delta) --> ulps(x, y) <= delta : x <= y
ala:
static double epsilon(double value)
{
int exponent;
double mantissa = frexp(value, &exponent);
double epsilon = ldexp(1.0, exponent - 53);
return epsilon;
}
static double distance(double x, double y)
{
return (x - y) / epsilon(x);
}
SQLITE_PRIVATE void _ulp(sqlite3_context *context, int argc, sqlite3_value
**argv)
{
sqlite3_result_double(context, epsilon(sqlite3_value_double(argv[0])));
}
SQLITE_PRIVATE void _ulps(sqlite3_context *context, int argc, sqlite3_value
**argv)
{
sqlite3_result_double(context, distance(sqlite3_value_double(argv[0]),
sqlite3_value_double(argv[1])));
}
SQLITE_PRIVATE void _feq(sqlite3_context *context, int argc, sqlite3_value
**argv)
{
double ulps = distance(sqlite3_value_double(argv[0]),
sqlite3_value_double(argv[1]));
if (fabs(ulps) <= fabs(sqlite3_value_double(argv[2])))
sqlite3_result_int(context, 1);
else
sqlite3_result_int(context, 0);
}
SQLITE_PRIVATE void _fne(sqlite3_context *context, int argc, sqlite3_value
**argv)
{
double ulps = distance(sqlite3_value_double(argv[0]),
sqlite3_value_double(argv[1]));
if (fabs(ulps) > fabs(sqlite3_value_double(argv[2])))
sqlite3_result_int(context, 1);
else
sqlite3_result_int(context, 0);
}
SQLITE_PRIVATE void _fgt(sqlite3_context *context, int argc, sqlite3_value
**argv)
{
double ulps = distance(sqlite3_value_double(argv[0]),
sqlite3_value_double(argv[1]));
if (ulps > fabs(sqlite3_value_double(argv[2])))
sqlite3_result_int(context, 1);
else
sqlite3_result_int(context, 0);
}
SQLITE_PRIVATE void _fge(sqlite3_context *context, int argc, sqlite3_value
**argv)
{
double ulps = distance(sqlite3_value_double(argv[0]),
sqlite3_value_double(argv[1]));
if (ulps >= -fabs(sqlite3_value_double(argv[2])))
sqlite3_result_int(context, 1);
else
sqlite3_result_int(context, 0);
}
SQLITE_PRIVATE void _flt(sqlite3_context *context, int argc, sqlite3_value
**argv)
{
double ulps = distance(sqlite3_value_double(argv[0]),
sqlite3_value_double(argv[1]));
if (ulps < -fabs(sqlite3_value_double(argv[2])))
sqlite3_result_int(context, 1);
else
sqlite3_result_int(context, 0);
}
SQLITE_PRIVATE void _fle(sqlite3_context *context, int argc, sqlite3_value
**argv)
{
double ulps = distance(sqlite3_value_double(argv[0]),
sqlite3_value_double(argv[1]));
if (ulps <= fabs(sqlite3_value_int(argv[2])))
sqlite3_result_int(context, 1);
else
sqlite3_result_int(context, 0);
}
SQLITE_PRIVATE void _nxtaft(sqlite3_context *context, int argc, sqlite3_value
**argv)
{
double x = sqlite3_value_double(argv[0]);
double y = sqlite3_value_double(argv[1]);
sqlite3_result_double(context, _nextafter(x, y));
}
nErr += sqlite3_create_function(db, "nextafter", 2,
SQLITE_ANY|SQLITE_DETERMINISTIC, 0, _nxtaft, 0, 0);
nErr += sqlite3_create_function(db, "ulp", 1,
SQLITE_ANY|SQLITE_DETERMINISTIC, 0, _ulp, 0, 0);
nErr += sqlite3_create_function(db, "ulps", 2,
SQLITE_ANY|SQLITE_DETERMINISTIC, 0, _ulps, 0, 0);
nErr += sqlite3_create_function(db, "feq", 3,
SQLITE_ANY|SQLITE_DETERMINISTIC, 0, _feq, 0, 0);
nErr += sqlite3_create_function(db, "fne", 3,
SQLITE_ANY|SQLITE_DETERMINISTIC, 0, _fne, 0, 0);
nErr += sqlite3_create_function(db, "fgt", 3,
SQLITE_ANY|SQLITE_DETERMINISTIC, 0, _fgt, 0, 0);
nErr += sqlite3_create_function(db, "fge", 3,
SQLITE_ANY|SQLITE_DETERMINISTIC, 0, _fge, 0, 0);
nErr += sqlite3_create_function(db, "flt", 3,
SQLITE_ANY|SQLITE_DETERMINISTIC, 0, _flt, 0, 0);
nErr += sqlite3_create_function(db, "fle", 3,
SQLITE_ANY|SQLITE_DETERMINISTIC, 0, _fle, 0, 0);
Giving:
>sqlite < ulps.sql
create table test (x real, y real);
insert into test values (9.2 + 7.9 + 0 + 4.0 + 2.6 + 1.3, 25.0);
insert into test values (9.2 + 7.8 + 0 + 3.0 + 1.3 + 1.7, 23.0);
insert into test values (9.2 + 7.9 + 0 + 1.0 + 1.3 + 1.6, 21.0);
insert into test values (30.0, 25.0);
insert into test values (25.0, 30.0);
select x, y, printf('%.3e', x - y) as diff, case when x == y then 'yes' else
'no' end as equal,
printf('%.3e', ulps(x, y)) as ulps,
feq(x,y,5),
fne(x,y,5),
fgt(x,y,5),
fge(x,y,5),
flt(x,y,5),
fle(x,y,5)
from test;
x y diff equal ulps feq(x,y,5)
fne(x,y,5) fgt(x,y,5) fge(x,y,5) flt(x,y,5) fle(x,y,5)
---------- ---------- ---------- ---------- ---------- ----------
---------- ---------- ---------- ---------- ----------
25.0 25.0 3.553e-15 no 1.000e+00 1 0
0 1 0 1
23.0 23.0 0.000e+00 yes 0.000e+00 1 0
0 1 0 1
21.0 21.0 3.553e-15 no 1.000e+00 1 0
0 1 0 1
30.0 25.0 5.000e+00 no 1.407e+15 0 1
1 1 0 0
25.0 30.0 -5.000e+00 no -1.407e+15 0 1
0 0 1 1
update test set x = x * 3984.2534, y = y * 3984.2534;
select x, y, printf('%.3e', x - y) as diff, case when x == y then 'yes' else
'no' end as equal,
printf('%.3e', ulps(x, y)) as ulps,
feq(x,y,5),
fne(x,y,5),
fgt(x,y,5),
fge(x,y,5),
flt(x,y,5),
fle(x,y,5)
from test;
x y diff equal ulps feq(x,y,5)
fne(x,y,5) fgt(x,y,5) fge(x,y,5) flt(x,y,5) fle(x,y,5)
---------- ---------- ---------- ---------- ---------- ----------
---------- ---------- ---------- ---------- ----------
99606.335 99606.335 1.455e-11 no 1.000e+00 1 0
0 1 0 1
91637.8282 91637.8282 0.000e+00 yes 0.000e+00 1 0
0 1 0 1
83669.3214 83669.3214 1.455e-11 no 1.000e+00 1 0
0 1 0 1
119527.602 99606.335 1.992e+04 no 1.369e+15 0 1
1 1 0 0
99606.335 119527.602 -1.992e+04 no -1.369e+15 0 1
0 0 1 1
update test set x = x * 3984.2534e60, y = y * 3984.2534e60;
select x, y, printf('%.3e', x - y) as diff, case when x == y then 'yes' else
'no' end as equal,
printf('%.3e', ulps(x, y)) as ulps,
feq(x,y,5),
fne(x,y,5),
fgt(x,y,5),
fge(x,y,5),
flt(x,y,5),
fle(x,y,5)
from test;
x y diff equal ulps
feq(x,y,5) fne(x,y,5) fgt(x,y,5) fge(x,y,5) flt(x,y,5) fle(x,y,5)
-------------------- -------------------- ---------- ---------- ----------
---------- ---------- ---------- ---------- ---------- ----------
3.96856878885289e+68 3.96856878885289e+68 4.789e+52 no 1.000e+00
1 0 0 1 0 1
3.65108328574466e+68 3.65108328574466e+68 0.000e+00 yes 0.000e+00
1 0 0 1 0 1
3.33359778263643e+68 3.33359778263643e+68 9.578e+52 no 2.000e+00
1 0 0 1 0 1
4.76228254662347e+68 3.96856878885289e+68 7.937e+67 no 8.287e+14
0 1 1 1 0 0
3.96856878885289e+68 4.76228254662347e+68 -7.937e+67 no -1.657e+15
0 1 0 0 1 1
The whole lot could probably be simplified into far fewer functions using the
context to pass the actual operation and checking the arg count rather than
having the function table specify the number of arguments -- and the "delta"
could also be optional with the tolerance defaulting to a standard error
tolerance of (2 * ulp(1) * fabs(x)) where ulp(1) gives the machine epsilon.
