bin/170206: complex arcsinh, log, etc.

[Bruce Evans]

On Sat, 28 Jul 2012, Stephen Montgomery-Smith wrote:

> On 07/28/2012 12:25 AM, Bruce Evans wrote:
>> 
>> #define    DE    DBL_EPSILON        // for clarity
>> 
>> (1)   1 + DE/2        = 1         (half way case rounded down to even)
>> (2)   1 + DE/2 + DE/2 = 1         (double rounding)
>> (3)   DE/2 + DE/2 + 1 = 1 + DE    (null rounding)
>> 
>> We want to add -1 to a value near 1 like the above.  Now a leading 1
>> in the above will cancel with the -1, and the the order in (3) becomes
>> the inaccurate one.
>
> Yes, but in my situation, I am rather sure that when I am adding highest to 
> lowest that this won't occur.  I am starting with -1, then adding something 
> close to 1, then adding lots of smaller terms.  And I find it very plausible 
> that the kind of situation you describe won't happen. x0*x0 is close to 1. 
> x0*x1 is at most sqrt(DE) times smaller.  And so on.  So I think the kind of 
> situation you describe should never happen.

Ahem.  FP^2 space is not nearly as large as the metaverse (only 2^256
cases even for sparc64), but it is large enough so that almost
everything that can happen in it does happen in it.  You are right
that problems are far away with the x* terms (x0*x0 had better not be
very close to 1 unless it is exactly 1, since it it is too close then
it will no longer be many times larger than x0*x1 after subtracting 1
from it; the other cases for x* are simpler).  The problem is with the
additional y* terms.  x and y are independent, so for many or most x,
there are many y's with bits that cause half-way cases when combined with x.
After splitting and squaring, the bits move around, so it hard to generate
or control the offending y's.

> As I said, I don't have a mathematical proof that the kind of thing you 
> describe can NEVER happen.  I just have never observed it happen.

There might be a measly 2^128 bad cases out of 2^256.  Then no one would
even observe them by chance :-).  But half-way cases are fairly common.

Bruce