usr/src/lib/libmvec/common/__vatan2f.c

   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 #ifdef __RESTRICT
  31 #define restrict _Restrict
  32 #else
  33 #define restrict
  34 #endif
  35
  36 extern const double __vlibm_TBL_atan1[];
  37
  38 static const double
  39 pio4    =  7.8539816339744827900e-01,
  40 pio2    =  1.5707963267948965580e+00,
  41 pi      =  3.1415926535897931160e+00;
  42
  43 static const float
  44 zero    =  0.0f,
  45 one     =  1.0f,
  46 q1      = -3.3333333333296428046e-01f,
  47 q2      =  1.9999999186853752618e-01f,
  48 twop24  =  16777216.0f;
  49
  50 void
  51 __vatan2f(int n, float * restrict y, int stridey, float * restrict x,
  52         int stridex, float * restrict z, int stridez)
  53 {
  54         float           x0, x1, x2, y0, y1, y2, *pz0 = 0, *pz1, *pz2;
  55         double          ah0, ah1, ah2;
  56         double          t0, t1, t2;
  57         double          sx0, sx1, sx2;
  58         double          sign0, sign1, sign2;
  59         int             i, k0 = 0, k1, k2, hx, sx, sy;
  60         int             hy0, hy1, hy2;
  61         float           base0 = 0.0, base1, base2;
  62         double          num0, num1, num2;
  63         double          den0, den1, den2;
  64         double          dx0, dx1, dx2;
  65         double          dy0, dy1, dy2;
  66         double          db0, db1, db2;
  67
  68         do
  69         {
  70 loop0:
  71                 hy0 = *(int*)y;
  72                 hx = *(int*)x;
  73                 sign0 = one;
  74                 sy = hy0 & 0x80000000;
  75                 hy0 &= ~0x80000000;
  76
  77                 sx = hx & 0x80000000;
  78                 hx &= ~0x80000000;
  79
  80                 if (hy0 > hx)
  81                 {
  82                         x0 = *y;
  83                         y0 = *x;
  84                         i = hx;
  85                         hx = hy0;
  86                         hy0 = i;
  87                         if (sy)
  88                         {
  89                                 x0 = -x0;
  90                                 sign0 = -sign0;
  91                         }
  92                         if (sx)
  93                         {
  94                                 y0 = -y0;
  95                                 ah0 = pio2;
  96                         }
  97                         else
  98                         {
  99                                 ah0 = -pio2;
 100                                 sign0 = -sign0;
 101                         }
 102                 }
 103                 else
 104                 {
 105                         y0 = *y;
 106                         x0 = *x;
 107                         if (sy)
 108                         {
 109                                 y0 = -y0;
 110                                 sign0 = -sign0;
 111                         }
 112                         if (sx)
 113                         {
 114                                 x0 = -x0;
 115                                 ah0 = -pi;
 116                                 sign0 = -sign0;
 117                         }
 118                         else
 119                                 ah0 = zero;
 120                 }
 121
 122                 if (hx >= 0x7f800000 || hx - hy0 >= 0x0c800000)
 123                 {
 124                         if (hx >= 0x7f800000)
 125                         {
 126                                 if (hx ^ 0x7f800000) /* nan */
 127                                         ah0 =  x0 + y0;
 128                                 else if (hy0 >= 0x7f800000)
 129                                         ah0 += pio4;
 130                         }
 131                         else if ((int) ah0 == 0)
 132                                 ah0 = y0 / x0;
 133                         *z = (sign0 == one) ? ah0 : -ah0;
 134 /* sign0*ah0 would change nan behavior relative to previous release */
 135                         x += stridex;
 136                         y += stridey;
 137                         z += stridez;
 138                         i = 0;
 139                         if (--n <= 0)
 140                                 break;
 141                         goto loop0;
 142                 }
 143                 if (hy0 < 0x00800000) {
 144                         if (hy0 == 0)
 145                         {
 146                                 *z = sign0 * (float) ah0;
 147                                 x += stridex;
 148                                 y += stridey;
 149                                 z += stridez;
 150                                 i = 0;
 151                                 if (--n <= 0)
 152                                         break;
 153                                 goto loop0;
 154                         }
 155                         y0 *= twop24; /* scale subnormal y */
 156                         x0 *= twop24; /* scale possibly subnormal x */
 157                         hy0 = *(int*)&y0;
 158                         hx = *(int*)&x0;
 159                 }
 160                 pz0 = z;
 161
 162                 k0 = (hy0 - hx + 0x3f800000) & 0xfff80000;
 163                 if (k0 >= 0x3C800000)          /* if |x| >= (1/64)... */
 164                 {
 165                         *(int*)&base0 = k0;
 166                         k0 = (k0 - 0x3C800000) >> 18; /* (index >> 19) << 1) */
 167                         k0 += 4;
 168                                 /* skip over 0,0,pi/2,pi/2 */
 169                 }
 170                 else                            /* |x| < 1/64 */
 171                 {
 172                         k0 = 0;
 173                         base0 = zero;
 174                 }
 175
 176                 x += stridex;
 177                 y += stridey;
 178                 z += stridez;
 179                 i = 1;
 180                 if (--n <= 0)
 181                         break;
 182
 183
 184 loop1:
 185                 hy1 = *(int*)y;
 186                 hx = *(int*)x;
 187                 sign1 = one;
 188                 sy = hy1 & 0x80000000;
 189                 hy1 &= ~0x80000000;
 190
 191                 sx = hx & 0x80000000;
 192                 hx &= ~0x80000000;
 193
 194                 if (hy1 > hx)
 195                 {
 196                         x1 = *y;
 197                         y1 = *x;
 198                         i = hx;
 199                         hx = hy1;
 200                         hy1 = i;
 201                         if (sy)
 202                         {
 203                                 x1 = -x1;
 204                                 sign1 = -sign1;
 205                         }
 206                         if (sx)
 207                         {
 208                                 y1 = -y1;
 209                                 ah1 = pio2;
 210                         }
 211                         else
 212                         {
 213                                 ah1 = -pio2;
 214                                 sign1 = -sign1;
 215                         }
 216                 }
 217                 else
 218                 {
 219                         y1 = *y;
 220                         x1 = *x;
 221                         if (sy)
 222                         {
 223                                 y1 = -y1;
 224                                 sign1 = -sign1;
 225                         }
 226                         if (sx)
 227                         {
 228                                 x1 = -x1;
 229                                 ah1 = -pi;
 230                                 sign1 = -sign1;
 231                         }
 232                         else
 233                                 ah1 = zero;
 234                 }
 235
 236                 if (hx >= 0x7f800000 || hx - hy1 >= 0x0c800000)
 237                 {
 238                         if (hx >= 0x7f800000)
 239                         {
 240                                 if (hx ^ 0x7f800000) /* nan */
 241                                         ah1 =  x1 + y1;
 242                                 else if (hy1 >= 0x7f800000)
 243                                         ah1 += pio4;
 244                         }
 245                         else if ((int) ah1 == 0)
 246                                 ah1 = y1 / x1;
 247                         *z = (sign1 == one)? ah1 : -ah1;
 248                         x += stridex;
 249                         y += stridey;
 250                         z += stridez;
 251                         i = 1;
 252                         if (--n <= 0)
 253                                 break;
 254                         goto loop1;
 255                 }
 256                 if (hy1 < 0x00800000) {
 257                         if (hy1 == 0)
 258                         {
 259                                 *z = sign1 * (float) ah1;
 260                                 x += stridex;
 261                                 y += stridey;
 262                                 z += stridez;
 263                                 i = 1;
 264                                 if (--n <= 0)
 265                                         break;
 266                                 goto loop1;
 267                         }
 268                         y1 *= twop24; /* scale subnormal y */
 269                         x1 *= twop24; /* scale possibly subnormal x */
 270                         hy1 = *(int*)&y1;
 271                         hx = *(int*)&x1;
 272                 }
 273                 pz1 = z;
 274
 275                 k1 = (hy1 - hx + 0x3f800000) & 0xfff80000;
 276                 if (k1 >= 0x3C800000)          /* if |x| >= (1/64)... */
 277                 {
 278                         *(int*)&base1 = k1;
 279                         k1 = (k1 - 0x3C800000) >> 18; /* (index >> 19) << 1) */
 280                         k1 += 4;
 281                                 /* skip over 0,0,pi/2,pi/2 */
 282                 }
 283                 else                            /* |x| < 1/64 */
 284                 {
 285                         k1 = 0;
 286                         base1 = zero;
 287                 }
 288
 289                 x += stridex;
 290                 y += stridey;
 291                 z += stridez;
 292                 i = 2;
 293                 if (--n <= 0)
 294                         break;
 295
 296 loop2:
 297                 hy2 = *(int*)y;
 298                 hx = *(int*)x;
 299                 sign2 = one;
 300                 sy = hy2 & 0x80000000;
 301                 hy2 &= ~0x80000000;
 302
 303                 sx = hx & 0x80000000;
 304                 hx &= ~0x80000000;
 305
 306                 if (hy2 > hx)
 307                 {
 308                         x2 = *y;
 309                         y2 = *x;
 310                         i = hx;
 311                         hx = hy2;
 312                         hy2 = i;
 313                         if (sy)
 314                         {
 315                                 x2 = -x2;
 316                                 sign2 = -sign2;
 317                         }
 318                         if (sx)
 319                         {
 320                                 y2 = -y2;
 321                                 ah2 = pio2;
 322                         }
 323                         else
 324                         {
 325                                 ah2 = -pio2;
 326                                 sign2 = -sign2;
 327                         }
 328                 }
 329                 else
 330                 {
 331                         y2 = *y;
 332                         x2 = *x;
 333                         if (sy)
 334                         {
 335                                 y2 = -y2;
 336                                 sign2 = -sign2;
 337                         }
 338                         if (sx)
 339                         {
 340                                 x2 = -x2;
 341                                 ah2 = -pi;
 342                                 sign2 = -sign2;
 343                         }
 344                         else
 345                                 ah2 = zero;
 346                 }
 347
 348                 if (hx >= 0x7f800000 || hx - hy2 >= 0x0c800000)
 349                 {
 350                         if (hx >= 0x7f800000)
 351                         {
 352                                 if (hx ^ 0x7f800000) /* nan */
 353                                         ah2 =  x2 + y2;
 354                                 else if (hy2 >= 0x7f800000)
 355                                         ah2 += pio4;
 356                         }
 357                         else if ((int) ah2 == 0)
 358                                 ah2 = y2 / x2;
 359                         *z = (sign2 == one)? ah2 : -ah2;
 360                         x += stridex;
 361                         y += stridey;
 362                         z += stridez;
 363                         i = 2;
 364                         if (--n <= 0)
 365                                 break;
 366                         goto loop2;
 367                 }
 368                 if (hy2 < 0x00800000) {
 369                         if (hy2 == 0)
 370                         {
 371                                 *z = sign2 * (float) ah2;
 372                                 x += stridex;
 373                                 y += stridey;
 374                                 z += stridez;
 375                                 i = 2;
 376                                 if (--n <= 0)
 377                                         break;
 378                                 goto loop2;
 379                         }
 380                         y2 *= twop24; /* scale subnormal y */
 381                         x2 *= twop24; /* scale possibly subnormal x */
 382                         hy2 = *(int*)&y2;
 383                         hx = *(int*)&x2;
 384                 }
 385
 386                 pz2 = z;
 387
 388                 k2 = (hy2 - hx + 0x3f800000) & 0xfff80000;
 389                 if (k2 >= 0x3C800000)          /* if |x| >= (1/64)... */
 390                 {
 391                         *(int*)&base2 = k2;
 392                         k2 = (k2 - 0x3C800000) >> 18; /* (index >> 19) << 1) */
 393                         k2 += 4;
 394                                 /* skip over 0,0,pi/2,pi/2 */
 395                 }
 396                 else                            /* |x| < 1/64 */
 397                 {
 398                         k2 = 0;
 399                         base2 = zero;
 400                 }
 401
 402                 goto endloop;
 403
 404 endloop:
 405
 406                 ah2 += __vlibm_TBL_atan1[k2];
 407                 ah1 += __vlibm_TBL_atan1[k1];
 408                 ah0 += __vlibm_TBL_atan1[k0];
 409
 410                 db2 = base2;
 411                 db1 = base1;
 412                 db0 = base0;
 413                 dy2 = y2;
 414                 dy1 = y1;
 415                 dy0 = y0;
 416                 dx2 = x2;
 417                 dx1 = x1;
 418                 dx0 = x0;
 419
 420                 num2 = dy2 - dx2 * db2;
 421                 den2 = dx2 + dy2 * db2;
 422
 423                 num1 = dy1 - dx1 * db1;
 424                 den1 = dx1 + dy1 * db1;
 425
 426                 num0 = dy0 - dx0 * db0;
 427                 den0 = dx0 + dy0 * db0;
 428
 429                 t2 = num2 / den2;
 430                 t1 = num1 / den1;
 431                 t0 = num0 / den0;
 432
 433                 sx2 = t2 * t2;
 434                 sx1 = t1 * t1;
 435                 sx0 = t0 * t0;
 436
 437                 t2 += t2 * sx2 * (q1 + sx2 * q2);
 438                 t1 += t1 * sx1 * (q1 + sx1 * q2);
 439                 t0 += t0 * sx0 * (q1 + sx0 * q2);
 440
 441                 t2 += ah2;
 442                 t1 += ah1;
 443                 t0 += ah0;
 444
 445                 *pz2 = sign2 * t2;
 446                 *pz1 = sign1 * t1;
 447                 *pz0 = sign0 * t0;
 448
 449                 x += stridex;
 450                 y += stridey;
 451                 z += stridez;
 452                 i = 0;
 453         } while (--n > 0);
 454
 455         if (i > 1)
 456         {
 457                 ah1 += __vlibm_TBL_atan1[k1];
 458                 t1 = (y1 - x1 * (double)base1) /
 459                         (x1 + y1 * (double)base1);
 460                 sx1 = t1 * t1;
 461                 t1 += t1 * sx1 * (q1 + sx1 * q2);
 462                 t1 += ah1;
 463                 *pz1 = sign1 * t1;
 464         }
 465
 466         if (i > 0)
 467         {
 468                 ah0 += __vlibm_TBL_atan1[k0];
 469                 t0 = (y0 - x0 * (double)base0) /
 470                         (x0 + y0 * (double)base0);
 471                 sx0 = t0 * t0;
 472                 t0 += t0 * sx0 * (q1 + sx0 * q2);
 473                 t0 += ah0;
 474                 *pz0 = sign0 * t0;
 475         }
 476 }