Date: Mon, 28 May 2012 16:30:35 -0700 From: Steve Kargl <sgk@troutmask.apl.washington.edu> To: Stephen Montgomery-Smith <stephen@missouri.edu> Cc: freebsd-current@freebsd.org Subject: Re: Use of C99 extra long double math functions after r236148 Message-ID: <20120528233035.GA77157@troutmask.apl.washington.edu> In-Reply-To: <4FC40449.3040602@missouri.edu> References: <4FC30090.4070003@gwdg.de> <4D8CF7D2-CBEE-438E-A9E7-9C47A8892622@FreeBSD.org> <4FC36FE1.9080908@gwdg.de> <DFE44442-A1B1-40B3-9A79-874BD2BE45E7@FreeBSD.org> <4FC38B81.6000302@gwdg.de> <4FC3A154.8030702@missouri.edu> <20120528203159.GA76340@troutmask.apl.washington.edu> <4FC3EBDA.2080502@missouri.edu> <20120528221731.GA76723@troutmask.apl.washington.edu> <4FC40449.3040602@missouri.edu>
next in thread | previous in thread | raw e-mail | index | archive | help
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
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20120528233035.GA77157>