usr/src/lib/libmvec/common/__vhypotf.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 #include "libm_inlines.h"
  31
  32 #ifdef __RESTRICT
  33 #define restrict _Restrict
  34 #else
  35 #define restrict
  36 #endif
  37
  38 extern double sqrt(double);
  39
  40 void
  41 __vhypotf(int n, float * restrict x, int stridex, float * restrict y,
  42         int stridey, float * restrict z, int stridez)
  43 {
  44         float           x0, x1, x2, y0, y1, y2, z0, z1, z2, *pz0, *pz1, *pz2;
  45         unsigned        hx0, hx1, hx2, hy0, hy1, hy2;
  46         int                     i, j0, j1, j2;
  47
  48         do
  49         {
  50 LOOP0:
  51                 hx0 = *(unsigned*)x & ~0x80000000;
  52                 hy0 = *(unsigned*)y & ~0x80000000;
  53                 *(unsigned*)&x0 = hx0;
  54                 *(unsigned*)&y0 = hy0;
  55                 if (hy0 > hx0)
  56                 {
  57                         i = hy0 - hx0;
  58                         j0 = hy0 & 0x7f800000;
  59                         if (hx0 == 0)
  60                                 i = 0x7f800000;
  61                 }
  62                 else
  63                 {
  64                         i = hx0 - hy0;
  65                         j0 = hx0 & 0x7f800000;
  66                         if (hy0 == 0)
  67                                 i = 0x7f800000;
  68                         else if (hx0 == 0)
  69                                 i = 0x7f800000;
  70                 }
  71                 if (i >= 0x0c800000 || j0 >= 0x7f800000)
  72                 {
  73                         z0 = x0 + y0;
  74                         if (hx0 == 0x7f800000)
  75                                 z0 = x0;
  76                         else if (hy0  == 0x7f800000)
  77                                 z0 = y0;
  78                         else if (hx0 > 0x7f800000 || hy0 > 0x7f800000)
  79                                 z0 = *x + *y;
  80                         *z = z0;
  81                         x += stridex;
  82                         y += stridey;
  83                         z += stridez;
  84                         i = 0;
  85                         if (--n <= 0)
  86                                 break;
  87                         goto LOOP0;
  88                 }
  89                 pz0 = z;
  90                 x += stridex;
  91                 y += stridey;
  92                 z += stridez;
  93                 i = 1;
  94                 if (--n <= 0)
  95                         break;
  96
  97 LOOP1:
  98                 hx1 = *(unsigned*)x & ~0x80000000;
  99                 hy1 = *(unsigned*)y & ~0x80000000;
 100                 *(unsigned*)&x1 = hx1;
 101                 *(unsigned*)&y1 = hy1;
 102                 if (hy1 > hx1)
 103                 {
 104                         i = hy1 - hx1;
 105                         j1 = hy1 & 0x7f800000;
 106                         if (hx1 == 0)
 107                                 i = 0x7f800000;
 108                 }
 109                 else
 110                 {
 111                         i = hx1 - hy1;
 112                         j1 = hx1 & 0x7f800000;
 113                         if (hy1 == 0)
 114                                 i = 0x7f800000;
 115                         else if (hx1 == 0)
 116                                 i = 0x7f800000;
 117                 }
 118                 if (i >= 0x0c800000 || j1 >= 0x7f800000)
 119                 {
 120                         z1 = x1 + y1;
 121                         if (hx1 == 0x7f800000)
 122                                 z1 = x1;
 123                         else if (hy1 == 0x7f800000)
 124                                 z1 = y1;
 125                         else if (hx1 > 0x7f800000 || hy1 > 0x7f800000)
 126                                 z1 = *x + *y;
 127                         *z = z1;
 128                         x += stridex;
 129                         y += stridey;
 130                         z += stridez;
 131                         i = 1;
 132                         if (--n <= 0)
 133                                 break;
 134                         goto LOOP1;
 135                 }
 136                 pz1 = z;
 137                 x += stridex;
 138                 y += stridey;
 139                 z += stridez;
 140                 i = 2;
 141                 if (--n <= 0)
 142                         break;
 143
 144 LOOP2:
 145                 hx2 = *(unsigned*)x & ~0x80000000;
 146                 hy2 = *(unsigned*)y & ~0x80000000;
 147                 *(unsigned*)&x2 = hx2;
 148                 *(unsigned*)&y2 = hy2;
 149                 if (hy2 > hx2)
 150                 {
 151                         i = hy2 - hx2;
 152                         j2 = hy2 & 0x7f800000;
 153                         if (hx2 == 0)
 154                                 i = 0x7f800000;
 155                 }
 156                 else
 157                 {
 158                         i = hx2 - hy2;
 159                         j2 = hx2 & 0x7f800000;
 160                         if (hy2 == 0)
 161                                 i = 0x7f800000;
 162                         else if (hx2 == 0)
 163                                 i = 0x7f800000;
 164                 }
 165                 if (i >= 0x0c800000 || j2 >= 0x7f800000)
 166                 {
 167                         z2 = x2 + y2;
 168                         if (hx2 == 0x7f800000)
 169                                 z2 = x2;
 170                         else if (hy2 == 0x7f800000)
 171                                 z2 = y2;
 172                         else if (hx2 > 0x7f800000 || hy2 > 0x7f800000)
 173                                 z2 = *x + *y;
 174                         *z = z2;
 175                         x += stridex;
 176                         y += stridey;
 177                         z += stridez;
 178                         i = 2;
 179                         if (--n <= 0)
 180                                 break;
 181                         goto LOOP2;
 182                 }
 183                 pz2 = z;
 184
 185                 z0 = sqrt(x0 * (double)x0 + y0 * (double)y0);
 186                 z1 = sqrt(x1 * (double)x1 + y1 * (double)y1);
 187                 z2 = sqrt(x2 * (double)x2 + y2 * (double)y2);
 188                 *pz0 = z0;
 189                 *pz1 = z1;
 190                 *pz2 = z2;
 191
 192                 x += stridex;
 193                 y += stridey;
 194                 z += stridez;
 195                 i = 0;
 196         } while (--n > 0);
 197
 198         if (i > 0)
 199         {
 200                 if (i > 1)
 201                 {
 202                         z1 = sqrt(x1 * (double)x1 + y1 * (double)y1);
 203                         *pz1 = z1;
 204                 }
 205                 z0 = sqrt(x0 * (double)x0 + y0 * (double)y0);
 206                 *pz0 = z0;
 207         }
 208 }