Date: Wed, 12 Oct 2005 17:23:28 +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: <20051012160109.I73531@delplex.bde.org> In-Reply-To: <20051010185153.GA55589@troutmask.apl.washington.edu> References: <20050930152734.GA20655@troutmask.apl.washington.edu> <20051001204005.N37704@delplex.bde.org> <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>
next in thread | previous in thread | raw e-mail | index | archive | help
On Mon, 10 Oct 2005, Steve Kargl wrote:
> (Attributions are a mess :-)
Now accumulating them again :-).
> I fixed this in my version of ccosh and used cpack(,) for
> constructing the complex values. I tried to use cpack(,)
> in your simplified version, but screwed something up in
> that the return value is the complex conjugate from that
> version.
>
> I've attached the latest version. Hopefully, I caught the
> rest of the style(9) problems. The cosh.3 man page has been
No :-).
> updated. I did not hook s_ccosh.c up to the Makefile
> because I wasn't sure were to put it in msun/src. If you have
> no further comments, can you commit this version (or your
> simplified version)? I'll look at the other trigometric and
> hyperbolic functions when we converge on s_ccosh.c.
I'm beginning to prefer the unsimplified version, but it needs some
more editing after the following changes. The simplified version
needed another case for ccosh(+-Inf + I NaN) and that already gives
it about half as many cases and the other version. It saves code
mainly by falling through to the code that works for finite args, but
eventually the case of finite args will need to be more special to
avoid overflows and then the non-finite cases won't be able to just
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.
Patch relative to your last version:
% --- s_ccosh.c~ Wed Oct 12 05:11:13 2005
% +++ s_ccosh.c Wed Oct 12 09:05:44 2005
% @@ -43,5 +43,5 @@
% ccosh(double complex z)
% {
% - double con, x, y;
% + double x, y;
% int32_t hx, hy, ix, iy, lx, ly;
%
% @@ -55,25 +55,10 @@
% iy = 0x7fffffff & hy; /* If iy >= 0x7ff00000, then inf or NaN. */
%
% - /* Even function: ccosh(-z) = ccosh(z). */
% - if (hx & 0x80000000 && !isnan(x)) {
% - x = -x;
% - hx = ix;
% - if (!isnan(y)) {
% - y = -y;
% - hy ^= 0x80000000;
% - }
% - }
% -
% - /* Required identity: ccosh(conj(z)) = conj(ccosh(z)). */
% - con = 1;
% - if (hy & 0x80000000 && !isnan(y)) {
% - con = -1;
% - y = -y;
% - hy = iy;
% - }
% -
Delete the explicit sign flipping since it mostly just complicates things.
% /* Handle the nearly-non-exceptional cases where x and y are finite. */
% - if (ix < 0x7ff00000 && iy < 0x7ff00000)
% - return (cpack(cosh(x) * cos(y), con * sinh(x) * sin(y)));
% + if (ix < 0x7ff00000 && iy < 0x7ff00000) {
% + if ((iy | ly) == 0)
% + return (cpack(cosh(x), x * y));
This special case was lost. Without this, cosh(x) * +-0 gives NaN if
cosh(x) is +-Inf.
x * y for the imaginary part is an efficient way to get +-0 with the
right sign.
% + return (cpack(cosh(x) * cos(y), sinh(x) * sin(y)));
% + }
%
% /*
% @@ -81,16 +66,7 @@
% * invalid exception.
% * cosh(+0 + I NaN) = NaN +- I 0, sign of 0 is unspecified.
% - *
% - * Use the sign of (sinh(original_x) * sin(original_y)) for
% - * the optional sign. This expression reduces to (y - y) for
% - * all cases checked (mainly i386's). The fdlibm sin(y) is
% - * certainly (y - y), but in the +-Inf cases we may have
% - * flipped the sign of original_y to get y, and we want to
% - * multiply by the sinh() term. These complications have no
% - * effect with normal NaN semantics (-Inf - -Inf) = same NaN
% - * as (Inf - Inf), and (finite * NaN) = same NaN).
% */
Delete this since it became out of date (x and y haven't changed) and
getting the sign right has fine details not just here. I don't want
comments in the code about all of them.
% if ((ix | lx) == 0 && iy >= 0x7ff00000)
% - return (cpack(y - y, copysign(0, y - y)));
% + return (cpack(y - y, copysign(0, x * (y - y))));
%
% /*
This change is a no-op on all machines that I know of. It is to get
as close as possible to what the hardware would do for sinh(x) *
sinh(y). x is +-0 so sinh(x) is the same as x. sinh(y) is (y - y)
(a NaN). We can't just muliply these since C99 specifies a result of
+-0, so we get as close as possible by taking the sign of the multiple.
Note that there is no additional invalid exception for this -- there
is already one for the real part. (Actually, there probably is a
second one since evaluating (y - y) only once would be an invalid
optimization, but these exceptions coalesce.) Multiplying by x has
no effect on all machines that I know of. Normally, if y is Inf then
(y - y) is some NaN and multiplying by x preserves this, and if y is
NaN then (y - y) is the same NaN and that is preseved.
% @@ -99,5 +75,5 @@
% */
% if (ix < 0x7ff00000 && iy >= 0x7ff00000)
% - return (cpack(y - y, y - y));
% + return (cpack(y - y, x * (y - y)));
% /*
% * cosh(+Inf + I 0) = +Inf + I 0
Similarly, except here we want a NaN imaginary part so we don't need a
copysign(). Cases like this can currently all use the formula.
% @@ -108,8 +84,8 @@
% if (ix >= 0x7ff00000 && ((hx & 0xfffff) | lx) == 0) {
% if ((iy | ly) == 0)
% - return (cpack(x * x, con * y));
% - else if (iy >= 0x7ff00000)
% - return (cpack(x * x, y - y));
% - return (cpack((x * x) * cos(y), (x * x) * con * sin(y)));
% + return (cpack(x * x, copysign(0, x) * y));
This is the only place that deleting the sign flipping made more complicated
later. The sign that was in "con" now needs to be determined.
% + if (iy >= 0x7ff00000)
Minor style fix (don't use "return; else").
% + return (cpack(x * x, x * (y - y)));
Let the sign of x possibly affect the NaN (y - y) as usual.
% + return (cpack((x * x) * cos(y), x * sin(y)));
"con" just goes away (it is in y), provided we use the right sign for the
x-related term here. (x * x) was logically wrong -- x is +-Inf and the
term is sinh(x) = +-Inf = x. (x * x) only worked because we normalized x
so it was +Inf.
% }
% /*
% @@ -119,5 +95,5 @@
% */
% if ((iy | ly) == 0)
% - return (cpack(x - x, copysign(0, x - x)));
% - return (cpack((x - x) + (y - y), (x - x) + (y - y)));
% + return (cpack(x * x, copysign(0, (x + x) * y)));
Now x is NaN. (x - x) was logically wrong for cosh(x) -- cosh(x) is
(x * x) for both Infs and NaNs -- it's a square so that it works for
+-Inf, and this gives a choice of sign and value for NaNs too, depending
on what the hardware does for NaN * NaN. Just use the same value here.
Similarly for the imaginary part:
- sinh(x) is (x + x), not (x - x), for Infs and NaNs, so as to preserve
the sign for Infs. Just use the same value.
- let the sign of y possibly affect the result like the sign of x did
in other cases.
% + return (cpack((x * x) * (y - y), (x + x) * (y - y)));
% }
Similarly. We want to let the sign of sin(y) possibly affect the real
part. sin(y) is (y - y), so just use that. For the imaginary part,
multiply instead of add to do the right thing with signs; just use
sinh(x) = (x + x) not (x - x); sin(y) = (y - y) was already correct.
% 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.
% diff -ruN /usr/src/lib/msun/src/complex/s_ccosh.c freebsd/src/lib/msun/src/complex/s_ccosh.c
% --- /usr/src/lib/msun/src/complex/s_ccosh.c Wed Dec 31 16:00:00 1969
% +++ freebsd/src/lib/msun/src/complex/s_ccosh.c Mon Oct 10 11:10:14 2005
Patch was relative to this.
% +
% +/*
% + * Hyperbolic cosine of a complex argument z = x + i y such that i = sqrt(-1).
Other changes: don't define "i" here...
% + ...
% + ix = 0x7fffffff & hx; /* If ix >= 0x7ff00000, then inf or NaN. */
% + iy = 0x7fffffff & hy; /* If iy >= 0x7ff00000, then inf or NaN. */
% +
% + /* Even function: ccosh(-z) = ccosh(z). */
% + if (hx & 0x80000000 && !isnan(x)) {
% + x = -x;
Still has lots of strange indents -- corrupt tab above, 4-chars only for
most second indents.
% + /*
% + * 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.
% + * cosh(+Inf + I y) = +Inf [cos(y) + I sin(y)], y != 0 and finite.
Still have the strange brackets here.
% diff -ruN /usr/src/lib/msun/src/complex/s_ccoshf.c freebsd/src/lib/msun/src/complex/s_ccoshf.c
% --- /usr/src/lib/msun/src/complex/s_ccoshf.c Wed Dec 31 16:00:00 1969
% +++ freebsd/src/lib/msun/src/complex/s_ccoshf.c Fri Oct 7 17:30:36 2005
% ...
% +float complex
% +ccoshf(float complex z)
% +{
% + return ((float complex) ccosh((double complex) z));
% +}
I prefer to use implicit conversions.
% 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)
%
% /*
% + * Inline functions that can be used to construct complex values.
% + */
I moved my comment about the gcc bug -- the reason for existence of these
functions -- to here.
% +static __inline float complex
% +cpackf(float __x, float __y)
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".
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.
Otherwise good. These interfaces seem to work well.
Here is my current "simpler" version. I plan to merge it into the
other version when the other version's indentationis fixed.
%%%
/*
* Hyperbolic cosine of a double complex argument.
*
* Most exceptional values are automatically correctly handled by the
* standard formula. See n1124.pdf for details.
*/
#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.
*/
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)));
}
Here is my program for comparing complex functions. It has been cleaned
up a bit (it now uses fixed-width types but not GET_FLOAT_WORD() etc.
to access them properly). Compile with -DTESTD for simplistic double
precision comparison. Oops, it needs to use cpack*().
%%%
#include <complex.h>
#include <math.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#ifndef FUNC
#define FUNC ccosh
#endif
double complex ccosh(double complex z);
double complex ccosh1(double complex z);
float complex ccoshf(float complex z);
float complex ccosh1f(float complex z);
int
main(void)
{
float complex vf0, vf1;
float complex z;
volatile float vf0i, vf0r, vf1i, vf1r;
float x, y;
uintmax_t na, tot_er;
uint32_t eri, err, max_er, stride, vf0iu, vf0ru, vf1iu, vf1ru, xu, yu;
max_er = 0;
na = 0;
stride = 0x80000;
tot_er = 0;
for (xu = 0;; xu += stride) {
for (yu = 0;; yu += stride) {
x = *(float *)&xu;
y = *(float *)&yu;
/* `z = x + I * y;', unbroken: */
crealf(z) = x;
cimagf(z) = y;
#ifdef TESTD
/*
* Hack for double precision functions. Testing them only
* on float args is OK, but we should report errors in
* double precision for them.
*/
vf0 = __CONCAT(FUNC, )(z);
vf1 = __CONCAT(FUNC, 1)(z);
#else
vf0 = __CONCAT(FUNC, f)(z);
vf1 = __CONCAT(FUNC, 1f)(z);
#endif
vf0r = crealf(vf0);
vf0ru = *(volatile uint32_t *)&vf0r;
vf0i = cimagf(vf0);
vf0iu = *(volatile uint32_t *)&vf0i;
vf1r = crealf(vf1);
vf1ru = *(volatile uint32_t *)&vf1r;
vf1i = cimagf(vf1);
vf1iu = *(volatile uint32_t *)&vf1i;
na++;
err = abs(vf1ru - vf0ru);
eri = abs(vf1iu - vf0iu);
tot_er += err + eri;
if (max_er < err)
max_er = err;
if (max_er < eri)
max_er = eri;
if (err > 4 || eri > 4) {
printf("z = %#10x+I*%-#10x %.15g+I*%-.15g\n",
xu, yu, x, y);
printf(
__XSTRING(FUNC) "f(z) = %#10x+I*%-#10x %.15g+I*%-.15g\n",
vf0ru, vf0iu, vf0r, vf0i);
printf(
__XSTRING(FUNC) "1f(z) = %#10x+I*%-#10x %.15g+I*%-.15g\n",
vf1ru, vf1iu, vf1r, vf1i);
printf("err = %s%#x+I*%s.%#x\n",
vf0ru <= vf1ru ? "+" : "-", err,
vf0iu <= vf1iu ? "+" : "-", eri);
fflush(stdout);
}
if (na % 10000000 == 1) {
printf(
"%ju: %#x+I*%#x: %g+I*%g: max_error = %#x, avg_error = %.3f\n",
na, xu, yu, x, y, max_er, (double)tot_er / na);
fflush(stdout);
}
if (yu == -stride)
break;
}
if (xu == -stride)
break;
}
printf("%ju cases: max_error = %#x, avg_error = %.3f\n",
na, max_er, (double)tot_er / na);
}
#include "s_ccosh.c"
#include "s_ccosh1.c"
#include "s_ccoshf.c"
#include "s_ccosh1f.c"
/* For ccosh(), must avoid using the broken hardware cos() and sin(): */
#undef rcsid
#define rcsid rcsid6
#include "../s_cos.c"
#undef rcsid
#define rcsid rcsid7
#include "../s_sin.c"
#ifdef DEBUG
#undef rcsid
#define rcsid rcsid0
#include "../e_coshf.c"
#undef one
#define one one0
#undef rcsid
#define rcsid rcsid1
#include "../e_sinhf.c"
#undef one
#define one one1
#undef rcsid
#define rcsid rcsid2
#include "../s_tanhf.c"
#endif
%%%
Bruce
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20051012160109.I73531>