Fix red in bootloaders
[maemo-rb.git] / apps / plugins / pdbox / pdbox-func.c
blobee4a8fd1661db5e2404bf93be38de22aaac3683f
1 /***************************************************************************
2 * __________ __ ___.
3 * Open \______ \ ____ ____ | | _\_ |__ _______ ___
4 * Source | _// _ \_/ ___\| |/ /| __ \ / _ \ \/ /
5 * Jukebox | | ( <_> ) \___| < | \_\ ( <_> > < <
6 * Firmware |____|_ /\____/ \___ >__|_ \|___ /\____/__/\_ \
7 * \/ \/ \/ \/ \/
8 * $Id$
10 * Copyright (C) 2009 Wincent Balin
12 * This program is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU General Public License
14 * as published by the Free Software Foundation; either version 2
15 * of the License, or (at your option) any later version.
17 * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY OF ANY
18 * KIND, either express or implied.
20 ****************************************************************************/
22 #include "plugin.h"
23 #include "pdbox.h"
24 #include "ctype.h"
26 #include "m_pd.h"
27 #include "s_stuff.h"
30 /* This implementation of strncat is taken from lua plug-in. */
32 /* gcc is broken and has a non-SUSv2 compliant internal prototype.
33 * This causes it to warn about a type mismatch here. Ignore it. */
34 char *rb_strncat(char *s, const char *t, size_t n)
36 char *dest = s;
37 register char *max;
38 s += strlen(s);
40 if((max = s + n) == s)
41 goto strncat_fini;
43 while(true)
45 if(!(*s = *t))
46 break;
47 if(++s == max)
48 break;
49 ++t;
51 #ifndef WANT_SMALL_STRING_ROUTINES
52 if(!(*s = *t))
53 break;
54 if(++s == max)
55 break;
56 ++t;
57 if(!(*s = *t))
58 break;
59 if(++s == max)
60 break;
61 ++t;
62 if(!(*s = *t))
63 break;
64 if(++s == max)
65 break;
66 ++t;
67 #endif
70 *s = 0;
72 strncat_fini:
73 return dest;
77 /* Implementation of floor, original. */
78 float rb_floor(float value)
80 /* If value is negative, decrement value to match function's definition. */
81 if(value < 0.0)
83 value -= 1.0;
86 /* Truncate fractional part (convert to integer)
87 and afterwards convert back to double. */
88 return (float) ((int) value);
92 /* Implementation of strtod() and atof(),
93 taken from SanOS (http://www.jbox.dk/sanos/). */
94 static int rb_errno = 0;
96 double rb_strtod(const char *str, char **endptr)
98 double number;
99 int exponent;
100 int negative;
101 char *p = (char *) str;
102 double p10;
103 int n;
104 int num_digits;
105 int num_decimals;
107 /* Reset Rockbox errno -- W.B. */
108 #ifdef ROCKBOX
109 rb_errno = 0;
110 #endif
112 // Skip leading whitespace
113 while (isspace(*p)) p++;
115 // Handle optional sign
116 negative = 0;
117 switch (*p)
119 case '-': negative = 1; // Fall through to increment position
120 case '+': p++;
123 number = 0.;
124 exponent = 0;
125 num_digits = 0;
126 num_decimals = 0;
128 // Process string of digits
129 while (isdigit(*p))
131 number = number * 10. + (*p - '0');
132 p++;
133 num_digits++;
136 // Process decimal part
137 if (*p == '.')
139 p++;
141 while (isdigit(*p))
143 number = number * 10. + (*p - '0');
144 p++;
145 num_digits++;
146 num_decimals++;
149 exponent -= num_decimals;
152 if (num_digits == 0)
154 #ifdef ROCKBOX
155 rb_errno = 1;
156 #else
157 errno = ERANGE;
158 #endif
159 return 0.0;
162 // Correct for sign
163 if (negative) number = -number;
165 // Process an exponent string
166 if (*p == 'e' || *p == 'E')
168 // Handle optional sign
169 negative = 0;
170 switch(*++p)
172 case '-': negative = 1; // Fall through to increment pos
173 case '+': p++;
176 // Process string of digits
177 n = 0;
178 while (isdigit(*p))
180 n = n * 10 + (*p - '0');
181 p++;
184 if (negative)
185 exponent -= n;
186 else
187 exponent += n;
190 #ifndef ROCKBOX
191 if (exponent < DBL_MIN_EXP || exponent > DBL_MAX_EXP)
193 errno = ERANGE;
194 return HUGE_VAL;
196 #endif
198 // Scale the result
199 p10 = 10.;
200 n = exponent;
201 if (n < 0) n = -n;
202 while (n)
204 if (n & 1)
206 if (exponent < 0)
207 number /= p10;
208 else
209 number *= p10;
211 n >>= 1;
212 p10 *= p10;
215 #ifndef ROCKBOX
216 if (number == HUGE_VAL) errno = ERANGE;
217 #endif
218 if (endptr) *endptr = p;
220 return number;
223 double rb_atof(const char *str)
225 return rb_strtod(str, NULL);
229 /* Implementation of ftoa(), original. */
230 void rb_ftoan(float f, char* out, int size)
232 #define SBUFSIZE 12
233 char sbuf[SBUFSIZE];
235 /* Zero out string. */
236 *out = '\0';
237 size--;
239 /* Handle negative numbers. */
240 if(f < 0.0)
242 f = -f;
243 strcat(out, "-");
244 size--;
247 /* Find and convert integer part. */
248 int int_part = (int) f;
249 snprintf(sbuf, SBUFSIZE-1, "%d", int_part);
250 int int_part_len = strlen(sbuf);
251 if(size < int_part_len)
252 return;
254 /* Append integral part to output string. */
255 strcat(out, sbuf);
256 size -= int_part_len;
258 /* Check whether further content is possible. */
259 if(size <= 0)
260 return;
262 /* Append decimal point. */
263 strcat(out, ".");
264 size--;
266 /* Calculate first rest and convert it. */
267 float rest1 = (f - (float) int_part) * 1000000000.0;
268 int irest1 = (int) rest1;
269 snprintf(sbuf, SBUFSIZE-1, "%09d", irest1);
271 /* Append first rest to output string. */
272 int rest1_len = strlen(sbuf);
273 int rest1_minlen = MIN(size, rest1_len);
274 strncat(out, sbuf, rest1_minlen);
275 size -= rest1_minlen;
277 /* Check whether output string still has enough space. */
278 if(size <= 0)
279 return;
281 /* Calculate second rest and convert it. */
282 float rest2 = (rest1 - (float) irest1) * 1000000000.0;
283 int irest2 = (int) rest2;
284 snprintf(sbuf, SBUFSIZE-1, "%09d", irest2);
286 /* Append second rest to the output string. */
287 int rest2_len = strlen(sbuf);
288 int rest2_minlen = MIN(size, rest2_len);
289 strncat(out, sbuf, rest2_minlen);
293 /* Implementation of atol(), adapted from
294 the atoi() implementation in Rockbox. */
295 long rb_atol(const char* str)
297 long value = 0L;
298 long sign = 1L;
300 while (isspace(*str))
302 str++;
305 if ('-' == *str)
307 sign = -1L;
308 str++;
310 else if ('+' == *str)
312 str++;
315 while ('0' == *str)
317 str++;
320 while (isdigit(*str))
322 value = (value * 10L) + (*str - '0');
323 str++;
326 return value * sign;
330 /* Implementation of sin() and cos(),
331 adapted from http://lab.polygonal.de/2007/07/18/fast-and-accurate-sinecosine-approximation/
334 float rb_sin(float rad)
336 /* Trim input value to -PI..PI interval. */
337 if(rad < -3.14159265)
338 rad += 6.28318531;
339 else if(rad > 3.14159265)
340 rad -= 6.28318531;
342 if(rad < 0)
343 return (1.27323954 * rad + 0.405284735 * rad * rad);
344 else
345 return (1.27323954 * rad - 0.405284735 * rad * rad);
348 float rb_cos(float rad)
350 /* Compute cosine: sin(x + PI/2) = cos(x) */
351 rad += 1.57079632;
352 if(rad > 3.14159265)
353 rad -= 6.28318531;
355 return rb_sin(rad);
359 /* Emulation of fscanf(fd, "%f", (float*) xxx);
360 Basically a reimplementation of rb_strtod() above. */
361 int rb_fscanf_f(int fd, float* f)
363 #define FSCANF_F_BUFSIZE 64
364 char buf[FSCANF_F_BUFSIZE];
366 /* Read line from file. */
367 int bytes_read = rb->read_line(fd, buf, FSCANF_F_BUFSIZE-1);
369 /* Terminate string. */
370 if(bytes_read >= FSCANF_F_BUFSIZE)
371 buf[FSCANF_F_BUFSIZE-1] = '\0';
372 else
373 buf[bytes_read-1] = '\0';
375 /* Convert buffer to float. */
376 *f = rb_atof(buf);
378 /* If there was an error, no float was read. */
379 if(rb_errno)
380 return 0;
382 return 1;
385 /* Emulation of fprintf(fd, "%f\n", (float*) xxx); */
386 int rb_fprintf_f(int fd, float f)
388 #define FPRINTF_F_BUFSIZE 64
389 char buf[FPRINTF_F_BUFSIZE];
390 const char* next_line = "\n";
392 /* Convert float to string. */
393 rb_ftoan(f, buf, sizeof(buf)-1);
395 /* Add next line character. */
396 strcat(buf, next_line);
398 /* Write string into file. */
399 return write(fd, buf, strlen(buf));
403 /* Natural logarithm.
404 Taken from glibc-2.8 */
405 static const float
406 ln2_hi = 6.9313812256e-01, /* 0x3f317180 */
407 ln2_lo = 9.0580006145e-06, /* 0x3717f7d1 */
408 two25 = 3.355443200e+07, /* 0x4c000000 */
409 Lg1 = 6.6666668653e-01, /* 3F2AAAAB */
410 Lg2 = 4.0000000596e-01, /* 3ECCCCCD */
411 Lg3 = 2.8571429849e-01, /* 3E924925 */
412 Lg4 = 2.2222198546e-01, /* 3E638E29 */
413 Lg5 = 1.8183572590e-01, /* 3E3A3325 */
414 Lg6 = 1.5313838422e-01, /* 3E1CD04F */
415 Lg7 = 1.4798198640e-01; /* 3E178897 */
417 static const float zero = 0.0;
419 /* A union which permits us to convert between a float and a 32 bit
420 int. */
422 typedef union
424 float value;
425 uint32_t word;
426 } ieee_float_shape_type;
428 /* Get a 32 bit int from a float. */
430 #define GET_FLOAT_WORD(i,d) \
431 do { \
432 ieee_float_shape_type gf_u; \
433 gf_u.value = (d); \
434 (i) = gf_u.word; \
435 } while (0)
437 /* Set a float from a 32 bit int. */
439 #define SET_FLOAT_WORD(d,i) \
440 do { \
441 ieee_float_shape_type sf_u; \
442 sf_u.word = (i); \
443 (d) = sf_u.value; \
444 } while (0)
447 float rb_log(float x)
449 float hfsq, f, s, z, R, w, t1, t2, dk;
450 int32_t k, ix, i, j;
452 GET_FLOAT_WORD(ix,x);
454 k=0;
455 if (ix < 0x00800000) { /* x < 2**-126 */
456 if ((ix&0x7fffffff)==0)
457 return -two25/(x-x); /* log(+-0)=-inf */
458 if (ix<0) return (x-x)/(x-x); /* log(-#) = NaN */
459 k -= 25; x *= two25; /* subnormal number, scale up x */
460 GET_FLOAT_WORD(ix,x);
462 if (ix >= 0x7f800000) return x+x;
463 k += (ix>>23)-127;
464 ix &= 0x007fffff;
465 i = (ix+(0x95f64<<3))&0x800000;
466 SET_FLOAT_WORD(x,ix|(i^0x3f800000)); /* normalize x or x/2 */
467 k += (i>>23);
468 f = x-(float)1.0;
469 if((0x007fffff&(15+ix))<16) { /* |f| < 2**-20 */
470 if(f==zero) {
471 if(k==0)
472 return zero;
473 else
475 dk=(float)k;
476 return dk*ln2_hi+dk*ln2_lo;
479 R = f*f*((float)0.5-(float)0.33333333333333333*f);
480 if(k==0)
481 return f-R;
482 else
484 dk=(float)k;
485 return dk*ln2_hi-((R-dk*ln2_lo)-f);
488 s = f/((float)2.0+f);
489 dk = (float)k;
490 z = s*s;
491 i = ix-(0x6147a<<3);
492 w = z*z;
493 j = (0x6b851<<3)-ix;
494 t1= w*(Lg2+w*(Lg4+w*Lg6));
495 t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7)));
496 i |= j;
497 R = t2+t1;
498 if(i>0) {
499 hfsq=(float)0.5*f*f;
500 if(k==0)
501 return f-(hfsq-s*(hfsq+R));
502 else
503 return dk*ln2_hi-((hfsq-(s*(hfsq+R)+dk*ln2_lo))-f);
504 } else {
505 if(k==0)
506 return f-s*(f-R);
507 else
508 return dk*ln2_hi-((s*(f-R)-dk*ln2_lo)-f);
513 /* Logarithm for 10th base,
514 taken from glibc-2.8 */
516 static const float
517 ivln10 = 4.3429449201e-01, /* 0x3ede5bd9 */
518 log10_2hi = 3.0102920532e-01, /* 0x3e9a2080 */
519 log10_2lo = 7.9034151668e-07; /* 0x355427db */
521 float rb_log10(float x)
523 float y,z;
524 int32_t i,k,hx;
526 GET_FLOAT_WORD(hx,x);
528 k=0;
529 if (hx < 0x00800000) { /* x < 2**-126 */
530 if ((hx&0x7fffffff)==0)
531 return -two25/(x-x); /* log(+-0)=-inf */
532 if (hx<0) return (x-x)/(x-x); /* log(-#) = NaN */
533 k -= 25; x *= two25; /* subnormal number, scale up x */
534 GET_FLOAT_WORD(hx,x);
536 if (hx >= 0x7f800000) return x+x;
537 k += (hx>>23)-127;
538 i = ((uint32_t)k&0x80000000)>>31;
539 hx = (hx&0x007fffff)|((0x7f-i)<<23);
540 y = (float)(k+i);
541 SET_FLOAT_WORD(x,hx);
542 z = y*log10_2lo + ivln10*rb_log(x);
543 return z+y*log10_2hi;
547 /* Power function,
548 Taken from glibc-2.8 */
550 int rb_isinf(float x)
552 int32_t ix, t;
553 GET_FLOAT_WORD(ix,x);
554 t = ix & 0x7fffffff;
555 t ^= 0x7f800000;
556 t |= -t;
557 return ~(t >> 31) & (ix >> 30);
560 float rb_copysignf(float x, float y)
562 uint32_t ix, iy;
563 GET_FLOAT_WORD(ix,x);
564 GET_FLOAT_WORD(iy,y);
565 SET_FLOAT_WORD(x,(ix&0x7fffffff)|(iy&0x80000000));
566 return x;
569 static const float
570 huge = 1.0e+30,
571 tiny = 1.0e-30,
572 twom25 = 2.9802322388e-08; /* 0x33000000 */
574 float rb_scalbnf(float x, int n)
576 int32_t k, ix;
577 GET_FLOAT_WORD(ix,x);
578 k = (ix&0x7f800000)>>23; /* extract exponent */
579 if (k==0) { /* 0 or subnormal x */
580 if ((ix&0x7fffffff)==0) return x; /* +-0 */
581 x *= two25;
582 GET_FLOAT_WORD(ix,x);
583 k = ((ix&0x7f800000)>>23) - 25;
585 if (k==0xff) return x+x; /* NaN or Inf */
586 k = k+n;
587 if (n> 50000 || k > 0xfe)
588 return huge*rb_copysignf(huge,x); /* overflow */
589 if (n< -50000)
590 return tiny*rb_copysignf(tiny,x); /*underflow*/
591 if (k > 0) /* normal result */
592 {SET_FLOAT_WORD(x,(ix&0x807fffff)|(k<<23)); return x;}
593 if (k <= -25)
594 return tiny*rb_copysignf(tiny,x); /*underflow*/
595 k += 25; /* subnormal result */
596 SET_FLOAT_WORD(x,(ix&0x807fffff)|(k<<23));
597 return x*twom25;
601 static const float
602 bp[] = {1.0, 1.5,},
603 dp_h[] = { 0.0, 5.84960938e-01,}, /* 0x3f15c000 */
604 dp_l[] = { 0.0, 1.56322085e-06,}, /* 0x35d1cfdc */
605 one = 1.0,
606 two = 2.0,
607 two24 = 16777216.0, /* 0x4b800000 */
608 /* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */
609 L1 = 6.0000002384e-01, /* 0x3f19999a */
610 L2 = 4.2857143283e-01, /* 0x3edb6db7 */
611 L3 = 3.3333334327e-01, /* 0x3eaaaaab */
612 L4 = 2.7272811532e-01, /* 0x3e8ba305 */
613 L5 = 2.3066075146e-01, /* 0x3e6c3255 */
614 L6 = 2.0697501302e-01, /* 0x3e53f142 */
615 P1 = 1.6666667163e-01, /* 0x3e2aaaab */
616 P2 = -2.7777778450e-03, /* 0xbb360b61 */
617 P3 = 6.6137559770e-05, /* 0x388ab355 */
618 P4 = -1.6533901999e-06, /* 0xb5ddea0e */
619 P5 = 4.1381369442e-08; /* 0x3331bb4c */
621 static const float
622 lg2 = 6.9314718246e-01, /* 0x3f317218 */
623 lg2_h = 6.93145752e-01, /* 0x3f317200 */
624 lg2_l = 1.42860654e-06, /* 0x35bfbe8c */
625 ovt = 4.2995665694e-08, /* -(128-log2(ovfl+.5ulp)) */
626 cp = 9.6179670095e-01, /* 0x3f76384f =2/(3ln2) */
627 cp_h = 9.6179199219e-01, /* 0x3f763800 =head of cp */
628 cp_l = 4.7017383622e-06, /* 0x369dc3a0 =tail of cp_h */
629 ivln2 = 1.4426950216e+00, /* 0x3fb8aa3b =1/ln2 */
630 ivln2_h = 1.4426879883e+00, /* 0x3fb8aa00 =16b 1/ln2*/
631 ivln2_l = 7.0526075433e-06; /* 0x36eca570 =1/ln2 tail*/
633 float rb_pow(float x, float y)
635 float z, ax, z_h, z_l, p_h, p_l;
636 float y1, t1, t2, r, s, t, u, v, w;
637 int32_t i, j, k, yisint, n;
638 int32_t hx, hy, ix, iy, is;
640 GET_FLOAT_WORD(hx,x);
641 GET_FLOAT_WORD(hy,y);
642 ix = hx&0x7fffffff;
643 iy = hy&0x7fffffff;
645 /* y==zero: x**0 = 1 */
646 if(iy==0) return one;
648 /* x==+-1 */
649 if(x == 1.0) return one;
650 if(x == -1.0 && rb_isinf(y)) return one;
652 /* +-NaN return x+y */
653 if(ix > 0x7f800000 || iy > 0x7f800000)
654 return x+y;
656 /* determine if y is an odd int when x < 0
657 * yisint = 0 ... y is not an integer
658 * yisint = 1 ... y is an odd int
659 * yisint = 2 ... y is an even int
661 yisint = 0;
662 if(hx<0) {
663 if(iy>=0x4b800000) yisint = 2; /* even integer y */
664 else if(iy>=0x3f800000) {
665 k = (iy>>23)-0x7f; /* exponent */
666 j = iy>>(23-k);
667 if((j<<(23-k))==iy) yisint = 2-(j&1);
671 /* special value of y */
672 if (iy==0x7f800000) { /* y is +-inf */
673 if (ix==0x3f800000)
674 return y - y; /* inf**+-1 is NaN */
675 else if (ix > 0x3f800000)/* (|x|>1)**+-inf = inf,0 */
676 return (hy>=0)? y: zero;
677 else /* (|x|<1)**-,+inf = inf,0 */
678 return (hy<0)?-y: zero;
680 if(iy==0x3f800000) { /* y is +-1 */
681 if(hy<0) return one/x; else return x;
683 if(hy==0x40000000) return x*x; /* y is 2 */
684 if(hy==0x3f000000) { /* y is 0.5 */
685 if(hx>=0) /* x >= +0 */
686 return rb_sqrt(x);
689 ax = rb_fabs(x);
690 /* special value of x */
691 if(ix==0x7f800000||ix==0||ix==0x3f800000){
692 z = ax; /*x is +-0,+-inf,+-1*/
693 if(hy<0) z = one/z; /* z = (1/|x|) */
694 if(hx<0) {
695 if(((ix-0x3f800000)|yisint)==0) {
696 z = (z-z)/(z-z); /* (-1)**non-int is NaN */
697 } else if(yisint==1)
698 z = -z; /* (x<0)**odd = -(|x|**odd) */
700 return z;
703 /* (x<0)**(non-int) is NaN */
704 if(((((uint32_t)hx>>31)-1)|yisint)==0) return (x-x)/(x-x);
706 /* |y| is huge */
707 if(iy>0x4d000000) { /* if |y| > 2**27 */
708 /* over/underflow if x is not close to one */
709 if(ix<0x3f7ffff8) return (hy<0)? huge*huge:tiny*tiny;
710 if(ix>0x3f800007) return (hy>0)? huge*huge:tiny*tiny;
711 /* now |1-x| is tiny <= 2**-20, suffice to compute
712 log(x) by x-x^2/2+x^3/3-x^4/4 */
713 t = x-1; /* t has 20 trailing zeros */
714 w = (t*t)*((float)0.5-t*((float)0.333333333333-t*(float)0.25));
715 u = ivln2_h*t; /* ivln2_h has 16 sig. bits */
716 v = t*ivln2_l-w*ivln2;
717 t1 = u+v;
718 GET_FLOAT_WORD(is,t1);
719 SET_FLOAT_WORD(t1,is&0xfffff000);
720 t2 = v-(t1-u);
721 } else {
722 float s2, s_h, s_l, t_h, t_l;
723 n = 0;
724 /* take care subnormal number */
725 if(ix<0x00800000)
726 {ax *= two24; n -= 24; GET_FLOAT_WORD(ix,ax); }
727 n += ((ix)>>23)-0x7f;
728 j = ix&0x007fffff;
729 /* determine interval */
730 ix = j|0x3f800000; /* normalize ix */
731 if(j<=0x1cc471) k=0; /* |x|<sqrt(3/2) */
732 else if(j<0x5db3d7) k=1; /* |x|<sqrt(3) */
733 else {k=0;n+=1;ix -= 0x00800000;}
734 SET_FLOAT_WORD(ax,ix);
736 /* compute s = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */
737 u = ax-bp[k]; /* bp[0]=1.0, bp[1]=1.5 */
738 v = one/(ax+bp[k]);
739 s = u*v;
740 s_h = s;
741 GET_FLOAT_WORD(is,s_h);
742 SET_FLOAT_WORD(s_h,is&0xfffff000);
743 /* t_h=ax+bp[k] High */
744 SET_FLOAT_WORD(t_h,((ix>>1)|0x20000000)+0x0040000+(k<<21));
745 t_l = ax - (t_h-bp[k]);
746 s_l = v*((u-s_h*t_h)-s_h*t_l);
747 /* compute log(ax) */
748 s2 = s*s;
749 r = s2*s2*(L1+s2*(L2+s2*(L3+s2*(L4+s2*(L5+s2*L6)))));
750 r += s_l*(s_h+s);
751 s2 = s_h*s_h;
752 t_h = (float)3.0+s2+r;
753 GET_FLOAT_WORD(is,t_h);
754 SET_FLOAT_WORD(t_h,is&0xfffff000);
755 t_l = r-((t_h-(float)3.0)-s2);
756 /* u+v = s*(1+...) */
757 u = s_h*t_h;
758 v = s_l*t_h+t_l*s;
759 /* 2/(3log2)*(s+...) */
760 p_h = u+v;
761 GET_FLOAT_WORD(is,p_h);
762 SET_FLOAT_WORD(p_h,is&0xfffff000);
763 p_l = v-(p_h-u);
764 z_h = cp_h*p_h; /* cp_h+cp_l = 2/(3*log2) */
765 z_l = cp_l*p_h+p_l*cp+dp_l[k];
766 /* log2(ax) = (s+..)*2/(3*log2) = n + dp_h + z_h + z_l */
767 t = (float)n;
768 t1 = (((z_h+z_l)+dp_h[k])+t);
769 GET_FLOAT_WORD(is,t1);
770 SET_FLOAT_WORD(t1,is&0xfffff000);
771 t2 = z_l-(((t1-t)-dp_h[k])-z_h);
774 s = one; /* s (sign of result -ve**odd) = -1 else = 1 */
775 if(((((uint32_t)hx>>31)-1)|(yisint-1))==0)
776 s = -one; /* (-ve)**(odd int) */
778 /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */
779 GET_FLOAT_WORD(is,y);
780 SET_FLOAT_WORD(y1,is&0xfffff000);
781 p_l = (y-y1)*t1+y*t2;
782 p_h = y1*t1;
783 z = p_l+p_h;
784 GET_FLOAT_WORD(j,z);
785 if (j>0x43000000) /* if z > 128 */
786 return s*huge*huge; /* overflow */
787 else if (j==0x43000000) { /* if z == 128 */
788 if(p_l+ovt>z-p_h) return s*huge*huge; /* overflow */
790 else if ((j&0x7fffffff)>0x43160000) /* z <= -150 */
791 return s*tiny*tiny; /* underflow */
792 else if ((uint32_t) j==0xc3160000){ /* z == -150 */
793 if(p_l<=z-p_h) return s*tiny*tiny; /* underflow */
796 * compute 2**(p_h+p_l)
798 i = j&0x7fffffff;
799 k = (i>>23)-0x7f;
800 n = 0;
801 if(i>0x3f000000) { /* if |z| > 0.5, set n = [z+0.5] */
802 n = j+(0x00800000>>(k+1));
803 k = ((n&0x7fffffff)>>23)-0x7f; /* new k for n */
804 SET_FLOAT_WORD(t,n&~(0x007fffff>>k));
805 n = ((n&0x007fffff)|0x00800000)>>(23-k);
806 if(j<0) n = -n;
807 p_h -= t;
809 t = p_l+p_h;
810 GET_FLOAT_WORD(is,t);
811 SET_FLOAT_WORD(t,is&0xfffff000);
812 u = t*lg2_h;
813 v = (p_l-(t-p_h))*lg2+t*lg2_l;
814 z = u+v;
815 w = v-(z-u);
816 t = z*z;
817 t1 = z - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
818 r = (z*t1)/(t1-two)-(w+z*w);
819 z = one-(r-z);
820 GET_FLOAT_WORD(j,z);
821 j += (n<<23);
822 if((j>>23)<=0) z = rb_scalbnf(z,n); /* subnormal output */
823 else SET_FLOAT_WORD(z,j);
824 return s*z;
829 /* Square root function, original. */
830 float rb_sqrt(float x)
832 float z;
833 int32_t sign = (int)0x80000000;
834 int32_t ix,s,q,m,t,i;
835 uint32_t r;
837 GET_FLOAT_WORD(ix,x);
839 /* take care of Inf and NaN */
840 if((ix&0x7f800000)==0x7f800000) {
841 return x*x+x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf
842 sqrt(-inf)=sNaN */
844 /* take care of zero */
845 if(ix<=0) {
846 if((ix&(~sign))==0) return x;/* sqrt(+-0) = +-0 */
847 else if(ix<0)
848 return (x-x)/(x-x); /* sqrt(-ve) = sNaN */
850 /* normalize x */
851 m = (ix>>23);
852 if(m==0) { /* subnormal x */
853 for(i=0;(ix&0x00800000)==0;i++) ix<<=1;
854 m -= i-1;
856 m -= 127; /* unbias exponent */
857 ix = (ix&0x007fffff)|0x00800000;
858 if(m&1) /* odd m, double x to make it even */
859 ix += ix;
860 m >>= 1; /* m = [m/2] */
862 /* generate sqrt(x) bit by bit */
863 ix += ix;
864 q = s = 0; /* q = sqrt(x) */
865 r = 0x01000000; /* r = moving bit from right to left */
867 while(r!=0) {
868 t = s+r;
869 if(t<=ix) {
870 s = t+r;
871 ix -= t;
872 q += r;
874 ix += ix;
875 r>>=1;
878 /* use floating add to find out rounding direction */
879 if(ix!=0) {
880 z = one-tiny; /* trigger inexact flag */
881 if (z>=one) {
882 z = one+tiny;
883 if (z>one)
884 q += 2;
885 else
886 q += (q&1);
889 ix = (q>>1)+0x3f000000;
890 ix += (m <<23);
891 SET_FLOAT_WORD(z,ix);
892 return z;
895 /* Absolute value,
896 taken from glibc-2.8 */
897 float rb_fabs(float x)
899 uint32_t ix;
900 GET_FLOAT_WORD(ix,x);
901 SET_FLOAT_WORD(x,ix&0x7fffffff);
902 return x;
905 /* Arc tangent,
906 taken from glibc-2.8. */
908 static const float atanhi[] = {
909 4.6364760399e-01, /* atan(0.5)hi 0x3eed6338 */
910 7.8539812565e-01, /* atan(1.0)hi 0x3f490fda */
911 9.8279368877e-01, /* atan(1.5)hi 0x3f7b985e */
912 1.5707962513e+00, /* atan(inf)hi 0x3fc90fda */
915 static const float atanlo[] = {
916 5.0121582440e-09, /* atan(0.5)lo 0x31ac3769 */
917 3.7748947079e-08, /* atan(1.0)lo 0x33222168 */
918 3.4473217170e-08, /* atan(1.5)lo 0x33140fb4 */
919 7.5497894159e-08, /* atan(inf)lo 0x33a22168 */
922 static const float aT[] = {
923 3.3333334327e-01, /* 0x3eaaaaaa */
924 -2.0000000298e-01, /* 0xbe4ccccd */
925 1.4285714924e-01, /* 0x3e124925 */
926 -1.1111110449e-01, /* 0xbde38e38 */
927 9.0908870101e-02, /* 0x3dba2e6e */
928 -7.6918758452e-02, /* 0xbd9d8795 */
929 6.6610731184e-02, /* 0x3d886b35 */
930 -5.8335702866e-02, /* 0xbd6ef16b */
931 4.9768779427e-02, /* 0x3d4bda59 */
932 -3.6531571299e-02, /* 0xbd15a221 */
933 1.6285819933e-02, /* 0x3c8569d7 */
937 float rb_atan(float x)
939 float w,s1,s2,z;
940 int32_t ix,hx,id;
942 GET_FLOAT_WORD(hx,x);
943 ix = hx&0x7fffffff;
944 if(ix>=0x50800000) { /* if |x| >= 2^34 */
945 if(ix>0x7f800000)
946 return x+x; /* NaN */
947 if(hx>0) return atanhi[3]+atanlo[3];
948 else return -atanhi[3]-atanlo[3];
949 } if (ix < 0x3ee00000) { /* |x| < 0.4375 */
950 if (ix < 0x31000000) { /* |x| < 2^-29 */
951 if(huge+x>one) return x; /* raise inexact */
953 id = -1;
954 } else {
955 x = rb_fabs(x);
956 if (ix < 0x3f980000) { /* |x| < 1.1875 */
957 if (ix < 0x3f300000) { /* 7/16 <=|x|<11/16 */
958 id = 0; x = ((float)2.0*x-one)/((float)2.0+x);
959 } else { /* 11/16<=|x|< 19/16 */
960 id = 1; x = (x-one)/(x+one);
962 } else {
963 if (ix < 0x401c0000) { /* |x| < 2.4375 */
964 id = 2; x = (x-(float)1.5)/(one+(float)1.5*x);
965 } else { /* 2.4375 <= |x| < 2^66 */
966 id = 3; x = -(float)1.0/x;
969 /* end of argument reduction */
970 z = x*x;
971 w = z*z;
972 /* break sum from i=0 to 10 aT[i]z**(i+1) into odd and even poly */
973 s1 = z*(aT[0]+w*(aT[2]+w*(aT[4]+w*(aT[6]+w*(aT[8]+w*aT[10])))));
974 s2 = w*(aT[1]+w*(aT[3]+w*(aT[5]+w*(aT[7]+w*aT[9]))));
975 if (id<0) return x - x*(s1+s2);
976 else {
977 z = atanhi[id] - ((x*(s1+s2) - atanlo[id]) - x);
978 return (hx<0)? -z:z;
982 /* Arc tangent from two variables, original. */
984 static const float
985 pi_o_4 = 7.8539818525e-01, /* 0x3f490fdb */
986 pi_o_2 = 1.5707963705e+00, /* 0x3fc90fdb */
987 pi = 3.1415927410e+00, /* 0x40490fdb */
988 pi_lo = -8.7422776573e-08; /* 0xb3bbbd2e */
990 float rb_atan2(float x, float y)
992 float z;
993 int32_t k,m,hx,hy,ix,iy;
995 GET_FLOAT_WORD(hx,x);
996 ix = hx&0x7fffffff;
997 GET_FLOAT_WORD(hy,y);
998 iy = hy&0x7fffffff;
999 if((ix>0x7f800000)||
1000 (iy>0x7f800000)) /* x or y is NaN */
1001 return x+y;
1002 if(hx==0x3f800000) return rb_atan(y); /* x=1.0 */
1003 m = ((hy>>31)&1)|((hx>>30)&2); /* 2*sign(x)+sign(y) */
1005 /* when y = 0 */
1006 if(iy==0) {
1007 switch(m) {
1008 case 0:
1009 case 1: return y; /* atan(+-0,+anything)=+-0 */
1010 case 2: return pi+tiny;/* atan(+0,-anything) = pi */
1011 case 3: return -pi-tiny;/* atan(-0,-anything) =-pi */
1014 /* when x = 0 */
1015 if(ix==0) return (hy<0)? -pi_o_2-tiny: pi_o_2+tiny;
1017 /* when x is INF */
1018 if(ix==0x7f800000) {
1019 if(iy==0x7f800000) {
1020 switch(m) {
1021 case 0: return pi_o_4+tiny;/* atan(+INF,+INF) */
1022 case 1: return -pi_o_4-tiny;/* atan(-INF,+INF) */
1023 case 2: return (float)3.0*pi_o_4+tiny;/*atan(+INF,-INF)*/
1024 case 3: return (float)-3.0*pi_o_4-tiny;/*atan(-INF,-INF)*/
1026 } else {
1027 switch(m) {
1028 case 0: return zero ; /* atan(+...,+INF) */
1029 case 1: return -zero ; /* atan(-...,+INF) */
1030 case 2: return pi+tiny ; /* atan(+...,-INF) */
1031 case 3: return -pi-tiny ; /* atan(-...,-INF) */
1035 /* when y is INF */
1036 if(iy==0x7f800000) return (hy<0)? -pi_o_2-tiny: pi_o_2+tiny;
1038 /* compute y/x */
1039 k = (iy-ix)>>23;
1040 if(k > 60) z=pi_o_2+(float)0.5*pi_lo; /* |y/x| > 2**60 */
1041 else if(hx<0&&k<-60) z=0.0; /* |y|/x < -2**60 */
1042 else z=rb_atan(rb_fabs(y/x)); /* safe to do y/x */
1043 switch (m) {
1044 case 0: return z ; /* atan(+,+) */
1045 case 1: {
1046 uint32_t zh;
1047 GET_FLOAT_WORD(zh,z);
1048 SET_FLOAT_WORD(z,zh ^ 0x80000000);
1050 return z ; /* atan(-,+) */
1051 case 2: return pi-(z-pi_lo);/* atan(+,-) */
1052 default: /* case 3 */
1053 return (z-pi_lo)-pi;/* atan(-,-) */
1058 /* Sine hyperbolic,
1059 taken from glibc-2.8 */
1061 static const float
1062 o_threshold = 8.8721679688e+01,/* 0x42b17180 */
1063 invln2 = 1.4426950216e+00,/* 0x3fb8aa3b */
1064 /* scaled coefficients related to expm1 */
1065 Q1 = -3.3333335072e-02, /* 0xbd088889 */
1066 Q2 = 1.5873016091e-03, /* 0x3ad00d01 */
1067 Q3 = -7.9365076090e-05, /* 0xb8a670cd */
1068 Q4 = 4.0082177293e-06, /* 0x36867e54 */
1069 Q5 = -2.0109921195e-07; /* 0xb457edbb */
1071 float rb_expm1(float x)
1073 float y,hi,lo,c=0,t,e,hxs,hfx,r1;
1074 int32_t k,xsb;
1075 uint32_t hx;
1077 GET_FLOAT_WORD(hx,x);
1078 xsb = hx&0x80000000; /* sign bit of x */
1079 if(xsb==0) y=x; else y= -x; /* y = |x| */
1080 hx &= 0x7fffffff; /* high word of |x| */
1082 /* filter out huge and non-finite argument */
1083 if(hx >= 0x4195b844) { /* if |x|>=27*ln2 */
1084 if(hx >= 0x42b17218) { /* if |x|>=88.721... */
1085 if(hx>0x7f800000)
1086 return x+x; /* NaN */
1087 if(hx==0x7f800000)
1088 return (xsb==0)? x:-1.0;/* exp(+-inf)={inf,-1} */
1089 if(x > o_threshold) return huge*huge; /* overflow */
1091 if(xsb!=0) { /* x < -27*ln2, return -1.0 with inexact */
1092 if(x+tiny<(float)0.0) /* raise inexact */
1093 return tiny-one; /* return -1 */
1097 /* argument reduction */
1098 if(hx > 0x3eb17218) { /* if |x| > 0.5 ln2 */
1099 if(hx < 0x3F851592) { /* and |x| < 1.5 ln2 */
1100 if(xsb==0)
1101 {hi = x - ln2_hi; lo = ln2_lo; k = 1;}
1102 else
1103 {hi = x + ln2_hi; lo = -ln2_lo; k = -1;}
1104 } else {
1105 k = invln2*x+((xsb==0)?(float)0.5:(float)-0.5);
1106 t = k;
1107 hi = x - t*ln2_hi; /* t*ln2_hi is exact here */
1108 lo = t*ln2_lo;
1110 x = hi - lo;
1111 c = (hi-x)-lo;
1113 else if(hx < 0x33000000) { /* when |x|<2**-25, return x */
1114 t = huge+x; /* return x with inexact flags when x!=0 */
1115 return x - (t-(huge+x));
1117 else k = 0;
1119 /* x is now in primary range */
1120 hfx = (float)0.5*x;
1121 hxs = x*hfx;
1122 r1 = one+hxs*(Q1+hxs*(Q2+hxs*(Q3+hxs*(Q4+hxs*Q5))));
1123 t = (float)3.0-r1*hfx;
1124 e = hxs*((r1-t)/((float)6.0 - x*t));
1125 if(k==0) return x - (x*e-hxs); /* c is 0 */
1126 else {
1127 e = (x*(e-c)-c);
1128 e -= hxs;
1129 if(k== -1) return (float)0.5*(x-e)-(float)0.5;
1130 if(k==1) {
1131 if(x < (float)-0.25) return -(float)2.0*(e-(x+(float)0.5));
1132 else return one+(float)2.0*(x-e);
1134 if (k <= -2 || k>56) { /* suffice to return exp(x)-1 */
1135 int32_t i;
1136 y = one-(e-x);
1137 GET_FLOAT_WORD(i,y);
1138 SET_FLOAT_WORD(y,i+(k<<23)); /* add k to y's exponent */
1139 return y-one;
1141 t = one;
1142 if(k<23) {
1143 int32_t i;
1144 SET_FLOAT_WORD(t,0x3f800000 - (0x1000000>>k)); /* t=1-2^-k */
1145 y = t-(e-x);
1146 GET_FLOAT_WORD(i,y);
1147 SET_FLOAT_WORD(y,i+(k<<23)); /* add k to y's exponent */
1148 } else {
1149 int32_t i;
1150 SET_FLOAT_WORD(t,((0x7f-k)<<23)); /* 2^-k */
1151 y = x-(e+t);
1152 y += one;
1153 GET_FLOAT_WORD(i,y);
1154 SET_FLOAT_WORD(y,i+(k<<23)); /* add k to y's exponent */
1157 return y;
1160 static const float shuge = 1.0e37;
1162 float rb_sinh(float x)
1164 float t,w,h;
1165 int32_t ix,jx;
1167 GET_FLOAT_WORD(jx,x);
1168 ix = jx&0x7fffffff;
1170 /* x is INF or NaN */
1171 if(ix>=0x7f800000) return x+x;
1173 h = 0.5;
1174 if (jx<0) h = -h;
1175 /* |x| in [0,22], return sign(x)*0.5*(E+E/(E+1))) */
1176 if (ix < 0x41b00000) { /* |x|<22 */
1177 if (ix<0x31800000) /* |x|<2**-28 */
1178 if(shuge+x>one) return x;/* sinh(tiny) = tiny with inexact */
1179 t = rb_expm1(rb_fabs(x));
1180 if(ix<0x3f800000) return h*((float)2.0*t-t*t/(t+one));
1181 return h*(t+t/(t+one));
1184 /* |x| in [22, log(maxdouble)] return 0.5*exp(|x|) */
1185 if (ix < 0x42b17180) return h*rb_exp(rb_fabs(x));
1187 /* |x| in [log(maxdouble), overflowthresold] */
1188 if (ix<=0x42b2d4fc) {
1189 w = rb_exp((float)0.5*rb_fabs(x));
1190 t = h*w;
1191 return t*w;
1194 /* |x| > overflowthresold, sinh(x) overflow */
1195 return x*shuge;
1199 /* Tangent,
1200 taken from glibc-2.8 */
1202 static const float
1203 pio4 = 7.8539812565e-01, /* 0x3f490fda */
1204 pio4lo= 3.7748947079e-08, /* 0x33222168 */
1205 T[] = {
1206 3.3333334327e-01, /* 0x3eaaaaab */
1207 1.3333334029e-01, /* 0x3e088889 */
1208 5.3968254477e-02, /* 0x3d5d0dd1 */
1209 2.1869488060e-02, /* 0x3cb327a4 */
1210 8.8632395491e-03, /* 0x3c11371f */
1211 3.5920790397e-03, /* 0x3b6b6916 */
1212 1.4562094584e-03, /* 0x3abede48 */
1213 5.8804126456e-04, /* 0x3a1a26c8 */
1214 2.4646313977e-04, /* 0x398137b9 */
1215 7.8179444245e-05, /* 0x38a3f445 */
1216 7.1407252108e-05, /* 0x3895c07a */
1217 -1.8558637748e-05, /* 0xb79bae5f */
1218 2.5907305826e-05, /* 0x37d95384 */
1221 float kernel_tan(float x, float y, int iy)
1223 float z,r,v,w,s;
1224 int32_t ix,hx;
1225 GET_FLOAT_WORD(hx,x);
1226 ix = hx&0x7fffffff; /* high word of |x| */
1227 if(ix<0x31800000) /* x < 2**-28 */
1228 {if((int)x==0) { /* generate inexact */
1229 if((ix|(iy+1))==0) return one/rb_fabs(x);
1230 else return (iy==1)? x: -one/x;
1233 if(ix>=0x3f2ca140) { /* |x|>=0.6744 */
1234 if(hx<0) {x = -x; y = -y;}
1235 z = pio4-x;
1236 w = pio4lo-y;
1237 x = z+w; y = 0.0;
1239 z = x*x;
1240 w = z*z;
1241 /* Break x^5*(T[1]+x^2*T[2]+...) into
1242 * x^5(T[1]+x^4*T[3]+...+x^20*T[11]) +
1243 * x^5(x^2*(T[2]+x^4*T[4]+...+x^22*[T12]))
1245 r = T[1]+w*(T[3]+w*(T[5]+w*(T[7]+w*(T[9]+w*T[11]))));
1246 v = z*(T[2]+w*(T[4]+w*(T[6]+w*(T[8]+w*(T[10]+w*T[12])))));
1247 s = z*x;
1248 r = y + z*(s*(r+v)+y);
1249 r += T[0]*s;
1250 w = x+r;
1251 if(ix>=0x3f2ca140) {
1252 v = (float)iy;
1253 return (float)(1-((hx>>30)&2))*(v-(float)2.0*(x-(w*w/(w+v)-r)));
1255 if(iy==1) return w;
1256 else { /* if allow error up to 2 ulp,
1257 simply return -1.0/(x+r) here */
1258 /* compute -1.0/(x+r) accurately */
1259 float a,t;
1260 int32_t i;
1261 z = w;
1262 GET_FLOAT_WORD(i,z);
1263 SET_FLOAT_WORD(z,i&0xfffff000);
1264 v = r-(z - x); /* z+v = r+x */
1265 t = a = -(float)1.0/w; /* a = -1.0/w */
1266 GET_FLOAT_WORD(i,t);
1267 SET_FLOAT_WORD(t,i&0xfffff000);
1268 s = (float)1.0+t*z;
1269 return t+a*(s+t*v);
1275 static const int init_jk[] = {4,7,9}; /* initial value for jk */
1277 static const float PIo2[] = {
1278 1.5703125000e+00, /* 0x3fc90000 */
1279 4.5776367188e-04, /* 0x39f00000 */
1280 2.5987625122e-05, /* 0x37da0000 */
1281 7.5437128544e-08, /* 0x33a20000 */
1282 6.0026650317e-11, /* 0x2e840000 */
1283 7.3896444519e-13, /* 0x2b500000 */
1284 5.3845816694e-15, /* 0x27c20000 */
1285 5.6378512969e-18, /* 0x22d00000 */
1286 8.3009228831e-20, /* 0x1fc40000 */
1287 3.2756352257e-22, /* 0x1bc60000 */
1288 6.3331015649e-25, /* 0x17440000 */
1291 static const float
1292 two8 = 2.5600000000e+02, /* 0x43800000 */
1293 twon8 = 3.9062500000e-03; /* 0x3b800000 */
1295 int kernel_rem_pio2(float *x, float *y, int e0, int nx, int prec, const int32_t *ipio2)
1297 int32_t jz,jx,jv,jp,jk,carry,n,iq[20],i,j,k,m,q0,ih;
1298 float z,fw,f[20],fq[20],q[20];
1300 /* initialize jk*/
1301 jk = init_jk[prec];
1302 jp = jk;
1304 /* determine jx,jv,q0, note that 3>q0 */
1305 jx = nx-1;
1306 jv = (e0-3)/8; if(jv<0) jv=0;
1307 q0 = e0-8*(jv+1);
1309 /* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */
1310 j = jv-jx; m = jx+jk;
1311 for(i=0;i<=m;i++,j++) f[i] = (j<0)? zero : (float) ipio2[j];
1313 /* compute q[0],q[1],...q[jk] */
1314 for (i=0;i<=jk;i++) {
1315 for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j]; q[i] = fw;
1318 jz = jk;
1319 recompute:
1320 /* distill q[] into iq[] reversingly */
1321 for(i=0,j=jz,z=q[jz];j>0;i++,j--) {
1322 fw = (float)((int32_t)(twon8* z));
1323 iq[i] = (int32_t)(z-two8*fw);
1324 z = q[j-1]+fw;
1327 /* compute n */
1328 z = rb_scalbnf(z,q0); /* actual value of z */
1329 z -= (float)8.0*rb_floor(z*(float)0.125); /* trim off integer >= 8 */
1330 n = (int32_t) z;
1331 z -= (float)n;
1332 ih = 0;
1333 if(q0>0) { /* need iq[jz-1] to determine n */
1334 i = (iq[jz-1]>>(8-q0)); n += i;
1335 iq[jz-1] -= i<<(8-q0);
1336 ih = iq[jz-1]>>(7-q0);
1338 else if(q0==0) ih = iq[jz-1]>>8;
1339 else if(z>=(float)0.5) ih=2;
1341 if(ih>0) { /* q > 0.5 */
1342 n += 1; carry = 0;
1343 for(i=0;i<jz ;i++) { /* compute 1-q */
1344 j = iq[i];
1345 if(carry==0) {
1346 if(j!=0) {
1347 carry = 1; iq[i] = 0x100- j;
1349 } else iq[i] = 0xff - j;
1351 if(q0>0) { /* rare case: chance is 1 in 12 */
1352 switch(q0) {
1353 case 1:
1354 iq[jz-1] &= 0x7f; break;
1355 case 2:
1356 iq[jz-1] &= 0x3f; break;
1359 if(ih==2) {
1360 z = one - z;
1361 if(carry!=0) z -= rb_scalbnf(one,q0);
1365 /* check if recomputation is needed */
1366 if(z==zero) {
1367 j = 0;
1368 for (i=jz-1;i>=jk;i--) j |= iq[i];
1369 if(j==0) { /* need recomputation */
1370 for(k=1;iq[jk-k]==0;k++); /* k = no. of terms needed */
1372 for(i=jz+1;i<=jz+k;i++) { /* add q[jz+1] to q[jz+k] */
1373 f[jx+i] = (float) ipio2[jv+i];
1374 for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j];
1375 q[i] = fw;
1377 jz += k;
1378 goto recompute;
1382 /* chop off zero terms */
1383 if(z==(float)0.0) {
1384 jz -= 1; q0 -= 8;
1385 while(iq[jz]==0) { jz--; q0-=8;}
1386 } else { /* break z into 8-bit if necessary */
1387 z = rb_scalbnf(z,-q0);
1388 if(z>=two8) {
1389 fw = (float)((int32_t)(twon8*z));
1390 iq[jz] = (int32_t)(z-two8*fw);
1391 jz += 1; q0 += 8;
1392 iq[jz] = (int32_t) fw;
1393 } else iq[jz] = (int32_t) z ;
1396 /* convert integer "bit" chunk to floating-point value */
1397 fw = rb_scalbnf(one,q0);
1398 for(i=jz;i>=0;i--) {
1399 q[i] = fw*(float)iq[i]; fw*=twon8;
1402 /* compute PIo2[0,...,jp]*q[jz,...,0] */
1403 for(i=jz;i>=0;i--) {
1404 for(fw=0.0,k=0;k<=jp&&k<=jz-i;k++) fw += PIo2[k]*q[i+k];
1405 fq[jz-i] = fw;
1408 /* compress fq[] into y[] */
1409 switch(prec) {
1410 case 0:
1411 fw = 0.0;
1412 for (i=jz;i>=0;i--) fw += fq[i];
1413 y[0] = (ih==0)? fw: -fw;
1414 break;
1415 case 1:
1416 case 2:
1417 fw = 0.0;
1418 for (i=jz;i>=0;i--) fw += fq[i];
1419 y[0] = (ih==0)? fw: -fw;
1420 fw = fq[0]-fw;
1421 for (i=1;i<=jz;i++) fw += fq[i];
1422 y[1] = (ih==0)? fw: -fw;
1423 break;
1424 case 3: /* painful */
1425 for (i=jz;i>0;i--) {
1426 fw = fq[i-1]+fq[i];
1427 fq[i] += fq[i-1]-fw;
1428 fq[i-1] = fw;
1430 for (i=jz;i>1;i--) {
1431 fw = fq[i-1]+fq[i];
1432 fq[i] += fq[i-1]-fw;
1433 fq[i-1] = fw;
1435 for (fw=0.0,i=jz;i>=2;i--) fw += fq[i];
1436 if(ih==0) {
1437 y[0] = fq[0]; y[1] = fq[1]; y[2] = fw;
1438 } else {
1439 y[0] = -fq[0]; y[1] = -fq[1]; y[2] = -fw;
1442 return n&7;
1446 static const int32_t two_over_pi[] = {
1447 0xA2, 0xF9, 0x83, 0x6E, 0x4E, 0x44, 0x15, 0x29, 0xFC,
1448 0x27, 0x57, 0xD1, 0xF5, 0x34, 0xDD, 0xC0, 0xDB, 0x62,
1449 0x95, 0x99, 0x3C, 0x43, 0x90, 0x41, 0xFE, 0x51, 0x63,
1450 0xAB, 0xDE, 0xBB, 0xC5, 0x61, 0xB7, 0x24, 0x6E, 0x3A,
1451 0x42, 0x4D, 0xD2, 0xE0, 0x06, 0x49, 0x2E, 0xEA, 0x09,
1452 0xD1, 0x92, 0x1C, 0xFE, 0x1D, 0xEB, 0x1C, 0xB1, 0x29,
1453 0xA7, 0x3E, 0xE8, 0x82, 0x35, 0xF5, 0x2E, 0xBB, 0x44,
1454 0x84, 0xE9, 0x9C, 0x70, 0x26, 0xB4, 0x5F, 0x7E, 0x41,
1455 0x39, 0x91, 0xD6, 0x39, 0x83, 0x53, 0x39, 0xF4, 0x9C,
1456 0x84, 0x5F, 0x8B, 0xBD, 0xF9, 0x28, 0x3B, 0x1F, 0xF8,
1457 0x97, 0xFF, 0xDE, 0x05, 0x98, 0x0F, 0xEF, 0x2F, 0x11,
1458 0x8B, 0x5A, 0x0A, 0x6D, 0x1F, 0x6D, 0x36, 0x7E, 0xCF,
1459 0x27, 0xCB, 0x09, 0xB7, 0x4F, 0x46, 0x3F, 0x66, 0x9E,
1460 0x5F, 0xEA, 0x2D, 0x75, 0x27, 0xBA, 0xC7, 0xEB, 0xE5,
1461 0xF1, 0x7B, 0x3D, 0x07, 0x39, 0xF7, 0x8A, 0x52, 0x92,
1462 0xEA, 0x6B, 0xFB, 0x5F, 0xB1, 0x1F, 0x8D, 0x5D, 0x08,
1463 0x56, 0x03, 0x30, 0x46, 0xFC, 0x7B, 0x6B, 0xAB, 0xF0,
1464 0xCF, 0xBC, 0x20, 0x9A, 0xF4, 0x36, 0x1D, 0xA9, 0xE3,
1465 0x91, 0x61, 0x5E, 0xE6, 0x1B, 0x08, 0x65, 0x99, 0x85,
1466 0x5F, 0x14, 0xA0, 0x68, 0x40, 0x8D, 0xFF, 0xD8, 0x80,
1467 0x4D, 0x73, 0x27, 0x31, 0x06, 0x06, 0x15, 0x56, 0xCA,
1468 0x73, 0xA8, 0xC9, 0x60, 0xE2, 0x7B, 0xC0, 0x8C, 0x6B,
1471 static const int32_t npio2_hw[] = {
1472 0x3fc90f00, 0x40490f00, 0x4096cb00, 0x40c90f00, 0x40fb5300, 0x4116cb00,
1473 0x412fed00, 0x41490f00, 0x41623100, 0x417b5300, 0x418a3a00, 0x4196cb00,
1474 0x41a35c00, 0x41afed00, 0x41bc7e00, 0x41c90f00, 0x41d5a000, 0x41e23100,
1475 0x41eec200, 0x41fb5300, 0x4203f200, 0x420a3a00, 0x42108300, 0x4216cb00,
1476 0x421d1400, 0x42235c00, 0x4229a500, 0x422fed00, 0x42363600, 0x423c7e00,
1477 0x4242c700, 0x42490f00
1481 * invpio2: 24 bits of 2/pi
1482 * pio2_1: first 17 bit of pi/2
1483 * pio2_1t: pi/2 - pio2_1
1484 * pio2_2: second 17 bit of pi/2
1485 * pio2_2t: pi/2 - (pio2_1+pio2_2)
1486 * pio2_3: third 17 bit of pi/2
1487 * pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3)
1490 static const float
1491 half = 5.0000000000e-01, /* 0x3f000000 */
1492 invpio2 = 6.3661980629e-01, /* 0x3f22f984 */
1493 pio2_1 = 1.5707855225e+00, /* 0x3fc90f80 */
1494 pio2_1t = 1.0804334124e-05, /* 0x37354443 */
1495 pio2_2 = 1.0804273188e-05, /* 0x37354400 */
1496 pio2_2t = 6.0770999344e-11, /* 0x2e85a308 */
1497 pio2_3 = 6.0770943833e-11, /* 0x2e85a300 */
1498 pio2_3t = 6.1232342629e-17; /* 0x248d3132 */
1500 int32_t rem_pio2(float x, float *y)
1502 float z,w,t,r,fn;
1503 float tx[3];
1504 int32_t e0,i,j,nx,n,ix,hx;
1506 GET_FLOAT_WORD(hx,x);
1507 ix = hx&0x7fffffff;
1508 if(ix<=0x3f490fd8) /* |x| ~<= pi/4 , no need for reduction */
1509 {y[0] = x; y[1] = 0; return 0;}
1510 if(ix<0x4016cbe4) { /* |x| < 3pi/4, special case with n=+-1 */
1511 if(hx>0) {
1512 z = x - pio2_1;
1513 if((ix&0xfffffff0)!=0x3fc90fd0) { /* 24+24 bit pi OK */
1514 y[0] = z - pio2_1t;
1515 y[1] = (z-y[0])-pio2_1t;
1516 } else { /* near pi/2, use 24+24+24 bit pi */
1517 z -= pio2_2;
1518 y[0] = z - pio2_2t;
1519 y[1] = (z-y[0])-pio2_2t;
1521 return 1;
1522 } else { /* negative x */
1523 z = x + pio2_1;
1524 if((ix&0xfffffff0)!=0x3fc90fd0) { /* 24+24 bit pi OK */
1525 y[0] = z + pio2_1t;
1526 y[1] = (z-y[0])+pio2_1t;
1527 } else { /* near pi/2, use 24+24+24 bit pi */
1528 z += pio2_2;
1529 y[0] = z + pio2_2t;
1530 y[1] = (z-y[0])+pio2_2t;
1532 return -1;
1535 if(ix<=0x43490f80) { /* |x| ~<= 2^7*(pi/2), medium size */
1536 t = rb_fabs(x);
1537 n = (int32_t) (t*invpio2+half);
1538 fn = (float)n;
1539 r = t-fn*pio2_1;
1540 w = fn*pio2_1t; /* 1st round good to 40 bit */
1541 if(n<32&&(int32_t)(ix&0xffffff00)!=npio2_hw[n-1]) {
1542 y[0] = r-w; /* quick check no cancellation */
1543 } else {
1544 uint32_t high;
1545 j = ix>>23;
1546 y[0] = r-w;
1547 GET_FLOAT_WORD(high,y[0]);
1548 i = j-((high>>23)&0xff);
1549 if(i>8) { /* 2nd iteration needed, good to 57 */
1550 t = r;
1551 w = fn*pio2_2;
1552 r = t-w;
1553 w = fn*pio2_2t-((t-r)-w);
1554 y[0] = r-w;
1555 GET_FLOAT_WORD(high,y[0]);
1556 i = j-((high>>23)&0xff);
1557 if(i>25) { /* 3rd iteration need, 74 bits acc */
1558 t = r; /* will cover all possible cases */
1559 w = fn*pio2_3;
1560 r = t-w;
1561 w = fn*pio2_3t-((t-r)-w);
1562 y[0] = r-w;
1566 y[1] = (r-y[0])-w;
1567 if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;}
1568 else return n;
1571 * all other (large) arguments
1573 if(ix>=0x7f800000) { /* x is inf or NaN */
1574 y[0]=y[1]=x-x; return 0;
1576 /* set z = scalbn(|x|,ilogb(x)-7) */
1577 e0 = (ix>>23)-134; /* e0 = ilogb(z)-7; */
1578 SET_FLOAT_WORD(z, ix - ((int32_t)(e0<<23)));
1579 for(i=0;i<2;i++) {
1580 tx[i] = (float)((int32_t)(z));
1581 z = (z-tx[i])*two8;
1583 tx[2] = z;
1584 nx = 3;
1585 while(tx[nx-1]==zero) nx--; /* skip zero term */
1586 n = kernel_rem_pio2(tx,y,e0,nx,2,two_over_pi);
1587 if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;}
1588 return n;
1591 float rb_tan(float x)
1593 float y[2],z=0.0;
1594 int32_t n, ix;
1596 GET_FLOAT_WORD(ix,x);
1598 /* |x| ~< pi/4 */
1599 ix &= 0x7fffffff;
1600 if(ix <= 0x3f490fda) return kernel_tan(x,z,1);
1602 /* tan(Inf or NaN) is NaN */
1603 else if (ix>=0x7f800000) return x-x; /* NaN */
1605 /* argument reduction needed */
1606 else {
1607 n = rem_pio2(x,y);
1608 return kernel_tan(y[0],y[1],1-((n&1)<<1)); /* 1 -- n even
1609 -1 -- n odd */
1615 /* Exponential function,
1616 taken from glibc-2.8
1617 As it uses double values and udefines some symbols,
1618 it was moved to the end of the source code */
1620 #define W52 (2.22044605e-16)
1621 #define W55 (2.77555756e-17)
1622 #define W58 (3.46944695e-18)
1623 #define W59 (1.73472348e-18)
1624 #define W60 (8.67361738e-19)
1625 const float __exp_deltatable[178] = {
1626 0*W60, 16558714*W60, -10672149*W59, 1441652*W60,
1627 -15787963*W55, 462888*W60, 7291806*W60, 1698880*W60,
1628 -14375103*W58, -2021016*W60, 728829*W60, -3759654*W60,
1629 3202123*W60, -10916019*W58, -251570*W60, -1043086*W60,
1630 8207536*W60, -409964*W60, -5993931*W60, -475500*W60,
1631 2237522*W60, 324170*W60, -244117*W60, 32077*W60,
1632 123907*W60, -1019734*W60, -143*W60, 813077*W60,
1633 743345*W60, 462461*W60, 629794*W60, 2125066*W60,
1634 -2339121*W60, -337951*W60, 9922067*W60, -648704*W60,
1635 149407*W60, -2687209*W60, -631608*W60, 2128280*W60,
1636 -4882082*W60, 2001360*W60, 175074*W60, 2923216*W60,
1637 -538947*W60, -1212193*W60, -1920926*W60, -1080577*W60,
1638 3690196*W60, 2643367*W60, 2911937*W60, 671455*W60,
1639 -1128674*W60, 593282*W60, -5219347*W60, -1941490*W60,
1640 11007953*W60, 239609*W60, -2969658*W60, -1183650*W60,
1641 942998*W60, 699063*W60, 450569*W60, -329250*W60,
1642 -7257875*W60, -312436*W60, 51626*W60, 555877*W60,
1643 -641761*W60, 1565666*W60, 884327*W60, -10960035*W60,
1644 -2004679*W60, -995793*W60, -2229051*W60, -146179*W60,
1645 -510327*W60, 1453482*W60, -3778852*W60, -2238056*W60,
1646 -4895983*W60, 3398883*W60, -252738*W60, 1230155*W60,
1647 346918*W60, 1109352*W60, 268941*W60, -2930483*W60,
1648 -1036263*W60, -1159280*W60, 1328176*W60, 2937642*W60,
1649 -9371420*W60, -6902650*W60, -1419134*W60, 1442904*W60,
1650 -1319056*W60, -16369*W60, 696555*W60, -279987*W60,
1651 -7919763*W60, 252741*W60, 459711*W60, -1709645*W60,
1652 354913*W60, 6025867*W60, -421460*W60, -853103*W60,
1653 -338649*W60, 962151*W60, 955965*W60, 784419*W60,
1654 -3633653*W60, 2277133*W60, -8847927*W52, 1223028*W60,
1655 5907079*W60, 623167*W60, 5142888*W60, 2599099*W60,
1656 1214280*W60, 4870359*W60, 593349*W60, -57705*W60,
1657 7761209*W60, -5564097*W60, 2051261*W60, 6216869*W60,
1658 4692163*W60, 601691*W60, -5264906*W60, 1077872*W60,
1659 -3205949*W60, 1833082*W60, 2081746*W60, -987363*W60,
1660 -1049535*W60, 2015244*W60, 874230*W60, 2168259*W60,
1661 -1740124*W60, -10068269*W60, -18242*W60, -3013583*W60,
1662 580601*W60, -2547161*W60, -535689*W60, 2220815*W60,
1663 1285067*W60, 2806933*W60, -983086*W60, -1729097*W60,
1664 -1162985*W60, -2561904*W60, 801988*W60, 244351*W60,
1665 1441893*W60, -7517981*W60, 271781*W60, -15021588*W60,
1666 -2341588*W60, -919198*W60, 1642232*W60, 4771771*W60,
1667 -1220099*W60, -3062372*W60, 628624*W60, 1278114*W60,
1668 13083513*W60, -10521925*W60, 3180310*W60, -1659307*W60,
1669 3543773*W60, 2501203*W60, 4151*W60, -340748*W60,
1670 -2285625*W60, 2495202*W60
1673 const double __exp_atable[355] /* __attribute__((mode(DF))) */ = {
1674 0.707722561055888932371, /* 0x0.b52d4e46605c27ffd */
1675 0.709106182438804188967, /* 0x0.b587fb96f75097ffb */
1676 0.710492508843861281234, /* 0x0.b5e2d649899167ffd */
1677 0.711881545564593931623, /* 0x0.b63dde74d36bdfffe */
1678 0.713273297897442870573, /* 0x0.b699142f945f87ffc */
1679 0.714667771153751463236, /* 0x0.b6f477909c4ea0001 */
1680 0.716064970655995725059, /* 0x0.b75008aec758f8004 */
1681 0.717464901723956938193, /* 0x0.b7abc7a0eea7e0002 */
1682 0.718867569715736398602, /* 0x0.b807b47e1586c7ff8 */
1683 0.720272979947266023271, /* 0x0.b863cf5d10e380003 */
1684 0.721681137825144314297, /* 0x0.b8c01855195c37ffb */
1685 0.723092048691992950199, /* 0x0.b91c8f7d213740004 */
1686 0.724505717938892290800, /* 0x0.b97934ec5002d0007 */
1687 0.725922150953176470431, /* 0x0.b9d608b9c92ea7ffc */
1688 0.727341353138962865022, /* 0x0.ba330afcc29e98003 */
1689 0.728763329918453162104, /* 0x0.ba903bcc8618b7ffc */
1690 0.730188086709957051568, /* 0x0.baed9b40591ba0000 */
1691 0.731615628948127705309, /* 0x0.bb4b296f931e30002 */
1692 0.733045962086486091436, /* 0x0.bba8e671a05617ff9 */
1693 0.734479091556371366251, /* 0x0.bc06d25dd49568001 */
1694 0.735915022857225542529, /* 0x0.bc64ed4bce8f6fff9 */
1695 0.737353761441304711410, /* 0x0.bcc33752f915d7ff9 */
1696 0.738795312814142124419, /* 0x0.bd21b08af98e78005 */
1697 0.740239682467211168593, /* 0x0.bd80590b65e9a8000 */
1698 0.741686875913991849885, /* 0x0.bddf30ebec4a10000 */
1699 0.743136898669507939299, /* 0x0.be3e38443c84e0007 */
1700 0.744589756269486091620, /* 0x0.be9d6f2c1d32a0002 */
1701 0.746045454254026796384, /* 0x0.befcd5bb59baf8004 */
1702 0.747503998175051087583, /* 0x0.bf5c6c09ca84c0003 */
1703 0.748965393601880857739, /* 0x0.bfbc322f5b18b7ff8 */
1704 0.750429646104262104698, /* 0x0.c01c2843f776fffff */
1705 0.751896761271877989160, /* 0x0.c07c4e5fa18b88002 */
1706 0.753366744698445112140, /* 0x0.c0dca49a5fb18fffd */
1707 0.754839601988627206827, /* 0x0.c13d2b0c444db0005 */
1708 0.756315338768691947122, /* 0x0.c19de1cd798578006 */
1709 0.757793960659406629066, /* 0x0.c1fec8f623723fffd */
1710 0.759275473314173443536, /* 0x0.c25fe09e8a0f47ff8 */
1711 0.760759882363831851927, /* 0x0.c2c128dedc88f8000 */
1712 0.762247193485956486805, /* 0x0.c322a1cf7d6e7fffa */
1713 0.763737412354726363781, /* 0x0.c3844b88cb9347ffc */
1714 0.765230544649828092739, /* 0x0.c3e626232bd8f7ffc */
1715 0.766726596071518051729, /* 0x0.c44831b719bf18002 */
1716 0.768225572321911687194, /* 0x0.c4aa6e5d12d078001 */
1717 0.769727479119219348810, /* 0x0.c50cdc2da64a37ffb */
1718 0.771232322196981678892, /* 0x0.c56f7b41744490001 */
1719 0.772740107296721268087, /* 0x0.c5d24bb1259e70004 */
1720 0.774250840160724651565, /* 0x0.c6354d95640dd0007 */
1721 0.775764526565368872643, /* 0x0.c6988106fec447fff */
1722 0.777281172269557396602, /* 0x0.c6fbe61eb1bd0ffff */
1723 0.778800783068235302750, /* 0x0.c75f7cf560942fffc */
1724 0.780323364758801041312, /* 0x0.c7c345a3f1983fffe */
1725 0.781848923151573727006, /* 0x0.c8274043594cb0002 */
1726 0.783377464064598849602, /* 0x0.c88b6cec94b3b7ff9 */
1727 0.784908993312207869935, /* 0x0.c8efcbb89cba27ffe */
1728 0.786443516765346961618, /* 0x0.c9545cc0a88c70003 */
1729 0.787981040257604625744, /* 0x0.c9b9201dc643bfffa */
1730 0.789521569657452682047, /* 0x0.ca1e15e92a5410007 */
1731 0.791065110849462849192, /* 0x0.ca833e3c1ae510005 */
1732 0.792611669712891875319, /* 0x0.cae8992fd84667ffd */
1733 0.794161252150049179450, /* 0x0.cb4e26ddbc207fff8 */
1734 0.795713864077794763584, /* 0x0.cbb3e75f301b60003 */
1735 0.797269511407239561694, /* 0x0.cc19dacd978cd8002 */
1736 0.798828200086368567220, /* 0x0.cc8001427e55d7ffb */
1737 0.800389937624300440456, /* 0x0.cce65ade24d360006 */
1738 0.801954725261124767840, /* 0x0.cd4ce7a5de839fffb */
1739 0.803522573691593189330, /* 0x0.cdb3a7c79a678fffd */
1740 0.805093487311204114563, /* 0x0.ce1a9b563965ffffc */
1741 0.806667472122675088819, /* 0x0.ce81c26b838db8000 */
1742 0.808244534127439906441, /* 0x0.cee91d213f8428002 */
1743 0.809824679342317166307, /* 0x0.cf50ab9144d92fff9 */
1744 0.811407913793616542005, /* 0x0.cfb86dd5758c2ffff */
1745 0.812994243520784198882, /* 0x0.d0206407c20e20005 */
1746 0.814583674571603966162, /* 0x0.d0888e4223facfff9 */
1747 0.816176213022088536960, /* 0x0.d0f0ec9eb3f7c8002 */
1748 0.817771864936188586101, /* 0x0.d1597f377d6768002 */
1749 0.819370636400374108252, /* 0x0.d1c24626a46eafff8 */
1750 0.820972533518165570298, /* 0x0.d22b41865ff1e7ff9 */
1751 0.822577562404315121269, /* 0x0.d2947170f32ec7ff9 */
1752 0.824185729164559344159, /* 0x0.d2fdd60097795fff8 */
1753 0.825797039949601741075, /* 0x0.d3676f4fb796d0001 */
1754 0.827411500902565544264, /* 0x0.d3d13d78b5f68fffb */
1755 0.829029118181348834154, /* 0x0.d43b40960546d8001 */
1756 0.830649897953322891022, /* 0x0.d4a578c222a058000 */
1757 0.832273846408250750368, /* 0x0.d50fe617a3ba78005 */
1758 0.833900969738858188772, /* 0x0.d57a88b1218e90002 */
1759 0.835531274148056613016, /* 0x0.d5e560a94048f8006 */
1760 0.837164765846411529371, /* 0x0.d6506e1aac8078003 */
1761 0.838801451086016225394, /* 0x0.d6bbb1204074e0001 */
1762 0.840441336100884561780, /* 0x0.d72729d4c28518004 */
1763 0.842084427144139224814, /* 0x0.d792d8530e12b0001 */
1764 0.843730730487052604790, /* 0x0.d7febcb61273e7fff */
1765 0.845380252404570153833, /* 0x0.d86ad718c308dfff9 */
1766 0.847032999194574087728, /* 0x0.d8d727962c69d7fff */
1767 0.848688977161248581090, /* 0x0.d943ae49621ce7ffb */
1768 0.850348192619261200615, /* 0x0.d9b06b4d832ef8005 */
1769 0.852010651900976245816, /* 0x0.da1d5ebdc22220005 */
1770 0.853676361342631029337, /* 0x0.da8a88b555baa0006 */
1771 0.855345327311054837175, /* 0x0.daf7e94f965f98004 */
1772 0.857017556155879489641, /* 0x0.db6580a7c98f7fff8 */
1773 0.858693054267390953857, /* 0x0.dbd34ed9617befff8 */
1774 0.860371828028939855647, /* 0x0.dc4153ffc8b65fff9 */
1775 0.862053883854957292436, /* 0x0.dcaf90368bfca8004 */
1776 0.863739228154875360306, /* 0x0.dd1e0399328d87ffe */
1777 0.865427867361348468455, /* 0x0.dd8cae435d303fff9 */
1778 0.867119807911702289458, /* 0x0.ddfb9050b1cee8006 */
1779 0.868815056264353846599, /* 0x0.de6aa9dced8448001 */
1780 0.870513618890481399881, /* 0x0.ded9fb03db7320006 */
1781 0.872215502247877139094, /* 0x0.df4983e1380657ff8 */
1782 0.873920712852848668986, /* 0x0.dfb94490ffff77ffd */
1783 0.875629257204025623884, /* 0x0.e0293d2f1cb01fff9 */
1784 0.877341141814212965880, /* 0x0.e0996dd786fff0007 */
1785 0.879056373217612985183, /* 0x0.e109d6a64f5d57ffc */
1786 0.880774957955916648615, /* 0x0.e17a77b78e72a7ffe */
1787 0.882496902590150900078, /* 0x0.e1eb5127722cc7ff8 */
1788 0.884222213673356738383, /* 0x0.e25c63121fb0c8006 */
1789 0.885950897802399772740, /* 0x0.e2cdad93ec5340003 */
1790 0.887682961567391237685, /* 0x0.e33f30c925fb97ffb */
1791 0.889418411575228162725, /* 0x0.e3b0ecce2d05ffff9 */
1792 0.891157254447957902797, /* 0x0.e422e1bf727718006 */
1793 0.892899496816652704641, /* 0x0.e4950fb9713fc7ffe */
1794 0.894645145323828439008, /* 0x0.e50776d8b0e60fff8 */
1795 0.896394206626591749641, /* 0x0.e57a1739c8fadfffc */
1796 0.898146687421414902124, /* 0x0.e5ecf0f97c5798007 */
1797 0.899902594367530173098, /* 0x0.e660043464e378005 */
1798 0.901661934163603406867, /* 0x0.e6d3510747e150006 */
1799 0.903424713533971135418, /* 0x0.e746d78f06cd97ffd */
1800 0.905190939194458810123, /* 0x0.e7ba97e879c91fffc */
1801 0.906960617885092856864, /* 0x0.e82e92309390b0007 */
1802 0.908733756358986566306, /* 0x0.e8a2c6845544afffa */
1803 0.910510361377119825629, /* 0x0.e9173500c8abc7ff8 */
1804 0.912290439722343249336, /* 0x0.e98bddc30f98b0002 */
1805 0.914073998177417412765, /* 0x0.ea00c0e84bc4c7fff */
1806 0.915861043547953501680, /* 0x0.ea75de8db8094fffe */
1807 0.917651582652244779397, /* 0x0.eaeb36d09d3137ffe */
1808 0.919445622318405764159, /* 0x0.eb60c9ce4ed3dffff */
1809 0.921243169397334638073, /* 0x0.ebd697a43995b0007 */
1810 0.923044230737526172328, /* 0x0.ec4ca06fc7768fffa */
1811 0.924848813220121135342, /* 0x0.ecc2e44e865b6fffb */
1812 0.926656923710931002014, /* 0x0.ed39635df34e70006 */
1813 0.928468569126343790092, /* 0x0.edb01dbbc2f5b7ffa */
1814 0.930283756368834757725, /* 0x0.ee2713859aab57ffa */
1815 0.932102492359406786818, /* 0x0.ee9e44d9342870004 */
1816 0.933924784042873379360, /* 0x0.ef15b1d4635438005 */
1817 0.935750638358567643520, /* 0x0.ef8d5a94f60f50007 */
1818 0.937580062297704630580, /* 0x0.f0053f38f345cffff */
1819 0.939413062815381727516, /* 0x0.f07d5fde3a2d98001 */
1820 0.941249646905368053689, /* 0x0.f0f5bca2d481a8004 */
1821 0.943089821583810716806, /* 0x0.f16e55a4e497d7ffe */
1822 0.944933593864477061592, /* 0x0.f1e72b028a2827ffb */
1823 0.946780970781518460559, /* 0x0.f2603cd9fb5430001 */
1824 0.948631959382661205081, /* 0x0.f2d98b497d2a87ff9 */
1825 0.950486566729423554277, /* 0x0.f353166f63e3dffff */
1826 0.952344799896018723290, /* 0x0.f3ccde6a11ae37ffe */
1827 0.954206665969085765512, /* 0x0.f446e357f66120000 */
1828 0.956072172053890279009, /* 0x0.f4c12557964f0fff9 */
1829 0.957941325265908139014, /* 0x0.f53ba48781046fffb */
1830 0.959814132734539637840, /* 0x0.f5b66106555d07ffa */
1831 0.961690601603558903308, /* 0x0.f6315af2c2027fffc */
1832 0.963570739036113010927, /* 0x0.f6ac926b8aeb80004 */
1833 0.965454552202857141381, /* 0x0.f728078f7c5008002 */
1834 0.967342048278315158608, /* 0x0.f7a3ba7d66a908001 */
1835 0.969233234469444204768, /* 0x0.f81fab543e1897ffb */
1836 0.971128118008140250896, /* 0x0.f89bda33122c78007 */
1837 0.973026706099345495256, /* 0x0.f9184738d4cf97ff8 */
1838 0.974929006031422851235, /* 0x0.f994f284d3a5c0008 */
1839 0.976835024947348973265, /* 0x0.fa11dc35bc7820002 */
1840 0.978744770239899142285, /* 0x0.fa8f046b4fb7f8007 */
1841 0.980658249138918636210, /* 0x0.fb0c6b449ab1cfff9 */
1842 0.982575468959622777535, /* 0x0.fb8a10e1088fb7ffa */
1843 0.984496437054508843888, /* 0x0.fc07f5602d79afffc */
1844 0.986421160608523028820, /* 0x0.fc8618e0e55e47ffb */
1845 0.988349647107594098099, /* 0x0.fd047b83571b1fffa */
1846 0.990281903873210800357, /* 0x0.fd831d66f4c018002 */
1847 0.992217938695037382475, /* 0x0.fe01fead3320bfff8 */
1848 0.994157757657894713987, /* 0x0.fe811f703491e8006 */
1849 0.996101369488558541238, /* 0x0.ff007fd5744490005 */
1850 0.998048781093141101932, /* 0x0.ff801ffa9b9280007 */
1851 1.000000000000000000000, /* 0x1.00000000000000000 */
1852 1.001955033605393285965, /* 0x1.0080200565d29ffff */
1853 1.003913889319761887310, /* 0x1.0100802aa0e80fff0 */
1854 1.005876574715736104818, /* 0x1.01812090377240007 */
1855 1.007843096764807100351, /* 0x1.020201541aad7fff6 */
1856 1.009813464316352327214, /* 0x1.0283229c4c9820007 */
1857 1.011787683565730677817, /* 0x1.030484836910a000e */
1858 1.013765762469146736174, /* 0x1.0386272b9c077fffe */
1859 1.015747708536026694351, /* 0x1.04080ab526304fff0 */
1860 1.017733529475172815584, /* 0x1.048a2f412375ffff0 */
1861 1.019723232714418781378, /* 0x1.050c94ef7ad5e000a */
1862 1.021716825883923762690, /* 0x1.058f3be0f1c2d0004 */
1863 1.023714316605201180057, /* 0x1.06122436442e2000e */
1864 1.025715712440059545995, /* 0x1.06954e0fec63afff2 */
1865 1.027721021151397406936, /* 0x1.0718b98f41c92fff6 */
1866 1.029730250269221158939, /* 0x1.079c66d49bb2ffff1 */
1867 1.031743407506447551857, /* 0x1.082056011a9230009 */
1868 1.033760500517691527387, /* 0x1.08a487359ebd50002 */
1869 1.035781537016238873464, /* 0x1.0928fa93490d4fff3 */
1870 1.037806524719013578963, /* 0x1.09adb03b3e5b3000d */
1871 1.039835471338248051878, /* 0x1.0a32a84e9e5760004 */
1872 1.041868384612101516848, /* 0x1.0ab7e2eea5340ffff */
1873 1.043905272300907460835, /* 0x1.0b3d603ca784f0009 */
1874 1.045946142174331239262, /* 0x1.0bc3205a042060000 */
1875 1.047991002016745332165, /* 0x1.0c4923682a086fffe */
1876 1.050039859627715177527, /* 0x1.0ccf698898f3a000d */
1877 1.052092722826109660856, /* 0x1.0d55f2dce5d1dfffb */
1878 1.054149599440827866881, /* 0x1.0ddcbf86b09a5fff6 */
1879 1.056210497317612961855, /* 0x1.0e63cfa7abc97fffd */
1880 1.058275424318780855142, /* 0x1.0eeb23619c146fffb */
1881 1.060344388322010722446, /* 0x1.0f72bad65714bffff */
1882 1.062417397220589476718, /* 0x1.0ffa9627c38d30004 */
1883 1.064494458915699715017, /* 0x1.1082b577d0eef0003 */
1884 1.066575581342167566880, /* 0x1.110b18e893a90000a */
1885 1.068660772440545025953, /* 0x1.1193c09c267610006 */
1886 1.070750040138235936705, /* 0x1.121cacb4959befff6 */
1887 1.072843392435016474095, /* 0x1.12a5dd543cf36ffff */
1888 1.074940837302467588937, /* 0x1.132f529d59552000b */
1889 1.077042382749654914030, /* 0x1.13b90cb250d08fff5 */
1890 1.079148036789447484528, /* 0x1.14430bb58da3dfff9 */
1891 1.081257807444460983297, /* 0x1.14cd4fc984c4a000e */
1892 1.083371702785017154417, /* 0x1.1557d910df9c7000e */
1893 1.085489730853784307038, /* 0x1.15e2a7ae292d30002 */
1894 1.087611899742884524772, /* 0x1.166dbbc422d8c0004 */
1895 1.089738217537583819804, /* 0x1.16f9157586772ffff */
1896 1.091868692357631731528, /* 0x1.1784b4e533cacfff0 */
1897 1.094003332327482702577, /* 0x1.18109a360fc23fff2 */
1898 1.096142145591650907149, /* 0x1.189cc58b155a70008 */
1899 1.098285140311341168136, /* 0x1.1929370751ea50002 */
1900 1.100432324652149906842, /* 0x1.19b5eecdd79cefff0 */
1901 1.102583706811727015711, /* 0x1.1a42ed01dbdba000e */
1902 1.104739294993289488947, /* 0x1.1ad031c69a2eafff0 */
1903 1.106899097422573863281, /* 0x1.1b5dbd3f66e120003 */
1904 1.109063122341542140286, /* 0x1.1beb8f8fa8150000b */
1905 1.111231377994659874592, /* 0x1.1c79a8dac6ad0fff4 */
1906 1.113403872669181282605, /* 0x1.1d0809445a97ffffc */
1907 1.115580614653132185460, /* 0x1.1d96b0effc9db000e */
1908 1.117761612217810673898, /* 0x1.1e25a001332190000 */
1909 1.119946873713312474002, /* 0x1.1eb4d69bdb2a9fff1 */
1910 1.122136407473298902480, /* 0x1.1f4454e3bfae00006 */
1911 1.124330221845670330058, /* 0x1.1fd41afcbb48bfff8 */
1912 1.126528325196519908506, /* 0x1.2064290abc98c0001 */
1913 1.128730725913251964394, /* 0x1.20f47f31c9aa7000f */
1914 1.130937432396844410880, /* 0x1.21851d95f776dfff0 */
1915 1.133148453059692917203, /* 0x1.2216045b6784efffa */
1916 1.135363796355857157764, /* 0x1.22a733a6692ae0004 */
1917 1.137583470716100553249, /* 0x1.2338ab9b3221a0004 */
1918 1.139807484614418608939, /* 0x1.23ca6c5e27aadfff7 */
1919 1.142035846532929888057, /* 0x1.245c7613b7f6c0004 */
1920 1.144268564977221958089, /* 0x1.24eec8e06b035000c */
1921 1.146505648458203463465, /* 0x1.258164e8cea85fff8 */
1922 1.148747105501412235671, /* 0x1.26144a5180d380009 */
1923 1.150992944689175123667, /* 0x1.26a7793f5de2efffa */
1924 1.153243174560058870217, /* 0x1.273af1d712179000d */
1925 1.155497803703682491111, /* 0x1.27ceb43d81d42fff1 */
1926 1.157756840726344771440, /* 0x1.2862c097a3d29000c */
1927 1.160020294239811677834, /* 0x1.28f7170a74cf4fff1 */
1928 1.162288172883275239058, /* 0x1.298bb7bb0faed0004 */
1929 1.164560485298402170388, /* 0x1.2a20a2ce920dffff4 */
1930 1.166837240167474476460, /* 0x1.2ab5d86a4631ffff6 */
1931 1.169118446164539637555, /* 0x1.2b4b58b36d5220009 */
1932 1.171404112007080167155, /* 0x1.2be123cf786790002 */
1933 1.173694246390975415341, /* 0x1.2c7739e3c0aac000d */
1934 1.175988858069749065617, /* 0x1.2d0d9b15deb58fff6 */
1935 1.178287955789017793514, /* 0x1.2da4478b627040002 */
1936 1.180591548323240091978, /* 0x1.2e3b3f69fb794fffc */
1937 1.182899644456603782686, /* 0x1.2ed282d76421d0004 */
1938 1.185212252993012693694, /* 0x1.2f6a11f96c685fff3 */
1939 1.187529382762033236513, /* 0x1.3001ecf60082ffffa */
1940 1.189851042595508889847, /* 0x1.309a13f30f28a0004 */
1941 1.192177241354644978669, /* 0x1.31328716a758cfff7 */
1942 1.194507987909589896687, /* 0x1.31cb4686e1e85fffb */
1943 1.196843291137896336843, /* 0x1.32645269dfd04000a */
1944 1.199183159977805113226, /* 0x1.32fdaae604c39000f */
1945 1.201527603343041317132, /* 0x1.339750219980dfff3 */
1946 1.203876630171082595692, /* 0x1.3431424300e480007 */
1947 1.206230249419600664189, /* 0x1.34cb8170b3fee000e */
1948 1.208588470077065268869, /* 0x1.35660dd14dbd4fffc */
1949 1.210951301134513435915, /* 0x1.3600e78b6bdfc0005 */
1950 1.213318751604272271958, /* 0x1.369c0ec5c38ebfff2 */
1951 1.215690830512196507537, /* 0x1.373783a718d29000f */
1952 1.218067546930756250870, /* 0x1.37d3465662f480007 */
1953 1.220448909901335365929, /* 0x1.386f56fa770fe0008 */
1954 1.222834928513994334780, /* 0x1.390bb5ba5fc540004 */
1955 1.225225611877684750397, /* 0x1.39a862bd3c7a8fff3 */
1956 1.227620969111500981433, /* 0x1.3a455e2a37bcafffd */
1957 1.230021009336254911271, /* 0x1.3ae2a8287dfbefff6 */
1958 1.232425741726685064472, /* 0x1.3b8040df76f39fffa */
1959 1.234835175450728295084, /* 0x1.3c1e287682e48fff1 */
1960 1.237249319699482263931, /* 0x1.3cbc5f151b86bfff8 */
1961 1.239668183679933477545, /* 0x1.3d5ae4e2cc0a8000f */
1962 1.242091776620540377629, /* 0x1.3df9ba07373bf0006 */
1963 1.244520107762172811399, /* 0x1.3e98deaa0d8cafffe */
1964 1.246953186383919165383, /* 0x1.3f3852f32973efff0 */
1965 1.249391019292643401078, /* 0x1.3fd816ffc72b90001 */
1966 1.251833623164381181797, /* 0x1.40782b17863250005 */
1967 1.254280999953110153911, /* 0x1.41188f42caf400000 */
1968 1.256733161434815393410, /* 0x1.41b943b42945bfffd */
1969 1.259190116985283935980, /* 0x1.425a4893e5f10000a */
1970 1.261651875958665236542, /* 0x1.42fb9e0a2df4c0009 */
1971 1.264118447754797758244, /* 0x1.439d443f608c4fff9 */
1972 1.266589841787181258708, /* 0x1.443f3b5bebf850008 */
1973 1.269066067469190262045, /* 0x1.44e183883e561fff7 */
1974 1.271547134259576328224, /* 0x1.45841cecf7a7a0001 */
1975 1.274033051628237434048, /* 0x1.462707b2c43020009 */
1976 1.276523829025464573684, /* 0x1.46ca44023aa410007 */
1977 1.279019475999373156531, /* 0x1.476dd2045d46ffff0 */
1978 1.281520002043128991825, /* 0x1.4811b1e1f1f19000b */
1979 1.284025416692967214122, /* 0x1.48b5e3c3edd74fff4 */
1980 1.286535729509738823464, /* 0x1.495a67d3613c8fff7 */
1981 1.289050950070396384145, /* 0x1.49ff3e396e19d000b */
1982 1.291571087985403654081, /* 0x1.4aa4671f5b401fff1 */
1983 1.294096152842774794011, /* 0x1.4b49e2ae56d19000d */
1984 1.296626154297237043484, /* 0x1.4befb10fd84a3fff4 */
1985 1.299161101984141142272, /* 0x1.4c95d26d41d84fff8 */
1986 1.301701005575179204100, /* 0x1.4d3c46f01d9f0fff3 */
1987 1.304245874766450485904, /* 0x1.4de30ec21097d0003 */
1988 1.306795719266019562007, /* 0x1.4e8a2a0ccce3d0002 */
1989 1.309350548792467483458, /* 0x1.4f3198fa10346fff5 */
1990 1.311910373099227200545, /* 0x1.4fd95bb3be8cffffd */
1991 1.314475201942565174546, /* 0x1.50817263bf0e5fffb */
1992 1.317045045107389400535, /* 0x1.5129dd3418575000e */
1993 1.319619912422941299109, /* 0x1.51d29c4f01c54ffff */
1994 1.322199813675649204855, /* 0x1.527bafde83a310009 */
1995 1.324784758729532718739, /* 0x1.5325180cfb8b3fffd */
1996 1.327374757430096474625, /* 0x1.53ced504b2bd0fff4 */
1997 1.329969819671041886272, /* 0x1.5478e6f02775e0001 */
1998 1.332569955346704748651, /* 0x1.55234df9d8a59fff8 */
1999 1.335175174370685002822, /* 0x1.55ce0a4c5a6a9fff6 */
2000 1.337785486688218616860, /* 0x1.56791c1263abefff7 */
2001 1.340400902247843806217, /* 0x1.57248376aef21fffa */
2002 1.343021431036279800211, /* 0x1.57d040a420c0bfff3 */
2003 1.345647083048053138662, /* 0x1.587c53c5a630f0002 */
2004 1.348277868295411074918, /* 0x1.5928bd063fd7bfff9 */
2005 1.350913796821875845231, /* 0x1.59d57c9110ad60006 */
2006 1.353554878672557082439, /* 0x1.5a8292913d68cfffc */
2007 1.356201123929036356254, /* 0x1.5b2fff3212db00007 */
2008 1.358852542671913132777, /* 0x1.5bddc29edcc06fff3 */
2009 1.361509145047255398051, /* 0x1.5c8bdd032ed16000f */
2010 1.364170941142184734180, /* 0x1.5d3a4e8a5bf61fff4 */
2011 1.366837941171020309735, /* 0x1.5de9176042f1effff */
2012 1.369510155261156381121, /* 0x1.5e9837b062f4e0005 */
2013 1.372187593620959988833, /* 0x1.5f47afa69436cfff1 */
2014 1.374870266463378287715, /* 0x1.5ff77f6eb3f8cfffd */
2015 1.377558184010425845733, /* 0x1.60a7a734a9742fff9 */
2016 1.380251356531521533853, /* 0x1.6158272490016000c */
2017 1.382949794301995272203, /* 0x1.6208ff6a8978a000f */
2018 1.385653507605306700170, /* 0x1.62ba3032c0a280004 */
2019 1.388362506772382154503, /* 0x1.636bb9a994784000f */
2020 1.391076802081129493127, /* 0x1.641d9bfb29a7bfff6 */
2021 1.393796403973427855412, /* 0x1.64cfd7545928b0002 */
2022 1.396521322756352656542, /* 0x1.65826be167badfff8 */
2023 1.399251568859207761660, /* 0x1.663559cf20826000c */
2024 1.401987152677323100733, /* 0x1.66e8a14a29486fffc */
2025 1.404728084651919228815, /* 0x1.679c427f5a4b6000b */
2026 1.407474375243217723560, /* 0x1.68503d9ba0add000f */
2027 1.410226034922914983815, /* 0x1.690492cbf6303fff9 */
2028 1.412983074197955213304, /* 0x1.69b9423d7b548fff6 */
2031 /* All floating-point numbers can be put in one of these categories. */
2032 enum
2034 FP_NAN,
2035 # define FP_NAN FP_NAN
2036 FP_INFINITE,
2037 # define FP_INFINITE FP_INFINITE
2038 FP_ZERO,
2039 # define FP_ZERO FP_ZERO
2040 FP_SUBNORMAL,
2041 # define FP_SUBNORMAL FP_SUBNORMAL
2042 FP_NORMAL
2043 # define FP_NORMAL FP_NORMAL
2048 __fpclassifyf (float x)
2050 uint32_t wx;
2051 int retval = FP_NORMAL;
2053 GET_FLOAT_WORD (wx, x);
2054 wx &= 0x7fffffff;
2055 if (wx == 0)
2056 retval = FP_ZERO;
2057 else if (wx < 0x800000)
2058 retval = FP_SUBNORMAL;
2059 else if (wx >= 0x7f800000)
2060 retval = wx > 0x7f800000 ? FP_NAN : FP_INFINITE;
2062 return retval;
2067 __isinff (float x)
2069 int32_t ix,t;
2070 GET_FLOAT_WORD(ix,x);
2071 t = ix & 0x7fffffff;
2072 t ^= 0x7f800000;
2073 t |= -t;
2074 return ~(t >> 31) & (ix >> 30);
2077 /* Return nonzero value if arguments are unordered. */
2078 # define fpclassify(x) \
2079 (sizeof (x) == sizeof (float) ? __fpclassifyf (x) : __fpclassifyf (x))
2081 # ifndef isunordered
2082 # define isunordered(u, v) \
2083 (__extension__ \
2084 ({ __typeof__(u) __u = (u); __typeof__(v) __v = (v); \
2085 fpclassify (__u) == FP_NAN || fpclassify (__v) == FP_NAN; }))
2086 # endif
2088 /* Return nonzero value if X is less than Y. */
2089 # ifndef isless
2090 # define isless(x, y) \
2091 (__extension__ \
2092 ({ __typeof__(x) __x = (x); __typeof__(y) __y = (y); \
2093 !isunordered (__x, __y) && __x < __y; }))
2094 # endif
2096 /* Return nonzero value if X is greater than Y. */
2097 # ifndef isgreater
2098 # define isgreater(x, y) \
2099 (__extension__ \
2100 ({ __typeof__(x) __x = (x); __typeof__(y) __y = (y); \
2101 !isunordered (__x, __y) && __x > __y; }))
2102 # endif
2104 union ieee754_double
2106 double d;
2108 /* This is the IEEE 754 double-precision format. */
2109 struct
2111 #if defined(ROCKBOX_BIG_ENDIAN)
2112 unsigned int negative:1;
2113 unsigned int exponent:11;
2114 /* Together these comprise the mantissa. */
2115 unsigned int mantissa0:20;
2116 unsigned int mantissa1:32;
2117 #else
2118 # if __FLOAT_WORD_ORDER == __BIG_ENDIAN
2119 unsigned int mantissa0:20;
2120 unsigned int exponent:11;
2121 unsigned int negative:1;
2122 unsigned int mantissa1:32;
2123 # else
2124 /* Together these comprise the mantissa. */
2125 unsigned int mantissa1:32;
2126 unsigned int mantissa0:20;
2127 unsigned int exponent:11;
2128 unsigned int negative:1;
2129 # endif
2130 #endif /* Little endian. */
2131 } ieee;
2133 /* This format makes it easier to see if a NaN is a signalling NaN. */
2134 struct
2136 #if defined(ROCKBOX_BIG_ENDIAN)
2137 unsigned int negative:1;
2138 unsigned int exponent:11;
2139 unsigned int quiet_nan:1;
2140 /* Together these comprise the mantissa. */
2141 unsigned int mantissa0:19;
2142 unsigned int mantissa1:32;
2143 #else
2144 # if __FLOAT_WORD_ORDER == __BIG_ENDIAN
2145 unsigned int mantissa0:19;
2146 unsigned int quiet_nan:1;
2147 unsigned int exponent:11;
2148 unsigned int negative:1;
2149 unsigned int mantissa1:32;
2150 # else
2151 /* Together these comprise the mantissa. */
2152 unsigned int mantissa1:32;
2153 unsigned int mantissa0:19;
2154 unsigned int quiet_nan:1;
2155 unsigned int exponent:11;
2156 unsigned int negative:1;
2157 # endif
2158 #endif
2159 } ieee_nan;
2163 static const volatile float TWOM100 = 7.88860905e-31;
2164 static const volatile float TWO127 = 1.7014118346e+38;
2166 float rb_exp(float x)
2168 static const float himark = 88.72283935546875;
2169 static const float lomark = -103.972084045410;
2170 /* Check for usual case. */
2171 if (isless (x, himark) && isgreater (x, lomark))
2173 static const float THREEp42 = 13194139533312.0;
2174 static const float THREEp22 = 12582912.0;
2175 /* 1/ln(2). */
2176 #undef M_1_LN2
2177 static const float M_1_LN2 = 1.44269502163f;
2178 /* ln(2) */
2179 #undef M_LN2
2180 static const double M_LN2 = .6931471805599452862;
2182 int tval;
2183 double x22, t, result, dx;
2184 float n, delta;
2185 union ieee754_double ex2_u;
2186 #ifndef ROCKBOX
2187 fenv_t oldenv;
2189 feholdexcept (&oldenv);
2190 #endif
2192 #ifdef FE_TONEAREST
2193 fesetround (FE_TONEAREST);
2194 #endif
2196 /* Calculate n. */
2197 n = x * M_1_LN2 + THREEp22;
2198 n -= THREEp22;
2199 dx = x - n*M_LN2;
2201 /* Calculate t/512. */
2202 t = dx + THREEp42;
2203 t -= THREEp42;
2204 dx -= t;
2206 /* Compute tval = t. */
2207 tval = (int) (t * 512.0);
2209 if (t >= 0)
2210 delta = - __exp_deltatable[tval];
2211 else
2212 delta = __exp_deltatable[-tval];
2214 /* Compute ex2 = 2^n e^(t/512+delta[t]). */
2215 ex2_u.d = __exp_atable[tval+177];
2216 ex2_u.ieee.exponent += (int) n;
2218 /* Approximate e^(dx+delta) - 1, using a second-degree polynomial,
2219 with maximum error in [-2^-10-2^-28,2^-10+2^-28]
2220 less than 5e-11. */
2221 x22 = (0.5000000496709180453 * dx + 1.0000001192102037084) * dx + delta;
2223 /* Return result. */
2224 #ifndef ROCKBOX
2225 fesetenv (&oldenv);
2226 #endif
2228 result = x22 * ex2_u.d + ex2_u.d;
2229 return (float) result;
2231 /* Exceptional cases: */
2232 else if (isless (x, himark))
2234 if (__isinff (x))
2235 /* e^-inf == 0, with no error. */
2236 return 0;
2237 else
2238 /* Underflow */
2239 return TWOM100 * TWOM100;
2241 else
2242 /* Return x, if x is a NaN or Inf; or overflow, otherwise. */
2243 return TWO127*x;
2247 /* Division with rest, original. */
2248 div_t div(int x, int y)
2250 div_t result;
2252 result.quot = x / y;
2253 result.rem = x % y;
2255 /* Addition from glibc-2.8: */
2256 if(x >= 0 && result.rem < 0)
2258 result.quot += 1;
2259 result.rem -= y;
2262 return result;
2266 /* Placeholder function. */
2267 /* Originally defined in s_midi.c */
2268 void sys_listmididevs(void)
2272 /* Placeholder function,
2273 as we do sleep in the core thread
2274 and not in the scheduler function. */
2275 /* Originally defined in s_inter.c */
2276 void sys_microsleep(int microsec)
2278 (void) microsec;
2281 /* Get running time in milliseconds. */
2282 /* Originally defined in s_inter.c */
2283 extern uint64_t runningtime;
2284 t_time sys_getrealtime(void)
2286 return runningtime;
2289 /* Place holder, as we do no IPC. */
2290 /* Originally defined in s_inter.c */
2291 void glob_ping(void* dummy)
2293 (void) dummy;
2296 /* Call to quit. */
2297 /* Originally defined in s_inter.c */
2298 extern bool quit;
2299 void glob_quit(void* dummy)
2301 (void) dummy;
2303 static bool reentered = false;
2305 if(!reentered)
2307 reentered = true;
2309 /* Close audio subsystem. */
2310 sys_close_audio();
2312 /* Stop main loop. */
2313 quit = true;
2317 /* Open file. Originally in s_main.c */
2318 void glob_evalfile(t_pd *ignore, t_symbol *name, t_symbol *dir);
2319 void openit(const char *dirname, const char *filename)
2321 char dirbuf[MAXPDSTRING], *nameptr;
2323 /* Workaround: If the file resides in the root directory,
2324 add a trailing slash to prevent directory part
2325 of the filename from being removed. */
2326 char ffilename[MAXPDSTRING];
2327 ffilename[0] = '/';
2328 ffilename[1] = '\0';
2329 strcat(ffilename, filename);
2331 int fd = open_via_path(dirname, ffilename, "", dirbuf, &nameptr,
2332 MAXPDSTRING, 0);
2333 if (fd)
2335 close (fd);
2336 glob_evalfile(0, gensym(nameptr), gensym(dirbuf));
2338 else
2339 error("%s: can't open", filename);
2343 /* Get current working directory. */
2344 extern char* filename;
2345 char* rb_getcwd(char* buf, ssize_t size)
2347 char* end_of_dir = strrchr(filename, '/');
2349 /* Check whether buffer may hold enough data. */
2350 if(size < end_of_dir - filename)
2351 return NULL;
2353 /* Copy current working directory to buffer. */
2354 strncpy(buf, filename, end_of_dir - filename);
2355 return buf;
2359 /* Execute instructions supplied on command-line.
2360 Basically a placeholder, because the only
2361 command line argument we get is the name of the .pd file. */
2362 /* Originally defined in s_main.c */
2363 extern t_namelist* sys_openlist;
2364 void glob_initfromgui(void *dummy, t_symbol *s, int argc, t_atom *argv)
2366 (void) dummy;
2367 (void) s;
2368 (void) argc;
2369 (void) argv;
2371 t_namelist *nl;
2372 char cwd[MAXPDSTRING];
2374 /* Get current working directory. */
2375 rb_getcwd(cwd, MAXPDSTRING);
2377 /* open patches specifies with "-open" args */
2378 for(nl = sys_openlist; nl; nl = nl->nl_next)
2379 openit(cwd, nl->nl_string);
2380 namelist_free(sys_openlist);
2381 sys_openlist = 0;
2384 /* Fake GUI start. Originally in s_inter.c */
2385 static int defaultfontshit[] = {
2386 8, 5, 9, 10, 6, 10, 12, 7, 13, 14, 9, 17, 16, 10, 19, 24, 15, 28,
2387 24, 15, 28};
2388 extern t_binbuf* inbinbuf;
2389 int sys_startgui(const char *guidir)
2391 unsigned int i;
2392 t_atom zz[23];
2393 char cmdbuf[4*MAXPDSTRING];
2395 (void) guidir;
2397 inbinbuf = binbuf_new();
2399 if(!rb_getcwd(cmdbuf, MAXPDSTRING))
2400 strcpy(cmdbuf, ".");
2401 SETSYMBOL(zz, gensym(cmdbuf));
2402 for (i = 1; i < 22; i++)
2403 SETFLOAT(zz + i, defaultfontshit[i-1]);
2404 SETFLOAT(zz+22,0);
2406 glob_initfromgui(0, 0, 23, zz);
2408 return 0;
2411 /* Return default DAC block size. */
2412 /* Originally defined in s_main.c */
2413 int sys_getblksize(void)
2415 return (DEFDACBLKSIZE);
2418 /* Find library directory and set it. */
2419 void sys_findlibdir(const char* filename)
2421 char sbuf[MAXPDSTRING];
2423 (void) filename;
2425 /* Make current working directory the system library directory. */
2426 rb_getcwd(sbuf, MAXPDSTRING);
2427 sys_libdir = gensym(sbuf);