diff options
Diffstat (limited to 'apps/plugins/pdbox/pdbox-func.c')
-rw-r--r-- | apps/plugins/pdbox/pdbox-func.c | 2428 |
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. */ | ||
34 | char *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 | |||
72 | strncat_fini: | ||
73 | return dest; | ||
74 | } | ||
75 | |||
76 | |||
77 | /* Implementation of floor, original. */ | ||
78 | float 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/). */ | ||
94 | static int rb_errno = 0; | ||
95 | |||
96 | double 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 | |||
223 | double rb_atof(const char *str) | ||
224 | { | ||
225 | return rb_strtod(str, NULL); | ||
226 | } | ||
227 | |||
228 | |||
229 | /* Implementation of ftoa(), original. */ | ||
230 | void 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. */ | ||
295 | long 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 | |||
334 | float 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 | |||
348 | float 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. */ | ||
361 | int 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); */ | ||
386 | int 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 */ | ||
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 */ | ||
416 | |||
417 | static const float zero = 0.0; | ||
418 | |||
419 | /* A union which permits us to convert between a float and a 32 bit | ||
420 | int. */ | ||
421 | |||
422 | typedef 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) \ | ||
431 | do { \ | ||
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) \ | ||
440 | do { \ | ||
441 | ieee_float_shape_type sf_u; \ | ||
442 | sf_u.word = (i); \ | ||
443 | (d) = sf_u.value; \ | ||
444 | } while (0) | ||
445 | |||
446 | |||
447 | float 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 | |||
516 | static const float | ||
517 | ivln10 = 4.3429449201e-01, /* 0x3ede5bd9 */ | ||
518 | log10_2hi = 3.0102920532e-01, /* 0x3e9a2080 */ | ||
519 | log10_2lo = 7.9034151668e-07; /* 0x355427db */ | ||
520 | |||
521 | float 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 | |||
550 | int 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 | |||
560 | float 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 | |||
569 | static const float | ||
570 | huge = 1.0e+30, | ||
571 | tiny = 1.0e-30, | ||
572 | twom25 = 2.9802322388e-08; /* 0x33000000 */ | ||
573 | |||
574 | float 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 | |||
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 */ | ||
620 | |||
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*/ | ||
632 | |||
633 | float 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. */ | ||
830 | float 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 */ | ||
897 | float 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 | |||
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 */ | ||
913 | }; | ||
914 | |||
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 */ | ||
920 | }; | ||
921 | |||
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 */ | ||
934 | }; | ||
935 | |||
936 | |||
937 | float 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 | |||
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 */ | ||
989 | |||
990 | float 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 | |||
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 */ | ||
1070 | |||
1071 | float 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 | |||
1160 | static const float shuge = 1.0e37; | ||
1161 | |||
1162 | float 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 | |||
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 */ | ||
1219 | }; | ||
1220 | |||
1221 | float 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 | |||
1275 | static const int init_jk[] = {4,7,9}; /* initial value for jk */ | ||
1276 | |||
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 */ | ||
1289 | }; | ||
1290 | |||
1291 | static const float | ||
1292 | two8 = 2.5600000000e+02, /* 0x43800000 */ | ||
1293 | twon8 = 3.9062500000e-03; /* 0x3b800000 */ | ||
1294 | |||
1295 | int 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; | ||
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; | ||
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 | |||
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, | ||
1469 | }; | ||
1470 | |||
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 | ||
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 | |||
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 */ | ||
1499 | |||
1500 | int32_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 | |||
1591 | float 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) | ||
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 | ||
1671 | }; | ||
1672 | |||
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 */ | ||
2029 | }; | ||
2030 | |||
2031 | /* All floating-point numbers can be put in one of these categories. */ | ||
2032 | enum | ||
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 | |||
2047 | int | ||
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 | |||
2066 | int | ||
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 | |||
2104 | union 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 | |||
2163 | static const volatile float TWOM100 = 7.88860905e-31; | ||
2164 | static const volatile float TWO127 = 1.7014118346e+38; | ||
2165 | |||
2166 | float 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. */ | ||
2248 | div_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 */ | ||
2268 | void 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 */ | ||
2276 | void sys_microsleep(int microsec) | ||
2277 | { | ||
2278 | (void) microsec; | ||
2279 | } | ||
2280 | |||
2281 | /* Get running time in milliseconds. */ | ||
2282 | /* Originally defined in s_inter.c */ | ||
2283 | extern uint64_t runningtime; | ||
2284 | t_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 */ | ||
2291 | void glob_ping(void* dummy) | ||
2292 | { | ||
2293 | (void) dummy; | ||
2294 | } | ||
2295 | |||
2296 | /* Call to quit. */ | ||
2297 | /* Originally defined in s_inter.c */ | ||
2298 | extern bool quit; | ||
2299 | void 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 */ | ||
2318 | void glob_evalfile(t_pd *ignore, t_symbol *name, t_symbol *dir); | ||
2319 | void 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. */ | ||
2344 | extern char* filename; | ||
2345 | char* 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 */ | ||
2363 | extern t_namelist* sys_openlist; | ||
2364 | void 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 */ | ||
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) | ||
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 */ | ||
2413 | int sys_getblksize(void) | ||
2414 | { | ||
2415 | return (DEFDACBLKSIZE); | ||
2416 | } | ||
2417 | |||
2418 | /* Find library directory and set it. */ | ||
2419 | void 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 | } | ||