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_synonyms.h"
  32 #include "libm_inlines.h"
  33 
  34 #ifdef _LITTLE_ENDIAN
  35 #define HI(x)   *(1+(int*)x)
  36 #define LO(x)   *(unsigned*)x
  37 #else
  38 #define HI(x)   *(int*)x
  39 #define LO(x)   *(1+(unsigned*)x)
  40 #endif
  41 
  42 #ifdef __RESTRICT
  43 #define restrict _Restrict
  44 #else
  45 #define restrict
  46 #endif
  47 
  48 /* double rhypot(double x, double y)
  49  *
  50  * Method :
  51  *      1. Special cases:
  52  *              x or  y = Inf                                   => 0
  53  *              x or  y = NaN                                   => QNaN
  54  *              x and y = 0                                     => Inf + divide-by-zero
  55  *      2. Computes rhypot(x,y):
  56  *              rhypot(x,y) = m * sqrt(1/(xnm * xnm + ynm * ynm))
  57  *      Where:
  58  *              m = 1/max(|x|,|y|)
  59  *              xnm = x * m
  60  *              ynm = y * m
  61  *
  62  *      Compute 1/(xnm * xnm + ynm * ynm) by simulating
  63  *      muti-precision arithmetic.
  64  *
  65  * Accuracy:
  66  *      Maximum error observed: less than 0.869 ulp after 1.000.000.000
  67  *      results.
  68  */
  69 
  70 #define sqrt __sqrt
  71 
  72 extern double sqrt( double );
  73 
  74 extern double fabs( double );
  75 
  76 static const int __vlibm_TBL_rhypot[] = {
  77 /* i = [0,127]
  78  * TBL[i] = 0x3ff00000 + *(int*)&(1.0 / *(double*)&(0x3ff0000000000000ULL + (i << 45))); */
  79  0x7fe00000,  0x7fdfc07f, 0x7fdf81f8,  0x7fdf4465,
  80  0x7fdf07c1,  0x7fdecc07, 0x7fde9131,  0x7fde573a,
  81  0x7fde1e1e,  0x7fdde5d6, 0x7fddae60,  0x7fdd77b6,
  82  0x7fdd41d4,  0x7fdd0cb5, 0x7fdcd856,  0x7fdca4b3,
  83  0x7fdc71c7,  0x7fdc3f8f, 0x7fdc0e07,  0x7fdbdd2b,
  84  0x7fdbacf9,  0x7fdb7d6c, 0x7fdb4e81,  0x7fdb2036,
  85  0x7fdaf286,  0x7fdac570, 0x7fda98ef,  0x7fda6d01,
  86  0x7fda41a4,  0x7fda16d3, 0x7fd9ec8e,  0x7fd9c2d1,
  87  0x7fd99999,  0x7fd970e4, 0x7fd948b0,  0x7fd920fb,
  88  0x7fd8f9c1,  0x7fd8d301, 0x7fd8acb9,  0x7fd886e5,
  89  0x7fd86186,  0x7fd83c97, 0x7fd81818,  0x7fd7f405,
  90  0x7fd7d05f,  0x7fd7ad22, 0x7fd78a4c,  0x7fd767dc,
  91  0x7fd745d1,  0x7fd72428, 0x7fd702e0,  0x7fd6e1f7,
  92  0x7fd6c16c,  0x7fd6a13c, 0x7fd68168,  0x7fd661ec,
  93  0x7fd642c8,  0x7fd623fa, 0x7fd60581,  0x7fd5e75b,
  94  0x7fd5c988,  0x7fd5ac05, 0x7fd58ed2,  0x7fd571ed,
  95  0x7fd55555,  0x7fd53909, 0x7fd51d07,  0x7fd50150,
  96  0x7fd4e5e0,  0x7fd4cab8, 0x7fd4afd6,  0x7fd49539,
  97  0x7fd47ae1,  0x7fd460cb, 0x7fd446f8,  0x7fd42d66,
  98  0x7fd41414,  0x7fd3fb01, 0x7fd3e22c,  0x7fd3c995,
  99  0x7fd3b13b,  0x7fd3991c, 0x7fd38138,  0x7fd3698d,
 100  0x7fd3521c,  0x7fd33ae4, 0x7fd323e3,  0x7fd30d19,
 101  0x7fd2f684,  0x7fd2e025, 0x7fd2c9fb,  0x7fd2b404,
 102  0x7fd29e41,  0x7fd288b0, 0x7fd27350,  0x7fd25e22,
 103  0x7fd24924,  0x7fd23456, 0x7fd21fb7,  0x7fd20b47,
 104  0x7fd1f704,  0x7fd1e2ef, 0x7fd1cf06,  0x7fd1bb4a,
 105  0x7fd1a7b9,  0x7fd19453, 0x7fd18118,  0x7fd16e06,
 106  0x7fd15b1e,  0x7fd1485f, 0x7fd135c8,  0x7fd12358,
 107  0x7fd11111,  0x7fd0fef0, 0x7fd0ecf5,  0x7fd0db20,
 108  0x7fd0c971,  0x7fd0b7e6, 0x7fd0a681,  0x7fd0953f,
 109  0x7fd08421,  0x7fd07326, 0x7fd0624d,  0x7fd05197,
 110  0x7fd04104,  0x7fd03091, 0x7fd02040,  0x7fd01010,
 111 };
 112 
 113 static const unsigned long long LCONST[] = {
 114 0x3ff0000000000000ULL,  /* DONE = 1.0           */
 115 0x4000000000000000ULL,  /* DTWO = 2.0           */
 116 0x4230000000000000ULL,  /* D2ON36 = 2**36       */
 117 0x7fd0000000000000ULL,  /* D2ON1022 = 2**1022   */
 118 0x3cb0000000000000ULL,  /* D2ONM52 = 2**-52     */
 119 };
 120 
 121 #define RET_SC(I)                                                                               \
 122         px += stridex;                                                                          \
 123         py += stridey;                                                                          \
 124         pz += stridez;                                                                          \
 125         if ( --n <= 0 )                                                                              \
 126                 break;                                                                          \
 127         goto start##I;
 128 
 129 #define RETURN(I, ret)                                                                          \
 130 {                                                                                               \
 131         pz[0] = (ret);                                                                          \
 132         RET_SC(I)                                                                               \
 133 }
 134 
 135 #define PREP(I)                                                                                 \
 136 hx##I = HI(px);                                                                         \
 137 hy##I = HI(py);                                                                         \
 138 hx##I &= 0x7fffffff;                                                                                \
 139 hy##I &= 0x7fffffff;                                                                                \
 140 pz##I = pz;                                                                                     \
 141 if ( hx##I >= 0x7ff00000 || hy##I >= 0x7ff00000 ) /* |X| or |Y| = Inf,NaN */              \
 142 {                                                                                               \
 143         lx = LO(px);                                                                    \
 144         ly = LO(py);                                                                    \
 145         x = *px;                                                                                \
 146         y = *py;                                                                                \
 147         if ( hx##I == 0x7ff00000 && lx == 0 ) res0 = 0.0;               /* |X| = Inf */         \
 148         else if ( hy##I == 0x7ff00000 && ly == 0 ) res0 = 0.0;  /* |Y| = Inf */         \
 149         else res0 = fabs(x) + fabs(y);                                                          \
 150                                                                                                 \
 151         RETURN (I, res0)                                                                        \
 152 }                                                                                               \
 153 x##I = *px;                                                                                     \
 154 y##I = *py;                                                                                     \
 155 diff0 = hy##I - hx##I;                                                                          \
 156 j0 = diff0 >> 31;                                                                         \
 157 if ( hx##I < 0x00100000 && hy##I < 0x00100000 )   /* |X| and |Y| = subnormal or zero */           \
 158 {                                                                                               \
 159         lx = LO(px);                                                                    \
 160         ly = LO(py);                                                                    \
 161         x = x##I;                                                                               \
 162         y = y##I;                                                                               \
 163                                                                                                 \
 164         if ( (hx##I | hy##I | lx | ly) == 0 )   /* |X| and |Y| = 0 */                           \
 165                 RETURN (I, DONE / 0.0)                                                  \
 166                                                                                                 \
 167         x = fabs(x);                                                                            \
 168         y = fabs(y);                                                                            \
 169                                                                                                 \
 170         x = *(long long*)&x;                                                                        \
 171         y = *(long long*)&y;                                                                        \
 172                                                                                                 \
 173         x *= D2ONM52;                                                                           \
 174         y *= D2ONM52;                                                                           \
 175                                                                                                 \
 176         x_hi0 = ( x + D2ON36 ) - D2ON36;                                                        \
 177         y_hi0 = ( y + D2ON36 ) - D2ON36;                                                        \
 178         x_lo0 = x - x_hi0;                                                                      \
 179         y_lo0 = y - y_hi0;                                                                      \
 180         res0_hi = (x_hi0 * x_hi0 + y_hi0 * y_hi0);                                              \
 181         res0_lo = ((x + x_hi0) * x_lo0 + (y + y_hi0) * y_lo0);                                  \
 182                                                                                                 \
 183         dres0 = res0_hi + res0_lo;                                                              \
 184                                                                                                 \
 185         iarr0 = HI(&dres0);                                                         \
 186         iexp0 = iarr0 & 0xfff00000;                                                         \
 187                                                                                                 \
 188         iarr0 = (iarr0 >> 11) & 0x1fc;                                                                \
 189         itbl0 = ((int*)((char*)__vlibm_TBL_rhypot + iarr0))[0];                                 \
 190         itbl0 -= iexp0;                                                                         \
 191         HI(&dd0) = itbl0;                                           \
 192         LO(&dd0) = 0;                                                               \
 193                                                                                                 \
 194         dd0 = dd0 * (DTWO - dd0 * dres0);                                                       \
 195         dd0 = dd0 * (DTWO - dd0 * dres0);                                                       \
 196         dres0 = dd0 * (DTWO - dd0 * dres0);                                                     \
 197                                                                                                 \
 198         HI(&res0) = HI(&dres0) & 0xffffff00;                                        \
 199         LO(&res0) = 0;                                                              \
 200         res0 += (DONE - res0_hi * res0 - res0_lo * res0) * dres0;                               \
 201         res0 = sqrt ( res0 );                                                                   \
 202                                                                                                 \
 203         res0 = D2ON1022 * res0;                                                                 \
 204         RETURN (I, res0)                                                                        \
 205 }                                                                                               \
 206 j0 = hy##I - (diff0 & j0);                                                                  \
 207 j0 &= 0x7ff00000;                                                                           \
 208 HI(&scl##I) = 0x7ff00000 - j0;
 209 
 210 void
 211 __vrhypot( int n, double * restrict px, int stridex, double * restrict py,
 212         int stridey, double * restrict pz, int stridez )
 213 {
 214         int             i = 0;
 215         double          x, y;
 216         double          x_hi0, x_lo0, y_hi0, y_lo0, scl0 = 0;
 217         double          x0, y0, res0, dd0;
 218         double          res0_hi,res0_lo, dres0;
 219         double          x_hi1, x_lo1, y_hi1, y_lo1, scl1 = 0;
 220         double          x1, y1, res1, dd1;
 221         double          res1_hi,res1_lo, dres1;
 222         double          x_hi2, x_lo2, y_hi2, y_lo2, scl2 = 0;
 223         double          x2, y2, res2, dd2;
 224         double          res2_hi,res2_lo, dres2;
 225 
 226         int             hx0, hy0, j0, diff0;
 227         int             iarr0, iexp0, itbl0;
 228         int             hx1, hy1;
 229         int             iarr1, iexp1, itbl1;
 230         int             hx2, hy2;
 231         int             iarr2, iexp2, itbl2;
 232 
 233         int             lx, ly;
 234 
 235         double          DONE = ((double*)LCONST)[0];
 236         double          DTWO = ((double*)LCONST)[1];
 237         double          D2ON36 = ((double*)LCONST)[2];
 238         double          D2ON1022 = ((double*)LCONST)[3];
 239         double          D2ONM52 = ((double*)LCONST)[4];
 240 
 241         double          *pz0, *pz1, *pz2;
 242 
 243         do
 244         {
 245 start0:
 246                 PREP(0)
 247                 px += stridex;
 248                 py += stridey;
 249                 pz += stridez;
 250                 i = 1;
 251                 if ( --n <= 0 )
 252                         break;
 253 
 254 start1:
 255                 PREP(1)
 256                 px += stridex;
 257                 py += stridey;
 258                 pz += stridez;
 259                 i = 2;
 260                 if ( --n <= 0 )
 261                         break;
 262 
 263 start2:
 264                 PREP(2)
 265 
 266                 x0 *= scl0;
 267                 y0 *= scl0;
 268                 x1 *= scl1;
 269                 y1 *= scl1;
 270                 x2 *= scl2;
 271                 y2 *= scl2;
 272 
 273                 x_hi0 = ( x0 + D2ON36 ) - D2ON36;
 274                 y_hi0 = ( y0 + D2ON36 ) - D2ON36;
 275                 x_hi1 = ( x1 + D2ON36 ) - D2ON36;
 276                 y_hi1 = ( y1 + D2ON36 ) - D2ON36;
 277                 x_hi2 = ( x2 + D2ON36 ) - D2ON36;
 278                 y_hi2 = ( y2 + D2ON36 ) - D2ON36;
 279                 x_lo0 = x0 - x_hi0;
 280                 y_lo0 = y0 - y_hi0;
 281                 x_lo1 = x1 - x_hi1;
 282                 y_lo1 = y1 - y_hi1;
 283                 x_lo2 = x2 - x_hi2;
 284                 y_lo2 = y2 - y_hi2;
 285                 res0_hi = (x_hi0 * x_hi0 + y_hi0 * y_hi0);
 286                 res1_hi = (x_hi1 * x_hi1 + y_hi1 * y_hi1);
 287                 res2_hi = (x_hi2 * x_hi2 + y_hi2 * y_hi2);
 288                 res0_lo = ((x0 + x_hi0) * x_lo0 + (y0 + y_hi0) * y_lo0);
 289                 res1_lo = ((x1 + x_hi1) * x_lo1 + (y1 + y_hi1) * y_lo1);
 290                 res2_lo = ((x2 + x_hi2) * x_lo2 + (y2 + y_hi2) * y_lo2);
 291 
 292                 dres0 = res0_hi + res0_lo;
 293                 dres1 = res1_hi + res1_lo;
 294                 dres2 = res2_hi + res2_lo;
 295 
 296                 iarr0 = HI(&dres0);
 297                 iarr1 = HI(&dres1);
 298                 iarr2 = HI(&dres2);
 299                 iexp0 = iarr0 & 0xfff00000;
 300                 iexp1 = iarr1 & 0xfff00000;
 301                 iexp2 = iarr2 & 0xfff00000;
 302 
 303                 iarr0 = (iarr0 >> 11) & 0x1fc;
 304                 iarr1 = (iarr1 >> 11) & 0x1fc;
 305                 iarr2 = (iarr2 >> 11) & 0x1fc;
 306                 itbl0 = ((int*)((char*)__vlibm_TBL_rhypot + iarr0))[0];
 307                 itbl1 = ((int*)((char*)__vlibm_TBL_rhypot + iarr1))[0];
 308                 itbl2 = ((int*)((char*)__vlibm_TBL_rhypot + iarr2))[0];
 309                 itbl0 -= iexp0;
 310                 itbl1 -= iexp1;
 311                 itbl2 -= iexp2;
 312                 HI(&dd0) = itbl0;
 313                 HI(&dd1) = itbl1;
 314                 HI(&dd2) = itbl2;
 315                 LO(&dd0) = 0;
 316                 LO(&dd1) = 0;
 317                 LO(&dd2) = 0;
 318 
 319                 dd0 = dd0 * (DTWO - dd0 * dres0);
 320                 dd1 = dd1 * (DTWO - dd1 * dres1);
 321                 dd2 = dd2 * (DTWO - dd2 * dres2);
 322                 dd0 = dd0 * (DTWO - dd0 * dres0);
 323                 dd1 = dd1 * (DTWO - dd1 * dres1);
 324                 dd2 = dd2 * (DTWO - dd2 * dres2);
 325                 dres0 = dd0 * (DTWO - dd0 * dres0);
 326                 dres1 = dd1 * (DTWO - dd1 * dres1);
 327                 dres2 = dd2 * (DTWO - dd2 * dres2);
 328 
 329                 HI(&res0) = HI(&dres0) & 0xffffff00;
 330                 HI(&res1) = HI(&dres1) & 0xffffff00;
 331                 HI(&res2) = HI(&dres2) & 0xffffff00;
 332                 LO(&res0) = 0;
 333                 LO(&res1) = 0;
 334                 LO(&res2) = 0;
 335                 res0 += (DONE - res0_hi * res0 - res0_lo * res0) * dres0;
 336                 res1 += (DONE - res1_hi * res1 - res1_lo * res1) * dres1;
 337                 res2 += (DONE - res2_hi * res2 - res2_lo * res2) * dres2;
 338                 res0 = sqrt ( res0 );
 339                 res1 = sqrt ( res1 );
 340                 res2 = sqrt ( res2 );
 341 
 342                 res0 = scl0 * res0;
 343                 res1 = scl1 * res1;
 344                 res2 = scl2 * res2;
 345 
 346                 *pz0 = res0;
 347                 *pz1 = res1;
 348                 *pz2 = res2;
 349 
 350                 px += stridex;
 351                 py += stridey;
 352                 pz += stridez;
 353                 i = 0;
 354 
 355         } while ( --n > 0 );
 356 
 357         if ( i > 0 )
 358         {
 359                 x0 *= scl0;
 360                 y0 *= scl0;
 361 
 362                 x_hi0 = ( x0 + D2ON36 ) - D2ON36;
 363                 y_hi0 = ( y0 + D2ON36 ) - D2ON36;
 364                 x_lo0 = x0 - x_hi0;
 365                 y_lo0 = y0 - y_hi0;
 366                 res0_hi = (x_hi0 * x_hi0 + y_hi0 * y_hi0);
 367                 res0_lo = ((x0 + x_hi0) * x_lo0 + (y0 + y_hi0) * y_lo0);
 368 
 369                 dres0 = res0_hi + res0_lo;
 370 
 371                 iarr0 = HI(&dres0);
 372                 iexp0 = iarr0 & 0xfff00000;
 373 
 374                 iarr0 = (iarr0 >> 11) & 0x1fc;
 375                 itbl0 = ((int*)((char*)__vlibm_TBL_rhypot + iarr0))[0];
 376                 itbl0 -= iexp0;
 377                 HI(&dd0) = itbl0;
 378                 LO(&dd0) = 0;
 379 
 380                 dd0 = dd0 * (DTWO - dd0 * dres0);
 381                 dd0 = dd0 * (DTWO - dd0 * dres0);
 382                 dres0 = dd0 * (DTWO - dd0 * dres0);
 383 
 384                 HI(&res0) = HI(&dres0) & 0xffffff00;
 385                 LO(&res0) = 0;
 386                 res0 += (DONE - res0_hi * res0 - res0_lo * res0) * dres0;
 387                 res0 = sqrt ( res0 );
 388 
 389                 res0 = scl0 * res0;
 390 
 391                 *pz0 = res0;
 392 
 393                 if ( i > 1 )
 394                 {
 395                         x1 *= scl1;
 396                         y1 *= scl1;
 397 
 398                         x_hi1 = ( x1 + D2ON36 ) - D2ON36;
 399                         y_hi1 = ( y1 + D2ON36 ) - D2ON36;
 400                         x_lo1 = x1 - x_hi1;
 401                         y_lo1 = y1 - y_hi1;
 402                         res1_hi = (x_hi1 * x_hi1 + y_hi1 * y_hi1);
 403                         res1_lo = ((x1 + x_hi1) * x_lo1 + (y1 + y_hi1) * y_lo1);
 404 
 405                         dres1 = res1_hi + res1_lo;
 406 
 407                         iarr1 = HI(&dres1);
 408                         iexp1 = iarr1 & 0xfff00000;
 409 
 410                         iarr1 = (iarr1 >> 11) & 0x1fc;
 411                         itbl1 = ((int*)((char*)__vlibm_TBL_rhypot + iarr1))[0];
 412                         itbl1 -= iexp1;
 413                         HI(&dd1) = itbl1;
 414                         LO(&dd1) = 0;
 415 
 416                         dd1 = dd1 * (DTWO - dd1 * dres1);
 417                         dd1 = dd1 * (DTWO - dd1 * dres1);
 418                         dres1 = dd1 * (DTWO - dd1 * dres1);
 419 
 420                         HI(&res1) = HI(&dres1) & 0xffffff00;
 421                         LO(&res1) = 0;
 422                         res1 += (DONE - res1_hi * res1 - res1_lo * res1) * dres1;
 423                         res1 = sqrt ( res1 );
 424 
 425                         res1 = scl1 * res1;
 426 
 427                         *pz1 = res1;
 428                 }
 429         }
 430 }
 431