Re: Slaying the HYPOTamus

Greg Stark wrote:
> We're implementing things like box_distance and point_distance which
> as it happens will already be doing earlier arithmetic on the doubles,
> so whatever HYPOT() does had better be consistent with that or the
> results will be just an inexplicable mishmash.
>
>
Let's look at what the HYPOT macro in PostgreSQL does right now:

#define HYPOT(A, B) sqrt((A)*(A)+(B)*(B))

And let's see how it is used in point_distance:

Datum
point_distance(PG_FUNCTION_ARGS)
{
Point *pt1 = PG_GETARG_POINT_P(0);
Point *pt2 = PG_GETARG_POINT_P(1);

PG_RETURN_FLOAT8(HYPOT(pt1->x - pt2->x, pt1->y - pt2->y));
}

Oh. Surprise. It does not look out for NaN's, InF's, overflows,
underflows or anything else. In addition, for the majority of it's input
space, it fails to work correctly at all. If x and y are both 1E200 the
hypotenuse should be 1.4142135e200, well within the range of the double.
However by naively squaring x yields 1E400, outside the range of a
double. Currently HYPOT() returns INFINITY, which is not the correct answer.

I am trying to introduce the hypot() function (where required) and
associated patch, in order to start correcting some of the more obvious
defects with the current geometry handling. Maybe my proposed hypot()
function is not perfect, but it's so far in front of current the HYPOT
macro is not funny.

> Incidentally, what on earth does it mean to multiply or divide a
> circle by a point?
>
>
Basically it's complex math. This comment, from my new geometry
implementation, might help.

/*-------------------------------------------------------------------------
* Additional Operators
*
* The +, -, * and / operators treat IPoint's as complex numbers.
* (This would indicate a requirement for a Complex type?)
*
* (a+bi)+(c+di) = ((a+c)+(b+d)i)
* (a+bi)-(c+di) = ((a-c)+(b-d)i)
* (a+bi)*(c+di) = ac + adi + bci + bdi^2
* = ac + (ad+bc)i - bd # as i^2 = -1
* = ((ac-bd)+(ad+bc)i)
* (a+bi)/(c+di) = (a+bi)(c-di) / (c+di)(c-di)
* = ((ac+bd) + (bc-ad)i) / (c^2+d^2)
* = ((ac + bd)/(c^2+d^2) + ((bc-ad)/(c^2+d^2))i)
*
* translation + IPoint_IPoint_add
* translation - IPoint_IPoint_sub
* multiplication * IPoint_IPoint_mul
* division / IPoint_IPoint_div
*-----------------------------------------------------------------------*/