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