/*
* 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
* 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.
*/
#include <sys/isa_defs.h>
#include "libm_inlines.h"
#ifdef _LITTLE_ENDIAN
#define LO(x) *(unsigned*)x
#else
#define HI(x) *(int*)x
#endif
#ifdef __RESTRICT
#define restrict _Restrict
#else
#define restrict
#endif
/* double rhypot(double x, double y)
*
* Method :
* 1. Special cases:
* x or y = Inf => 0
* x or y = NaN => QNaN
* x and y = 0 => Inf + divide-by-zero
* 2. Computes rhypot(x,y):
* rhypot(x,y) = m * sqrt(1/(xnm * xnm + ynm * ynm))
* Where:
* m = 1/max(|x|,|y|)
* xnm = x * m
* ynm = y * m
*
* Compute 1/(xnm * xnm + ynm * ynm) by simulating
* muti-precision arithmetic.
*
* Accuracy:
* Maximum error observed: less than 0.869 ulp after 1.000.000.000
* results.
*/
extern double sqrt(double);
extern double fabs(double);
static const int __vlibm_TBL_rhypot[] = {
/* i = [0,127]
* TBL[i] = 0x3ff00000 + *(int*)&(1.0 / *(double*)&(0x3ff0000000000000ULL + (i << 45))); */
0x7fe00000, 0x7fdfc07f, 0x7fdf81f8, 0x7fdf4465,
0x7fdf07c1, 0x7fdecc07, 0x7fde9131, 0x7fde573a,
0x7fde1e1e, 0x7fdde5d6, 0x7fddae60, 0x7fdd77b6,
0x7fdd41d4, 0x7fdd0cb5, 0x7fdcd856, 0x7fdca4b3,
0x7fdc71c7, 0x7fdc3f8f, 0x7fdc0e07, 0x7fdbdd2b,
0x7fdbacf9, 0x7fdb7d6c, 0x7fdb4e81, 0x7fdb2036,
0x7fdaf286, 0x7fdac570, 0x7fda98ef, 0x7fda6d01,
0x7fda41a4, 0x7fda16d3, 0x7fd9ec8e, 0x7fd9c2d1,
0x7fd99999, 0x7fd970e4, 0x7fd948b0, 0x7fd920fb,
0x7fd8f9c1, 0x7fd8d301, 0x7fd8acb9, 0x7fd886e5,
0x7fd86186, 0x7fd83c97, 0x7fd81818, 0x7fd7f405,
0x7fd7d05f, 0x7fd7ad22, 0x7fd78a4c, 0x7fd767dc,
0x7fd745d1, 0x7fd72428, 0x7fd702e0, 0x7fd6e1f7,
0x7fd6c16c, 0x7fd6a13c, 0x7fd68168, 0x7fd661ec,
0x7fd642c8, 0x7fd623fa, 0x7fd60581, 0x7fd5e75b,
0x7fd5c988, 0x7fd5ac05, 0x7fd58ed2, 0x7fd571ed,
0x7fd55555, 0x7fd53909, 0x7fd51d07, 0x7fd50150,
0x7fd4e5e0, 0x7fd4cab8, 0x7fd4afd6, 0x7fd49539,
0x7fd47ae1, 0x7fd460cb, 0x7fd446f8, 0x7fd42d66,
0x7fd41414, 0x7fd3fb01, 0x7fd3e22c, 0x7fd3c995,
0x7fd3b13b, 0x7fd3991c, 0x7fd38138, 0x7fd3698d,
0x7fd3521c, 0x7fd33ae4, 0x7fd323e3, 0x7fd30d19,
0x7fd2f684, 0x7fd2e025, 0x7fd2c9fb, 0x7fd2b404,
0x7fd29e41, 0x7fd288b0, 0x7fd27350, 0x7fd25e22,
0x7fd24924, 0x7fd23456, 0x7fd21fb7, 0x7fd20b47,
0x7fd1f704, 0x7fd1e2ef, 0x7fd1cf06, 0x7fd1bb4a,
0x7fd1a7b9, 0x7fd19453, 0x7fd18118, 0x7fd16e06,
0x7fd15b1e, 0x7fd1485f, 0x7fd135c8, 0x7fd12358,
0x7fd11111, 0x7fd0fef0, 0x7fd0ecf5, 0x7fd0db20,
0x7fd0c971, 0x7fd0b7e6, 0x7fd0a681, 0x7fd0953f,
0x7fd08421, 0x7fd07326, 0x7fd0624d, 0x7fd05197,
0x7fd04104, 0x7fd03091, 0x7fd02040, 0x7fd01010,
};
static const unsigned long long LCONST[] = {
0x3ff0000000000000ULL, /* DONE = 1.0 */
0x4000000000000000ULL, /* DTWO = 2.0 */
0x4230000000000000ULL, /* D2ON36 = 2**36 */
0x7fd0000000000000ULL, /* D2ON1022 = 2**1022 */
0x3cb0000000000000ULL, /* D2ONM52 = 2**-52 */
};
#define RET_SC(I) \
if (--n <= 0) \
break; \
goto start##I;
{ \
RET_SC(I) \
}
#define PREP(I) \
hx##I &= 0x7fffffff; \
hy##I &= 0x7fffffff; \
{ \
x = *px; \
y = *py; \
\
} \
x##I = *px; \
y##I = *py; \
{ \
x = x##I; \
y = y##I; \
\
\
x = fabs(x); \
y = fabs(y); \
\
x = *(long long*)&x; \
y = *(long long*)&y; \
\
x *= D2ONM52; \
y *= D2ONM52; \
\
\
\
\
\
\
\
} \
j0 &= 0x7ff00000; \
void
{
int i = 0;
double x, y;
do
{
PREP(0)
i = 1;
if (--n <= 0)
break;
PREP(1)
i = 2;
if (--n <= 0)
break;
PREP(2)
i = 0;
} while (--n > 0);
if (i > 0)
{
if (i > 1)
{
}
}
}