1 /*
2 * This file is part of the MicroPython project, http://micropython.org/
3 *
4 * The MIT License (MIT)
5 *
6 * Copyright (c) 2013, 2014 Damien P. George
7 *
8 * Permission is hereby granted, free of charge, to any person obtaining a copy
9 * of this software and associated documentation files (the "Software"), to deal
10 * in the Software without restriction, including without limitation the rights
11 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
12 * copies of the Software, and to permit persons to whom the Software is
13 * furnished to do so, subject to the following conditions:
14 *
15 * The above copyright notice and this permission notice shall be included in
16 * all copies or substantial portions of the Software.
17 *
18 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
19 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
20 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
21 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
22 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
23 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
24 * THE SOFTWARE.
25 */
26
27 #include "libm.h"
28
29 typedef float float_t;
30 typedef union {
31 float f;
32 struct {
33 uint32_t m : 23;
34 uint32_t e : 8;
35 uint32_t s : 1;
36 };
37 } float_s_t;
38
__signbitf(float f)39 int __signbitf(float f) {
40 float_s_t u = {.f = f};
41 return u.s;
42 }
43
44 #ifndef NDEBUG
copysignf(float x,float y)45 float copysignf(float x, float y) {
46 float_s_t fx={.f = x};
47 float_s_t fy={.f = y};
48
49 // copy sign bit;
50 fx.s = fy.s;
51
52 return fx.f;
53 }
54 #endif
55
56 static const float _M_LN10 = 2.30258509299404f; // 0x40135d8e
log10f(float x)57 float log10f(float x) { return logf(x) / (float)_M_LN10; }
58
tanhf(float x)59 float tanhf(float x) {
60 int sign = 0;
61 if (x < 0) {
62 sign = 1;
63 x = -x;
64 }
65 x = expm1f(-2 * x);
66 x = x / (x + 2);
67 return sign ? x : -x;
68 }
69
70 /*****************************************************************************/
71 /*****************************************************************************/
72 // __fpclassifyf from musl-0.9.15
73 /*****************************************************************************/
74 /*****************************************************************************/
75
__fpclassifyf(float x)76 int __fpclassifyf(float x)
77 {
78 union {float f; uint32_t i;} u = {x};
79 int e = u.i>>23 & 0xff;
80 if (!e) return u.i<<1 ? FP_SUBNORMAL : FP_ZERO;
81 if (e==0xff) return u.i<<9 ? FP_NAN : FP_INFINITE;
82 return FP_NORMAL;
83 }
84
85 /*****************************************************************************/
86 /*****************************************************************************/
87 // scalbnf from musl-0.9.15
88 /*****************************************************************************/
89 /*****************************************************************************/
90
scalbnf(float x,int n)91 float scalbnf(float x, int n)
92 {
93 union {float f; uint32_t i;} u;
94 float_t y = x;
95
96 if (n > 127) {
97 y *= 0x1p127f;
98 n -= 127;
99 if (n > 127) {
100 y *= 0x1p127f;
101 n -= 127;
102 if (n > 127)
103 n = 127;
104 }
105 } else if (n < -126) {
106 y *= 0x1p-126f;
107 n += 126;
108 if (n < -126) {
109 y *= 0x1p-126f;
110 n += 126;
111 if (n < -126)
112 n = -126;
113 }
114 }
115 u.i = (uint32_t)(0x7f+n)<<23;
116 x = y * u.f;
117 return x;
118 }
119
120 /*****************************************************************************/
121 /*****************************************************************************/
122 // powf from musl-0.9.15
123 /*****************************************************************************/
124 /*****************************************************************************/
125
126 /* origin: FreeBSD /usr/src/lib/msun/src/e_powf.c */
127 /*
128 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
129 */
130 /*
131 * ====================================================
132 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
133 *
134 * Developed at SunPro, a Sun Microsystems, Inc. business.
135 * Permission to use, copy, modify, and distribute this
136 * software is freely granted, provided that this notice
137 * is preserved.
138 * ====================================================
139 */
140
141 static const float
142 bp[] = {1.0f, 1.5f,},
143 dp_h[] = { 0.0f, 5.84960938e-01f,}, /* 0x3f15c000 */
144 dp_l[] = { 0.0f, 1.56322085e-06f,}, /* 0x35d1cfdc */
145 two24 = 16777216.0f, /* 0x4b800000 */
146 huge = 1.0e30f,
147 tiny = 1.0e-30f,
148 /* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */
149 L1 = 6.0000002384e-01f, /* 0x3f19999a */
150 L2 = 4.2857143283e-01f, /* 0x3edb6db7 */
151 L3 = 3.3333334327e-01f, /* 0x3eaaaaab */
152 L4 = 2.7272811532e-01f, /* 0x3e8ba305 */
153 L5 = 2.3066075146e-01f, /* 0x3e6c3255 */
154 L6 = 2.0697501302e-01f, /* 0x3e53f142 */
155 P1 = 1.6666667163e-01f, /* 0x3e2aaaab */
156 P2 = -2.7777778450e-03f, /* 0xbb360b61 */
157 P3 = 6.6137559770e-05f, /* 0x388ab355 */
158 P4 = -1.6533901999e-06f, /* 0xb5ddea0e */
159 P5 = 4.1381369442e-08f, /* 0x3331bb4c */
160 lg2 = 6.9314718246e-01f, /* 0x3f317218 */
161 lg2_h = 6.93145752e-01f, /* 0x3f317200 */
162 lg2_l = 1.42860654e-06f, /* 0x35bfbe8c */
163 ovt = 4.2995665694e-08f, /* -(128-log2(ovfl+.5ulp)) */
164 cp = 9.6179670095e-01f, /* 0x3f76384f =2/(3ln2) */
165 cp_h = 9.6191406250e-01f, /* 0x3f764000 =12b cp */
166 cp_l = -1.1736857402e-04f, /* 0xb8f623c6 =tail of cp_h */
167 ivln2 = 1.4426950216e+00f, /* 0x3fb8aa3b =1/ln2 */
168 ivln2_h = 1.4426879883e+00f, /* 0x3fb8aa00 =16b 1/ln2*/
169 ivln2_l = 7.0526075433e-06f; /* 0x36eca570 =1/ln2 tail*/
170
powf(float x,float y)171 float powf(float x, float y)
172 {
173 float z,ax,z_h,z_l,p_h,p_l;
174 float y1,t1,t2,r,s,sn,t,u,v,w;
175 int32_t i,j,k,yisint,n;
176 int32_t hx,hy,ix,iy,is;
177
178 GET_FLOAT_WORD(hx, x);
179 GET_FLOAT_WORD(hy, y);
180 ix = hx & 0x7fffffff;
181 iy = hy & 0x7fffffff;
182
183 /* x**0 = 1, even if x is NaN */
184 if (iy == 0)
185 return 1.0f;
186 /* 1**y = 1, even if y is NaN */
187 if (hx == 0x3f800000)
188 return 1.0f;
189 /* NaN if either arg is NaN */
190 if (ix > 0x7f800000 || iy > 0x7f800000)
191 return x + y;
192
193 /* determine if y is an odd int when x < 0
194 * yisint = 0 ... y is not an integer
195 * yisint = 1 ... y is an odd int
196 * yisint = 2 ... y is an even int
197 */
198 yisint = 0;
199 if (hx < 0) {
200 if (iy >= 0x4b800000)
201 yisint = 2; /* even integer y */
202 else if (iy >= 0x3f800000) {
203 k = (iy>>23) - 0x7f; /* exponent */
204 j = iy>>(23-k);
205 if ((j<<(23-k)) == iy)
206 yisint = 2 - (j & 1);
207 }
208 }
209
210 /* special value of y */
211 if (iy == 0x7f800000) { /* y is +-inf */
212 if (ix == 0x3f800000) /* (-1)**+-inf is 1 */
213 return 1.0f;
214 else if (ix > 0x3f800000) /* (|x|>1)**+-inf = inf,0 */
215 return hy >= 0 ? y : 0.0f;
216 else if (ix != 0) /* (|x|<1)**+-inf = 0,inf if x!=0 */
217 return hy >= 0 ? 0.0f: -y;
218 }
219 if (iy == 0x3f800000) /* y is +-1 */
220 return hy >= 0 ? x : 1.0f/x;
221 if (hy == 0x40000000) /* y is 2 */
222 return x*x;
223 if (hy == 0x3f000000) { /* y is 0.5 */
224 if (hx >= 0) /* x >= +0 */
225 return sqrtf(x);
226 }
227
228 ax = fabsf(x);
229 /* special value of x */
230 if (ix == 0x7f800000 || ix == 0 || ix == 0x3f800000) { /* x is +-0,+-inf,+-1 */
231 z = ax;
232 if (hy < 0) /* z = (1/|x|) */
233 z = 1.0f/z;
234 if (hx < 0) {
235 if (((ix-0x3f800000)|yisint) == 0) {
236 z = (z-z)/(z-z); /* (-1)**non-int is NaN */
237 } else if (yisint == 1)
238 z = -z; /* (x<0)**odd = -(|x|**odd) */
239 }
240 return z;
241 }
242
243 sn = 1.0f; /* sign of result */
244 if (hx < 0) {
245 if (yisint == 0) /* (x<0)**(non-int) is NaN */
246 return (x-x)/(x-x);
247 if (yisint == 1) /* (x<0)**(odd int) */
248 sn = -1.0f;
249 }
250
251 /* |y| is huge */
252 if (iy > 0x4d000000) { /* if |y| > 2**27 */
253 /* over/underflow if x is not close to one */
254 if (ix < 0x3f7ffff8)
255 return hy < 0 ? sn*huge*huge : sn*tiny*tiny;
256 if (ix > 0x3f800007)
257 return hy > 0 ? sn*huge*huge : sn*tiny*tiny;
258 /* now |1-x| is tiny <= 2**-20, suffice to compute
259 log(x) by x-x^2/2+x^3/3-x^4/4 */
260 t = ax - 1; /* t has 20 trailing zeros */
261 w = (t*t)*(0.5f - t*(0.333333333333f - t*0.25f));
262 u = ivln2_h*t; /* ivln2_h has 16 sig. bits */
263 v = t*ivln2_l - w*ivln2;
264 t1 = u + v;
265 GET_FLOAT_WORD(is, t1);
266 SET_FLOAT_WORD(t1, is & 0xfffff000);
267 t2 = v - (t1-u);
268 } else {
269 float s2,s_h,s_l,t_h,t_l;
270 n = 0;
271 /* take care subnormal number */
272 if (ix < 0x00800000) {
273 ax *= two24;
274 n -= 24;
275 GET_FLOAT_WORD(ix, ax);
276 }
277 n += ((ix)>>23) - 0x7f;
278 j = ix & 0x007fffff;
279 /* determine interval */
280 ix = j | 0x3f800000; /* normalize ix */
281 if (j <= 0x1cc471) /* |x|<sqrt(3/2) */
282 k = 0;
283 else if (j < 0x5db3d7) /* |x|<sqrt(3) */
284 k = 1;
285 else {
286 k = 0;
287 n += 1;
288 ix -= 0x00800000;
289 }
290 SET_FLOAT_WORD(ax, ix);
291
292 /* compute s = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */
293 u = ax - bp[k]; /* bp[0]=1.0, bp[1]=1.5 */
294 v = 1.0f/(ax+bp[k]);
295 s = u*v;
296 s_h = s;
297 GET_FLOAT_WORD(is, s_h);
298 SET_FLOAT_WORD(s_h, is & 0xfffff000);
299 /* t_h=ax+bp[k] High */
300 is = ((ix>>1) & 0xfffff000) | 0x20000000;
301 SET_FLOAT_WORD(t_h, is + 0x00400000 + (k<<21));
302 t_l = ax - (t_h - bp[k]);
303 s_l = v*((u - s_h*t_h) - s_h*t_l);
304 /* compute log(ax) */
305 s2 = s*s;
306 r = s2*s2*(L1+s2*(L2+s2*(L3+s2*(L4+s2*(L5+s2*L6)))));
307 r += s_l*(s_h+s);
308 s2 = s_h*s_h;
309 t_h = 3.0f + s2 + r;
310 GET_FLOAT_WORD(is, t_h);
311 SET_FLOAT_WORD(t_h, is & 0xfffff000);
312 t_l = r - ((t_h - 3.0f) - s2);
313 /* u+v = s*(1+...) */
314 u = s_h*t_h;
315 v = s_l*t_h + t_l*s;
316 /* 2/(3log2)*(s+...) */
317 p_h = u + v;
318 GET_FLOAT_WORD(is, p_h);
319 SET_FLOAT_WORD(p_h, is & 0xfffff000);
320 p_l = v - (p_h - u);
321 z_h = cp_h*p_h; /* cp_h+cp_l = 2/(3*log2) */
322 z_l = cp_l*p_h + p_l*cp+dp_l[k];
323 /* log2(ax) = (s+..)*2/(3*log2) = n + dp_h + z_h + z_l */
324 t = (float)n;
325 t1 = (((z_h + z_l) + dp_h[k]) + t);
326 GET_FLOAT_WORD(is, t1);
327 SET_FLOAT_WORD(t1, is & 0xfffff000);
328 t2 = z_l - (((t1 - t) - dp_h[k]) - z_h);
329 }
330
331 /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */
332 GET_FLOAT_WORD(is, y);
333 SET_FLOAT_WORD(y1, is & 0xfffff000);
334 p_l = (y-y1)*t1 + y*t2;
335 p_h = y1*t1;
336 z = p_l + p_h;
337 GET_FLOAT_WORD(j, z);
338 if (j > 0x43000000) /* if z > 128 */
339 return sn*huge*huge; /* overflow */
340 else if (j == 0x43000000) { /* if z == 128 */
341 if (p_l + ovt > z - p_h)
342 return sn*huge*huge; /* overflow */
343 } else if ((j&0x7fffffff) > 0x43160000) /* z < -150 */ // FIXME: check should be (uint32_t)j > 0xc3160000
344 return sn*tiny*tiny; /* underflow */
345 else if (j == 0xc3160000) { /* z == -150 */
346 if (p_l <= z-p_h)
347 return sn*tiny*tiny; /* underflow */
348 }
349 /*
350 * compute 2**(p_h+p_l)
351 */
352 i = j & 0x7fffffff;
353 k = (i>>23) - 0x7f;
354 n = 0;
355 if (i > 0x3f000000) { /* if |z| > 0.5, set n = [z+0.5] */
356 n = j + (0x00800000>>(k+1));
357 k = ((n&0x7fffffff)>>23) - 0x7f; /* new k for n */
358 SET_FLOAT_WORD(t, n & ~(0x007fffff>>k));
359 n = ((n&0x007fffff)|0x00800000)>>(23-k);
360 if (j < 0)
361 n = -n;
362 p_h -= t;
363 }
364 t = p_l + p_h;
365 GET_FLOAT_WORD(is, t);
366 SET_FLOAT_WORD(t, is & 0xffff8000);
367 u = t*lg2_h;
368 v = (p_l-(t-p_h))*lg2 + t*lg2_l;
369 z = u + v;
370 w = v - (z - u);
371 t = z*z;
372 t1 = z - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
373 r = (z*t1)/(t1-2.0f) - (w+z*w);
374 z = 1.0f - (r - z);
375 GET_FLOAT_WORD(j, z);
376 j += n<<23;
377 if ((j>>23) <= 0) /* subnormal output */
378 z = scalbnf(z, n);
379 else
380 SET_FLOAT_WORD(z, j);
381 return sn*z;
382 }
383
384 /*****************************************************************************/
385 /*****************************************************************************/
386 // expf from musl-0.9.15
387 /*****************************************************************************/
388 /*****************************************************************************/
389
390 /* origin: FreeBSD /usr/src/lib/msun/src/e_expf.c */
391 /*
392 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
393 */
394 /*
395 * ====================================================
396 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
397 *
398 * Developed at SunPro, a Sun Microsystems, Inc. business.
399 * Permission to use, copy, modify, and distribute this
400 * software is freely granted, provided that this notice
401 * is preserved.
402 * ====================================================
403 */
404
405 static const float
406 half[2] = {0.5f,-0.5f},
407 ln2hi = 6.9314575195e-1f, /* 0x3f317200 */
408 ln2lo = 1.4286067653e-6f, /* 0x35bfbe8e */
409 invln2 = 1.4426950216e+0f, /* 0x3fb8aa3b */
410 /*
411 * Domain [-0.34568, 0.34568], range ~[-4.278e-9, 4.447e-9]:
412 * |x*(exp(x)+1)/(exp(x)-1) - p(x)| < 2**-27.74
413 */
414 expf_P1 = 1.6666625440e-1f, /* 0xaaaa8f.0p-26 */
415 expf_P2 = -2.7667332906e-3f; /* -0xb55215.0p-32 */
416
expf(float x)417 float expf(float x)
418 {
419 float_t hi, lo, c, xx, y;
420 int k, sign;
421 uint32_t hx;
422
423 GET_FLOAT_WORD(hx, x);
424 sign = hx >> 31; /* sign bit of x */
425 hx &= 0x7fffffff; /* high word of |x| */
426
427 /* special cases */
428 if (hx >= 0x42aeac50) { /* if |x| >= -87.33655f or NaN */
429 if (hx >= 0x42b17218 && !sign) { /* x >= 88.722839f */
430 /* overflow */
431 x *= 0x1p127f;
432 return x;
433 }
434 if (sign) {
435 /* underflow */
436 FORCE_EVAL(-0x1p-149f/x);
437 if (hx >= 0x42cff1b5) /* x <= -103.972084f */
438 return 0;
439 }
440 }
441
442 /* argument reduction */
443 if (hx > 0x3eb17218) { /* if |x| > 0.5 ln2 */
444 if (hx > 0x3f851592) /* if |x| > 1.5 ln2 */
445 k = (int)(invln2*x + half[sign]);
446 else
447 k = 1 - sign - sign;
448 hi = x - k*ln2hi; /* k*ln2hi is exact here */
449 lo = k*ln2lo;
450 x = hi - lo;
451 } else if (hx > 0x39000000) { /* |x| > 2**-14 */
452 k = 0;
453 hi = x;
454 lo = 0;
455 } else {
456 /* raise inexact */
457 FORCE_EVAL(0x1p127f + x);
458 return 1 + x;
459 }
460
461 /* x is now in primary range */
462 xx = x*x;
463 c = x - xx*(expf_P1+xx*expf_P2);
464 y = 1 + (x*c/(2-c) - lo + hi);
465 if (k == 0)
466 return y;
467 return scalbnf(y, k);
468 }
469
470 /*****************************************************************************/
471 /*****************************************************************************/
472 // expm1f from musl-0.9.15
473 /*****************************************************************************/
474 /*****************************************************************************/
475
476 /* origin: FreeBSD /usr/src/lib/msun/src/s_expm1f.c */
477 /*
478 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
479 */
480 /*
481 * ====================================================
482 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
483 *
484 * Developed at SunPro, a Sun Microsystems, Inc. business.
485 * Permission to use, copy, modify, and distribute this
486 * software is freely granted, provided that this notice
487 * is preserved.
488 * ====================================================
489 */
490
491 static const float
492 o_threshold = 8.8721679688e+01f, /* 0x42b17180 */
493 ln2_hi = 6.9313812256e-01f, /* 0x3f317180 */
494 ln2_lo = 9.0580006145e-06f, /* 0x3717f7d1 */
495 //invln2 = 1.4426950216e+00, /* 0x3fb8aa3b */
496 /*
497 * Domain [-0.34568, 0.34568], range ~[-6.694e-10, 6.696e-10]:
498 * |6 / x * (1 + 2 * (1 / (exp(x) - 1) - 1 / x)) - q(x)| < 2**-30.04
499 * Scaled coefficients: Qn_here = 2**n * Qn_for_q (see s_expm1.c):
500 */
501 Q1 = -3.3333212137e-2f, /* -0x888868.0p-28 */
502 Q2 = 1.5807170421e-3f; /* 0xcf3010.0p-33 */
503
expm1f(float x)504 float expm1f(float x)
505 {
506 float_t y,hi,lo,c,t,e,hxs,hfx,r1,twopk;
507 union {float f; uint32_t i;} u = {x};
508 uint32_t hx = u.i & 0x7fffffff;
509 int k, sign = u.i >> 31;
510
511 /* filter out huge and non-finite argument */
512 if (hx >= 0x4195b844) { /* if |x|>=27*ln2 */
513 if (hx > 0x7f800000) /* NaN */
514 return x;
515 if (sign)
516 return -1;
517 if (x > o_threshold) {
518 x *= 0x1p127f;
519 return x;
520 }
521 }
522
523 /* argument reduction */
524 if (hx > 0x3eb17218) { /* if |x| > 0.5 ln2 */
525 if (hx < 0x3F851592) { /* and |x| < 1.5 ln2 */
526 if (!sign) {
527 hi = x - ln2_hi;
528 lo = ln2_lo;
529 k = 1;
530 } else {
531 hi = x + ln2_hi;
532 lo = -ln2_lo;
533 k = -1;
534 }
535 } else {
536 k = (int)(invln2*x + (sign ? -0.5f : 0.5f));
537 t = k;
538 hi = x - t*ln2_hi; /* t*ln2_hi is exact here */
539 lo = t*ln2_lo;
540 }
541 x = hi-lo;
542 c = (hi-x)-lo;
543 } else if (hx < 0x33000000) { /* when |x|<2**-25, return x */
544 if (hx < 0x00800000)
545 FORCE_EVAL(x*x);
546 return x;
547 } else
548 k = 0;
549
550 /* x is now in primary range */
551 hfx = 0.5f*x;
552 hxs = x*hfx;
553 r1 = 1.0f+hxs*(Q1+hxs*Q2);
554 t = 3.0f - r1*hfx;
555 e = hxs*((r1-t)/(6.0f - x*t));
556 if (k == 0) /* c is 0 */
557 return x - (x*e-hxs);
558 e = x*(e-c) - c;
559 e -= hxs;
560 /* exp(x) ~ 2^k (x_reduced - e + 1) */
561 if (k == -1)
562 return 0.5f*(x-e) - 0.5f;
563 if (k == 1) {
564 if (x < -0.25f)
565 return -2.0f*(e-(x+0.5f));
566 return 1.0f + 2.0f*(x-e);
567 }
568 u.i = (0x7f+k)<<23; /* 2^k */
569 twopk = u.f;
570 if (k < 0 || k > 56) { /* suffice to return exp(x)-1 */
571 y = x - e + 1.0f;
572 if (k == 128)
573 y = y*2.0f*0x1p127f;
574 else
575 y = y*twopk;
576 return y - 1.0f;
577 }
578 u.i = (0x7f-k)<<23; /* 2^-k */
579 if (k < 23)
580 y = (x-e+(1-u.f))*twopk;
581 else
582 y = (x-(e+u.f)+1)*twopk;
583 return y;
584 }
585
586 /*****************************************************************************/
587 /*****************************************************************************/
588 // __expo2f from musl-0.9.15
589 /*****************************************************************************/
590 /*****************************************************************************/
591
592 /* k is such that k*ln2 has minimal relative error and x - kln2 > log(FLT_MIN) */
593 static const int k = 235;
594 static const float kln2 = 0x1.45c778p+7f;
595
596 /* expf(x)/2 for x >= log(FLT_MAX), slightly better than 0.5f*expf(x/2)*expf(x/2) */
__expo2f(float x)597 float __expo2f(float x)
598 {
599 float scale;
600
601 /* note that k is odd and scale*scale overflows */
602 SET_FLOAT_WORD(scale, (uint32_t)(0x7f + k/2) << 23);
603 /* exp(x - k ln2) * 2**(k-1) */
604 return expf(x - kln2) * scale * scale;
605 }
606
607 /*****************************************************************************/
608 /*****************************************************************************/
609 // logf from musl-0.9.15
610 /*****************************************************************************/
611 /*****************************************************************************/
612
613 /* origin: FreeBSD /usr/src/lib/msun/src/e_logf.c */
614 /*
615 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
616 */
617 /*
618 * ====================================================
619 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
620 *
621 * Developed at SunPro, a Sun Microsystems, Inc. business.
622 * Permission to use, copy, modify, and distribute this
623 * software is freely granted, provided that this notice
624 * is preserved.
625 * ====================================================
626 */
627
628 static const float
629 /* |(log(1+s)-log(1-s))/s - Lg(s)| < 2**-34.24 (~[-4.95e-11, 4.97e-11]). */
630 Lg1 = 0xaaaaaa.0p-24, /* 0.66666662693 */
631 Lg2 = 0xccce13.0p-25, /* 0.40000972152 */
632 Lg3 = 0x91e9ee.0p-25, /* 0.28498786688 */
633 Lg4 = 0xf89e26.0p-26; /* 0.24279078841 */
634
logf(float x)635 float logf(float x)
636 {
637 union {float f; uint32_t i;} u = {x};
638 float_t hfsq,f,s,z,R,w,t1,t2,dk;
639 uint32_t ix;
640 int k;
641
642 ix = u.i;
643 k = 0;
644 if (ix < 0x00800000 || ix>>31) { /* x < 2**-126 */
645 if (ix<<1 == 0)
646 return -1/(x*x); /* log(+-0)=-inf */
647 if (ix>>31)
648 return (x-x)/0.0f; /* log(-#) = NaN */
649 /* subnormal number, scale up x */
650 k -= 25;
651 x *= 0x1p25f;
652 u.f = x;
653 ix = u.i;
654 } else if (ix >= 0x7f800000) {
655 return x;
656 } else if (ix == 0x3f800000)
657 return 0;
658
659 /* reduce x into [sqrt(2)/2, sqrt(2)] */
660 ix += 0x3f800000 - 0x3f3504f3;
661 k += (int)(ix>>23) - 0x7f;
662 ix = (ix&0x007fffff) + 0x3f3504f3;
663 u.i = ix;
664 x = u.f;
665
666 f = x - 1.0f;
667 s = f/(2.0f + f);
668 z = s*s;
669 w = z*z;
670 t1= w*(Lg2+w*Lg4);
671 t2= z*(Lg1+w*Lg3);
672 R = t2 + t1;
673 hfsq = 0.5f*f*f;
674 dk = k;
675 return s*(hfsq+R) + dk*ln2_lo - hfsq + f + dk*ln2_hi;
676 }
677
678 /*****************************************************************************/
679 /*****************************************************************************/
680 // coshf from musl-0.9.15
681 /*****************************************************************************/
682 /*****************************************************************************/
683
coshf(float x)684 float coshf(float x)
685 {
686 union {float f; uint32_t i;} u = {.f = x};
687 uint32_t w;
688 float t;
689
690 /* |x| */
691 u.i &= 0x7fffffff;
692 x = u.f;
693 w = u.i;
694
695 /* |x| < log(2) */
696 if (w < 0x3f317217) {
697 if (w < 0x3f800000 - (12<<23)) {
698 FORCE_EVAL(x + 0x1p120f);
699 return 1;
700 }
701 t = expm1f(x);
702 return 1 + t*t/(2*(1+t));
703 }
704
705 /* |x| < log(FLT_MAX) */
706 if (w < 0x42b17217) {
707 t = expf(x);
708 return 0.5f*(t + 1/t);
709 }
710
711 /* |x| > log(FLT_MAX) or nan */
712 t = __expo2f(x);
713 return t;
714 }
715
716 /*****************************************************************************/
717 /*****************************************************************************/
718 // sinhf from musl-0.9.15
719 /*****************************************************************************/
720 /*****************************************************************************/
721
sinhf(float x)722 float sinhf(float x)
723 {
724 union {float f; uint32_t i;} u = {.f = x};
725 uint32_t w;
726 float t, h, absx;
727
728 h = 0.5;
729 if (u.i >> 31)
730 h = -h;
731 /* |x| */
732 u.i &= 0x7fffffff;
733 absx = u.f;
734 w = u.i;
735
736 /* |x| < log(FLT_MAX) */
737 if (w < 0x42b17217) {
738 t = expm1f(absx);
739 if (w < 0x3f800000) {
740 if (w < 0x3f800000 - (12<<23))
741 return x;
742 return h*(2*t - t*t/(t+1));
743 }
744 return h*(t + t/(t+1));
745 }
746
747 /* |x| > logf(FLT_MAX) or nan */
748 t = 2*h*__expo2f(absx);
749 return t;
750 }
751
752 /*****************************************************************************/
753 /*****************************************************************************/
754 // ceilf, floorf and truncf from musl-0.9.15
755 /*****************************************************************************/
756 /*****************************************************************************/
757
ceilf(float x)758 float ceilf(float x)
759 {
760 union {float f; uint32_t i;} u = {x};
761 int e = (int)(u.i >> 23 & 0xff) - 0x7f;
762 uint32_t m;
763
764 if (e >= 23)
765 return x;
766 if (e >= 0) {
767 m = 0x007fffff >> e;
768 if ((u.i & m) == 0)
769 return x;
770 FORCE_EVAL(x + 0x1p120f);
771 if (u.i >> 31 == 0)
772 u.i += m;
773 u.i &= ~m;
774 } else {
775 FORCE_EVAL(x + 0x1p120f);
776 if (u.i >> 31)
777 u.f = -0.0;
778 else if (u.i << 1)
779 u.f = 1.0;
780 }
781 return u.f;
782 }
783
floorf(float x)784 float floorf(float x)
785 {
786 union {float f; uint32_t i;} u = {x};
787 int e = (int)(u.i >> 23 & 0xff) - 0x7f;
788 uint32_t m;
789
790 if (e >= 23)
791 return x;
792 if (e >= 0) {
793 m = 0x007fffff >> e;
794 if ((u.i & m) == 0)
795 return x;
796 FORCE_EVAL(x + 0x1p120f);
797 if (u.i >> 31)
798 u.i += m;
799 u.i &= ~m;
800 } else {
801 FORCE_EVAL(x + 0x1p120f);
802 if (u.i >> 31 == 0)
803 u.i = 0;
804 else if (u.i << 1)
805 u.f = -1.0;
806 }
807 return u.f;
808 }
809
truncf(float x)810 float truncf(float x)
811 {
812 union {float f; uint32_t i;} u = {x};
813 int e = (int)(u.i >> 23 & 0xff) - 0x7f + 9;
814 uint32_t m;
815
816 if (e >= 23 + 9)
817 return x;
818 if (e < 9)
819 e = 1;
820 m = -1U >> e;
821 if ((u.i & m) == 0)
822 return x;
823 FORCE_EVAL(x + 0x1p120f);
824 u.i &= ~m;
825 return u.f;
826 }
827