Use of C99 extra long double math functions after r236148

[Peter Jeremy]

On 2012-Jul-27 17:27:45 +1000, Bruce Evans <brde at optusnet.com.au> wrote:
>On Fri, 27 Jul 2012, Peter Jeremy wrote:
>
>> On 2012-Jul-26 22:02:05 -0500, Stephen Montgomery-Smith <stephen at missouri.edu> wrote:
>>> I am not getting many responses to the programs I submitted.  I
>>> understand that you may be all very busy.
>
>I'm still working on testing and fixing clog.  Haven't got near the more
>complex functions.

I'd suggest that clog() is the most difficult complex function.  The
others mostly build from other real and complex functions.

>For clog, the worst case that I've found so far has x^2+y^2-1 ~= 1e-47:
...
>so it needs more than tripled double precision for a brute force
>evaluation, and the general case is probably worse.  I'm working
>on a rearrangement so that doubled double precision works in the
>general case.

Whilst FreeBSD includes ld128 soft-float primitives, clogl() on sparc
needs to provide a ld128 result - which implies ld240 in required for
intermediate calculations (and, even on x86, ld128 won't provide
sufficient additional precision for clogl()).

How easy is it to determine which cases need higher precision?  Can
the majority of the cases be done in double precision or does
everything need double double precision?

Would the following approach work (all variables are double):
    x1 = (double)(float)x;  /* x1 is top 24 bits of fraction */
    x2 = x - x1;            /* x2 is bottom 29 bits of fraction */
    x3 = 2.0 * x1 * x2;     /* 24b * 29b = 53b therefore this is exact */
    x1 *= x1;               /* x1 now has 48 bits */
    x2 *= x2;               /* 58 bits rounded to 53 bits */

    y1 = (double)(float)y;
    y2 = y - y1;
    y3 = 2.0 * y1 * y2;
    y1 *= y1;
    y2 *= y2;

    /* add low to high to minimise cancellation */
    result = x2 + y2;
    result += x3 + y3;
    result += x1 + y1 - 1.0;

The difficulty is that the final result summation needs to be
rearranged to if the magnitudes of x and y are sufficiently different.
For the case you presented, you need:
    result = y2;
    result += x2 + y3;
    result += x3 + y1;
    result += x1 - 1.0;
However this gives a result 1.094764425253763337e-47 without needing
any more than double precision.

>> just Appendix G.6 of WG14/N1256 turned into a C array, plus code to
>> actually run the tests & interpret the results.  So far, it's about
>> 1100 lines of which about 1/3 is the test cases and is intended to run
>
>Yikes.  My basic test program is getting too large and complex at 412
>lines.  It basically just compares with a known good or different bad
>function (with zillions of parameters to select the function and args).
>It tests exceptional cases as a side effect.

Portably supporting arm, sparc & x86 means copying structures from
fpmath.h as well as various _fpmath.h and each family needs different
code to display long double values.

>Do you do tests for fenv (i.e, that the specified exceptions and no
>others are raised)?  This might be a problem with external libraries.

I have included details of expected exceptions in the data array but
haven't written code to actually check for them.

>The can deliver values and NaNs for comparison, but nothing requires
>them to have the same fenv behaviour as the C library.

Currently, I just check that an expected NaN result is a NaN.  I don't
bother to verify that the NaN returned is tho "correct" NaN.

On 2012-Jul-27 08:21:13 -0500, Stephen Montgomery-Smith <stephen at missouri.edu> wrote:
>If I remember correctly, the C99 specification is very liberal in its 
>requirements for cpow, so that
>cpow(x,y) = cexp(y*clog(x))
>is compliant.

The main reason I skipped cpow() is that it's dyadic, which means it
needs a different test harness.  And the lack of explicit requirements
makes very difficult to determine when it should fail.

One problem with cexp(y*clog(x)) is that you automatically lose about
exponent-bits of precision.  Ideally, you need to decompose it in
much the same way as pow() - though it's messier.

-- 
Peter Jeremy
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/pgp-signature
Size: 196 bytes
Desc: not available
Url : http://lists.freebsd.org/pipermail/freebsd-numerics/attachments/20120812/36cd6e3c/attachment.pgp