1 /*
   2  * CDDL HEADER START
   3  *
   4  * The contents of this file are subject to the terms of the
   5  * Common Development and Distribution License (the "License").
   6  * You may not use this file except in compliance with the License.
   7  *
   8  * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
   9  * or http://www.opensolaris.org/os/licensing.
  10  * See the License for the specific language governing permissions
  11  * and limitations under the License.
  12  *
  13  * When distributing Covered Code, include this CDDL HEADER in each
  14  * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
  15  * If applicable, add the following below this CDDL HEADER, with the
  16  * fields enclosed by brackets "[]" replaced with your own identifying
  17  * information: Portions Copyright [yyyy] [name of copyright owner]
  18  *
  19  * CDDL HEADER END
  20  */
  21 
  22 /*
  23  * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
  24  */
  25 /*
  26  * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
  27  * Use is subject to license terms.
  28  */
  29 
  30 #include <sys/isa_defs.h>
  31 #include "libm_inlines.h"
  32 
  33 #ifdef _LITTLE_ENDIAN
  34 #define HI(x)   *(1+(int*)x)
  35 #define LO(x)   *(unsigned*)x
  36 #else
  37 #define HI(x)   *(int*)x
  38 #define LO(x)   *(1+(unsigned*)x)
  39 #endif
  40 
  41 #ifdef __RESTRICT
  42 #define restrict _Restrict
  43 #else
  44 #define restrict
  45 #endif
  46 
  47 /* double rsqrt(double x)
  48  *
  49  * Method :
  50  *      1. Special cases:
  51  *              for x = NaN                             => QNaN;
  52  *              for x = +Inf                            => 0;
  53  *              for x is negative, -Inf                 => QNaN + invalid;
  54  *              for x = +0                              => +Inf + divide-by-zero;
  55  *              for x = -0                              => -Inf + divide-by-zero.
  56  *      2. Computes reciprocal square root from:
  57  *              x = m * 2**n
  58  *      Where:
  59  *              m = [0.5, 2),
  60  *              n = ((exponent + 1) & ~1).
  61  *      Then:
  62  *              rsqrt(x) = 1/sqrt( m * 2**n ) = (2 ** (-n/2)) * (1/sqrt(m))
  63  *      2. Computes 1/sqrt(m) from:
  64  *              1/sqrt(m) = (1/sqrt(m0)) * (1/sqrt(1 + (1/m0)*dm))
  65  *      Where:
  66  *              m = m0 + dm,
  67  *              m0 = 0.5 * (1 + k/64) for m = [0.5,         0.5+127/256), k = [0, 63];
  68  *              m0 = 1.0 * (0 + k/64) for m = [0.5+127/256, 1.0+127/128), k = [64, 127];
  69  *              m0 = 2.0              for m = [1.0+127/128, 2.0),         k = 128.
  70  *      Then:
  71  *              1/sqrt(m0) is looked up in a table,
  72  *              1/m0 is computed as (1/sqrt(m0)) * (1/sqrt(m0)).
  73  *              1/sqrt(1 + (1/m0)*dm) is computed using approximation:
  74  *                      1/sqrt(1 + z) = (((((a6 * z + a5) * z + a4) * z + a3)
  75  *                                              * z + a2) * z + a1) * z + a0
  76  *                      where z = [-1/128, 1/128].
  77  *
  78  * Accuracy:
  79  *      The maximum relative error for the approximating
  80  *      polynomial is 2**(-56.26).
  81  *      Maximum error observed: less than 0.563 ulp after 1.500.000.000
  82  *      results.
  83  */
  84 
  85 extern double sqrt (double);
  86 extern const double __vlibm_TBL_rsqrt[];
  87 
  88 static void
  89 __vrsqrt_n(int n, double * restrict px, int stridex, double * restrict py, int stridey);
  90 
  91 #pragma no_inline(__vrsqrt_n)
  92 
  93 #define RETURN(ret)                                             \
  94 {                                                               \
  95         *py = (ret);                                            \
  96         py += stridey;                                          \
  97         if (n_n == 0)                                           \
  98         {                                                       \
  99                 spx = px; spy = py;                             \
 100                 hx = HI(px);                                    \
 101                 continue;                                       \
 102         }                                                       \
 103         n--;                                                    \
 104         break;                                                  \
 105 }
 106 
 107 static const double
 108         DONE = 1.0,
 109         K1 = -5.00000000000005209867e-01,
 110         K2 =  3.75000000000004884257e-01,
 111         K3 = -3.12499999317136886551e-01,
 112         K4 =  2.73437499359815081532e-01,
 113         K5 = -2.46116125605037803130e-01,
 114         K6 =  2.25606914648617522896e-01;
 115 
 116 void
 117 __vrsqrt(int n, double * restrict px, int stridex, double * restrict py, int stridey)
 118 {
 119         double          *spx, *spy;
 120         int             ax, lx, hx, n_n;
 121         double          res;
 122 
 123         while (n > 1)
 124         {
 125                 n_n = 0;
 126                 spx = px;
 127                 spy = py;
 128                 hx = HI(px);
 129                 for (; n > 1 ; n--)
 130                 {
 131                         px += stridex;
 132                         if (hx >= 0x7ff00000)                /* X = NaN or Inf       */
 133                         {
 134                                 res = *(px - stridex);
 135                                 RETURN (DONE / res)
 136                         }
 137 
 138                         py += stridey;
 139 
 140                         if (hx < 0x00100000)         /* X = denormal, zero or negative       */
 141                         {
 142                                 py -= stridey;
 143                                 ax = hx & 0x7fffffff;
 144                                 lx = LO((px - stridex));
 145                                 res = *(px - stridex);
 146 
 147                                 if ((ax | lx) == 0)     /* |X| = zero   */
 148                                 {
 149                                         RETURN (DONE / res)
 150                                 }
 151                                 else if (hx >= 0)    /* X = denormal */
 152                                 {
 153                                         double          res_c0, dsqrt_exp0;
 154                                         int             ind0, sqrt_exp0;
 155                                         double          xx0, dexp_hi0, dexp_lo0;
 156                                         int             hx0, resh0, res_ch0;
 157 
 158                                         res = *(long long*)&res;
 159 
 160                                         hx0 = HI(&res);
 161                                         sqrt_exp0 = (0x817 - (hx0 >> 21)) << 20;
 162                                         ind0 = (((hx0 >> 10) & 0x7f8) + 8) & -16;
 163 
 164                                         resh0 = (hx0 & 0x001fffff) | 0x3fe00000;
 165                                         res_ch0 = (resh0 + 0x00002000) & 0x7fffc000;
 166                                         HI(&res) = resh0;
 167                                         HI(&res_c0) = res_ch0;
 168                                         LO(&res_c0) = 0;
 169 
 170                                         dexp_hi0 = ((double*)((char*)__vlibm_TBL_rsqrt + ind0))[0];
 171                                         dexp_lo0 = ((double*)((char*)__vlibm_TBL_rsqrt + ind0))[1];
 172                                         xx0 = dexp_hi0 * dexp_hi0;
 173                                         xx0 = (res - res_c0) * xx0;
 174                                         res = (((((K6 * xx0 + K5) * xx0 + K4) * xx0 + K3) * xx0 + K2) * xx0 + K1) * xx0;
 175 
 176                                         res = dexp_hi0 * res + dexp_lo0 + dexp_hi0;
 177 
 178                                         HI(&dsqrt_exp0) = sqrt_exp0;
 179                                         LO(&dsqrt_exp0) = 0;
 180                                         res *= dsqrt_exp0;
 181 
 182                                         RETURN (res)
 183                                 }
 184                                 else    /* X = negative */
 185                                 {
 186                                         RETURN (sqrt(res))
 187                                 }
 188                         }
 189                         n_n++;
 190                         hx = HI(px);
 191                 }
 192                 if (n_n > 0)
 193                         __vrsqrt_n(n_n, spx, stridex, spy, stridey);
 194         }
 195         if (n > 0)
 196         {
 197                 hx = HI(px);
 198 
 199                 if (hx >= 0x7ff00000)                /* X = NaN or Inf       */
 200                 {
 201                         res = *px;
 202                         *py = DONE / res;
 203                 }
 204                 else if (hx < 0x00100000)    /* X = denormal, zero or negative       */
 205                 {
 206                         ax = hx & 0x7fffffff;
 207                         lx = LO(px);
 208                         res = *px;
 209 
 210                         if ((ax | lx) == 0)     /* |X| = zero   */
 211                         {
 212                                 *py = DONE / res;
 213                         }
 214                         else if (hx >= 0)    /* X = denormal */
 215                         {
 216                                 double          res_c0, dsqrt_exp0;
 217                                 int             ind0, sqrt_exp0;
 218                                 double          xx0, dexp_hi0, dexp_lo0;
 219                                 int             hx0, resh0, res_ch0;
 220 
 221                                 res = *(long long*)&res;
 222 
 223                                 hx0 = HI(&res);
 224                                 sqrt_exp0 = (0x817 - (hx0 >> 21)) << 20;
 225                                 ind0 = (((hx0 >> 10) & 0x7f8) + 8) & -16;
 226 
 227                                 resh0 = (hx0 & 0x001fffff) | 0x3fe00000;
 228                                 res_ch0 = (resh0 + 0x00002000) & 0x7fffc000;
 229                                 HI(&res) = resh0;
 230                                 HI(&res_c0) = res_ch0;
 231                                 LO(&res_c0) = 0;
 232 
 233                                 dexp_hi0 = ((double*)((char*)__vlibm_TBL_rsqrt + ind0))[0];
 234                                 dexp_lo0 = ((double*)((char*)__vlibm_TBL_rsqrt + ind0))[1];
 235                                 xx0 = dexp_hi0 * dexp_hi0;
 236                                 xx0 = (res - res_c0) * xx0;
 237                                 res = (((((K6 * xx0 + K5) * xx0 + K4) * xx0 + K3) * xx0 + K2) * xx0 + K1) * xx0;
 238 
 239                                 res = dexp_hi0 * res + dexp_lo0 + dexp_hi0;
 240 
 241                                 HI(&dsqrt_exp0) = sqrt_exp0;
 242                                 LO(&dsqrt_exp0) = 0;
 243                                 res *= dsqrt_exp0;
 244 
 245                                 *py = res;
 246                         }
 247                         else    /* X = negative */
 248                         {
 249                                 *py = sqrt(res);
 250                         }
 251                 }
 252                 else
 253                 {
 254                         double          res_c0, dsqrt_exp0;
 255                         int             ind0, sqrt_exp0;
 256                         double          xx0, dexp_hi0, dexp_lo0;
 257                         int             resh0, res_ch0;
 258 
 259                         sqrt_exp0 = (0x5fe - (hx >> 21)) << 20;
 260                         ind0 = (((hx >> 10) & 0x7f8) + 8) & -16;
 261 
 262                         resh0 = (hx & 0x001fffff) | 0x3fe00000;
 263                         res_ch0 = (resh0 + 0x00002000) & 0x7fffc000;
 264                         HI(&res) = resh0;
 265                         LO(&res) = LO(px);
 266                         HI(&res_c0) = res_ch0;
 267                         LO(&res_c0) = 0;
 268 
 269                         dexp_hi0 = ((double*)((char*)__vlibm_TBL_rsqrt + ind0))[0];
 270                         dexp_lo0 = ((double*)((char*)__vlibm_TBL_rsqrt + ind0))[1];
 271                         xx0 = dexp_hi0 * dexp_hi0;
 272                         xx0 = (res - res_c0) * xx0;
 273                         res = (((((K6 * xx0 + K5) * xx0 + K4) * xx0 + K3) * xx0 + K2) * xx0 + K1) * xx0;
 274 
 275                         res = dexp_hi0 * res + dexp_lo0 + dexp_hi0;
 276 
 277                         HI(&dsqrt_exp0) = sqrt_exp0;
 278                         LO(&dsqrt_exp0) = 0;
 279                         res *= dsqrt_exp0;
 280 
 281                         *py = res;
 282                 }
 283         }
 284 }
 285 
 286 static void
 287 __vrsqrt_n(int n, double * restrict px, int stridex, double * restrict py, int stridey)
 288 {
 289         double          res0, res_c0, dsqrt_exp0;
 290         double          res1, res_c1, dsqrt_exp1;
 291         double          res2, res_c2, dsqrt_exp2;
 292         int             ind0, sqrt_exp0;
 293         int             ind1, sqrt_exp1;
 294         int             ind2, sqrt_exp2;
 295         double          xx0, dexp_hi0, dexp_lo0;
 296         double          xx1, dexp_hi1, dexp_lo1;
 297         double          xx2, dexp_hi2, dexp_lo2;
 298         int             hx0, resh0, res_ch0;
 299         int             hx1, resh1, res_ch1;
 300         int             hx2, resh2, res_ch2;
 301 
 302         LO(&dsqrt_exp0) = 0;
 303         LO(&dsqrt_exp1) = 0;
 304         LO(&dsqrt_exp2) = 0;
 305         LO(&res_c0) = 0;
 306         LO(&res_c1) = 0;
 307         LO(&res_c2) = 0;
 308 
 309         for(; n > 2 ; n -= 3)
 310         {
 311                 hx0 = HI(px);
 312                 LO(&res0) = LO(px);
 313                 px += stridex;
 314 
 315                 hx1 = HI(px);
 316                 LO(&res1) = LO(px);
 317                 px += stridex;
 318 
 319                 hx2 = HI(px);
 320                 LO(&res2) = LO(px);
 321                 px += stridex;
 322 
 323                 sqrt_exp0 = (0x5fe - (hx0 >> 21)) << 20;
 324                 sqrt_exp1 = (0x5fe - (hx1 >> 21)) << 20;
 325                 sqrt_exp2 = (0x5fe - (hx2 >> 21)) << 20;
 326                 ind0 = (((hx0 >> 10) & 0x7f8) + 8) & -16;
 327                 ind1 = (((hx1 >> 10) & 0x7f8) + 8) & -16;
 328                 ind2 = (((hx2 >> 10) & 0x7f8) + 8) & -16;
 329 
 330                 resh0 = (hx0 & 0x001fffff) | 0x3fe00000;
 331                 resh1 = (hx1 & 0x001fffff) | 0x3fe00000;
 332                 resh2 = (hx2 & 0x001fffff) | 0x3fe00000;
 333                 res_ch0 = (resh0 + 0x00002000) & 0x7fffc000;
 334                 res_ch1 = (resh1 + 0x00002000) & 0x7fffc000;
 335                 res_ch2 = (resh2 + 0x00002000) & 0x7fffc000;
 336                 HI(&res0) = resh0;
 337                 HI(&res1) = resh1;
 338                 HI(&res2) = resh2;
 339                 HI(&res_c0) = res_ch0;
 340                 HI(&res_c1) = res_ch1;
 341                 HI(&res_c2) = res_ch2;
 342 
 343                 dexp_hi0 = ((double*)((char*)__vlibm_TBL_rsqrt + ind0))[0];
 344                 dexp_hi1 = ((double*)((char*)__vlibm_TBL_rsqrt + ind1))[0];
 345                 dexp_hi2 = ((double*)((char*)__vlibm_TBL_rsqrt + ind2))[0];
 346                 dexp_lo0 = ((double*)((char*)__vlibm_TBL_rsqrt + ind0))[1];
 347                 dexp_lo1 = ((double*)((char*)__vlibm_TBL_rsqrt + ind1))[1];
 348                 dexp_lo2 = ((double*)((char*)__vlibm_TBL_rsqrt + ind2))[1];
 349                 xx0 = dexp_hi0 * dexp_hi0;
 350                 xx1 = dexp_hi1 * dexp_hi1;
 351                 xx2 = dexp_hi2 * dexp_hi2;
 352                 xx0 = (res0 - res_c0) * xx0;
 353                 xx1 = (res1 - res_c1) * xx1;
 354                 xx2 = (res2 - res_c2) * xx2;
 355                 res0 = (((((K6 * xx0 + K5) * xx0 + K4) * xx0 + K3) * xx0 + K2) * xx0 + K1) * xx0;
 356                 res1 = (((((K6 * xx1 + K5) * xx1 + K4) * xx1 + K3) * xx1 + K2) * xx1 + K1) * xx1;
 357                 res2 = (((((K6 * xx2 + K5) * xx2 + K4) * xx2 + K3) * xx2 + K2) * xx2 + K1) * xx2;
 358 
 359                 res0 = dexp_hi0 * res0 + dexp_lo0 + dexp_hi0;
 360                 res1 = dexp_hi1 * res1 + dexp_lo1 + dexp_hi1;
 361                 res2 = dexp_hi2 * res2 + dexp_lo2 + dexp_hi2;
 362 
 363                 HI(&dsqrt_exp0) = sqrt_exp0;
 364                 HI(&dsqrt_exp1) = sqrt_exp1;
 365                 HI(&dsqrt_exp2) = sqrt_exp2;
 366                 res0 *= dsqrt_exp0;
 367                 res1 *= dsqrt_exp1;
 368                 res2 *= dsqrt_exp2;
 369 
 370                 *py = res0;
 371                 py += stridey;
 372 
 373                 *py = res1;
 374                 py += stridey;
 375 
 376                 *py = res2;
 377                 py += stridey;
 378         }
 379 
 380         for(; n > 0 ; n--)
 381         {
 382                 hx0 = HI(px);
 383 
 384                 sqrt_exp0 = (0x5fe - (hx0 >> 21)) << 20;
 385                 ind0 = (((hx0 >> 10) & 0x7f8) + 8) & -16;
 386 
 387                 resh0 = (hx0 & 0x001fffff) | 0x3fe00000;
 388                 res_ch0 = (resh0 + 0x00002000) & 0x7fffc000;
 389                 HI(&res0) = resh0;
 390                 LO(&res0) = LO(px);
 391                 HI(&res_c0) = res_ch0;
 392                 LO(&res_c0) = 0;
 393 
 394                 px += stridex;
 395 
 396                 dexp_hi0 = ((double*)((char*)__vlibm_TBL_rsqrt + ind0))[0];
 397                 dexp_lo0 = ((double*)((char*)__vlibm_TBL_rsqrt + ind0))[1];
 398                 xx0 = dexp_hi0 * dexp_hi0;
 399                 xx0 = (res0 - res_c0) * xx0;
 400                 res0 = (((((K6 * xx0 + K5) * xx0 + K4) * xx0 + K3) * xx0 + K2) * xx0 + K1) * xx0;
 401 
 402                 res0 = dexp_hi0 * res0 + dexp_lo0 + dexp_hi0;
 403 
 404                 HI(&dsqrt_exp0) = sqrt_exp0;
 405                 LO(&dsqrt_exp0) = 0;
 406                 res0 *= dsqrt_exp0;
 407 
 408                 *py = res0;
 409                 py += stridey;
 410         }
 411 }