common/Q/sqrtl.c

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