? sqrtl.diff Index: Makefile =================================================================== RCS file: /home/ncvs/src/lib/msun/Makefile,v retrieving revision 1.78 diff -u -p -r1.78 Makefile --- Makefile 21 May 2007 02:49:08 -0000 1.78 +++ Makefile 6 Aug 2007 02:30:58 -0000 @@ -68,6 +68,10 @@ 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 +.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 Index: Symbol.map =================================================================== RCS file: /home/ncvs/src/lib/msun/Symbol.map,v retrieving revision 1.4 diff -u -p -r1.4 Symbol.map --- Symbol.map 29 Apr 2007 14:05:20 -0000 1.4 +++ Symbol.map 6 Aug 2007 02:30:58 -0000 @@ -56,6 +56,7 @@ FBSD_1.0 { sinhf; sqrt; sqrtf; + sqrtl; asinh; asinhf; atan; Index: 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 --- man/sqrt.3 9 Jan 2007 01:02:06 -0000 1.13 +++ man/sqrt.3 6 Aug 2007 02:30:58 -0000 @@ -28,14 +28,15 @@ .\" 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 sqrt , -.Nm sqrtf +.Nm sqrtf , +.Nm sqrtl .Nd cube root and square root functions .Sh LIBRARY .Lb libm @@ -50,6 +51,9 @@ .Ft float .Fn sqrtf "float x" .Sh DESCRIPTION +.Ft "long double" +.Fn sqrtl "long double x" +.Sh DESCRIPTION The .Fn cbrt and the Index: src/e_sqrtl.c =================================================================== RCS file: src/e_sqrtl.c diff -N src/e_sqrtl.c --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ src/e_sqrtl.c 6 Aug 2007 02:30:58 -0000 @@ -0,0 +1,118 @@ +/*- + * 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 +#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, z; + fexcept_t flagp; + + r = fegetround(); + fegetexceptflag(&flagp, FE_INEXACT); + + u.e = x; + + /* If x = +-0, then sqrt(x) = +-0. */ + if ((u.bits.manh | u.bits.manl) == 0) + return (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 (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. */ + z = x / u.e; /* Chopped quotient (inexact?). */ + + if (!fetestexcept(FE_INEXACT)) { /* Quotient is exact. */ + if (z == u.e) { + /* Restore inexact and rounding mode. */ + fesetexceptflag(&flagp, FE_INEXACT); + fesetround(r); + return (u.e); + } + z = nexttowardl(z, 0.); /* z = z - ulp. */ + } + + fegetexceptflag(&flagp, FE_INEXACT); /* sqrt(x) is inexact. */ + if (r == FE_TONEAREST) + z = nexttowardl(z, x); /* z = z + ulp. */ + if (r == FE_UPWARD) { + u.e = nexttowardl(u.e, x); /* u.e = u.e + ulp. */ + z = nexttowardl(z, x); /* z = z + ulp. */ + } + u.e = u.e + z; /* Chopped sum. */ + u.bits.exp--; /* Correctly rounded. */ + /* Restore inexact and rounding mode. */ + fesetexceptflag(&flagp, FE_INEXACT); + fesetround(r); + return (u.e); +} Index: 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 --- src/math.h 7 Jan 2007 07:54:21 -0000 1.62 +++ src/math.h 6 Aug 2007 02:30:58 -0000 @@ -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);