Date: Fri, 21 Oct 2005 20:50:53 +1000 (EST) From: Bruce Evans <bde@zeta.org.au> To: Steve Kargl <sgk@troutmask.apl.washington.edu> Cc: freebsd-standards@freebsd.org Subject: Re: complex.h math functions Message-ID: <20051021194624.H598@epsplex.bde.org> In-Reply-To: <20051016184129.GA24651@troutmask.apl.washington.edu> References: <20051002191841.GA40568@troutmask.apl.washington.edu> <20051004164618.U47262@delplex.bde.org> <20051005032400.GA6736@troutmask.apl.washington.edu> <20051005191936.N51036@delplex.bde.org> <20051006212602.GA40609@troutmask.apl.washington.edu> <20051007231439.F58005@delplex.bde.org> <20051007190239.GA78674@troutmask.apl.washington.edu> <20051008052850.S59139@delplex.bde.org> <20051010185153.GA55589@troutmask.apl.washington.edu> <20051012160109.I73531@delplex.bde.org> <20051016184129.GA24651@troutmask.apl.washington.edu>
next in thread | previous in thread | raw e-mail | index | archive | help
On Sun, 16 Oct 2005, Steve Kargl wrote:
> On Wed, Oct 12, 2005 at 05:23:28PM +1000, Bruce Evans wrote:
> Hopefully, this new version is closer to KNF.
It still has only 4 chars for all secondary indents
>> use it. Listing all the special cases also serves as documentation.
>> I think the current order of special cases is not quite the best,
>> however.
>
> I don't follow why you think the special cases need to be re-order.
> Without knowing how the user base will abuse ccosh, there should
> be no prefered order.
The current order is almost that in the standard. This order is not
bad but it might not give the simplest or shortest classification.
In particular, I think it might be better to group by 0's before
grouping by Infs and NaNs.
> This version incorporates your changes, fixes
> a few more whitespace problems, and expands the
> comments on the exceptional cases.
I modified to to try to document all the sign combinations, and fixed a
couple of bugs:
% --- s_ccosh.c~ Mon Oct 17 12:11:18 2005
% +++ s_ccosh.c Tue Oct 18 18:11:22 2005
% @@ -63,9 +63,17 @@
%
% /*
% - * cosh(+0 + I Inf) = NaN +- I 0. The sign of 0 in the result
% - * is unspecified. Raise the invalid floating-point exception.
% + * cosh(+-0 +- I Inf) = dNaN + I sign(d(+-0, dNaN))0.
% + * The sign of 0 in the result is unspecified. Choice = normally
% + * the same as dNaN. Raise the invalid floating-point exception.
% *
% - * cosh(+0 + I NaN) = NaN +- I 0. The sign of 0 in the result
% - * is unspecified.
% + * cosh(+-0 +- I NaN) = d(NaN) + I sign(d(+-0, NaN))0.
% + * The sign of 0 in the result is unspecified. Choice = normally
% + * the same as d(NaN).
% + *
% + * Say something about our sign choices generally but not specifically?
% + * The above is specific; dNaN = default NaN and d(NaN) = input NaN
% + * with default transformation under operations (e.g., on i386's,
% + * convert signaling NaNs to quiet NaNs by setting the quiet bit).
% + * XXX indentation of all of these.
% */
% if ((ix | lx) == 0 && iy >= 0x7ff00000)
This also says something about NaN not always being the same value. I'm
not sure where this should be documented.
Ther is another indentation problem in these comments. 2 spaces after the
"*" is unusual, and is not useful since the same indentation is used for
the formulas and the descriptions. I only changed the formatting to start
a new line after the formulas.
% @@ -73,9 +81,11 @@
%
% /*
% - * cosh(x + I Inf) = NaN + I NaN. Raise the invalid
% + * cosh(x +- I Inf) = dNaN + I dNaN.
% + * Raise the invalid
% * floating-point exception for finite nonzero x.
% *
% - * cosh(x + I NaN) = NaN + I NaN. Optionally raise the invalid
% - * floating-point exception for finite nonzero x.
% + * cosh(x + I NaN) = d(NaN) + I d(NaN).
% + * Optionally raises the invalid
% + * floating-point exception for finite nonzero x. Choice = raise.
% */
% if (ix < 0x7ff00000 && iy >= 0x7ff00000)
An earlier version tried to document the sign of the resulting NaN. I
decided that this was too much detail. But the detail that the sign of
0's is copied from the sign of a NaN is documented for other cases.
This also improves the wording for options.
% @@ -83,10 +93,13 @@
%
% /*
% - * cosh(+Inf + I 0) = +Inf + I 0.
% - * cosh(+Inf + I NaN) = +Inf + I NaN.
% - * cosh(+Inf + I y) = +Inf cis(y) for finite nonzero y.
% - * cis(y) = cos(y) + I sin(y).
% - * cosh(+Inf + I Inf) = +Inf + I NaN. The sign of Inf in the result
% - * is unspecified. Raise the invalid floating-point exception.
% + * cosh(+-Inf +- I 0) = +Inf + I (+-)(+-)0.
% + *
% + * 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)
% */
% if (ix >= 0x7ff00000 && ((hx & 0xfffff) | lx) == 0) {
This also adds newlines and drops cis(). cis() doesn't work in the
precence of signs.
% @@ -99,11 +112,14 @@
%
% /*
% - * cosh(NaN + I 0) = NaN +- I 0. The sign of 0 in the result
% + * cosh(NaN +- I 0) = d(NaN) + I sign(d(NaN, +-0))0.
% + * The sign of 0 in the result
% * is unspecified.
% *
% - * cosh(NaN + I y) = NaN + I NaN. Raise the invalid
% - * floating-point exception for finite nonzero y.
% + * cosh(NaN + I y) = d(NaN) + I d(NaN).
% + * Optionally raises the invalid
% + * floating-point exception for finite or infinite nonzero y.
% + * Choice = raise for infinite y only.
% *
% - * cosh(NaN + I NaN) = NaN + I NaN.
% + * cosh(NaN + I NaN) = d(NaN) + I d(NaN).
% */
% if ((iy | ly) == 0)
This also fixes the description of the exception for cos(NaN + I y). I
checked all your decriptions against n1124.pdf and got this fix from there.
Note that the case cosh(NaN +- I Inf) is together with cosh(NaN + I y)
for finite y since finite and infinite y's act the same here because the
NaN dominates. Hmm, it seems wrong to even optionally raise an exception
for finite y.
>
>> % diff -ruN /usr/src/lib/msun/man/cosh.3 freebsd/src/lib/msun/man/cosh.3
>> % --- /usr/src/lib/msun/man/cosh.3 Fri Jan 14 15:28:28 2005
>> % +++ freebsd/src/lib/msun/man/cosh.3 Sat Oct 8 11:46:29 2005
>> % ...
>>
>> OK; could be more detailed.
>
> Do you want a section that describes the range of inputs that
> return a non-infinite results? Do you want a section that
> describes the return value and the exceptional cases? I'll
> have to do this later because I'll need to learning some
> additional mdoc features.
Exceptional cases are mostly not described at all for the real functions,
so I wouldn't describe them better for the complex functions. I was
thinking mainly of doumenting the intention that the error is < 2 ulps
and noting cases where this isn't implemented in the BUGS section.
>> % + * cosh(x + I Inf) = NaN + I NaN, finite x != 0, invalid exception.
>> % + * cosh(x + I NaN) = NaN + I NaN, finite x != 0, opt. inval. except.
>> % + */
>>
>> Not quite enough space to describe it all on 1 line; the inputs and outputs
>> get mixed cryptically; might use more formal abbreviations.
>
> I've expanded all comments to essentially the language in n1124.pdf.
I'd still like to find a more concise version. Have a look at e_pow.c.
> I've updated [c_coshf.c] also. Note, the forward declaration of
> ccosh() is not needed because complex.h should contain a
> prototype for it.
My version with these and some other cosmetic changes. Everything except
the copyright:
% /*
% * Hyperbolic cosine of a float complex argument.
% */
%
% #include <complex.h>
%
% float complex
% ccoshf(float complex z)
% {
%
% return (ccosh(z));
% }
>> % diff -ruN /usr/src/lib/msun/src/math_private.h
>> freebsd/src/lib/msun/src/math_private.h
>> % --- /usr/src/lib/msun/src/math_private.h Fri Feb 4 12:05:39 2005
>> % +++ freebsd/src/lib/msun/src/math_private.h Fri Oct 7 17:30:40 2005
>> % @@ -155,6 +155,36 @@
>> % } while (0)
>> %
>>
>> I needed a namespace hack to make this compile. Clients that don't use
>> this shouldn't need to include <complex.h>, and this file shouldn't
>> include it. I used "#ifdef I".
>
> A single letter #ifdef makes me nervous, and it isn't very descriptive.
> I would suggest using #ifdef _COMPLEX_H. OTOH, math_private.h isn't
> an installed header file, so only the person or persons writing the
> complex math functions need to be bothered by this.
I actually used `complex'. _COMPLEX_H is much better. The only problem
with it is to misspell it consistently (it should be spelled _COMPLEX_H_;
tail -1 *.h in /usr/include shows similar misspellings and other style
bugs like backwards comments in about half of the headers).
>> I think the underscores shouldn't be used here (not even for __inline).
>> This is not a public interface so we don't need to be very careful with
>> the namespace.
>
> I removed the underscores in my local tree. There, of course, were no problems.
>
>> Otherwise good. These interfaces seem to work well.
>
> You did not provide a patch to your local math_private.h.
Here it is, after fixing the ifdef and removing most undercores:
% Index: math_private.h
% ===================================================================
% RCS file: /home/ncvs/src/lib/msun/src/math_private.h,v
% retrieving revision 1.17
% diff -u -2 -r1.17 math_private.h
% --- math_private.h 4 Feb 2005 20:05:39 -0000 1.17
% +++ math_private.h 21 Oct 2005 10:43:22 -0000
% @@ -155,4 +162,45 @@
% } while (0)
%
% +#ifdef _COMPLEX_H
% +/*
% + * Inline functions that can be used to construct complex values.
% + *
% + * The C99 standard intends x+I*y to work for this, but x+I*y is
% + * currently unusable in general since gcc introduces many overflow,
% + * underflow, sign and efficiency bugs by rewriting I*y as (0+I)*(y+0*I)
% + * and laboriously computing the full complex product. In particular,
% + * I*Inf is corrupted to NaN+I*Inf, and I*-0 is corrupted to 0+I*0.
% + */
% +static __inline float complex
% +cpackf(float x, float y)
% +{
% + float complex z;
% +
% + __real__ z = x;
% + __imag__ z = y;
% + return (z);
% +}
% +
% +static __inline double complex
% +cpack(double x, double y)
% +{
% + double complex z;
% +
% + __real__ z = x;
% + __imag__ z = y;
% + return (z);
% +}
% +
% +static __inline long double complex
% +cpackl(long double x, long double y)
% +{
% + long double complex z;
% +
% + __real__ z = x;
% + __imag__ z = y;
% + return (z);
% +}
% +#endif /* _COMPLEX_H */
% +
% /*
% * ieee style elementary functions
> BTW, the exceptional cases for ccos(z) are defined in n1124.pdf
> in terms of ccos(z) = ccosh(I z). We can do
>
> double complex
> ccos(double complex z)
> {
> return (ccosh(-cimag(z),creal(z));
> }
>
> or if you prefer we can copy s_ccosh.c and change
>
> x = creal(z);
> y = cimag(z);
> to
> x = -cimag(z);
> y = creal(x);
Use the version with the extra function call for now.
Bruce
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20051021194624.H598>