summaryrefslogtreecommitdiff
path: root/usr/src/lib/libm/common/Q/sqrtl.c
diff options
context:
space:
mode:
Diffstat (limited to 'usr/src/lib/libm/common/Q/sqrtl.c')
-rw-r--r--usr/src/lib/libm/common/Q/sqrtl.c479
1 files changed, 479 insertions, 0 deletions
diff --git a/usr/src/lib/libm/common/Q/sqrtl.c b/usr/src/lib/libm/common/Q/sqrtl.c
new file mode 100644
index 0000000000..30c8f5e097
--- /dev/null
+++ b/usr/src/lib/libm/common/Q/sqrtl.c
@@ -0,0 +1,479 @@
+/*
+ * CDDL HEADER START
+ *
+ * The contents of this file are subject to the terms of the
+ * Common Development and Distribution License (the "License").
+ * You may not use this file except in compliance with the License.
+ *
+ * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
+ * or http://www.opensolaris.org/os/licensing.
+ * See the License for the specific language governing permissions
+ * and limitations under the License.
+ *
+ * When distributing Covered Code, include this CDDL HEADER in each
+ * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
+ * If applicable, add the following below this CDDL HEADER, with the
+ * fields enclosed by brackets "[]" replaced with your own identifying
+ * information: Portions Copyright [yyyy] [name of copyright owner]
+ *
+ * CDDL HEADER END
+ */
+
+/*
+ * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
+ */
+/*
+ * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
+ * Use is subject to license terms.
+ */
+
+#pragma weak sqrtl = __sqrtl
+
+#include "libm.h"
+#include "longdouble.h"
+
+extern int __swapTE(int);
+extern int __swapEX(int);
+extern enum fp_direction_type __swapRD(enum fp_direction_type);
+
+/*
+ * in struct longdouble, msw consists of
+ * unsigned short sgn:1;
+ * unsigned short exp:15;
+ * unsigned short frac1:16;
+ */
+
+#ifdef __LITTLE_ENDIAN
+
+/* array indices used to access words within a double */
+#define HIWORD 1
+#define LOWORD 0
+
+/* structure used to access words within a quad */
+union longdouble {
+ struct {
+ unsigned int frac4;
+ unsigned int frac3;
+ unsigned int frac2;
+ unsigned int msw;
+ } l;
+ long double d;
+};
+
+/* default NaN returned for sqrt(neg) */
+static const union longdouble
+ qnan = { 0xffffffff, 0xffffffff, 0xffffffff, 0x7fffffff };
+
+/* signalling NaN used to raise invalid */
+static const union {
+ unsigned u[2];
+ double d;
+} snan = { 0, 0x7ff00001 };
+
+#else
+
+/* array indices used to access words within a double */
+#define HIWORD 0
+#define LOWORD 1
+
+/* structure used to access words within a quad */
+union longdouble {
+ struct {
+ unsigned int msw;
+ unsigned int frac2;
+ unsigned int frac3;
+ unsigned int frac4;
+ } l;
+ long double d;
+};
+
+/* default NaN returned for sqrt(neg) */
+static const union longdouble
+ qnan = { 0x7fffffff, 0xffffffff, 0xffffffff, 0xffffffff };
+
+/* signalling NaN used to raise invalid */
+static const union {
+ unsigned u[2];
+ double d;
+} snan = { 0x7ff00001, 0 };
+
+#endif /* __LITTLE_ENDIAN */
+
+
+static const double
+ zero = 0.0,
+ half = 0.5,
+ one = 1.0,
+ huge = 1.0e300,
+ tiny = 1.0e-300,
+ two36 = 6.87194767360000000000e+10,
+ two30 = 1.07374182400000000000e+09,
+ two6 = 6.40000000000000000000e+01,
+ two4 = 1.60000000000000000000e+01,
+ twom18 = 3.81469726562500000000e-06,
+ twom28 = 3.72529029846191406250e-09,
+ twom42 = 2.27373675443232059479e-13,
+ twom60 = 8.67361737988403547206e-19,
+ twom62 = 2.16840434497100886801e-19,
+ twom66 = 1.35525271560688054251e-20,
+ twom90 = 8.07793566946316088742e-28,
+ twom113 = 9.62964972193617926528e-35,
+ twom124 = 4.70197740328915003187e-38;
+
+
+/*
+* Extract the exponent and normalized significand (represented as
+* an array of five doubles) from a finite, nonzero quad.
+*/
+static int
+__q_unpack(const union longdouble *x, double *s)
+{
+ union {
+ double d;
+ unsigned int l[2];
+ } u;
+ double b;
+ unsigned int lx, w[3];
+ int ex;
+
+ /* get the normalized significand and exponent */
+ ex = (int) ((x->l.msw & 0x7fffffff) >> 16);
+ lx = x->l.msw & 0xffff;
+ if (ex)
+ {
+ lx |= 0x10000;
+ w[0] = x->l.frac2;
+ w[1] = x->l.frac3;
+ w[2] = x->l.frac4;
+ }
+ else
+ {
+ if (lx | (x->l.frac2 & 0xfffe0000))
+ {
+ w[0] = x->l.frac2;
+ w[1] = x->l.frac3;
+ w[2] = x->l.frac4;
+ ex = 1;
+ }
+ else if (x->l.frac2 | (x->l.frac3 & 0xfffe0000))
+ {
+ lx = x->l.frac2;
+ w[0] = x->l.frac3;
+ w[1] = x->l.frac4;
+ w[2] = 0;
+ ex = -31;
+ }
+ else if (x->l.frac3 | (x->l.frac4 & 0xfffe0000))
+ {
+ lx = x->l.frac3;
+ w[0] = x->l.frac4;
+ w[1] = w[2] = 0;
+ ex = -63;
+ }
+ else
+ {
+ lx = x->l.frac4;
+ w[0] = w[1] = w[2] = 0;
+ ex = -95;
+ }
+ while ((lx & 0x10000) == 0)
+ {
+ lx = (lx << 1) | (w[0] >> 31);
+ w[0] = (w[0] << 1) | (w[1] >> 31);
+ w[1] = (w[1] << 1) | (w[2] >> 31);
+ w[2] <<= 1;
+ ex--;
+ }
+ }
+
+ /* extract the significand into five doubles */
+ u.l[HIWORD] = 0x42300000;
+ u.l[LOWORD] = 0;
+ b = u.d;
+ u.l[LOWORD] = lx;
+ s[0] = u.d - b;
+
+ u.l[HIWORD] = 0x40300000;
+ u.l[LOWORD] = 0;
+ b = u.d;
+ u.l[LOWORD] = w[0] & 0xffffff00;
+ s[1] = u.d - b;
+
+ u.l[HIWORD] = 0x3e300000;
+ u.l[LOWORD] = 0;
+ b = u.d;
+ u.l[HIWORD] |= w[0] & 0xff;
+ u.l[LOWORD] = w[1] & 0xffff0000;
+ s[2] = u.d - b;
+
+ u.l[HIWORD] = 0x3c300000;
+ u.l[LOWORD] = 0;
+ b = u.d;
+ u.l[HIWORD] |= w[1] & 0xffff;
+ u.l[LOWORD] = w[2] & 0xff000000;
+ s[3] = u.d - b;
+
+ u.l[HIWORD] = 0x3c300000;
+ u.l[LOWORD] = 0;
+ b = u.d;
+ u.l[LOWORD] = w[2] & 0xffffff;
+ s[4] = u.d - b;
+
+ return ex - 0x3fff;
+}
+
+
+/*
+* Pack an exponent and array of three doubles representing a finite,
+* nonzero number into a quad. Assume the sign is already there and
+* the rounding mode has been fudged accordingly.
+*/
+static void
+__q_pack(const double *z, int exp, enum fp_direction_type rm,
+ union longdouble *x, int *inexact)
+{
+ union {
+ double d;
+ unsigned int l[2];
+ } u;
+ double s[3], t, t2;
+ unsigned int msw, frac2, frac3, frac4;
+
+ /* bias exponent and strip off integer bit */
+ exp += 0x3fff;
+ s[0] = z[0] - one;
+ s[1] = z[1];
+ s[2] = z[2];
+
+ /*
+ * chop the significand to obtain the fraction;
+ * use round-to-minus-infinity to ensure chopping
+ */
+ (void) __swapRD(fp_negative);
+
+ /* extract the first eighty bits of fraction */
+ t = s[1] + s[2];
+ u.d = two36 + (s[0] + t);
+ msw = u.l[LOWORD];
+ s[0] -= (u.d - two36);
+
+ u.d = two4 + (s[0] + t);
+ frac2 = u.l[LOWORD];
+ s[0] -= (u.d - two4);
+
+ u.d = twom28 + (s[0] + t);
+ frac3 = u.l[LOWORD];
+ s[0] -= (u.d - twom28);
+
+ /* condense the remaining fraction; errors here won't matter */
+ t = s[0] + s[1];
+ s[1] = ((s[0] - t) + s[1]) + s[2];
+ s[0] = t;
+
+ /* get the last word of fraction */
+ u.d = twom60 + (s[0] + s[1]);
+ frac4 = u.l[LOWORD];
+ s[0] -= (u.d - twom60);
+
+ /*
+ * keep track of what's left for rounding; note that
+ * t2 will be non-negative due to rounding mode
+ */
+ t = s[0] + s[1];
+ t2 = (s[0] - t) + s[1];
+
+ if (t != zero)
+ {
+ *inexact = 1;
+
+ /* decide whether to round the fraction up */
+ if (rm == fp_positive || (rm == fp_nearest && (t > twom113 ||
+ (t == twom113 && (t2 != zero || frac4 & 1)))))
+ {
+ /* round up and renormalize if necessary */
+ if (++frac4 == 0)
+ if (++frac3 == 0)
+ if (++frac2 == 0)
+ if (++msw == 0x10000)
+ {
+ msw = 0;
+ exp++;
+ }
+ }
+ }
+
+ /* assemble the result */
+ x->l.msw |= msw | (exp << 16);
+ x->l.frac2 = frac2;
+ x->l.frac3 = frac3;
+ x->l.frac4 = frac4;
+}
+
+
+/*
+* Compute the square root of x and place the TP result in s.
+*/
+static void
+__q_tp_sqrt(const double *x, double *s)
+{
+ double c, rr, r[3], tt[3], t[5];
+
+ /* approximate the divisor for the Newton iteration */
+ c = sqrt((x[0] + x[1]) + x[2]);
+ rr = half / c;
+
+ /* compute the first five "digits" of the square root */
+ t[0] = (c + two30) - two30;
+ tt[0] = t[0] + t[0];
+ r[0] = ((x[0] - t[0] * t[0]) + x[1]) + x[2];
+
+ t[1] = (rr * (r[0] + x[3]) + two6) - two6;
+ tt[1] = t[1] + t[1];
+ r[0] -= tt[0] * t[1];
+ r[1] = x[3] - t[1] * t[1];
+ c = (r[1] + twom18) - twom18;
+ r[0] += c;
+ r[1] = (r[1] - c) + x[4];
+
+ t[2] = (rr * (r[0] + r[1]) + twom18) - twom18;
+ tt[2] = t[2] + t[2];
+ r[0] -= tt[0] * t[2];
+ r[1] -= tt[1] * t[2];
+ c = (r[1] + twom42) - twom42;
+ r[0] += c;
+ r[1] = (r[1] - c) - t[2] * t[2];
+
+ t[3] = (rr * (r[0] + r[1]) + twom42) - twom42;
+ r[0] = ((r[0] - tt[0] * t[3]) + r[1]) - tt[1] * t[3];
+ r[1] = -tt[2] * t[3];
+ c = (r[1] + twom90) - twom90;
+ r[0] += c;
+ r[1] = (r[1] - c) - t[3] * t[3];
+
+ t[4] = (rr * (r[0] + r[1]) + twom66) - twom66;
+
+ /* here we just need to get the sign of the remainder */
+ c = (((((r[0] - tt[0] * t[4]) - tt[1] * t[4]) + r[1])
+ - tt[2] * t[4]) - (t[3] + t[3]) * t[4]) - t[4] * t[4];
+
+ /* reduce to three doubles */
+ t[0] += t[1];
+ t[1] = t[2] + t[3];
+ t[2] = t[4];
+
+ /* if the third term might lie on a rounding boundary, perturb it */
+ if (c != zero && t[2] == (twom62 + t[2]) - twom62)
+ {
+ if (c < zero)
+ t[2] -= twom124;
+ else
+ t[2] += twom124;
+ }
+
+ /* condense the square root */
+ c = t[1] + t[2];
+ t[2] += (t[1] - c);
+ t[1] = c;
+ c = t[0] + t[1];
+ s[1] = t[1] + (t[0] - c);
+ s[0] = c;
+ if (s[1] == zero)
+ {
+ c = s[0] + t[2];
+ s[1] = t[2] + (s[0] - c);
+ s[0] = c;
+ s[2] = zero;
+ }
+ else
+ {
+ c = s[1] + t[2];
+ s[2] = t[2] + (s[1] - c);
+ s[1] = c;
+ }
+}
+
+
+long double
+sqrtl(long double ldx)
+{
+ union longdouble x;
+ volatile double t;
+ double xx[5], zz[3];
+ enum fp_direction_type rm;
+ int ex, inexact, exc, traps;
+
+ /* clear cexc */
+ t = zero;
+ t -= zero;
+
+ /* check for zero operand */
+ x.d = ldx;
+ if (!((x.l.msw & 0x7fffffff) | x.l.frac2 | x.l.frac3 | x.l.frac4))
+ return ldx;
+
+ /* handle nan and inf cases */
+ if ((x.l.msw & 0x7fffffff) >= 0x7fff0000)
+ {
+ if ((x.l.msw & 0xffff) | x.l.frac2 | x.l.frac3 | x.l.frac4)
+ {
+ if (!(x.l.msw & 0x8000))
+ {
+ /* snan, signal invalid */
+ t += snan.d;
+ }
+ x.l.msw |= 0x8000;
+ return x.d;
+ }
+ if (x.l.msw & 0x80000000)
+ {
+ /* sqrt(-inf), signal invalid */
+ t = -one;
+ t = sqrt(t);
+ return qnan.d;
+ }
+ /* sqrt(inf), return inf */
+ return x.d;
+ }
+
+ /* handle negative numbers */
+ if (x.l.msw & 0x80000000)
+ {
+ t = -one;
+ t = sqrt(t);
+ return qnan.d;
+ }
+
+ /* now x is finite, positive */
+
+ traps = __swapTE(0);
+ exc = __swapEX(0);
+ rm = __swapRD(fp_nearest);
+
+ ex = __q_unpack(&x, xx);
+ if (ex & 1)
+ {
+ /* make exponent even */
+ xx[0] += xx[0];
+ xx[1] += xx[1];
+ xx[2] += xx[2];
+ xx[3] += xx[3];
+ xx[4] += xx[4];
+ ex--;
+ }
+ __q_tp_sqrt(xx, zz);
+
+ /* put everything together */
+ x.l.msw = 0;
+ inexact = 0;
+ __q_pack(zz, ex >> 1, rm, &x, &inexact);
+
+ (void) __swapRD(rm);
+ (void) __swapEX(exc);
+ (void) __swapTE(traps);
+ if (inexact)
+ {
+ t = huge;
+ t += tiny;
+ }
+ return x.d;
+}
generated by cgit v1.2.3 (git 2.39.1) at 2025-04-27 13:55:45 +0000