Use of C99 extra long double math functions after r236148

Stephen Montgomery-Smith stephen at
Mon Jul 30 00:21:07 UTC 2012

On 07/29/2012 06:53 PM, Peter Jeremy wrote:
> [Pruning CC list to keep mailman happy]]

I had mailman complain to me in my last post.  But it only took the 
moderator a few minutes to let it through.

> On 2012-Jul-29 17:39:27 -0500, Stephen Montgomery-Smith <stephen at> wrote:
>> On 07/29/2012 05:27 PM, Peter Jeremy wrote:
>>> WG14/N1256 G.6.  I hadn't considered extending that to verifying that
>>> purely real or imaginary inputs give purely real or imaginary outputs,
>>> with the appropriately signed zero.  This might be reasonable but it's
>>> not completely trivial to implement in general since the domains of
>>> the real part can be different.
>> Maybe this should be a different program, since its logical structure
>> would be quite different.  In particular, you wouldn't be checking the
>> value of the non-zero parts.
> Adding code to skip checks on the real or imaginary part of the
> result is quite easy.
>> Also I forgot that the real part of casinh(0+I*x) isn't always 0.  If
>> |x|>1, it is something non-zero.  And so you need to check that
>> creal(casinh(0+I*x)) and creal(casinh(-0+I*x)) have opposite signs in
>> this case.
> This is related to my "domains can be different" comment.  Adding code
> to restrict the domain of the argument to be compatible with the real
> function isn't too hard (off the top of my head, I think the domains
> are all one of [-1,1], [0,Inf] or (0,Inf]).  Handling behaviour
> outside that domain requires more special-casing because the behaviour
> is less consistent.

I have a rather good handle on what that behavior will be for the 
complex arc-functions (since I have been working hard on them recently).

casinh(x+I*y) and casin(x+I*y) have positive real and imaginary parts if 
x and y are positive.  (Positive in this context includes 0, but not -0.)

In particular:
the sign of creal(casinh(z)) is the same as the sign of x;
the sign of cimag(casinh(z)) is the same as the sign of y;
the sign of creal(casin(z)) is the same as the sign of x;
the sign of cimag(casin(z)) is the same as the sign of y;

cacosh(z) and cacos(z) always have positive real part.

The imaginary part of cacos(x+I*y) has the opposite sign to x (since it 
is PI/2 - casin(x+I*y)

The imaginary part of cacosh(x+I*y) has the same sign as x.

The signs for catanh and catan are exactly the same as for casinh and casin.

For clog:

The imaginary part of clog(x+I*y) has the same sign as y.
The real part of clog will never be -0, and this doesn't have to be checked.

I am fairly sure I got these correct.  If your program starts spitting 
out huge numbers of errors, then I am wrong.  But it won't take me long 
to figure out which ones I got wrong.

>>> I'm less sure of the next logical
>>> step, which is to check things like
>>>     casinh(x + I*0) = asinh(x) + I*0
>> Does C99 mandate this?
> Nope.  They are just mathematical equivalences (at least within the
> domains supported by the real function).  POLA implies that they
> should be true but unless they are special-cased, the complex variant
> probably has less accuracy as a result of the additional calculations
> to support the imaginary component.
>>   My programs probably won't satisfy this, because
>> I realized that the computation works in these cases anyway.  Of course,
>> it would be easy to make it happen.
> It's probably up to the implementation - special casing pure real or
> imaginary arguments should give those cases a shorter and simpler (and
> therefore faster and more accurate) calculation but it's a matter of
> whether this case in common enough to justify the additional test(s)
> in all cases.
> It's also just occurred to me that doing so may result in unexpected
> output discontinuities between cfoo(x-I*tiny), cfoo(x-I*0), cfoo(x+I*0)
> and cfoo(x+I*tiny).

Since the functions are correct within a few ULP, you shouldn't see any 
discontinuities like this.  The algorithms I coded are meant to be 
rather good at handling these situations.

More information about the freebsd-numerics mailing list