At Mon, 2 Oct 2017 08:23:49 -0400, Robert Haas <robertmh...@gmail.com> wrote in
<ca+tgmoysgw0tcjjq1ce_6vdoxgehxyqkfnx93mfwx23wolm...@mail.gmail.com>
> On Mon, Oct 2, 2017 at 4:23 AM, Kyotaro HORIGUCHI
> <horiguchi.kyot...@lab.ntt.co.jp> wrote:
> > For other potential reviewers:
> >
> > I found the origin of the function here.
> >
> > https://www.postgresql.org/message-id/4a90bd76.7070...@netspace.net.au
> > https://www.postgresql.org/message-id/AANLkTim4cHELcGPf5w7Zd43_dQi_2RJ_b5_F_idSSbZI%40mail.gmail.com
> >
> > And the reason for pg_hypot is seen here.
> >
> > https://www.postgresql.org/message-id/407d949e0908222139t35ad3ad2q3e6b15646a27d...@mail.gmail.com
> >
> > I think the replacement is now acceptable according to the discussion.
> > ======
>
> I think if we're going to do this it should be separated out as its
> own patch.
+1
> Also, I think someone should explain what the reasoning
> behind the change is. Do we, for example, foresee that the built-in
> code might be faster or less likely to overflow? Because we're
> clearly taking a risk -- most trivially, that the BF will break, or
> more seriously, that some machines will have versions of this function
> that don't actually behave quite the same.
>
> That brings up a related point. How good is our test case coverage
> for hypot(), especially in strange corner cases, like this one
> mentioned in pg_hypot()'s comment:
>
> * This implementation conforms to IEEE Std 1003.1 and GLIBC, in that the
> * case of hypot(inf,nan) results in INF, and not NAN.
I'm not sure how precise we practically need them to be
identical. FWIW as a rough confirmation on my box, I compared
hypot and pg_hypot for the following arbitrary choosed pairs of
parameters.
{2.2e-308, 2.2e-308},
{2.2e-308, 1.7e307},
{1.7e307, 1.7e307},
{1.7e308, 1.7e308},
{2.2e-308, DBL_MAX},
{1.7e308, DBL_MAX},
{DBL_MIN, DBL_MAX},
{DBL_MAX, DBL_MAX},
{1.7e307, INFINITY},
{2.2e-308, INFINITY},
{0, INFINITY},
{DBL_MIN, INFINITY},
{INFINITY, INFINITY},
{1, NAN},
{INFINITY, NAN},
{NAN, NAN},
Only the first pair gave slightly not-exactly-equal results but
it seems to do no harm. hypot set underflow flag.
0: hypot=3.111269837220809e-308 (== 0.0 is 0, < DBL_MIN is 0)
pg_hypot=3.11126983722081e-308 (== 0.0 is 0, < DBL_MIN is 0)
equal=0,
hypot(errno:0, inval:0, div0:0, of=0, uf=1),
pg_hypot(errno:0, inval:0, div0:0, of=0, uf=0)
But not sure how the both functions behave on other platforms.
> I'm potentially willing to commit a patch that just makes the
> pg_hypot() -> hypot() change and does nothing else, if there are not
> objections to that change, but I want to be sure that we'll know right
> away if that turns out to break.
--
Kyotaro Horiguchi
NTT Open Source Software Center
#include <stdio.h>
#include <float.h>
#include <math.h>
#include <fenv.h>
#include <errno.h>
double
pg_hypot(double x, double y)
{
double yx;
/* Handle INF and NaN properly */
if (isinf(x) || isinf(y))
return (double) INFINITY;
if (isnan(x) || isnan(y))
return (double) NAN;
/* Else, drop any minus signs */
x = fabs(x);
y = fabs(y);
/* Swap x and y if needed to make x the larger one */
if (x < y)
{
double temp = x;
x = y;
y = temp;
}
/*
* If y is zero, the hypotenuse is x. This test saves a few cycles in
* such cases, but more importantly it also protects against
* divide-by-zero errors, since now x >= y.
*/
if (y == 0.0)
return x;
/* Determine the hypotenuse */
yx = y / x;
return x * sqrt(1.0 + (yx * yx));
}
void
setfeflags(int *invalid, int *divzero, int *overflow, int *underflow)
{
int err = fetestexcept(FE_INVALID | FE_DIVBYZERO |
FE_OVERFLOW | FE_UNDERFLOW);
*invalid = ((err & FE_INVALID) != 0);
*divzero = ((err & FE_DIVBYZERO) != 0);
*overflow = ((err & FE_OVERFLOW) != 0);
*underflow = ((err & FE_UNDERFLOW) != 0);
}
int
main(void)
{
double x;
double y;
double p[][2] =
{
{2.2e-308, 2.2e-308},
{2.2e-308, 1.7e307},
{1.7e307, 1.7e307},
{1.7e308, 1.7e308},
{2.2e-308, DBL_MAX},
{1.7e308, DBL_MAX},
{DBL_MIN, DBL_MAX},
{DBL_MAX, DBL_MAX},
{1.7e307, INFINITY},
{2.2e-308, INFINITY},
{0, INFINITY},
{DBL_MIN, INFINITY},
{INFINITY, INFINITY},
{1, NAN},
{INFINITY, NAN},
{NAN, NAN},
{0.0, 0.0}
};
int i;
for (i = 0 ; p[i][0] != 0.0 || p[i][1] != 0.0 ; i++)
{
int errno1, errno2;
int invalid1 = 0, divzero1 = 0, overflow1 = 0, underflow1 = 0;
int invalid2 = 0, divzero2 = 0, overflow2 = 0, underflow2 = 0;
int cmp;
errno = 0;
feclearexcept(FE_ALL_EXCEPT);
x = hypot(p[i][0], p[i][1]);
errno1 = errno;
setfeflags(&invalid1, &divzero1, &overflow1, &underflow1);
errno = 0;
feclearexcept(FE_ALL_EXCEPT);
y = pg_hypot(p[i][0], p[i][1]);
errno2 = errno;
setfeflags(&invalid2, &divzero2, &overflow2, &underflow2);
/* assume NaN == NaN here */
if (isnan(x) && isnan(y))
cmp = 1;
else
cmp = (x == y);
printf("%d: hypot=%.16lg (== 0.0 is %d, < DBL_MIN is %d)\n
pg_hypot=%.16lg (== 0.0 is %d, < DBL_MIN is %d)\n equal=%d, hypot(errno:%d,
inval:%d, div0:%d, of=%d, uf=%d), pg_hypot(errno:%d, inval:%d, div0:%d, of=%d,
uf=%d)\n",
i, x, x == 0.0, x < DBL_MIN, y, y == 0.0, y <
DBL_MIN, cmp,
errno1, invalid1, divzero1, overflow1, underflow1,
errno2, invalid2, divzero2, overflow2, underflow2);
}
return 0;
}
--
Sent via pgsql-hackers mailing list (pgsql-hackers@postgresql.org)
To make changes to your subscription:
http://www.postgresql.org/mailpref/pgsql-hackers
