Re: [HACKERS] [PATCH] Improve geometric types - Mailing list pgsql-hackers

At Mon, 2 Oct 2017 08:23:49 -0400, Robert Haas <robertmhaas@gmail.com> wrote in
<CA+TgmoYsgw0TcjJQ1CE_6vDOxgEhxYQkfNx93Mfwx23WOLM0NA@mail.gmail.com>
> On Mon, Oct 2, 2017 at 4:23 AM, Kyotaro HORIGUCHI
> <horiguchi.kyotaro@lab.ntt.co.jp> wrote:
> > For other potential reviewers:
> >
> > I found the origin of the function here.
> >
> > https://www.postgresql.org/message-id/4A90BD76.7070804@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/407d949e0908222139t35ad3ad2q3e6b15646a27dd64@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
is0)  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
protectsagainst * 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}
};inti;
 

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