Use of C99 extra long double math functions after r236148
Steve Kargl
sgk at troutmask.apl.washington.edu
Mon May 28 23:30:36 UTC 2012
On Mon, May 28, 2012 at 06:03:37PM -0500, Stephen Montgomery-Smith wrote:
> On 05/28/2012 05:17 PM, Steve Kargl wrote:
> >On Mon, May 28, 2012 at 04:19:22PM -0500, Stephen Montgomery-Smith wrote:
> >>On 05/28/2012 03:31 PM, Steve Kargl wrote:
> >>>On Mon, May 28, 2012 at 11:01:24AM -0500, Stephen Montgomery-Smith wrote:
> >>>>One thing that could be done is to have a "math/cephes" port that adds
> >>>>the extra C99 math functions. This is already done in the math/sage
> >>>>port, using a rather clever patch due to Peter Jeremy, that applies to
> >>>>the cephes code.
> >>>>
> >>>>What it would do is to create a /usr/local/lib/libm.so that would
> >>>>provide the extra functions not currently included in /lib/libm.so, and
> >>>>then link in /lib/libm.so as well. It would also create its own
> >>>>/usr/local/include/math.h and /usr/local/include/complex.h as well.
> >>>>
> >>>>What do you guys think? Do you want someone to start experimenting with
> >>>>this idea? I could do it, but probably not for a little while.
> >>>>
> >>>
> >>>This is a horrible, horrible, horrible idea. Have you
> >>>looked at the cephes code, particularly the complex.h
> >>>functions?
> >>
> >>I have only taken a very cursory look. What should I specifically look
> >>for in seeing that the code is bad?
> >
> >Well, to start with, the extra C99 math functions that
> >are missing in libm include a few for the long double type
> >and many (if not most) of the complex functions for any
> >type. Cephes at best will give you float, double, and ld80, but
> >not ld128 versions of the functions. Then, there is fun little fact
> >that neither (base) gcc nor clang nor gcc less than 4.6 does
> >complex arithematic correctly; so, one needs to jump through
> >hoops to get the correct answer (See Annex G in n1256.pdf).
> >One item of particular importance in Annex G is the behavior
> >of the functions for Re(z) and/or Im(z) being +-0, +-Inf, and
> >NaN.
> >
>
> IMHO these problems are definitely not bad enough to avoid making a
> math/cephes port. I for one would definitely like to have some kind of
> implementation of ccosh, even if it gets a few things wrong when it is
> fed Nan or I*Inf or such like as its input. I mean, if clang or gcc46
> doesn't even get this right, how can we expect better from libm?
Search the llvm and gcc bug databases, I'm the person that informed
both that they can't deal with complex arithematic correctly.
Interesting that you mention ccosh.
>From clog.c in http://www.netlib.org/cephes/cmath.tgz
void
ccosh (z, w)
cmplx *z, *w;
{
double x, y;
x = z->r;
y = z->i;
w->r = cosh (x) * cos (y);
w->i = sinh (x) * sin (y);
}
>From clog.c in http://www.netlib.org/cephes/c9x-complex
double complex
ccosh (z)
double complex z;
{
double complex w;
double x, y;
x = creal(z);
y = cimag(z);
w = cosh (x) * cos (y) + (sinh (x) * sin (y)) * I;
return (w);
}
See math_private.h about the above.
And, finally,
:
/*-
* Copyright (c) 2005 Bruce D. Evans and Steven G. Kargl
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* 1. Redistributions of source code must retain the above copyright
* notice unmodified, this list of conditions, and the following
* disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
* IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
* OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
* IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
* NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
* DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
* THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
* THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
/*
* Hyperbolic cosine of a complex argument z = x + i y.
*
* cosh(z) = cosh(x+iy)
* = cosh(x) cos(y) + i sinh(x) sin(y).
*
* Exceptional values are noted in the comments within the source code.
* These values and the return value were taken from n1124.pdf.
*/
#include <sys/cdefs.h>
__FBSDID("$FreeBSD: head/lib/msun/src/s_ccosh.c 226599 2011-10-21 06:29:32Z das $");
#include <complex.h>
#include <math.h>
#include "math_private.h"
static const double huge = 0x1p1023;
double complex
ccosh(double complex z)
{
double x, y, h;
int32_t hx, hy, ix, iy, lx, ly;
x = creal(z);
y = cimag(z);
EXTRACT_WORDS(hx, lx, x);
EXTRACT_WORDS(hy, ly, y);
ix = 0x7fffffff & hx;
iy = 0x7fffffff & hy;
/* Handle the nearly-non-exceptional cases where x and y are finite. */
if (ix < 0x7ff00000 && iy < 0x7ff00000) {
if ((iy | ly) == 0)
return (cpack(cosh(x), x * y));
if (ix < 0x40360000) /* small x: normal case */
return (cpack(cosh(x) * cos(y), sinh(x) * sin(y)));
/* |x| >= 22, so cosh(x) ~= exp(|x|) */
if (ix < 0x40862e42) {
/* x < 710: exp(|x|) won't overflow */
h = exp(fabs(x)) * 0.5;
return (cpack(h * cos(y), copysign(h, x) * sin(y)));
} else if (ix < 0x4096bbaa) {
/* x < 1455: scale to avoid overflow */
z = __ldexp_cexp(cpack(fabs(x), y), -1);
return (cpack(creal(z), cimag(z) * copysign(1, x)));
} else {
/* x >= 1455: the result always overflows */
h = huge * x;
return (cpack(h * h * cos(y), h * sin(y)));
}
}
/*
* 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) = d(NaN) + I sign(d(+-0, NaN))0.
* The sign of 0 in the result is unspecified. Choice = normally
* the same as d(NaN).
*/
if ((ix | lx) == 0 && iy >= 0x7ff00000)
return (cpack(y - y, copysign(0, x * (y - y))));
/*
* cosh(+-Inf +- I 0) = +Inf + I (+-)(+-)0.
*
* cosh(NaN +- I 0) = d(NaN) + I sign(d(NaN, +-0))0.
* The sign of 0 in the result is unspecified.
*/
if ((iy | ly) == 0 && ix >= 0x7ff00000) {
if (((hx & 0xfffff) | lx) == 0)
return (cpack(x * x, copysign(0, x) * y));
return (cpack(x * x, copysign(0, (x + x) * y)));
}
/*
* cosh(x +- I Inf) = dNaN + I dNaN.
* 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 = don't raise (except for signaling NaNs).
*/
if (ix < 0x7ff00000 && iy >= 0x7ff00000)
return (cpack(y - y, x * (y - y)));
/*
* 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) {
if (iy >= 0x7ff00000)
return (cpack(x * x, x * (y - y)));
return (cpack((x * x) * cos(y), x * sin(y)));
}
/*
* cosh(NaN + I NaN) = d(NaN) + I d(NaN).
*
* cosh(NaN +- I Inf) = d(NaN) + I d(NaN).
* Optionally raises the invalid floating-point exception.
* Choice = raise.
*
* cosh(NaN + I y) = d(NaN) + I d(NaN).
* Optionally raises the invalid floating-point exception for finite
* nonzero y. Choice = don't raise (except for signaling NaNs).
*/
return (cpack((x * x) * (y - y), (x + x) * (y - y)));
}
And, one gets ccos for free.
double complex
ccos(double complex z)
{
/* ccos(z) = ccosh(I * z) */
return (ccosh(cpack(-cimag(z), creal(z))));
}
> Also, a really nice thing about Jeremy's patch is that it automatically
> detects which functions are already included in the base libm, and only
> adds functions not already in libm.
>
> And ld80 is better than nothing if ld128 isn't available.
>
> I would avoid cephes only if it got some of the answers horribly wrong
> (as in erf(100) being something like -14232 as the openoffice
> implementation of the erf function used to report).
Who's writing the code to test the implementations? That is
better much the problem. Without testing, one might get an
implementation that appears to work until it doesn't! It took
me 3+ years to get sqrtl() into libm, but bde and das (and
myself) wanted to make sure the code worked.
> By the way, do you have an opinion on the complex functions used in
> Linux? (I tried reading the glibc code, but I could only find chains of
> macros and functions calling each other, and I could never find where
> the actual work was performed.)
I haven't looked at glibc code in years, because I hack on libm
when I can. I do not want to run into questions about whether
my code is tainted by the gpl.
--
Steve
More information about the freebsd-current
mailing list