summaryrefslogtreecommitdiff
path: root/apps/plugins/pdbox/pdbox-func.c
diff options
context:
space:
mode:
Diffstat (limited to 'apps/plugins/pdbox/pdbox-func.c')
-rw-r--r--apps/plugins/pdbox/pdbox-func.c2428
1 files changed, 2428 insertions, 0 deletions
diff --git a/apps/plugins/pdbox/pdbox-func.c b/apps/plugins/pdbox/pdbox-func.c
new file mode 100644
index 0000000000..ee4a8fd166
--- /dev/null
+++ b/apps/plugins/pdbox/pdbox-func.c
@@ -0,0 +1,2428 @@
1/***************************************************************************
2 * __________ __ ___.
3 * Open \______ \ ____ ____ | | _\_ |__ _______ ___
4 * Source | _// _ \_/ ___\| |/ /| __ \ / _ \ \/ /
5 * Jukebox | | ( <_> ) \___| < | \_\ ( <_> > < <
6 * Firmware |____|_ /\____/ \___ >__|_ \|___ /\____/__/\_ \
7 * \/ \/ \/ \/ \/
8 * $Id$
9 *
10 * Copyright (C) 2009 Wincent Balin
11 *
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.
16 *
17 * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY OF ANY
18 * KIND, either express or implied.
19 *
20 ****************************************************************************/
21
22#include "plugin.h"
23#include "pdbox.h"
24#include "ctype.h"
25
26#include "m_pd.h"
27#include "s_stuff.h"
28
29
30/* This implementation of strncat is taken from lua plug-in. */
31
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. */
34char *rb_strncat(char *s, const char *t, size_t n)
35{
36 char *dest = s;
37 register char *max;
38 s += strlen(s);
39
40 if((max = s + n) == s)
41 goto strncat_fini;
42
43 while(true)
44 {
45 if(!(*s = *t))
46 break;
47 if(++s == max)
48 break;
49 ++t;
50
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
68 }
69
70 *s = 0;
71
72strncat_fini:
73 return dest;
74}
75
76
77/* Implementation of floor, original. */
78float rb_floor(float value)
79{
80 /* If value is negative, decrement value to match function's definition. */
81 if(value < 0.0)
82 {
83 value -= 1.0;
84 }
85
86 /* Truncate fractional part (convert to integer)
87 and afterwards convert back to double. */
88 return (float) ((int) value);
89}
90
91
92/* Implementation of strtod() and atof(),
93 taken from SanOS (http://www.jbox.dk/sanos/). */
94static int rb_errno = 0;
95
96double rb_strtod(const char *str, char **endptr)
97{
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;
106
107 /* Reset Rockbox errno -- W.B. */
108#ifdef ROCKBOX
109 rb_errno = 0;
110#endif
111
112 // Skip leading whitespace
113 while (isspace(*p)) p++;
114
115 // Handle optional sign
116 negative = 0;
117 switch (*p)
118 {
119 case '-': negative = 1; // Fall through to increment position
120 case '+': p++;
121 }
122
123 number = 0.;
124 exponent = 0;
125 num_digits = 0;
126 num_decimals = 0;
127
128 // Process string of digits
129 while (isdigit(*p))
130 {
131 number = number * 10. + (*p - '0');
132 p++;
133 num_digits++;
134 }
135
136 // Process decimal part
137 if (*p == '.')
138 {
139 p++;
140
141 while (isdigit(*p))
142 {
143 number = number * 10. + (*p - '0');
144 p++;
145 num_digits++;
146 num_decimals++;
147 }
148
149 exponent -= num_decimals;
150 }
151
152 if (num_digits == 0)
153 {
154#ifdef ROCKBOX
155 rb_errno = 1;
156#else
157 errno = ERANGE;
158#endif
159 return 0.0;
160 }
161
162 // Correct for sign
163 if (negative) number = -number;
164
165 // Process an exponent string
166 if (*p == 'e' || *p == 'E')
167 {
168 // Handle optional sign
169 negative = 0;
170 switch(*++p)
171 {
172 case '-': negative = 1; // Fall through to increment pos
173 case '+': p++;
174 }
175
176 // Process string of digits
177 n = 0;
178 while (isdigit(*p))
179 {
180 n = n * 10 + (*p - '0');
181 p++;
182 }
183
184 if (negative)
185 exponent -= n;
186 else
187 exponent += n;
188 }
189
190#ifndef ROCKBOX
191 if (exponent < DBL_MIN_EXP || exponent > DBL_MAX_EXP)
192 {
193 errno = ERANGE;
194 return HUGE_VAL;
195 }
196#endif
197
198 // Scale the result
199 p10 = 10.;
200 n = exponent;
201 if (n < 0) n = -n;
202 while (n)
203 {
204 if (n & 1)
205 {
206 if (exponent < 0)
207 number /= p10;
208 else
209 number *= p10;
210 }
211 n >>= 1;
212 p10 *= p10;
213 }
214
215#ifndef ROCKBOX
216 if (number == HUGE_VAL) errno = ERANGE;
217#endif
218 if (endptr) *endptr = p;
219
220 return number;
221}
222
223double rb_atof(const char *str)
224{
225 return rb_strtod(str, NULL);
226}
227
228
229/* Implementation of ftoa(), original. */
230void rb_ftoan(float f, char* out, int size)
231{
232 #define SBUFSIZE 12
233 char sbuf[SBUFSIZE];
234
235 /* Zero out string. */
236 *out = '\0';
237 size--;
238
239 /* Handle negative numbers. */
240 if(f < 0.0)
241 {
242 f = -f;
243 strcat(out, "-");
244 size--;
245 }
246
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;
253
254 /* Append integral part to output string. */
255 strcat(out, sbuf);
256 size -= int_part_len;
257
258 /* Check whether further content is possible. */
259 if(size <= 0)
260 return;
261
262 /* Append decimal point. */
263 strcat(out, ".");
264 size--;
265
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);
270
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;
276
277 /* Check whether output string still has enough space. */
278 if(size <= 0)
279 return;
280
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);
285
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);
290}
291
292
293/* Implementation of atol(), adapted from
294 the atoi() implementation in Rockbox. */
295long rb_atol(const char* str)
296{
297 long value = 0L;
298 long sign = 1L;
299
300 while (isspace(*str))
301 {
302 str++;
303 }
304
305 if ('-' == *str)
306 {
307 sign = -1L;
308 str++;
309 }
310 else if ('+' == *str)
311 {
312 str++;
313 }
314
315 while ('0' == *str)
316 {
317 str++;
318 }
319
320 while (isdigit(*str))
321 {
322 value = (value * 10L) + (*str - '0');
323 str++;
324 }
325
326 return value * sign;
327}
328
329
330/* Implementation of sin() and cos(),
331 adapted from http://lab.polygonal.de/2007/07/18/fast-and-accurate-sinecosine-approximation/
332*/
333
334float rb_sin(float rad)
335{
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;
341
342 if(rad < 0)
343 return (1.27323954 * rad + 0.405284735 * rad * rad);
344 else
345 return (1.27323954 * rad - 0.405284735 * rad * rad);
346}
347
348float rb_cos(float rad)
349{
350 /* Compute cosine: sin(x + PI/2) = cos(x) */
351 rad += 1.57079632;
352 if(rad > 3.14159265)
353 rad -= 6.28318531;
354
355 return rb_sin(rad);
356}
357
358
359/* Emulation of fscanf(fd, "%f", (float*) xxx);
360 Basically a reimplementation of rb_strtod() above. */
361int rb_fscanf_f(int fd, float* f)
362{
363 #define FSCANF_F_BUFSIZE 64
364 char buf[FSCANF_F_BUFSIZE];
365
366 /* Read line from file. */
367 int bytes_read = rb->read_line(fd, buf, FSCANF_F_BUFSIZE-1);
368
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';
374
375 /* Convert buffer to float. */
376 *f = rb_atof(buf);
377
378 /* If there was an error, no float was read. */
379 if(rb_errno)
380 return 0;
381
382 return 1;
383}
384
385/* Emulation of fprintf(fd, "%f\n", (float*) xxx); */
386int rb_fprintf_f(int fd, float f)
387{
388 #define FPRINTF_F_BUFSIZE 64
389 char buf[FPRINTF_F_BUFSIZE];
390 const char* next_line = "\n";
391
392 /* Convert float to string. */
393 rb_ftoan(f, buf, sizeof(buf)-1);
394
395 /* Add next line character. */
396 strcat(buf, next_line);
397
398 /* Write string into file. */
399 return write(fd, buf, strlen(buf));
400}
401
402
403/* Natural logarithm.
404 Taken from glibc-2.8 */
405static const float
406ln2_hi = 6.9313812256e-01, /* 0x3f317180 */
407ln2_lo = 9.0580006145e-06, /* 0x3717f7d1 */
408two25 = 3.355443200e+07, /* 0x4c000000 */
409Lg1 = 6.6666668653e-01, /* 3F2AAAAB */
410Lg2 = 4.0000000596e-01, /* 3ECCCCCD */
411Lg3 = 2.8571429849e-01, /* 3E924925 */
412Lg4 = 2.2222198546e-01, /* 3E638E29 */
413Lg5 = 1.8183572590e-01, /* 3E3A3325 */
414Lg6 = 1.5313838422e-01, /* 3E1CD04F */
415Lg7 = 1.4798198640e-01; /* 3E178897 */
416
417static const float zero = 0.0;
418
419/* A union which permits us to convert between a float and a 32 bit
420 int. */
421
422typedef union
423{
424 float value;
425 uint32_t word;
426} ieee_float_shape_type;
427
428/* Get a 32 bit int from a float. */
429
430#define GET_FLOAT_WORD(i,d) \
431do { \
432 ieee_float_shape_type gf_u; \
433 gf_u.value = (d); \
434 (i) = gf_u.word; \
435} while (0)
436
437/* Set a float from a 32 bit int. */
438
439#define SET_FLOAT_WORD(d,i) \
440do { \
441 ieee_float_shape_type sf_u; \
442 sf_u.word = (i); \
443 (d) = sf_u.value; \
444} while (0)
445
446
447float rb_log(float x)
448{
449 float hfsq, f, s, z, R, w, t1, t2, dk;
450 int32_t k, ix, i, j;
451
452 GET_FLOAT_WORD(ix,x);
453
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);
461 }
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
474 {
475 dk=(float)k;
476 return dk*ln2_hi+dk*ln2_lo;
477 }
478 }
479 R = f*f*((float)0.5-(float)0.33333333333333333*f);
480 if(k==0)
481 return f-R;
482 else
483 {
484 dk=(float)k;
485 return dk*ln2_hi-((R-dk*ln2_lo)-f);
486 }
487 }
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);
509 }
510}
511
512
513/* Logarithm for 10th base,
514 taken from glibc-2.8 */
515
516static const float
517ivln10 = 4.3429449201e-01, /* 0x3ede5bd9 */
518log10_2hi = 3.0102920532e-01, /* 0x3e9a2080 */
519log10_2lo = 7.9034151668e-07; /* 0x355427db */
520
521float rb_log10(float x)
522{
523 float y,z;
524 int32_t i,k,hx;
525
526 GET_FLOAT_WORD(hx,x);
527
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);
535 }
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;
544}
545
546
547/* Power function,
548 Taken from glibc-2.8 */
549
550int rb_isinf(float x)
551{
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);
558}
559
560float rb_copysignf(float x, float y)
561{
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;
567}
568
569static const float
570huge = 1.0e+30,
571tiny = 1.0e-30,
572twom25 = 2.9802322388e-08; /* 0x33000000 */
573
574float rb_scalbnf(float x, int n)
575{
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;
584 }
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;
598}
599
600
601static const float
602bp[] = {1.0, 1.5,},
603dp_h[] = { 0.0, 5.84960938e-01,}, /* 0x3f15c000 */
604dp_l[] = { 0.0, 1.56322085e-06,}, /* 0x35d1cfdc */
605one = 1.0,
606two = 2.0,
607two24 = 16777216.0, /* 0x4b800000 */
608 /* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */
609L1 = 6.0000002384e-01, /* 0x3f19999a */
610L2 = 4.2857143283e-01, /* 0x3edb6db7 */
611L3 = 3.3333334327e-01, /* 0x3eaaaaab */
612L4 = 2.7272811532e-01, /* 0x3e8ba305 */
613L5 = 2.3066075146e-01, /* 0x3e6c3255 */
614L6 = 2.0697501302e-01, /* 0x3e53f142 */
615P1 = 1.6666667163e-01, /* 0x3e2aaaab */
616P2 = -2.7777778450e-03, /* 0xbb360b61 */
617P3 = 6.6137559770e-05, /* 0x388ab355 */
618P4 = -1.6533901999e-06, /* 0xb5ddea0e */
619P5 = 4.1381369442e-08; /* 0x3331bb4c */
620
621static const float
622lg2 = 6.9314718246e-01, /* 0x3f317218 */
623lg2_h = 6.93145752e-01, /* 0x3f317200 */
624lg2_l = 1.42860654e-06, /* 0x35bfbe8c */
625ovt = 4.2995665694e-08, /* -(128-log2(ovfl+.5ulp)) */
626cp = 9.6179670095e-01, /* 0x3f76384f =2/(3ln2) */
627cp_h = 9.6179199219e-01, /* 0x3f763800 =head of cp */
628cp_l = 4.7017383622e-06, /* 0x369dc3a0 =tail of cp_h */
629ivln2 = 1.4426950216e+00, /* 0x3fb8aa3b =1/ln2 */
630ivln2_h = 1.4426879883e+00, /* 0x3fb8aa00 =16b 1/ln2*/
631ivln2_l = 7.0526075433e-06; /* 0x36eca570 =1/ln2 tail*/
632
633float rb_pow(float x, float y)
634{
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;
639
640 GET_FLOAT_WORD(hx,x);
641 GET_FLOAT_WORD(hy,y);
642 ix = hx&0x7fffffff;
643 iy = hy&0x7fffffff;
644
645 /* y==zero: x**0 = 1 */
646 if(iy==0) return one;
647
648 /* x==+-1 */
649 if(x == 1.0) return one;
650 if(x == -1.0 && rb_isinf(y)) return one;
651
652 /* +-NaN return x+y */
653 if(ix > 0x7f800000 || iy > 0x7f800000)
654 return x+y;
655
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
660 */
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);
668 }
669 }
670
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;
679 }
680 if(iy==0x3f800000) { /* y is +-1 */
681 if(hy<0) return one/x; else return x;
682 }
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);
687 }
688
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) */
699 }
700 return z;
701 }
702
703 /* (x<0)**(non-int) is NaN */
704 if(((((uint32_t)hx>>31)-1)|yisint)==0) return (x-x)/(x-x);
705
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);
735
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);
772 }
773
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) */
777
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 */
789 }
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 */
794 }
795 /*
796 * compute 2**(p_h+p_l)
797 */
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;
808 }
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;
825}
826
827
828
829/* Square root function, original. */
830float rb_sqrt(float x)
831{
832 float z;
833 int32_t sign = (int)0x80000000;
834 int32_t ix,s,q,m,t,i;
835 uint32_t r;
836
837 GET_FLOAT_WORD(ix,x);
838
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 */
843 }
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 */
849 }
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;
855 }
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] */
861
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 */
866
867 while(r!=0) {
868 t = s+r;
869 if(t<=ix) {
870 s = t+r;
871 ix -= t;
872 q += r;
873 }
874 ix += ix;
875 r>>=1;
876 }
877
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);
887 }
888 }
889 ix = (q>>1)+0x3f000000;
890 ix += (m <<23);
891 SET_FLOAT_WORD(z,ix);
892 return z;
893}
894
895/* Absolute value,
896 taken from glibc-2.8 */
897float rb_fabs(float x)
898{
899 uint32_t ix;
900 GET_FLOAT_WORD(ix,x);
901 SET_FLOAT_WORD(x,ix&0x7fffffff);
902 return x;
903}
904
905/* Arc tangent,
906 taken from glibc-2.8. */
907
908static 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 */
913};
914
915static 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 */
920};
921
922static 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 */
934};
935
936
937float rb_atan(float x)
938{
939 float w,s1,s2,z;
940 int32_t ix,hx,id;
941
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 */
952 }
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);
961 }
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;
967 }
968 }}
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;
979 }
980}
981
982/* Arc tangent from two variables, original. */
983
984static const float
985pi_o_4 = 7.8539818525e-01, /* 0x3f490fdb */
986pi_o_2 = 1.5707963705e+00, /* 0x3fc90fdb */
987pi = 3.1415927410e+00, /* 0x40490fdb */
988pi_lo = -8.7422776573e-08; /* 0xb3bbbd2e */
989
990float rb_atan2(float x, float y)
991{
992 float z;
993 int32_t k,m,hx,hy,ix,iy;
994
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) */
1004
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 */
1012 }
1013 }
1014 /* when x = 0 */
1015 if(ix==0) return (hy<0)? -pi_o_2-tiny: pi_o_2+tiny;
1016
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)*/
1025 }
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) */
1032 }
1033 }
1034 }
1035 /* when y is INF */
1036 if(iy==0x7f800000) return (hy<0)? -pi_o_2-tiny: pi_o_2+tiny;
1037
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);
1049 }
1050 return z ; /* atan(-,+) */
1051 case 2: return pi-(z-pi_lo);/* atan(+,-) */
1052 default: /* case 3 */
1053 return (z-pi_lo)-pi;/* atan(-,-) */
1054 }
1055}
1056
1057
1058/* Sine hyperbolic,
1059 taken from glibc-2.8 */
1060
1061static const float
1062o_threshold = 8.8721679688e+01,/* 0x42b17180 */
1063invln2 = 1.4426950216e+00,/* 0x3fb8aa3b */
1064 /* scaled coefficients related to expm1 */
1065Q1 = -3.3333335072e-02, /* 0xbd088889 */
1066Q2 = 1.5873016091e-03, /* 0x3ad00d01 */
1067Q3 = -7.9365076090e-05, /* 0xb8a670cd */
1068Q4 = 4.0082177293e-06, /* 0x36867e54 */
1069Q5 = -2.0109921195e-07; /* 0xb457edbb */
1070
1071float rb_expm1(float x)
1072{
1073 float y,hi,lo,c=0,t,e,hxs,hfx,r1;
1074 int32_t k,xsb;
1075 uint32_t hx;
1076
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| */
1081
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 */
1090 }
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 */
1094 }
1095 }
1096
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;
1109 }
1110 x = hi - lo;
1111 c = (hi-x)-lo;
1112 }
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));
1116 }
1117 else k = 0;
1118
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);
1133 }
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;
1140 }
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 */
1155 }
1156 }
1157 return y;
1158}
1159
1160static const float shuge = 1.0e37;
1161
1162float rb_sinh(float x)
1163{
1164 float t,w,h;
1165 int32_t ix,jx;
1166
1167 GET_FLOAT_WORD(jx,x);
1168 ix = jx&0x7fffffff;
1169
1170 /* x is INF or NaN */
1171 if(ix>=0x7f800000) return x+x;
1172
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));
1182 }
1183
1184 /* |x| in [22, log(maxdouble)] return 0.5*exp(|x|) */
1185 if (ix < 0x42b17180) return h*rb_exp(rb_fabs(x));
1186
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;
1192 }
1193
1194 /* |x| > overflowthresold, sinh(x) overflow */
1195 return x*shuge;
1196}
1197
1198
1199/* Tangent,
1200 taken from glibc-2.8 */
1201
1202static const float
1203pio4 = 7.8539812565e-01, /* 0x3f490fda */
1204pio4lo= 3.7748947079e-08, /* 0x33222168 */
1205T[] = {
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 */
1219};
1220
1221float kernel_tan(float x, float y, int iy)
1222{
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;
1231 }
1232 }
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;
1238 }
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]))
1244 */
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)));
1254 }
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);
1270 }
1271}
1272
1273
1274
1275static const int init_jk[] = {4,7,9}; /* initial value for jk */
1276
1277static 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 */
1289};
1290
1291static const float
1292two8 = 2.5600000000e+02, /* 0x43800000 */
1293twon8 = 3.9062500000e-03; /* 0x3b800000 */
1294
1295int kernel_rem_pio2(float *x, float *y, int e0, int nx, int prec, const int32_t *ipio2)
1296{
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];
1299
1300 /* initialize jk*/
1301 jk = init_jk[prec];
1302 jp = jk;
1303
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);
1308
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];
1312
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;
1316 }
1317
1318 jz = jk;
1319recompute:
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;
1325 }
1326
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);
1337 }
1338 else if(q0==0) ih = iq[jz-1]>>8;
1339 else if(z>=(float)0.5) ih=2;
1340
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;
1348 }
1349 } else iq[i] = 0xff - j;
1350 }
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;
1357 }
1358 }
1359 if(ih==2) {
1360 z = one - z;
1361 if(carry!=0) z -= rb_scalbnf(one,q0);
1362 }
1363 }
1364
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 */
1371
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;
1376 }
1377 jz += k;
1378 goto recompute;
1379 }
1380 }
1381
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 ;
1394 }
1395
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;
1400 }
1401
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;
1406 }
1407
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;
1429 }
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;
1434 }
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;
1440 }
1441 }
1442 return n&7;
1443}
1444
1445
1446static const int32_t two_over_pi[] = {
14470xA2, 0xF9, 0x83, 0x6E, 0x4E, 0x44, 0x15, 0x29, 0xFC,
14480x27, 0x57, 0xD1, 0xF5, 0x34, 0xDD, 0xC0, 0xDB, 0x62,
14490x95, 0x99, 0x3C, 0x43, 0x90, 0x41, 0xFE, 0x51, 0x63,
14500xAB, 0xDE, 0xBB, 0xC5, 0x61, 0xB7, 0x24, 0x6E, 0x3A,
14510x42, 0x4D, 0xD2, 0xE0, 0x06, 0x49, 0x2E, 0xEA, 0x09,
14520xD1, 0x92, 0x1C, 0xFE, 0x1D, 0xEB, 0x1C, 0xB1, 0x29,
14530xA7, 0x3E, 0xE8, 0x82, 0x35, 0xF5, 0x2E, 0xBB, 0x44,
14540x84, 0xE9, 0x9C, 0x70, 0x26, 0xB4, 0x5F, 0x7E, 0x41,
14550x39, 0x91, 0xD6, 0x39, 0x83, 0x53, 0x39, 0xF4, 0x9C,
14560x84, 0x5F, 0x8B, 0xBD, 0xF9, 0x28, 0x3B, 0x1F, 0xF8,
14570x97, 0xFF, 0xDE, 0x05, 0x98, 0x0F, 0xEF, 0x2F, 0x11,
14580x8B, 0x5A, 0x0A, 0x6D, 0x1F, 0x6D, 0x36, 0x7E, 0xCF,
14590x27, 0xCB, 0x09, 0xB7, 0x4F, 0x46, 0x3F, 0x66, 0x9E,
14600x5F, 0xEA, 0x2D, 0x75, 0x27, 0xBA, 0xC7, 0xEB, 0xE5,
14610xF1, 0x7B, 0x3D, 0x07, 0x39, 0xF7, 0x8A, 0x52, 0x92,
14620xEA, 0x6B, 0xFB, 0x5F, 0xB1, 0x1F, 0x8D, 0x5D, 0x08,
14630x56, 0x03, 0x30, 0x46, 0xFC, 0x7B, 0x6B, 0xAB, 0xF0,
14640xCF, 0xBC, 0x20, 0x9A, 0xF4, 0x36, 0x1D, 0xA9, 0xE3,
14650x91, 0x61, 0x5E, 0xE6, 0x1B, 0x08, 0x65, 0x99, 0x85,
14660x5F, 0x14, 0xA0, 0x68, 0x40, 0x8D, 0xFF, 0xD8, 0x80,
14670x4D, 0x73, 0x27, 0x31, 0x06, 0x06, 0x15, 0x56, 0xCA,
14680x73, 0xA8, 0xC9, 0x60, 0xE2, 0x7B, 0xC0, 0x8C, 0x6B,
1469};
1470
1471static const int32_t npio2_hw[] = {
14720x3fc90f00, 0x40490f00, 0x4096cb00, 0x40c90f00, 0x40fb5300, 0x4116cb00,
14730x412fed00, 0x41490f00, 0x41623100, 0x417b5300, 0x418a3a00, 0x4196cb00,
14740x41a35c00, 0x41afed00, 0x41bc7e00, 0x41c90f00, 0x41d5a000, 0x41e23100,
14750x41eec200, 0x41fb5300, 0x4203f200, 0x420a3a00, 0x42108300, 0x4216cb00,
14760x421d1400, 0x42235c00, 0x4229a500, 0x422fed00, 0x42363600, 0x423c7e00,
14770x4242c700, 0x42490f00
1478};
1479
1480/*
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)
1488 */
1489
1490static const float
1491half = 5.0000000000e-01, /* 0x3f000000 */
1492invpio2 = 6.3661980629e-01, /* 0x3f22f984 */
1493pio2_1 = 1.5707855225e+00, /* 0x3fc90f80 */
1494pio2_1t = 1.0804334124e-05, /* 0x37354443 */
1495pio2_2 = 1.0804273188e-05, /* 0x37354400 */
1496pio2_2t = 6.0770999344e-11, /* 0x2e85a308 */
1497pio2_3 = 6.0770943833e-11, /* 0x2e85a300 */
1498pio2_3t = 6.1232342629e-17; /* 0x248d3132 */
1499
1500int32_t rem_pio2(float x, float *y)
1501{
1502 float z,w,t,r,fn;
1503 float tx[3];
1504 int32_t e0,i,j,nx,n,ix,hx;
1505
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;
1520 }
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;
1531 }
1532 return -1;
1533 }
1534 }
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;
1563 }
1564 }
1565 }
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;
1569 }
1570 /*
1571 * all other (large) arguments
1572 */
1573 if(ix>=0x7f800000) { /* x is inf or NaN */
1574 y[0]=y[1]=x-x; return 0;
1575 }
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;
1582 }
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;
1589}
1590
1591float rb_tan(float x)
1592{
1593 float y[2],z=0.0;
1594 int32_t n, ix;
1595
1596 GET_FLOAT_WORD(ix,x);
1597
1598 /* |x| ~< pi/4 */
1599 ix &= 0x7fffffff;
1600 if(ix <= 0x3f490fda) return kernel_tan(x,z,1);
1601
1602 /* tan(Inf or NaN) is NaN */
1603 else if (ix>=0x7f800000) return x-x; /* NaN */
1604
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 */
1610 }
1611}
1612
1613
1614
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 */
1619
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)
1625const 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
1671};
1672
1673const 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 */
2029};
2030
2031/* All floating-point numbers can be put in one of these categories. */
2032enum
2033 {
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
2044 };
2045
2046
2047int
2048__fpclassifyf (float x)
2049{
2050 uint32_t wx;
2051 int retval = FP_NORMAL;
2052
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;
2061
2062 return retval;
2063}
2064
2065
2066int
2067__isinff (float x)
2068{
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);
2075}
2076
2077/* Return nonzero value if arguments are unordered. */
2078# define fpclassify(x) \
2079 (sizeof (x) == sizeof (float) ? __fpclassifyf (x) : __fpclassifyf (x))
2080
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
2087
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
2095
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
2103
2104union ieee754_double
2105 {
2106 double d;
2107
2108 /* This is the IEEE 754 double-precision format. */
2109 struct
2110 {
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;
2132
2133 /* This format makes it easier to see if a NaN is a signalling NaN. */
2134 struct
2135 {
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;
2160 };
2161
2162
2163static const volatile float TWOM100 = 7.88860905e-31;
2164static const volatile float TWO127 = 1.7014118346e+38;
2165
2166float rb_exp(float x)
2167{
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))
2172 {
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;
2181
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;
2188
2189 feholdexcept (&oldenv);
2190#endif
2191
2192#ifdef FE_TONEAREST
2193 fesetround (FE_TONEAREST);
2194#endif
2195
2196 /* Calculate n. */
2197 n = x * M_1_LN2 + THREEp22;
2198 n -= THREEp22;
2199 dx = x - n*M_LN2;
2200
2201 /* Calculate t/512. */
2202 t = dx + THREEp42;
2203 t -= THREEp42;
2204 dx -= t;
2205
2206 /* Compute tval = t. */
2207 tval = (int) (t * 512.0);
2208
2209 if (t >= 0)
2210 delta = - __exp_deltatable[tval];
2211 else
2212 delta = __exp_deltatable[-tval];
2213
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;
2217
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;
2222
2223 /* Return result. */
2224#ifndef ROCKBOX
2225 fesetenv (&oldenv);
2226#endif
2227
2228 result = x22 * ex2_u.d + ex2_u.d;
2229 return (float) result;
2230 }
2231 /* Exceptional cases: */
2232 else if (isless (x, himark))
2233 {
2234 if (__isinff (x))
2235 /* e^-inf == 0, with no error. */
2236 return 0;
2237 else
2238 /* Underflow */
2239 return TWOM100 * TWOM100;
2240 }
2241 else
2242 /* Return x, if x is a NaN or Inf; or overflow, otherwise. */
2243 return TWO127*x;
2244}
2245
2246
2247/* Division with rest, original. */
2248div_t div(int x, int y)
2249{
2250 div_t result;
2251
2252 result.quot = x / y;
2253 result.rem = x % y;
2254
2255 /* Addition from glibc-2.8: */
2256 if(x >= 0 && result.rem < 0)
2257 {
2258 result.quot += 1;
2259 result.rem -= y;
2260 }
2261
2262 return result;
2263}
2264
2265
2266/* Placeholder function. */
2267/* Originally defined in s_midi.c */
2268void sys_listmididevs(void)
2269{
2270}
2271
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 */
2276void sys_microsleep(int microsec)
2277{
2278 (void) microsec;
2279}
2280
2281/* Get running time in milliseconds. */
2282/* Originally defined in s_inter.c */
2283extern uint64_t runningtime;
2284t_time sys_getrealtime(void)
2285{
2286 return runningtime;
2287}
2288
2289/* Place holder, as we do no IPC. */
2290/* Originally defined in s_inter.c */
2291void glob_ping(void* dummy)
2292{
2293 (void) dummy;
2294}
2295
2296/* Call to quit. */
2297/* Originally defined in s_inter.c */
2298extern bool quit;
2299void glob_quit(void* dummy)
2300{
2301 (void) dummy;
2302
2303 static bool reentered = false;
2304
2305 if(!reentered)
2306 {
2307 reentered = true;
2308
2309 /* Close audio subsystem. */
2310 sys_close_audio();
2311
2312 /* Stop main loop. */
2313 quit = true;
2314 }
2315}
2316
2317/* Open file. Originally in s_main.c */
2318void glob_evalfile(t_pd *ignore, t_symbol *name, t_symbol *dir);
2319void openit(const char *dirname, const char *filename)
2320{
2321 char dirbuf[MAXPDSTRING], *nameptr;
2322
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);
2330
2331 int fd = open_via_path(dirname, ffilename, "", dirbuf, &nameptr,
2332 MAXPDSTRING, 0);
2333 if (fd)
2334 {
2335 close (fd);
2336 glob_evalfile(0, gensym(nameptr), gensym(dirbuf));
2337 }
2338 else
2339 error("%s: can't open", filename);
2340}
2341
2342
2343/* Get current working directory. */
2344extern char* filename;
2345char* rb_getcwd(char* buf, ssize_t size)
2346{
2347 char* end_of_dir = strrchr(filename, '/');
2348
2349 /* Check whether buffer may hold enough data. */
2350 if(size < end_of_dir - filename)
2351 return NULL;
2352
2353 /* Copy current working directory to buffer. */
2354 strncpy(buf, filename, end_of_dir - filename);
2355 return buf;
2356}
2357
2358
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 */
2363extern t_namelist* sys_openlist;
2364void glob_initfromgui(void *dummy, t_symbol *s, int argc, t_atom *argv)
2365{
2366 (void) dummy;
2367 (void) s;
2368 (void) argc;
2369 (void) argv;
2370
2371 t_namelist *nl;
2372 char cwd[MAXPDSTRING];
2373
2374 /* Get current working directory. */
2375 rb_getcwd(cwd, MAXPDSTRING);
2376
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;
2382}
2383
2384/* Fake GUI start. Originally in s_inter.c */
2385static 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};
2388extern t_binbuf* inbinbuf;
2389int sys_startgui(const char *guidir)
2390{
2391 unsigned int i;
2392 t_atom zz[23];
2393 char cmdbuf[4*MAXPDSTRING];
2394
2395 (void) guidir;
2396
2397 inbinbuf = binbuf_new();
2398
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);
2405
2406 glob_initfromgui(0, 0, 23, zz);
2407
2408 return 0;
2409}
2410
2411/* Return default DAC block size. */
2412/* Originally defined in s_main.c */
2413int sys_getblksize(void)
2414{
2415 return (DEFDACBLKSIZE);
2416}
2417
2418/* Find library directory and set it. */
2419void sys_findlibdir(const char* filename)
2420{
2421 char sbuf[MAXPDSTRING];
2422
2423 (void) filename;
2424
2425 /* Make current working directory the system library directory. */
2426 rb_getcwd(sbuf, MAXPDSTRING);
2427 sys_libdir = gensym(sbuf);
2428}
generated by cgit v1.2.3 (git 2.39.1) at 2024-09-12 18:26:27 +0000