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