Thread: Fixing geometic calculation

Fixing geometic calculation

From
Paul Matthews
Date:
Let us consider the ordering of real numbers in postgres. As you can see
from
the results below it has clearly returned the correct results.
 select( 1.0000000000000000 = 1.0000000000000002 ); => f select( 1.0000000000000000 < 1.0000000000000002 ); => t
select(1.0000000000000000 > 1.0000000000000002 ); => f
 

Imagine the situation however where postgres returned the following
values to
simple numerical inequalities. In such a case postgresql would be clearly
defective and unfit for purpose.
 select( 1.000000 = 1.000001 ); => f select( 1.000000 < 1.000001 ); => f select( 1.000000 > 1.000001 ); => f

If such a situation is unacceptable for the real number line, then in
what way
can it be acceptable for the real number plain.
 select( point(1.00000,0)  <> point(1.00001,0) ); => f select( point(1.00000,0)  << point(1.00001,0) ); => f select(
point(1.00000,0) >> point(1.00001,0) ); => f select( point(1.00000,0) <-> point(1.00001,0) ); => 1.00000000000655e-05
 

We have two points with a finite separation in the x axis. Postgres
thinks they
are not the same point, nor one left of the other, nor to the right. This is
clearly a both a physical and logical impossibility.

The cause of this is the ill conceived FP* macros. They seem represent a
solution to a problem that simply does not exist.

The first effect of these macros is to reduce the accuracy of all geometric
comparisons from double precision, to less than single precision. The
following
program correctly prints the correct answer. Whereas as we have seen above,
postgres falls in a heap.
 int main() {   float f = 1.00000;   float g = 1.00001;   if( f==g ) { printf( "f=g\n" ); }   if( f<g )  { printf(
"f<g\n"); }   if( f>g )  { printf( "f>g\n" ); }   return 0; }
 

The second effect is to take operations that would of worked correctly
even in
single precision, and to cause them to produce nonsensical result. For
example
points that can be both inside and outside a polygon at the same time.

Simple analysis of the postgres source code shows that the only places
where the
FPzero, FPeq, FPne, FPlt, FPle FPgt and FPge macros are defined and used
are in
the src/backend/utils/adt/geo_ops.c and src/include/utils/geo_decls.h files.

What is the justification for these macros? Why do they only affect
geometric
calculations, and not all numeric calculations? Why should these macro's
not be
abandoned?

Does anyone any any objections to me:
1) removing these macros, or at least disabling EPSILON by default.
2) adding in the obviously missing operators (ie: box @> point)




Re: Fixing geometic calculation

From
Kenneth Marshall
Date:
On Fri, Aug 07, 2009 at 11:29:47PM +1000, Paul Matthews wrote:
> Let us consider the ordering of real numbers in postgres. As you can see
> from
> the results below it has clearly returned the correct results.
> 
>   select( 1.0000000000000000 = 1.0000000000000002 ); => f
>   select( 1.0000000000000000 < 1.0000000000000002 ); => t
>   select( 1.0000000000000000 > 1.0000000000000002 ); => f
> 
> Imagine the situation however where postgres returned the following
> values to
> simple numerical inequalities. In such a case postgresql would be clearly
> defective and unfit for purpose.
> 
>   select( 1.000000 = 1.000001 ); => f
>   select( 1.000000 < 1.000001 ); => f
>   select( 1.000000 > 1.000001 ); => f
> 
> If such a situation is unacceptable for the real number line, then in
> what way
> can it be acceptable for the real number plain.
> 
>   select( point(1.00000,0)  <> point(1.00001,0) ); => f
>   select( point(1.00000,0)  << point(1.00001,0) ); => f
>   select( point(1.00000,0)  >> point(1.00001,0) ); => f
>   select( point(1.00000,0) <-> point(1.00001,0) ); => 1.00000000000655e-05
> 
> We have two points with a finite separation in the x axis. Postgres
> thinks they
> are not the same point, nor one left of the other, nor to the right. This is
> clearly a both a physical and logical impossibility.
> 
> The cause of this is the ill conceived FP* macros. They seem represent a
> solution to a problem that simply does not exist.
> 
> The first effect of these macros is to reduce the accuracy of all geometric
> comparisons from double precision, to less than single precision. The
> following
> program correctly prints the correct answer. Whereas as we have seen above,
> postgres falls in a heap.
> 
>   int main() {
>     float f = 1.00000;
>     float g = 1.00001;
>     if( f==g ) { printf( "f=g\n" ); }
>     if( f<g )  { printf( "f<g\n" ); }
>     if( f>g )  { printf( "f>g\n" ); }
>     return 0;
>   }
> 
> The second effect is to take operations that would of worked correctly
> even in
> single precision, and to cause them to produce nonsensical result. For
> example
> points that can be both inside and outside a polygon at the same time.
> 
> Simple analysis of the postgres source code shows that the only places
> where the
> FPzero, FPeq, FPne, FPlt, FPle FPgt and FPge macros are defined and used
> are in
> the src/backend/utils/adt/geo_ops.c and src/include/utils/geo_decls.h files.
> 
> What is the justification for these macros? Why do they only affect
> geometric
> calculations, and not all numeric calculations? Why should these macro's
> not be
> abandoned?
> 
> Does anyone any any objections to me:
> 1) removing these macros, or at least disabling EPSILON by default.
> 2) adding in the obviously missing operators (ie: box @> point)
> 

Hi Paul,

Floating point calculations always have a bit of inaccuracy
because at the very minimum some values do not have exact
floating point representations and the results can be
implimentation dependent. I think disabling EPLSILON by
default is a bad idea. In my work with numeric methods,
we actually calculated EPSILON for the system we where using
at runtime. Maybe postgresql could do the same on startup.

Regards,
Ken


Re: Fixing geometic calculation

From
Kenneth Marshall
Date:
On Fri, Aug 07, 2009 at 09:12:34AM -0500, Kenneth Marshall wrote:
> On Fri, Aug 07, 2009 at 11:29:47PM +1000, Paul Matthews wrote:
> > Let us consider the ordering of real numbers in postgres. As you can see
> > from
> > the results below it has clearly returned the correct results.
> > 
> >   select( 1.0000000000000000 = 1.0000000000000002 ); => f
> >   select( 1.0000000000000000 < 1.0000000000000002 ); => t
> >   select( 1.0000000000000000 > 1.0000000000000002 ); => f
> > 
> > Imagine the situation however where postgres returned the following
> > values to
> > simple numerical inequalities. In such a case postgresql would be clearly
> > defective and unfit for purpose.
> > 
> >   select( 1.000000 = 1.000001 ); => f
> >   select( 1.000000 < 1.000001 ); => f
> >   select( 1.000000 > 1.000001 ); => f
> > 
> > If such a situation is unacceptable for the real number line, then in
> > what way
> > can it be acceptable for the real number plain.
> > 
> >   select( point(1.00000,0)  <> point(1.00001,0) ); => f
> >   select( point(1.00000,0)  << point(1.00001,0) ); => f
> >   select( point(1.00000,0)  >> point(1.00001,0) ); => f
> >   select( point(1.00000,0) <-> point(1.00001,0) ); => 1.00000000000655e-05
> > 
> > We have two points with a finite separation in the x axis. Postgres
> > thinks they
> > are not the same point, nor one left of the other, nor to the right. This is
> > clearly a both a physical and logical impossibility.

Actually, quantum theory will allow this to happen. :)

> > 
> > The cause of this is the ill conceived FP* macros. They seem represent a
> > solution to a problem that simply does not exist.
> > 
> > The first effect of these macros is to reduce the accuracy of all geometric
> > comparisons from double precision, to less than single precision. The
> > following
> > program correctly prints the correct answer. Whereas as we have seen above,
> > postgres falls in a heap.
> > 
> >   int main() {
> >     float f = 1.00000;
> >     float g = 1.00001;
> >     if( f==g ) { printf( "f=g\n" ); }
> >     if( f<g )  { printf( "f<g\n" ); }
> >     if( f>g )  { printf( "f>g\n" ); }
> >     return 0;
> >   }
> > 
> > The second effect is to take operations that would of worked correctly
> > even in
> > single precision, and to cause them to produce nonsensical result. For
> > example
> > points that can be both inside and outside a polygon at the same time.
> > 
> > Simple analysis of the postgres source code shows that the only places
> > where the
> > FPzero, FPeq, FPne, FPlt, FPle FPgt and FPge macros are defined and used
> > are in
> > the src/backend/utils/adt/geo_ops.c and src/include/utils/geo_decls.h files.
> > 
> > What is the justification for these macros? Why do they only affect
> > geometric
> > calculations, and not all numeric calculations? Why should these macro's
> > not be
> > abandoned?
> > 
> > Does anyone any any objections to me:
> > 1) removing these macros, or at least disabling EPSILON by default.
> > 2) adding in the obviously missing operators (ie: box @> point)
> > 
> 
> Hi Paul,
> 
> Floating point calculations always have a bit of inaccuracy
> because at the very minimum some values do not have exact
> floating point representations and the results can be
> implimentation dependent. I think disabling EPLSILON by
> default is a bad idea. In my work with numeric methods,
> we actually calculated EPSILON for the system we where using
> at runtime. Maybe postgresql could do the same on startup.
> 
> Regards,
> Ken


Re: Fixing geometic calculation

From
Sam Mason
Date:
On Fri, Aug 07, 2009 at 09:49:41AM -0500, Kenneth Marshall wrote:
> On Fri, Aug 07, 2009 at 09:12:34AM -0500, Kenneth Marshall wrote:
> > On Fri, Aug 07, 2009 at 11:29:47PM +1000, Paul Matthews wrote:
> > > We have two points with a finite separation in the x axis.
> > > Postgres thinks they are not the same point, nor one left of the
> > > other, nor to the right. This is clearly a both a physical and
> > > logical impossibility.
>
> Actually, quantum theory will allow this to happen. :)

I'm not a physicist, but I don't think it does.  QM defines the
probability distribution within which the particle will be found.  Once
you've actually observed both "points" you will know their physical
relation--you'll also have given them energy them so next time you look
they'll be somewhere else, but the act of observation causes the above
distribution to be collapsed.  This sidesteps the whole issue of the
fact that points in PG are defined in euclidean space and do indeed
have a definite location and can be compared at all times---they don't
arbitrarily go jumping off millions of miles away or being annihilated
by their anti-particle just because it's possible.

I would agree with Paul that EPSILON is a hack and probably should be
removed.  However it will cause user visible changes so it's not quite
as simple as that to change.  I don't have anything really very useful
to add apart from saying that maybe the default should be the other way
around?

--  Sam  http://samason.me.uk/


Re: Fixing geometic calculation

From
Kenneth Marshall
Date:
On Fri, Aug 07, 2009 at 04:16:56PM +0100, Sam Mason wrote:
> On Fri, Aug 07, 2009 at 09:49:41AM -0500, Kenneth Marshall wrote:
> > On Fri, Aug 07, 2009 at 09:12:34AM -0500, Kenneth Marshall wrote:
> > > On Fri, Aug 07, 2009 at 11:29:47PM +1000, Paul Matthews wrote:
> > > > We have two points with a finite separation in the x axis.
> > > > Postgres thinks they are not the same point, nor one left of the
> > > > other, nor to the right. This is clearly a both a physical and
> > > > logical impossibility.
> >
> > Actually, quantum theory will allow this to happen. :)
> 
> I'm not a physicist, but I don't think it does.  QM defines the
> probability distribution within which the particle will be found.  Once
> you've actually observed both "points" you will know their physical
> relation--you'll also have given them energy them so next time you look
> they'll be somewhere else, but the act of observation causes the above
> distribution to be collapsed.  This sidesteps the whole issue of the
> fact that points in PG are defined in euclidean space and do indeed
> have a definite location and can be compared at all times---they don't
> arbitrarily go jumping off millions of miles away or being annihilated
> by their anti-particle just because it's possible.
> 
> I would agree with Paul that EPSILON is a hack and probably should be
> removed.  However it will cause user visible changes so it's not quite
> as simple as that to change.  I don't have anything really very useful
> to add apart from saying that maybe the default should be the other way
> around?
> 
> -- 
>   Sam  http://samason.me.uk/
> 

It was definitely a tongue-in-cheek response since QT is not really
a topic for this mailing list. However, removing EPSILON completely
is not a good idea for the exact reason it was included originally.
Floating point numbers are approximations and since their precision
is neccessarily limited this fact must be included in any calculation
using them. I do agree that hard-coding it to a value that does not
reflect the reality of the calculation is not good. It would be
better to have a GUC to allow it to be specified than to have it
be zero. Maybe one setting would allow the system to calculate the
appropriate value for EPSILON based on the hardward. One way to
address the duplicity issue is to define for yourself what it means
if a point is both inside and outside, i.e. in this case the point
is always defined to be inside or the point is always defined to
be outside.

Regards,
Ken


Re: Fixing geometic calculation

From
Tom Lane
Date:
Sam Mason <sam@samason.me.uk> writes:
> I would agree with Paul that EPSILON is a hack and probably should be
> removed.

It's a hack but it's dealing with an extremely real problem, namely
the built-in inaccuracy of floating-point arithmetic.  You can't just
close your eyes to that and hope that everything will be okay.

A quick look through the geometry sources says that we might not be
critically dependent on anything except the assumption that two values
that aren't FPeq() will have a nonzero difference.  (If you think this
is a tautology, you don't know enough about floating point arithmetic
to be qualified to offer an opinion here...)  We might be able to base
a tighter comparison procedure on that rule.  It would take a lot more
investigation to be sure though.
        regards, tom lane


Re: Fixing geometic calculation

From
Sam Mason
Date:
On Fri, Aug 07, 2009 at 10:29:27AM -0500, Kenneth Marshall wrote:
> On Fri, Aug 07, 2009 at 04:16:56PM +0100, Sam Mason wrote:
> > points in PG [..] don't
> > arbitrarily go jumping off millions of miles away or being annihilated
> > by their anti-particle just because it's possible.

> It was definitely a tongue-in-cheek response since QT is not really
> a topic for this mailing list.

Yup, I know.  Hence my somewhat over the top examples.

> > I would agree with Paul that EPSILON is a hack and probably should be
> > removed.  However it will cause user visible changes so it's not quite
> > as simple as that to change.  I don't have anything really very useful
> > to add apart from saying that maybe the default should be the other way
> > around?
>
> However, removing EPSILON completely
> is not a good idea for the exact reason it was included originally.

Hum, I think it's good in some limited situations but not by default.  I
personally think that PG should be exposing rawer access here, mainly
because FP math is hard to get right and the more we fiddle trying to
make it easier to appear to do the right thing in the common case the
more general cases become impossible.  It's similar to the auto TEXT
casting thing that was changed in 8.3, but at least you get a nice error
when things aren't automatically cast to TEXT.

There are also much more reliable ways of solving the inaccuracies than
what's done now by just relying on a simple test, interval arithmetic
is my favorite at the moment but it is slower and can make thins more
complicated.

--  Sam  http://samason.me.uk/


Re: Fixing geometic calculation

From
Sam Mason
Date:
On Fri, Aug 07, 2009 at 11:40:58AM -0400, Tom Lane wrote:
> Sam Mason <sam@samason.me.uk> writes:
> > I would agree with Paul that EPSILON is a hack and probably should be
> > removed.
> 
> It's a hack but it's dealing with an extremely real problem, namely
> the built-in inaccuracy of floating-point arithmetic.  You can't just
> close your eyes to that and hope that everything will be okay.

Yes, I know it's a fiddle to get right.  Choosing the right primitives
is generally the most difficult part.

> A quick look through the geometry sources says that we might not be
> critically dependent on anything except the assumption that two values
> that aren't FPeq() will have a nonzero difference.

Sorry, I'm struggling to parse that.  I think it's all the double
negatives.  Are you saying that HYPOT() should really return zero when
it's currently giving back would be FPzero?

> (If you think this
> is a tautology, you don't know enough about floating point arithmetic
> to be qualified to offer an opinion here...)

I think I have a reasonable idea about FP arithmetic, I have had to
worry about rounding modes and such like before.  Never tried to write a
FP emulator though so I'm sure I could know more.

--  Sam  http://samason.me.uk/


Re: Fixing geometic calculation

From
Tom Lane
Date:
Sam Mason <sam@samason.me.uk> writes:
> Sorry, I'm struggling to parse that.  I think it's all the double
> negatives.  Are you saying that HYPOT() should really return zero when
> it's currently giving back would be FPzero?

No, I'm worried about code that supposes that it can divide by (x - y)
after testing that FPeq(x,y) is not true.  point_sl() for instance.

We could perhaps fix those specific issues by testing the difference
explicitly instead of doing it like that.  But there's still the overall
problem of error accumulation ...
        regards, tom lane


Re: Fixing geometic calculation

From
Sam Mason
Date:
On Fri, Aug 07, 2009 at 12:50:39PM -0400, Tom Lane wrote:
> Sam Mason <sam@samason.me.uk> writes:
> > Sorry, I'm struggling to parse that.  I think it's all the double
> > negatives.  Are you saying that HYPOT() should really return zero when
> > it's currently giving back would be FPzero?
> 
> No, I'm worried about code that supposes that it can divide by (x - y)
> after testing that FPeq(x,y) is not true.  point_sl() for instance.

OK, but I'm still not sure what you're getting at.  If it's infinities
and NaNs then they shouldn't matter and will be taken care of by the
normal FP rules anyway.

> We could perhaps fix those specific issues by testing the difference
> explicitly instead of doing it like that.  But there's still the overall
> problem of error accumulation ...

Errors will accumulate whatever happens, that's why things like interval
arithmetic exist that usefully track those errors and why I said testing
EPSILON isn't a useful.

--  Sam  http://samason.me.uk/


Re: Fixing geometic calculation

From
Tom Lane
Date:
Sam Mason <sam@samason.me.uk> writes:
> On Fri, Aug 07, 2009 at 12:50:39PM -0400, Tom Lane wrote:
>> No, I'm worried about code that supposes that it can divide by (x - y)
>> after testing that FPeq(x,y) is not true.  point_sl() for instance.

> OK, but I'm still not sure what you're getting at.

Underflow.  x!=y does not imply (x-y) != 0, if x and y are sufficiently
small and close together.  The difference could underflow to zero.
        regards, tom lane


Re: Fixing geometic calculation

From
Sam Mason
Date:
On Fri, Aug 07, 2009 at 02:13:26PM -0400, Tom Lane wrote:
> Sam Mason <sam@samason.me.uk> writes:
> > On Fri, Aug 07, 2009 at 12:50:39PM -0400, Tom Lane wrote:
> >> No, I'm worried about code that supposes that it can divide by (x - y)
> >> after testing that FPeq(x,y) is not true.  point_sl() for instance.
> 
> > OK, but I'm still not sure what you're getting at.
> 
> Underflow.  x!=y does not imply (x-y) != 0, if x and y are sufficiently
> small and close together.  The difference could underflow to zero.

I've just realized why this discussion hasn't been making any sense.
I thought you were talking about correctness of the code with EPSILON
still there and not about what would happen if EPSILON was removed.
Thanks for the patience.

If EPSILON is indeed removed then yes, this will become a problem and
the easiest fix would seem to be to calculate the difference first and
test it explicitly.

The "error accumulation" comment also makes sense now!  Does anyone
know the original use case for using the EPSILON (need some shorthand
for that, a mail client that supports Unicode?) based comparisons so
liberally?  It only makes sense to me if they're done right at the end
of all the calculations, not all the way though.  What defines the "end"
seems up to the user as well, or am I missing something.

--  Sam  http://samason.me.uk/


Re: Fixing geometic calculation

From
Greg Stark
Date:
On Fri, Aug 7, 2009 at 7:13 PM, Tom Lane<tgl@sss.pgh.pa.us> wrote:
> Sam Mason <sam@samason.me.uk> writes:
>> On Fri, Aug 07, 2009 at 12:50:39PM -0400, Tom Lane wrote:
>>> No, I'm worried about code that supposes that it can divide by (x - y)
>>> after testing that FPeq(x,y) is not true.  point_sl() for instance.
>
>> OK, but I'm still not sure what you're getting at.
>
> Underflow.  x!=y does not imply (x-y) != 0, if x and y are sufficiently
> small and close together.  The difference could underflow to zero.


Actually I don't think subtraction can underflow with IEEE floats but
I don't think we want to count on IEEE floats everywhere. Even if we
did there's the risk on intel that FPeq() gets called on values which
have just been calculated and are still in registers but then get
spilled to RAM and lose precision before the division happens.

--
greg
http://mit.edu/~gsstark/resume.pdf


Re: Fixing geometic calculation

From
Sam Mason
Date:
On Fri, Aug 07, 2009 at 07:48:15PM +0100, Greg Stark wrote:
> On Fri, Aug 7, 2009 at 7:13 PM, Tom Lane<tgl@sss.pgh.pa.us> wrote:
> > Underflow.  x!=y does not imply (x-y) != 0, if x and y are sufficiently
> > small and close together.  The difference could underflow to zero.
> 
> Actually I don't think subtraction can underflow with IEEE floats but
> I don't think we want to count on IEEE floats everywhere. Even if we
> did there's the risk on intel that FPeq() gets called on values which
> have just been calculated and are still in registers but then get
> spilled to RAM and lose precision before the division happens.

If it does one subtraction in registers you can be reasonably certain
the other will be, either way just doing the subtraction and explicitly
testing if it's zero will do the right thing--the semantics of C are bad
but not that bad.

--  Sam  http://samason.me.uk/


Re: Fixing geometic calculation

From
Paul Matthews
Date:
Tom Lane wrote:
> It's a hack but it's dealing with an extremely real problem, namely
> the built-in inaccuracy of floating-point arithmetic.  You can't just
> close your eyes to that and hope that everything will be okay.
>   
If the above statement was true, then the FP* macros should be extended
to all numerical calculations in postgres. And I don't think anyone here
would suggest that doing that, as the results, as per my previous email,
would be immediately and clearly ludicrous.

Yes, floating point arithmetic has built in inaccuracy. However the FP*
macros produce results that are even less accurate than single point
arithmetic, and less accurate than my 25 year old calculator. And if I
could find my slide rule, it could probably do better.

At best, the EPSILON value might have be appropriate for single
precision arithmetic, but in no way is it  appropriate for double
precision arithmetic. Visual C++ defines EPSILON as "The difference
between 1 and smallest value greater than 1". The postgres EPSILON is 10
orders of magnitude greater than the precision supported by the
underlying hardware. This is clearly preposterous.
 EPSILON in postgres is 1.0E-06 EPSILON for floats is 1.19209e-07 EPSILON for doubles is 2.22045E-016




Re: Fixing geometic calculation

From
Paul Matthews
Date:
Tom Lane wrote:
> No, I'm worried about code that supposes that it can divide by (x - y)
> after testing that FPeq(x,y) is not true.  point_sl() for instance.
>
> We could perhaps fix those specific issues by testing the difference
> explicitly instead of doing it like that.  But there's still the overall
> problem of error accumulation ...
>
>             regards, tom lane
>   
IEEE754 does not allow two number X and Y, such that X!=Y and (X-Y)==0.
And since IEEE754 has been around since the 70's or 80's I think we can
start relying on its existence and behavior by now.

I don't see why it should be is postgres job to worry about error
accumulation. And then why only worry about it in geometric
calculations. Where in normal postgres code, or in perl, python, lisp,
or any in any other general purpose tool, is it where you ask it to tell
you 1.0+1.0 it responds with 2.0 +or- 0.0000000000000001?


Re: Fixing geometic calculation

From
Paul Matthews
Date:
Paul Matthews wrote:
>   EPSILON in postgres is 1.0E-06
>   EPSILON for floats is 1.19209e-07
>   EPSILON for doubles is 2.22045E-016
>   
Bad form to reply to my own post and all. If EPSILON for double
represented 1mm, then postgres is rounding to the nearest 10,000,000 km.
Since its only about 380,000 km to the moon it's a good thing <diety of
choice> does not use postgres.


Re: Fixing geometic calculation

From
marcin mank
Date:
On Sat, Aug 8, 2009 at 3:07 AM, Paul Matthews<plm@netspace.net.au> wrote:

> IEEE754 does not allow two number X and Y, such that X!=Y and (X-Y)==0.
> And since IEEE754 has been around since the 70's or 80's I think we can
> start relying on its existence and behavior by now.
>

You are correct, I think, though this does not solve the division problem:

$ cat t.c
#include <stdio.h>
int main(){   double a=1.000001e-307, b=1e-307, c=a-b;   printf("a=%le, b=%le, c=%le, c==0:%d, a==b:%d
1/c=%le\n",a,b,c,c==0,a==b,1.0/c);   return 0;
}
$ gcc -Wall -O2 t.c
$ ./a.out
a=1.000001e-307, b=1.000000e-307, c=1.000000e-313, c==0:0, a==b:0 1/c=inf


Greetings
Marcin Mańk


Re: Fixing geometic calculation

From
Paul Matthews
Date:
marcin mank wrote:
> You are correct, I think, though this does not solve the division problem:
>   
As a first goal I'm just attempting to reduce the EPSILON from 1.0E-6
down to 1.0E-015 (give or take). The current regression test suite works
fine down to 1.0E-09. At 1.0E-10 errors appear, not in the geometry
sections, but in the select_view test of all things. This is most likely
due to postgresql now giving the more correct (hence different) answers. 

A real test suite is needed for this. Setting up PostGIS + MySQL +
OtherCommerical for comparison purposes. The other problem is many of
the basic geometric operators in postgres, such a left of, above, etc,
are so incorrectly implemented, they are not even wrong.