```/* <![CDATA[ */
function get_sym_list(){return [["Macro","xm",[["INF",103],["PREPROCESS",116],["PROCESS",105],["restrict",31],["restrict",33]]],["Variable","xv",[["K1",96],["K2",95],["K3",94],["__TBL_logf",66],["inf",101]]],["Function","xf",[["__vlogf",144]]]];} /* ]]> */1 /*
3  *
4  * The contents of this file are subject to the terms of the
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
15  * If applicable, add the following below this CDDL HEADER, with the
16  * fields enclosed by brackets "[]" replaced with your own identifying
18  *
20  */
21
22 /*
24  */
25 /*
27  * Use is subject to license terms.
28  */
29
30 #ifdef __RESTRICT
31 #define restrict _Restrict
32 #else
33 #define restrict
34 #endif
35
36 /* float logf(float x)
37  *
38  * Method :
39  *	1. Special cases:
40  *		for x is negative, -Inf => QNaN + invalid;
41  *		for x = 0		=> -Inf + divide-by-zero;
42  *		for x = +Inf		=> Inf;
43  *		for x = NaN		=> QNaN.
44  *	2. Computes logarithm from:
45  *		x = m * 2**n => log(x) = n * log(2) + log(m),
46  *		m = [1, 2).
47  *	Let m = m0 + dm, where m0 = 1 + k / 32,
48  *		k = [0, 32],
49  *		dm = [-1/64, 1/64].
50  *	Then log(m) = log(m0 + dm) = log(m0) + log(1+y),
51  *		where y = dm*(1/m0), y = [-1/66, 1/64].
52  *	Then
53  *		1/m0 is looked up in a table of 1, 1/(1+1/32), ..., 1/(1+32/32);
54  *		log(m0) is looked up in a table of log(1), log(1+1/32),
55  *		..., log(1+32/32).
56  *		log(1+y) is computed using approximation:
57  *		log(1+y) = ((a3*y + a2)*y + a1)*y*y + y.
58  * Accuracy:
59  *	The maximum relative error for the approximating
60  *	polynomial is 2**(-28.41).  All calculations are of
61  *	double precision.
62  *	Maximum error observed: less than 0.545 ulp for the
63  *	whole float type range.
64  */
65
66 static const double __TBL_logf[] = {
67 	/* __TBL_logf[2*i] = log(1+i/32), i = [0, 32] */
68 	/* __TBL_logf[2*i+1] = 2**(-23)/(1+i/32), i = [0, 32] */
69 0.000000000000000000e+00, 1.192092895507812500e-07, 3.077165866675368733e-02,
70 1.155968868371212153e-07, 6.062462181643483994e-02, 1.121969784007352926e-07,
71 8.961215868968713805e-02, 1.089913504464285680e-07, 1.177830356563834557e-01,
72 1.059638129340277719e-07, 1.451820098444978890e-01, 1.030999260979729787e-07,
73 1.718502569266592284e-01, 1.003867701480263102e-07, 1.978257433299198675e-01,
74 9.781275040064102225e-08, 2.231435513142097649e-01, 9.536743164062500529e-08,
75 2.478361639045812692e-01, 9.304139672256097884e-08, 2.719337154836417580e-01,
76 9.082612537202380448e-08, 2.954642128938358980e-01, 8.871388989825581272e-08,
77 3.184537311185345887e-01, 8.669766512784091150e-08, 3.409265869705931928e-01,
78 8.477105034722222546e-08, 3.629054936893684746e-01, 8.292820142663043248e-08,
79 3.844116989103320559e-01, 8.116377160904255122e-08, 4.054651081081643849e-01,
80 7.947285970052082892e-08, 4.260843953109000881e-01, 7.785096460459183052e-08,
81 4.462871026284195297e-01, 7.629394531250000159e-08, 4.660897299245992387e-01,
82 7.479798560049019504e-08, 4.855078157817008244e-01, 7.335956280048077330e-08,
83 5.045560107523953119e-01, 7.197542010613207272e-08, 5.232481437645478684e-01,
84 7.064254195601851460e-08, 5.415972824327444091e-01, 6.935813210227272390e-08,
85 5.596157879354226594e-01, 6.811959402901785336e-08, 5.773153650348236132e-01,
86 6.692451343201754014e-08, 5.947071077466927758e-01, 6.577064251077586116e-08,
87 6.118015411059929409e-01, 6.465588585805084723e-08, 6.286086594223740942e-01,
88 6.357828776041666578e-08, 6.451379613735847007e-01, 6.253602074795082293e-08,
89 6.613984822453650159e-01, 6.152737525201612732e-08, 6.773988235918061429e-01,
90 6.055075024801586965e-08, 6.931471805599452862e-01, 5.960464477539062500e-08
91 };
92
93 static const double
94 	K3 = -2.49887584306188944706e-01,
95 	K2 =  3.33368809981254554946e-01,
96 	K1 = -5.00000008402474976565e-01;
97
98 static const union {
99 	int	i;
100 	float	f;
101 } inf = { 0x7f800000 };
102
103 #define INF	inf.f
104
105 #define PROCESS(N)								\
106 	iy##N = ival##N & 0x007fffff;						\
107 	ival##N = (iy##N + 0x20000) & 0xfffc0000;				\
108 	i##N  = ival##N >> 17;							\
109 	iy##N = iy##N - ival##N;						\
110 	ty##N = LN2 * (double) exp##N + __TBL_logf[i##N];			\
111 	yy##N = (double) iy##N * __TBL_logf[i##N + 1];				\
112 	yy##N = ((K3 * yy##N + K2) * yy##N + K1) * yy##N * yy##N + yy##N;	\
113 	y[0] = (float)(yy##N + ty##N);						\
114 	y += stridey;
115
116 #define PREPROCESS(N, index, label)						\
117 	ival##N = *(int*)x;							\
118 	value = x[0];								\
119 	x += stridex;								\
120 	exp##N = (ival##N >> 23) - 127;						\
121 	if ((ival##N & 0x7fffffff) >= 0x7f800000) /* X = NaN or Inf */	\
122 	{									\
123 		y[index] = value + INF;						\
124 		goto label;							\
125 	}									\
126 	if (ival##N < 0x00800000)						\
127 	{									\
128 		if (ival##N > 0)	/* X = denormal */			\
129 		{								\
130 			value = (float) ival##N;				\
131 			ival##N = *(int*) &value;				\
132 			exp##N = (ival##N >> 23) - (127 + 149);			\
133 		}								\
134 		else								\
135 		{								\
136 			value = 0.0f;						\
137 			y[index] = ((ival##N & 0x7fffffff) == 0) ?		\
138 				-1.0f / value : value / value;			\
139 			goto label;						\
140 		}								\
141 	}
142
143 void
__vlogf(int n,float * restrict x,int stridex,float * restrict y,int stridey)144 __vlogf(int n, float * restrict x, int stridex, float * restrict y,
145 	int stridey)
146 {
147 	double	LN2 = __TBL_logf[64];		/* log(2) = 0.6931471805599453094 	*/
148 	double	yy0, yy1, yy2, yy3, yy4;
149 	double	ty0, ty1, ty2, ty3, ty4;
150 	float	value;
151 	int	i0, i1, i2, i3, i4;
152 	int	ival0, ival1, ival2, ival3, ival4;
153 	int	exp0, exp1, exp2, exp3, exp4;
154 	int	iy0, iy1, iy2, iy3, iy4;
155
156 	y -= stridey;
157
158 	for (; ;)
159 	{
160 begin:
161 		y += stridey;
162
163 		if (--n < 0)
164 			break;
165
166 		PREPROCESS(0, 0, begin)
167
168 		if (--n < 0)
169 			goto process1;
170
171 		PREPROCESS(1, stridey, process1)
172
173 		if (--n < 0)
174 			goto process2;
175
176 		PREPROCESS(2, (stridey << 1), process2)
177
178 		if (--n < 0)
179 			goto process3;
180
181 		PREPROCESS(3, (stridey << 1) + stridey, process3)
182
183 		if (--n < 0)
184 			goto process4;
185
186 		PREPROCESS(4, (stridey << 2), process4)
187
188 		iy0 = ival0 & 0x007fffff;
189 		iy1 = ival1 & 0x007fffff;
190 		iy2 = ival2 & 0x007fffff;
191 		iy3 = ival3 & 0x007fffff;
192 		iy4 = ival4 & 0x007fffff;
193
194 		ival0 = (iy0 + 0x20000) & 0xfffc0000;
195 		ival1 = (iy1 + 0x20000) & 0xfffc0000;
196 		ival2 = (iy2 + 0x20000) & 0xfffc0000;
197 		ival3 = (iy3 + 0x20000) & 0xfffc0000;
198 		ival4 = (iy4 + 0x20000) & 0xfffc0000;
199
200 		i0 = ival0 >> 17;
201 		i1 = ival1 >> 17;
202 		i2 = ival2 >> 17;
203 		i3 = ival3 >> 17;
204 		i4 = ival4 >> 17;
205
206 		iy0 = iy0 - ival0;
207 		iy1 = iy1 - ival1;
208 		iy2 = iy2 - ival2;
209 		iy3 = iy3 - ival3;
210 		iy4 = iy4 - ival4;
211
212 		ty0 = LN2 * (double) exp0 + __TBL_logf[i0];
213 		ty1 = LN2 * (double) exp1 + __TBL_logf[i1];
214 		ty2 = LN2 * (double) exp2 + __TBL_logf[i2];
215 		ty3 = LN2 * (double) exp3 + __TBL_logf[i3];
216 		ty4 = LN2 * (double) exp4 + __TBL_logf[i4];
217
218 		yy0 = (double) iy0 * __TBL_logf[i0 + 1];
219 		yy1 = (double) iy1 * __TBL_logf[i1 + 1];
220 		yy2 = (double) iy2 * __TBL_logf[i2 + 1];
221 		yy3 = (double) iy3 * __TBL_logf[i3 + 1];
222 		yy4 = (double) iy4 * __TBL_logf[i4 + 1];
223
224 		yy0 = ((K3 * yy0 + K2) * yy0 + K1) * yy0 * yy0 + yy0;
225 		yy1 = ((K3 * yy1 + K2) * yy1 + K1) * yy1 * yy1 + yy1;
226 		yy2 = ((K3 * yy2 + K2) * yy2 + K1) * yy2 * yy2 + yy2;
227 		yy3 = ((K3 * yy3 + K2) * yy3 + K1) * yy3 * yy3 + yy3;
228 		yy4 = ((K3 * yy4 + K2) * yy4 + K1) * yy4 * yy4 + yy4;
229
230 		y[0] = (float)(yy0 + ty0);
231 		y += stridey;
232 		y[0] = (float)(yy1 + ty1);
233 		y += stridey;
234 		y[0] = (float)(yy2 + ty2);
235 		y += stridey;
236 		y[0] = (float)(yy3 + ty3);
237 		y += stridey;
238 		y[0] = (float)(yy4 + ty4);
239 		continue;
240
241 process1:
242 		PROCESS(0)
243 		continue;
244
245 process2:
246 		PROCESS(0)
247 		PROCESS(1)
248 		continue;
249
250 process3:
251 		PROCESS(0)
252 		PROCESS(1)
253 		PROCESS(2)
254 		continue;
255
256 process4:
257 		PROCESS(0)
258 		PROCESS(1)
259 		PROCESS(2)
260 		PROCESS(3)
261 	}
262 }
263 ```