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 }