__vrsqrt.c revision 25c28e83beb90e7c80452a7c818c5e6f73a07dc8
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * CDDL HEADER START
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 * 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 * 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 * CDDL HEADER END
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Use is subject to license terms.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#define LO(x) *(unsigned*)x
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#define HI(x) *(int*)x
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#define restrict _Restrict
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#define restrict
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis/* double rsqrt(double x)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * 1. Special cases:
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * for x = NaN => QNaN;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * for x = +Inf => 0;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * for x is negative, -Inf => QNaN + invalid;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * for x = +0 => +Inf + divide-by-zero;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * for x = -0 => -Inf + divide-by-zero.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * 2. Computes reciprocal square root from:
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * x = m * 2**n
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * m = [0.5, 2),
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * n = ((exponent + 1) & ~1).
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * rsqrt(x) = 1/sqrt( m * 2**n ) = (2 ** (-n/2)) * (1/sqrt(m))
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * 2. Computes 1/sqrt(m) from:
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * 1/sqrt(m) = (1/sqrt(m0)) * (1/sqrt(1 + (1/m0)*dm))
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * m = m0 + dm,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * m0 = 0.5 * (1 + k/64) for m = [0.5, 0.5+127/256), k = [0, 63];
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * m0 = 1.0 * (0 + k/64) for m = [0.5+127/256, 1.0+127/128), k = [64, 127];
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * m0 = 2.0 for m = [1.0+127/128, 2.0), k = 128.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * 1/sqrt(m0) is looked up in a table,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * 1/m0 is computed as (1/sqrt(m0)) * (1/sqrt(m0)).
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * 1/sqrt(1 + (1/m0)*dm) is computed using approximation:
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * 1/sqrt(1 + z) = (((((a6 * z + a5) * z + a4) * z + a3)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * * z + a2) * z + a1) * z + a0
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * where z = [-1/128, 1/128].
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * The maximum relative error for the approximating
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * polynomial is 2**(-56.26).
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Maximum error observed: less than 0.563 ulp after 1.500.000.000
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisextern double sqrt (double);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisextern const double __vlibm_TBL_rsqrt[];
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis__vrsqrt_n(int n, double * restrict px, int stridex, double * restrict py, int stridey);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisstatic const double
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis__vrsqrt(int n, double * restrict px, int stridex, double * restrict py, int stridey)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis for (; n > 1 ; n--)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (hx >= 0x7ff00000) /* X = NaN or Inf */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (hx < 0x00100000) /* X = denormal, zero or negative */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis res_ch0 = (resh0 + 0x00002000) & 0x7fffc000;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis dexp_hi0 = ((double*)((char*)__vlibm_TBL_rsqrt + ind0))[0];
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis dexp_lo0 = ((double*)((char*)__vlibm_TBL_rsqrt + ind0))[1];
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis res = (((((K6 * xx0 + K5) * xx0 + K4) * xx0 + K3) * xx0 + K2) * xx0 + K1) * xx0;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis res = dexp_hi0 * res + dexp_lo0 + dexp_hi0;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis else /* X = negative */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis __vrsqrt_n(n_n, spx, stridex, spy, stridey);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (hx >= 0x7ff00000) /* X = NaN or Inf */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis else if (hx < 0x00100000) /* X = denormal, zero or negative */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis res_ch0 = (resh0 + 0x00002000) & 0x7fffc000;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis dexp_hi0 = ((double*)((char*)__vlibm_TBL_rsqrt + ind0))[0];
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis dexp_lo0 = ((double*)((char*)__vlibm_TBL_rsqrt + ind0))[1];
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis res = (((((K6 * xx0 + K5) * xx0 + K4) * xx0 + K3) * xx0 + K2) * xx0 + K1) * xx0;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis res = dexp_hi0 * res + dexp_lo0 + dexp_hi0;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis else /* X = negative */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis res_ch0 = (resh0 + 0x00002000) & 0x7fffc000;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis dexp_hi0 = ((double*)((char*)__vlibm_TBL_rsqrt + ind0))[0];
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis dexp_lo0 = ((double*)((char*)__vlibm_TBL_rsqrt + ind0))[1];
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis res = (((((K6 * xx0 + K5) * xx0 + K4) * xx0 + K3) * xx0 + K2) * xx0 + K1) * xx0;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis res = dexp_hi0 * res + dexp_lo0 + dexp_hi0;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis__vrsqrt_n(int n, double * restrict px, int stridex, double * restrict py, int stridey)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis res_ch0 = (resh0 + 0x00002000) & 0x7fffc000;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis res_ch1 = (resh1 + 0x00002000) & 0x7fffc000;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis res_ch2 = (resh2 + 0x00002000) & 0x7fffc000;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis dexp_hi0 = ((double*)((char*)__vlibm_TBL_rsqrt + ind0))[0];
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis dexp_hi1 = ((double*)((char*)__vlibm_TBL_rsqrt + ind1))[0];
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis dexp_hi2 = ((double*)((char*)__vlibm_TBL_rsqrt + ind2))[0];
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis dexp_lo0 = ((double*)((char*)__vlibm_TBL_rsqrt + ind0))[1];
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis dexp_lo1 = ((double*)((char*)__vlibm_TBL_rsqrt + ind1))[1];
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis dexp_lo2 = ((double*)((char*)__vlibm_TBL_rsqrt + ind2))[1];
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis res0 = (((((K6 * xx0 + K5) * xx0 + K4) * xx0 + K3) * xx0 + K2) * xx0 + K1) * xx0;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis res1 = (((((K6 * xx1 + K5) * xx1 + K4) * xx1 + K3) * xx1 + K2) * xx1 + K1) * xx1;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis res2 = (((((K6 * xx2 + K5) * xx2 + K4) * xx2 + K3) * xx2 + K2) * xx2 + K1) * xx2;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis res0 = dexp_hi0 * res0 + dexp_lo0 + dexp_hi0;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis res1 = dexp_hi1 * res1 + dexp_lo1 + dexp_hi1;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis res2 = dexp_hi2 * res2 + dexp_lo2 + dexp_hi2;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis for(; n > 0 ; n--)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis res_ch0 = (resh0 + 0x00002000) & 0x7fffc000;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis dexp_hi0 = ((double*)((char*)__vlibm_TBL_rsqrt + ind0))[0];
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis dexp_lo0 = ((double*)((char*)__vlibm_TBL_rsqrt + ind0))[1];
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis res0 = (((((K6 * xx0 + K5) * xx0 + K4) * xx0 + K3) * xx0 + K2) * xx0 + K1) * xx0;