C99 Complex Math Functions

[Peter Jeremy]

Hi David,

I've recently had a requirement for the C99 complex math functions.
I know you have been doing some work on implementing functions in
complex.h and wonder what (if any) plans you have for implementing
the rest of the C99 functions.

--
Peter Jeremy

[David Schultz]

On Sat, Aug 15, 2009, Peter Jeremy wrote:
> Hi David,
>
> I've recently had a requirement for the C99 complex math functions.
> I know you have been doing some work on implementing functions in
> complex.h and wonder what (if any) plans you have for implementing
> the rest of the C99 functions.

Which ones do you need?

They are on my todo list, but they're not a priority, and I'm
swamped right now, so don't hold your breath!  Steve Kargle and
Bruce Evans also have some interest in this.

Many of the functions can be fudged with some simple formulas that
you can look up.  The simple formulas work great in most cases,
but often fail near asymptotes and for special cases (infinity,
NaN, underflow, etc.)
_______________________________________________
freebsd-arch@... mailing list
http://lists.freebsd.org/mailman/listinfo/freebsd-arch
To unsubscribe, send any mail to "freebsd-arch-unsubscribe@..."

[Steve Kargl]

On Mon, Aug 17, 2009 at 04:13:19PM -0400, David Schultz wrote:

> On Sat, Aug 15, 2009, Peter Jeremy wrote:
> > Hi David,
> >
> > I've recently had a requirement for the C99 complex math functions.
> > I know you have been doing some work on implementing functions in
> > complex.h and wonder what (if any) plans you have for implementing
> > the rest of the C99 functions.
>
> Which ones do you need?
>
> They are on my todo list, but they're not a priority, and I'm
> swamped right now, so don't hold your breath!  Steve Kargle and
> Bruce Evans also have some interest in this.
>
> Many of the functions can be fudged with some simple formulas that
> you can look up.  The simple formulas work great in most cases,
> but often fail near asymptotes and for special cases (infinity,
> NaN, underflow, etc.)

In a private response to Peter, I directed him to NetBSD.  NetBSD
contains some/most(?)/all(?) of the simple formula implementations.
I sent Peter the guts of ccosh as implemented by NetBSD (the simple
formula) and the current version of ccosh that Bruce and I had written.
The netbsd version is 15 lines (excluding the license).  The bde+me
version is 91 lines (excluding license and initial comment).  To be
fair, a large portion of those 91 lines are comments of the form

        /*
         * cosh(+-Inf + I NaN)  = +Inf + I d(NaN).
         *
         * cosh(+-Inf +- I Inf) = +Inf + I dNaN.
         * The sign of Inf in the result is unspecified.  Choice = always +.
         * Raise the invalid floating-point exception.
         *
         * cosh(+-Inf + I y)   = +Inf cos(y) +- I Inf sin(y)
         */

To be more fair, most of the comments with the above form are due
to Bruce as well as the details of the implementation.

PS: yes, I'll get around to completing the commit bit stuff one of
these days.  But, for now, I need to return to a report that is
6 weeks late and 2 proposals.

--
Steve
_______________________________________________
freebsd-arch@... mailing list
http://lists.freebsd.org/mailman/listinfo/freebsd-arch
To unsubscribe, send any mail to "freebsd-arch-unsubscribe@..."

[Peter Jeremy]

On 2009-Aug-20 16:04:00 -0700, Steve Kargl <sgk@...> wrote:
>On Mon, Aug 17, 2009 at 04:13:19PM -0400, David Schultz wrote:
>> On Sat, Aug 15, 2009, Peter Jeremy wrote:
>> > I've recently had a requirement for the C99 complex math functions.
>> > I know you have been doing some work on implementing functions in
>> > complex.h and wonder what (if any) plans you have for implementing
>> > the rest of the C99 functions.
>>
>> Which ones do you need?

cacos, cacosh, casin, casinh, catan, catanh, ccos, ccosh, cexp, clog,
cpow, csin, csinh, ctan, ctanh :-)

>> Many of the functions can be fudged with some simple formulas that
>> you can look up.

I have all the simple formulae.

>>  The simple formulas work great in most cases,
>> but often fail near asymptotes and for special cases (infinity,
>> NaN, underflow, etc.)

And simplistic implementations can lose significant amounts of
precision as well.

>In a private response to Peter, I directed him to NetBSD.

Thanks for that.  I've also discovered that math/gsl also implements
the required functions (though with non-standard names and argument
types).  Since I already have gsl available, I might use it instead.
(As a warning for anyone else interested in math/gsl - it's complex
functions are implemented very simplisticly - no better than NetBSD).

--
Peter Jeremy

[Bruce Evans-4]

On Thu, 20 Aug 2009, Steve Kargl wrote:

> On Mon, Aug 17, 2009 at 04:13:19PM -0400, David Schultz wrote:
>
> In a private response to Peter, I directed him to NetBSD.  NetBSD
> contains some/most(?)/all(?) of the simple formula implementations.
> I sent Peter the guts of ccosh as implemented by NetBSD (the simple
> formula) and the current version of ccosh that Bruce and I had written.
> The netbsd version is 15 lines (excluding the license).  The bde+me
> version is 91 lines (excluding license and initial comment).  To be
> fair, a large portion of those 91 lines are comments of the form
>
> /*
> * cosh(+-Inf + I NaN)  = +Inf + I d(NaN).
> *
> * cosh(+-Inf +- I Inf) = +Inf + I dNaN.
> * The sign of Inf in the result is unspecified.  Choice = always +.
> * Raise the invalid floating-point exception.
> *
> * cosh(+-Inf + I y)   = +Inf cos(y) +- I Inf sin(y)
> */
> ...

We also have a ~15 line version for ccosh[f]() only (after removing
comments that are less detailed than the above).  Howver, the error handling
has only been checked for the float version, and the overflow handling is
mostly nonexistent for all versions.  Here is a simple version:

% /*
%  * Hyperbolic cosine of a double complex argument.
%  *
%  * Most exceptional values are automatically correctly handled by the
%  * standard formula.  See n1124.pdf for details.
                             ^^^^^^^^^

This (draft) C99+ standard gives details like the above, except of course
our choice of sign and our details of NaN handling).

%  */
%
% #include <complex.h>
% #include <math.h>
%
% #include "../math_private.h"
%
% double complex
% ccosh1(double complex z)
% {
% double x, y;
%
% x = creal(z);
% y = cimag(z);
%
% /*
% * This is subtler than it looks.  We handle the cases x == 0 and
% * y == 0 directly not for efficiency, but to avoid multiplications
% * that don't work like we need.  In these cases, the result must
% * be almost pure real or a pure imaginary, except it has type
% * float complex and its impure part may be -0.  Multiplication of
% * +-0 by an infinity or a NaN would give a NaN for the impure part,
% * and would raise an unwanted exception.
% *
% * We depend on cos(y) and sin(y) never being precisely +-0 except
% * when y is +-0 to prevent getting NaNs from other cases of
% * +-Inf*+-0.  This is true in infinite precision (since pi is
% * irrational), and checking shows that it is also true after
% * rounding to float precision.

This is impossible to check in a reasonable time for double and higher
precisions, but our implementation of accuracy for trig functions already
depends on more (on certain bounding of the results away from 0).

% */
% if (x == 0 && !isfinite(y))
% return (cpack(y - y, copysign(0, x * (y - y))));
% if (y == 0)
% return (cpack(cosh(x), isnan(x) ? copysign(0, (x + x) * y) :
%    copysign(0, x) * y));
% if (isinf(x) && !isfinite(y))
% return (cpack(x * x, x * (y - y)));
% return (cpack(cosh(x) * cos(y), sinh(x) * sin(y)));
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Standard formulas like this almost work for this and for all elementary
transcendental functions.  There are obvious accuracy and overflow
problems in such formulas.  Here cosh(x) may overflow although all of
the infinitely precise results are less than DBL_MAX.  All these problems
can be reduced by evaluating expressions in extra precision, but that
would be slow, and no extra precision is availble for the long double
versions.

% }

Bruce
_______________________________________________
freebsd-arch@... mailing list
http://lists.freebsd.org/mailman/listinfo/freebsd-arch
To unsubscribe, send any mail to "freebsd-arch-unsubscribe@..."