xref: /titanic_44/usr/src/lib/libmvec/common/__vhypotf.c (revision 48a4016cae8aa2b8b3d8b258eb22e0c781912bed)
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 }
209