Index: include/complex.h =================================================================== RCS file: /home/ncvs/src/include/complex.h,v retrieving revision 1.6 diff -u -p -r1.6 complex.h --- include/complex.h 14 Aug 2004 18:03:21 -0000 1.6 +++ include/complex.h 8 Aug 2007 00:56:41 -0000 @@ -45,6 +45,7 @@ __BEGIN_DECLS double cabs(double complex); float cabsf(float complex); +double complex ccosh(double complex); double cimag(double complex); float cimagf(float complex); long double cimagl(long double complex); Index: lib/msun/Makefile =================================================================== RCS file: /home/ncvs/src/lib/msun/Makefile,v retrieving revision 1.78 diff -u -p -r1.78 Makefile --- lib/msun/Makefile 21 May 2007 02:49:08 -0000 1.78 +++ lib/msun/Makefile 8 Aug 2007 00:56:59 -0000 @@ -68,13 +68,18 @@ SYMBOL_MAPS= ${SYM_MAPS} # C99 long double functions COMMON_SRCS+= s_copysignl.c s_fabsl.c s_modfl.c +.if ${MACHINE_ARCH} == "i386" || ${MACHINE_ARCH} == "amd64" +COMMON_SRCS+= e_sqrtl.c s_cbrtl.c +.endif + .if ${LDBL_PREC} != 53 # If long double != double use these; otherwise, we alias the double versions. COMMON_SRCS+= s_fmal.c s_frexpl.c s_nextafterl.c s_nexttoward.c s_scalbnl.c .endif # C99 complex functions -COMMON_SRCS+= s_cimag.c s_cimagf.c s_cimagl.c s_conj.c s_conjf.c s_conjl.c \ +COMMON_SRCS+= s_ccosh.c \ + s_cimag.c s_cimagf.c s_cimagl.c s_conj.c s_conjf.c s_conjl.c \ s_creal.c s_crealf.c s_creall.c # FreeBSD's C library supplies these functions: Index: lib/msun/Symbol.map =================================================================== RCS file: /home/ncvs/src/lib/msun/Symbol.map,v retrieving revision 1.4 diff -u -p -r1.4 Symbol.map --- lib/msun/Symbol.map 29 Apr 2007 14:05:20 -0000 1.4 +++ lib/msun/Symbol.map 8 Aug 2007 00:56:59 -0000 @@ -56,12 +56,15 @@ FBSD_1.0 { sinhf; sqrt; sqrtf; + sqrtl; asinh; asinhf; atan; atanf; cbrt; cbrtf; + cbrtl; + ccosh; ceil; ceilf; ceill; Index: lib/msun/man/sqrt.3 =================================================================== RCS file: /home/ncvs/src/lib/msun/man/sqrt.3,v retrieving revision 1.13 diff -u -p -r1.13 sqrt.3 --- lib/msun/man/sqrt.3 9 Jan 2007 01:02:06 -0000 1.13 +++ lib/msun/man/sqrt.3 8 Aug 2007 00:56:59 -0000 @@ -28,14 +28,16 @@ .\" from: @(#)sqrt.3 6.4 (Berkeley) 5/6/91 .\" $FreeBSD: src/lib/msun/man/sqrt.3,v 1.13 2007/01/09 01:02:06 imp Exp $ .\" -.Dd May 6, 1991 +.Dd August 5, 2007 .Dt SQRT 3 .Os .Sh NAME .Nm cbrt , .Nm cbrtf , +.Nm cbrtl , .Nm sqrt , -.Nm sqrtf +.Nm sqrtf , +.Nm sqrtl .Nd cube root and square root functions .Sh LIBRARY .Lb libm @@ -45,35 +47,43 @@ .Fn cbrt "double x" .Ft float .Fn cbrtf "float x" +.Ft "long double" +.Fn cbrtl "long double x" .Ft double .Fn sqrt "double x" .Ft float .Fn sqrtf "float x" +.Ft "long double" +.Fn sqrtl "long double x" .Sh DESCRIPTION The -.Fn cbrt -and the -.Fn cbrtf -functions compute -the cube root of +.Fn cbrt , +.Fn cbrtf , +and +.Fn cbrtl +functions compute the cube root of .Ar x . .Pp The .Fn sqrt -and the .Fn sqrtf -functions compute the -non-negative square root of x. +and +.Fn sqrtl +functions compute the non-negative square root of x except +.Fn sqrt "-0" += -0. .Sh RETURN VALUES The -.Fn cbrt +.Fn cbrt , +.Fn cbrtf , and the -.Fn cbrtf +.Fn cbrtl functions return the requested cube root. The -.Fn sqrt +.Fn sqrt , +.Fn sqrtf , and the -.Fn sqrtf +.Fn sqrtl functions return the requested square root unless an error occurs. An attempt to take the Index: lib/msun/src/e_sqrtl.c =================================================================== RCS file: lib/msun/src/e_sqrtl.c diff -N lib/msun/src/e_sqrtl.c --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ lib/msun/src/e_sqrtl.c 8 Aug 2007 00:56:59 -0000 @@ -0,0 +1,92 @@ +#include +#include + +#include +#pragma STDC FENV_ACCESS ON + +#include "fpmath.h" + +long double +sqrtl(long double x) +{ + union IEEEl2bits u; + int k, r, i; + long double e, y; + fexcept_t flagp; + + r = fegetround(); + fegetexceptflag(&flagp, FE_INEXACT); + + u.e = x; + + /* If x < 0, then sqrt(x) = sNaN */ + if (u.bits.sign) + return ((x - x) / (x - x)); + + /* If x = NaN, then sqrt(x) = NaN. */ + /* If x = Inf, then sqrt(x) = Inf. */ + if (u.bits.exp == 32767) + return (x * x + x); + + /* If x = +-0, then sqrt(x) = +-0. */ + if ((u.bits.manh | u.bits.manl) == 0) + return (x); + + if (u.bits.exp == 0) { + /* Adjust subnormal numbers. */ + u.e *= 0x1.0p514; + k = u.bits.exp - 0x4200; + u.bits.exp = 0x3ffe; + } else { + /* + * x is a normal number, so break it into x = e*2^n where + * x = (2*e)*2^2k for odd n and x = (4*e)*2^2k for even n. + */ + if ((u.bits.exp - 0x3ffe) & 1) { /* n is odd. */ + k = u.bits.exp - 0x3fff; /* 2k = n - 1. */ + u.bits.exp = 0x3fff; /* u.e in [1,2). */ + } else { + k = u.bits.exp - 0x4000; /* 2k = n - 2. */ + u.bits.exp = 0x4000; /* u.e in [2,4). */ + } + } + /* + * Newton's iteration. Split u.e into a high and low part + * to achieve additional precision. + */ + y = sqrt(u.e); + e = u.e; + u.bits.manl = 0; + e = (e - u.e) / y; + y = (y + u.e / y); + u.e = y + e; + u.bits.exp += (k >> 1) - 1; + + feclearexcept(FE_INEXACT); /* Reset INEXACT flag. */ + fesetround(FE_TOWARDZERO); /* Set to round-toward-zero. */ + y = x / u.e; /* Chopped quotient (inexact?). */ + + if (!fetestexcept(FE_INEXACT)) { /* Quotient is exact. */ + if (y == u.e) { + /* Restore inexact and rounding mode. */ + fesetexceptflag(&flagp, FE_INEXACT); + fesetround(r); + return (u.e); + } + y = nexttowardl(y, 0.); /* y = y - ulp. */ + } + + fegetexceptflag(&flagp, FE_INEXACT); /* sqrt(x) is inexact. */ + if (r == FE_TONEAREST) + y = nexttowardl(y, x); /* y = y + ulp. */ + if (r == FE_UPWARD) { + u.e = nexttowardl(u.e, x); /* u.e = u.e + ulp. */ + y = nexttowardl(y, x); /* y = y + ulp. */ + } + u.e = u.e + y; /* Chopped sum. */ + u.bits.exp--; /* Correctly rounded. */ + /* Restore inexact and rounding mode. */ + fesetexceptflag(&flagp, FE_INEXACT); + fesetround(r); + return (u.e); +} Index: lib/msun/src/math.h =================================================================== RCS file: /home/ncvs/src/lib/msun/src/math.h,v retrieving revision 1.62 diff -u -p -r1.62 math.h --- lib/msun/src/math.h 7 Jan 2007 07:54:21 -0000 1.62 +++ lib/msun/src/math.h 8 Aug 2007 00:56:59 -0000 @@ -397,8 +397,8 @@ long double asinl(long double); long double atan2l(long double, long double); long double atanhl(long double); long double atanl(long double); -long double cbrtl(long double); #endif +long double cbrtl(long double); long double ceill(long double); long double copysignl(long double, long double) __pure2; #if 0 @@ -460,7 +460,9 @@ long double scalbnl(long double, int); #if 0 long double sinhl(long double); long double sinl(long double); +#endif long double sqrtl(long double); +#if 0 long double tanhl(long double); long double tanl(long double); long double tgammal(long double); Index: lib/msun/src/s_cbrtl.c =================================================================== RCS file: lib/msun/src/s_cbrtl.c diff -N lib/msun/src/s_cbrtl.c --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ lib/msun/src/s_cbrtl.c 8 Aug 2007 00:56:59 -0000 @@ -0,0 +1,90 @@ +/*- + * Copyright (c) 2007, 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. + */ + +#include +#include + +#include "fpmath.h" + +long double +cbrtl(long double x) +{ + union IEEEl2bits u; + int k; + long double e = 1., y, z; + + u.e = x; + + /* If x = +-Inf, then cbrt(x) = +-Inf. */ + /* If x = NaN, then sqrt(x) = NaN. */ + if (u.bits.exp == 32767) + return (x + x); + + /* If x = +-0, then sqrt(x) = +-0. */ + if ((u.bits.manh | u.bits.manl) == 0) + return (x); + + if (u.bits.exp == 0) { + /* Adjust subnormal numbers. */ + u.e *= 0x1.0p514; + k = u.bits.exp - 0x4200; + } else /* x is a normal number, remove the bias. */ + k = u.bits.exp - 0x3ffe; + + u.bits.exp = 0x3ffe; + + /* x = u.e * 2^n. */ + switch (k - 3 * (k / 3)) { + case 0: /* mod(n,3) == 0. */ + k -= 3; /* k = n - 3. */ + u.bits.exp += 3; /* u.e in [4,8). */ + break; + case -2: + case 1: /* mod(n,3) == 1. */ + k--; /* k = n - 1. */ + u.bits.exp++; /* u.e in [1,2). */ + break; + case -1: + case 2: /* mod(n,3) == 2. */ + k -= 2; /* k = n - 2. */ + u.bits.exp +=2; /* u.e in [2,4). */ + break; + } + /* + * Newton's iteration. Split u.e into a high and low part + * to achieve additional precision. + */ + y = cbrt(u.e); + z = 1 / (y * y); + e = u.e; + u.bits.manl = 0; + e = y + (e - u.e) * z; + y += u.e * z; + u.e = (y + e) / 3; + u.bits.exp += (k / 3); + + return (u.e); +} Index: lib/msun/src/s_ccosh.c =================================================================== RCS file: lib/msun/src/s_ccosh.c diff -N lib/msun/src/s_ccosh.c --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ lib/msun/src/s_ccosh.c 8 Aug 2007 00:56:59 -0000 @@ -0,0 +1,127 @@ +/*- + * 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 +#include + +#include "math_private.h" + +double complex +ccosh(double complex z) +{ + double x, y; + 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; /* If ix >= 0x7ff00000, then inf or NaN. */ + iy = 0x7fffffff & hy; /* If iy >= 0x7ff00000, then inf or NaN. */ + + /* 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)); + return (cpack(cosh(x) * cos(y), sinh(x) * 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). + * + * 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). + */ + 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) = NaN +- I 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 = raise. + */ + 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) = NaN + I NaN. + * + * cosh(NaN + I y) = NaN + I NaN. + * Raise the invalid floating-point exception for finite nonzero y. + * + */ + return (cpack((x * x) * (y - y), (x + x) * (y - y))); +}