Use of C99 extra long double math functions after r236148
Stephen Montgomery-Smith
stephen at missouri.edu
Sun Aug 12 23:06:31 UTC 2012
On 07/17/2012 06:13 AM, Bruce Evans wrote:
> On Mon, 16 Jul 2012, Stephen Montgomery-Smith wrote:
>
>> On 07/16/2012 06:37 PM, Stephen Montgomery-Smith wrote:
>>> ...
>>> The difficulty is that catanh, cacosh, and cacosh do not have unique
>>> answers unless one makes some kind of accepted convention as to what
>>> they should be.
>>>
>>> I did find some definitions at
>>> http://publib.boulder.ibm.com/infocenter/zos/v1r11/index.jsp?topic=/com.ibm.zos.r11.bpxbd00/catan.htm.
>>>
>>> I'm going to guess that there is a typo in that the range of the
>>> imaginary part should be [0,i pi], because the usual convention is that
>>> acos of a real number is a real number between 0 and pi.
>
> I think C99 specificies the branch cus and boundary behaviour completely.
>
>>> We might get lucky, and find that the definitions of csqrt and clog in
>>> the C99 standard are already set up so that the naive formulas for
>>> cacosh, etc, just work. But whether they do or whether they don't, I
>>> think I can do it. (As a first guess, I think that catanh and casinh
>>> will work "out of the box" but cacosh is going to take a bit more work.)
>
> See below what happened for naive formulars for ccosh.
>
>>> Also casin, catan and cacos are essentially the same as casinh, etc,
>>> using formulas like sin(ix) = i sinh(x). (The hardest part is to avoid
>>> making a sign error.)
>>
>> I came up with pseudo code that looks a bit like this.
>>
>> complex casinh(complex z) {
>> double x = z.re, y = z.im;
>>
>> if (y==0)
>> return asinh(x);
>> if (x==0) {
>> if (fabs(y)<=1) return I*asin(y);
>> else return signum(y)* (
>> log(fabs(y)+sqrt(y*y-1))
>> + I*PI/2);
>> }
>> if (x>0)
>> return clog(z+csqrt(z*z+1));
>> else
>> return -clog(-z+csqrt(z*z+1));
>> }
>
> I translated this to pari. There was a sign error for the log() term
> (assuming that pari asinh() has the same branch cuts as C99):
>
> % \p 100
> % % PI = Pi;
> % clog(z) = log(z);
> % csqrt(z) = sqrt(z);
> % fabs(x) = abs(x);
> % signum(x) = sign(x);
> % % casinh(z) =
> % {
> % local(x, y);
> % % x = real(z);
> % y = imag(z);
> % if (y == 0,
> % return (asinh(x));
> % );
> % if (x == 0,
> % if (fabs(y) <= 1,
> % return (I * asin(y));
> % ,
> % return (signum(y) *
> % (-log(fabs(y) + sqrt(y * y - 1)) + I * PI / 2));
> % );
> % );
> % if (x > 0,
> % return (clog(z + csqrt(z * z + 1)));
> % ,
> % return (-clog(-z + csqrt(z * z + 1)));
> % );
> % }
> % % {
> % forstep (x = -10, 10, 0.1,
> % forstep (y = -10, 10, 0.1,
> % z = x + I*y;
> % r = casinh(z) - asinh(z);
> % if (abs(r) > 1e-30,
> % print("z = " z);
> % print("casinh(z) = " casinh(z));
> % print(" asinh(z) = " asinh(z));
> % print("diff = " r);
> % );
> % );
> % );
> % }
>
> (No differences found after fixing the sign.)
>
> Pari of course does all the calculations almost perfectly accurately (I
> told it to provide 100 decimal digits). Most multi-precision calculators
> can do the same. So pari can be used to develop the logic, but it is hard
> to use it develop accurate routines for limited precision.
>
> The most obvious immediate difficulty in translating the above into C is
> that y*y and z*z may overflow when the result shouldn't. hypot() and
> cabs() handle similar problems, with remarkably large complications.
> C99 with IEEE754 FP actually handles some aspects of overflow better
> than pari. E.g., exp(10^9) gives oveflow in pari, with no way of handling
> it AFAIK, while in C99 + IEEE754 it gives +Inf with FE_OVERFLOW.
>
> Bruce
>
>
Excellent. I think I will use pari to write the test code to check the
ULP of the result.
And I'll look into using hypot, or its logic, to compute sqrt(y*y-1).
More information about the freebsd-numerics
mailing list