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 }