Re: [PATCH] Perform constant folding of math builtins

[Brad Lucier <lucier at math dot purdue dot edu>]

> 
> 
> > After reading Geoff Keating's post giving, for example, a/b/c -> a/(b*c)
> > as a valid transformation for -funsafe-math-optimizations, I realize it's
> > OK to accept sqrt(x)*sqrt(y)->sqrt(x*y) as OK under -ffast-math.
> >
> > I guess I'm just someone for whom this option does not hold a whole
> > lot of promise.
> 
> Hi Brad,
> 
> I think the situation is far less gloomy than you indicate above.  These
> optimizations should benefit much more code than they're likely to break.
> I doubt there's any major body of code that would be adversely affected,
> and those users that rely on long doubles probably understand the issues
> with "-ffast-math" and "-ffloat-store".

The issue isn't about long doubles, etc., it arises even with, e.g., PowerPC.
If I really wanted to be safe, i would write

#define len(v) sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2])

as (in Scheme notation)

(define (len v)

  (define (square x)
    (* x x))

  (let ((max-v (max (abs (vector-ref v 0)) 
                    (abs (vector-ref v 1)) 
                    (abs (vector-ref v 2)))))
    (* max-v (sqrt (+ (square (/ (vector-ref v 0) max-v))
                      (square (/ (vector-ref v 0) max-v))
                      (square (/ (vector-ref v 0) max-v)))))

etc.  This avoids overflow, underflow, etc.  I figure that if you
don't go to this trouble, then transforming (C notation)

sqrt(a[0]*a[0] + a[1]*a[1] + a[2]*a[2]) * sqrt(b[0]*b[0] + b[1]*b[1] + b[2]*b[2])

to

sqrt((a[0]*a[0] + a[1]*a[1] + a[2]*a[2]) * (b[0]*b[0] + b[1]*b[1] + b[2]*b[2]))

probably isn't so bad.  And if you do go to all that trouble, you're probably
going to set -fno-unsafe-math-optimizations anyway.

Of course, if you could ensure that all these temporary values were kept
in extended double precision on x86, you wouldn't have any problems either.
Perhaps that's why you mentioned long doubles.

Brad