1 /*
2 * This file is part of DGD, https://github.com/dworkin/dgd
3 * Copyright (C) 1993-2010 Dworkin B.V.
4 * Copyright (C) 2010-2012 DGD Authors (see the commit log for details)
5 *
6 * This program is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU Affero General Public License as
8 * published by the Free Software Foundation, either version 3 of the
9 * License, or (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU Affero General Public License for more details.
15 *
16 * You should have received a copy of the GNU Affero General Public License
17 * along with this program. If not, see <http://www.gnu.org/licenses/>.
18 */
19
20 # define INCLUDE_CTYPE
21 # include "dgd.h"
22 # include "xfloat.h"
23
24 typedef struct {
25 unsigned short sign; /* 0: positive, 0x8000: negative */
26 unsigned short exp; /* bias: 32767 */
27 unsigned short high; /* 0 / 1 / 14 bits */
28 Uint low; /* 0 / 29 bits / 0 / 0 */
29 } flt;
30
31 # define NBITS 44
32 # define BIAS 0x7fff
33
34
35 /* constants */
36
37 xfloat max_int = { 0x41df, 0xffffffc0L }; /* 0x7fffffff */
38 xfloat thousand = { 0x408f, 0x40000000L }; /* 1e3 */
39 xfloat thousandth = { 0x3f50, 0x624dd2f2L }; /* 1e-3 */
40
41 static flt half = { 0x0000, 0x7ffe, 0x4000, 0x00000000L };
42 static flt one = { 0x0000, 0x7fff, 0x4000, 0x00000000L };
43 static flt maxlog = { 0x0000, 0x8008, 0x58b9, 0x05fdf474L };
44 static flt minlog = { 0x8000, 0x8008, 0x58b9, 0x05fdf474L };
45 static flt sqrth = { 0x0000, 0x7ffe, 0x5a82, 0x3cccfe78L };
46 static flt pi = { 0x0000, 0x8000, 0x6487, 0x76a8885cL };
47 static flt pio2 = { 0x0000, 0x7fff, 0x6487, 0x76a8885cL };
48 static flt pio4 = { 0x0000, 0x7ffe, 0x6487, 0x76a8885cL };
49
50
51 /*
52 * NAME: f_edom()
53 * DESCRIPTION: Domain error
54 */
f_edom()55 static void f_edom()
56 {
57 error("Math argument");
58 }
59
60 /*
61 * NAME: f_erange()
62 * DESCRIPTION: Out of range
63 */
f_erange()64 static void f_erange()
65 {
66 error("Result too large");
67 }
68
69 static void f_sub (flt*, flt*);
70
71 /*
72 * NAME: f_add()
73 * DESCRIPTION: a = a + b. The result is normalized, but not guaranteed to
74 * be in range.
75 */
f_add(flt * a,flt * b)76 static void f_add(flt *a, flt *b)
77 {
78 unsigned short h, n;
79 Uint l;
80 flt tmp;
81
82 if (b->exp == 0) {
83 /* b is 0 */
84 return;
85 }
86 if (a->exp == 0) {
87 /* a is 0 */
88 *a = *b;
89 return;
90 }
91 if (a->sign != b->sign) {
92 a->sign ^= 0x8000;
93 f_sub(a, b);
94 a->sign ^= 0x8000;
95 return;
96 }
97 if (a->exp < b->exp) {
98 /* b is the largest; exchange a and b */
99 tmp = *a;
100 *a = *b;
101 b = &tmp;
102 }
103
104 n = a->exp - b->exp;
105 if (n <= NBITS) {
106 h = a->high;
107 l = a->low;
108
109 /*
110 * perform addition
111 */
112 if (n < 31) {
113 h += (Uint) b->high >> n;
114 l += (((Uint) b->high << (31 - n)) & 0x7fffffffL) | (b->low >> n);
115 } else {
116 l += b->high >> (n - 31);
117 }
118 if ((Int) l < 0) {
119 /* carry */
120 l &= 0x7fffffffL;
121 h++;
122 }
123 if ((short) h < 0) {
124 /* too large */
125 l |= (Uint) h << 31;
126 l >>= 1;
127 h >>= 1;
128 a->exp++;
129 }
130
131 /*
132 * rounding off
133 */
134 if ((Int) (l += 2) < 0 && (short) ++h < 0) {
135 h >>= 1;
136 a->exp++;
137 }
138 l &= 0x7ffffffcL;
139
140 a->high = h;
141 a->low = l;
142 }
143 }
144
145 /*
146 * NAME: f_sub()
147 * DESCRIPTION: a = a - b. The result is normalized, but not guaranteed to be
148 * in range.
149 */
f_sub(flt * a,flt * b)150 static void f_sub(flt *a, flt *b)
151 {
152 unsigned short h, n;
153 Uint l;
154 flt tmp;
155
156 if (b->exp == 0) {
157 /* b is 0 */
158 return;
159 }
160 if (a->exp == 0) {
161 *a = *b;
162 a->sign ^= 0x8000;
163 return;
164 }
165 if (a->sign != b->sign) {
166 a->sign ^= 0x8000;
167 f_add(a, b);
168 a->sign ^= 0x8000;
169 return;
170 }
171 if (a->exp <= b->exp &&
172 (a->exp < b->exp || (a->high <= b->high &&
173 (a->high < b->high || a->low < b->low)))) {
174 /* b is the largest; exchange a and b */
175 tmp = *a;
176 *a = *b;
177 b = &tmp;
178 a->sign ^= 0x8000;
179 }
180
181 n = a->exp - b->exp;
182 if (n <= NBITS) {
183 h = a->high;
184 l = a->low;
185
186 /*
187 * perform subtraction
188 */
189 if (n < 31) {
190 h -= (Uint) b->high >> n;
191 l -= (((Uint) b->high << (31 - n)) & 0x7fffffffL) | (b->low >> n);
192 if (b->low & ((1 << n) - 1)) {
193 --l;
194 }
195 } else {
196 n -= 31;
197 l -= b->high >> n;
198 if (b->low != 0 || (b->high & ((1 << n) - 1))) {
199 --l;
200 }
201 }
202 if ((Int) l < 0) {
203 /* borrow */
204 l &= 0x7fffffffL;
205 --h;
206 }
207
208 /*
209 * normalize
210 */
211 if (h == 0) {
212 if (l == 0) {
213 a->exp = 0;
214 return;
215 }
216 n = 15;
217 if ((l & 0xffff0000L) == 0) {
218 l <<= 15;
219 n += 15;
220 }
221 h = l >> 16;
222 l <<= 15;
223 l &= 0x7fffffffL;
224 a->exp -= n;
225 }
226 if (h < 0x4000) {
227 n = 0;
228 if ((h & 0xff00) == 0) {
229 h <<= 7;
230 n += 7;
231 }
232 while (h < 0x4000) {
233 h <<= 1;
234 n++;
235 }
236 h |= l >> (31 - n);
237 l <<= n;
238 l &= 0x7fffffffL;
239 a->exp -= n;
240 }
241
242 /*
243 * rounding off
244 */
245 if ((Int) (l += 2) < 0 && (short) ++h < 0) {
246 h >>= 1;
247 a->exp++;
248 }
249 l &= 0x7ffffffcL;
250
251 a->high = h;
252 a->low = l;
253 }
254 }
255
256 /*
257 * NAME: f_mult()
258 * DESCRIPTION: a = a * b. The result is normalized, but may be out of range.
259 */
f_mult(flt * a,flt * b)260 static void f_mult(flt *a, flt *b)
261 {
262 Uint m, l, albl, ambm, ahbh;
263 short al, am, ah, bl, bm, bh;
264
265 if (a->exp == 0) {
266 /* a is 0 */
267 return;
268 }
269 if (b->exp == 0) {
270 /* b is 0 */
271 a->exp = 0;
272 return;
273 }
274
275 al = ((unsigned short) a->low) >> 1;
276 bl = ((unsigned short) b->low) >> 1;
277 am = a->low >> 16;
278 bm = b->low >> 16;
279 ah = a->high;
280 bh = b->high;
281
282 albl = (Uint) al * bl;
283 ambm = (Uint) am * bm;
284 ahbh = (Uint) ah * bh;
285 m = albl;
286 m >>= 15;
287 m += albl + ambm + (Int) (al - am) * (bm - bl);
288 m >>= 15;
289 m += albl + ambm + ahbh + (Int) (al - ah) * (bh - bl);
290 m >>= 13;
291 l = m & 0x03;
292 m >>= 2;
293 m += ambm + ahbh + (Int) (am - ah) * (bh - bm);
294 l |= (m & 0x7fff) << 2;
295 m >>= 15;
296 m += ahbh;
297 l |= m << 17;
298 ah = m >> 14;
299
300 a->sign ^= b->sign;
301 a->exp += b->exp - BIAS;
302 if (ah < 0) {
303 ah = (unsigned short) ah >> 1;
304 l >>= 1;
305 a->exp++;
306 }
307 l &= 0x7fffffffL;
308
309 /*
310 * rounding off
311 */
312 if ((Int) (l += 2) < 0 && ++ah < 0) {
313 ah = (unsigned short) ah >> 1;
314 a->exp++;
315 }
316 l &= 0x7ffffffcL;
317
318 a->high = ah;
319 a->low = l;
320 }
321
322 /*
323 * NAME: f_div()
324 * DESCRIPTION: a = a / b. b must be non-zero. The result is normalized,
325 * but may be out of range.
326 */
f_div(flt * a,flt * b)327 static void f_div(flt *a, flt *b)
328 {
329 unsigned short n[3];
330 Uint numh, numl, divl, high, low, q;
331 unsigned short divh, i;
332
333 if (b->exp == 0) {
334 error("Division by zero");
335 }
336 if (a->exp == 0) {
337 /* a is 0 */
338 return;
339 }
340
341 numh = ((Uint) a->high << 16) | (a->low >> 15);
342 numl = a->low << 17;
343
344 divh = (b->high << 1) | (b->low >> 30);
345 divl = b->low << 2;
346
347 n[0] = 0;
348 n[1] = 0;
349 i = 3;
350 do {
351 /* estimate the high word of the quotient */
352 q = numh / divh;
353 /* highlow = div * q */
354 low = (unsigned short) (high = q * (unsigned short) divl);
355 high >>= 16;
356 high += q * (divl >> 16);
357 low |= high << 16;
358 high >>= 16;
359 high += q * divh;
360
361 /* the estimated quotient may be 2 off; correct it if needed */
362 if (high >= numh && (high > numh || low > numl)) {
363 high -= divh;
364 if (low < divl) {
365 --high;
366 }
367 low -= divl;
368 --q;
369 if (high >= numh && (high > numh || low > numl)) {
370 high -= divh;
371 if (low < divl) {
372 --high;
373 }
374 low -= divl;
375 --q;
376 }
377 }
378
379 n[--i] = q;
380 if (i == 0) {
381 break;
382 }
383
384 /* subtract highlow */
385 numh -= high;
386 if (numl < low) {
387 --numh;
388 }
389 numl -= low;
390 numh <<= 16;
391 numh |= numl >> 16;
392 numl <<= 16;
393
394 } while (numh != 0 || numl != 0);
395
396 a->sign ^= b->sign;
397 a->exp -= b->exp - BIAS + 1;
398 high = n[2];
399 low = ((Uint) n[1] << 15) | (n[0] >> 1);
400 if ((short) high < 0) {
401 low |= high << 31;
402 low >>= 1;
403 high >>= 1;
404 a->exp++;
405 }
406
407 /*
408 * rounding off
409 */
410 if ((Int) (low += 2) < 0 && (short) ++high < 0) {
411 high >>= 1;
412 a->exp++;
413 }
414 low &= 0x7ffffffcL;
415
416 a->high = high;
417 a->low = low;
418 }
419
420 /*
421 * NAME: f_trunc()
422 * DESCRIPTION: truncate a flt
423 */
f_trunc(flt * a)424 static void f_trunc(flt *a)
425 {
426 static unsigned short maskh[] = {
427 0x4000, 0x6000, 0x7000, 0x7800, 0x7c00, 0x7e00, 0x7f00,
428 0x7f80, 0x7fc0, 0x7fe0, 0x7ff0, 0x7ff8, 0x7ffc, 0x7ffe, 0x7fff,
429 0x7fff, 0x7fff, 0x7fff, 0x7fff, 0x7fff, 0x7fff, 0x7fff,
430 0x7fff, 0x7fff, 0x7fff, 0x7fff, 0x7fff, 0x7fff, 0x7fff, 0x7fff,
431 0x7fff, 0x7fff, 0x7fff, 0x7fff, 0x7fff, 0x7fff, 0x7fff, 0x7fff,
432 0x7fff, 0x7fff, 0x7fff, 0x7fff, 0x7fff
433 };
434 static Uint maskl[] = {
435 0x00000000L, 0x00000000L, 0x00000000L,
436 0x00000000L, 0x00000000L, 0x00000000L, 0x00000000L,
437 0x00000000L, 0x00000000L, 0x00000000L, 0x00000000L,
438 0x00000000L, 0x00000000L, 0x00000000L, 0x00000000L,
439 0x40000000L, 0x60000000L, 0x70000000L,
440 0x78000000L, 0x7c000000L, 0x7e000000L, 0x7f000000L,
441 0x7f800000L, 0x7fc00000L, 0x7fe00000L, 0x7ff00000L,
442 0x7ff80000L, 0x7ffc0000L, 0x7ffe0000L, 0x7fff0000L,
443 0x7fff8000L, 0x7fffc000L, 0x7fffe000L, 0x7ffff000L,
444 0x7ffff800L, 0x7ffffc00L, 0x7ffffe00L, 0x7fffff00L,
445 0x7fffff80L, 0x7fffffc0L, 0x7fffffe0L, 0x7ffffff0L,
446 0x7ffffff8L
447 };
448
449 if (a->exp < BIAS) {
450 a->exp = 0;
451 } else if (a->exp < BIAS + NBITS - 1) {
452 a->high &= maskh[a->exp - BIAS];
453 a->low &= maskl[a->exp - BIAS];
454 }
455 }
456
457 /*
458 * NAME: f_37bits()
459 * DESCRIPTION: round a flt to 37 binary digits of precision
460 */
f_37bits(flt * a)461 static void f_37bits(flt *a)
462 {
463 if ((Int) (a->low += 0x100) < 0) {
464 a->low = 0;
465 if ((short) ++(a->high) < 0) {
466 a->high >>= 1;
467 a->exp++;
468 }
469 } else {
470 a->low &= 0xfffffe00L;
471 }
472 }
473
474 /*
475 * NAME: f_round()
476 * DESCRIPTION: round off a flt
477 */
f_round(flt * a)478 static void f_round(flt *a)
479 {
480 static flt half = { 0, 0x7ffe, 0x4000, 0x00000000L };
481
482 half.sign = a->sign;
483 f_add(a, &half);
484 f_trunc(a);
485 }
486
487 /*
488 * NAME: f_cmp()
489 * DESCRIPTION: compate two flts
490 */
f_cmp(flt * a,flt * b)491 static int f_cmp(flt *a, flt *b)
492 {
493 if (a->exp == 0) {
494 if (b->exp == 0) {
495 return 0;
496 }
497 return (b->sign == 0) ? -1 : 1;
498 } else if (b->exp == 0) {
499 if (a->exp == 0) {
500 return 0;
501 }
502 return (a->sign == 0) ? 1 : -1;
503 } else if (a->sign != b->sign) {
504 return (a->sign == 0) ? 1 : -1;
505 } else {
506 if (a->exp == b->exp && a->high == b->high && a->low == b->low) {
507 return 0;
508 }
509 if (a->exp <= b->exp &&
510 (a->exp < b->exp || (a->high <= b->high &&
511 (a->high < b->high || a->low < b->low)))) {
512 return (a->sign == 0) ? -1 : 1;
513 }
514 return (a->sign == 0) ? 1 : -1;
515 }
516 }
517
518 /*
519 * NAME: f_itof()
520 * DESCRIPTION: convert an integer to a flt
521 */
f_itof(Int i,flt * a)522 static void f_itof(Int i, flt *a)
523 {
524 Uint n;
525 unsigned short shift;
526
527 /* deal with zero and sign */
528 if (i == 0) {
529 a->exp = 0;
530 return;
531 } else if (i < 0) {
532 a->sign = 0x8000;
533 n = -i;
534 } else {
535 a->sign = 0;
536 n = i;
537 }
538
539 shift = 0;
540 while ((n & 0xff000000L) == 0) {
541 n <<= 8;
542 shift += 8;
543 }
544 while ((Int) n >= 0) {
545 n <<= 1;
546 shift++;
547 }
548 a->exp = BIAS + 31 - shift;
549 a->high = n >> 17;
550 a->low = (n << 14) & 0x7fffffff;
551 }
552
553 /*
554 * NAME: f_ftoi()
555 * DESCRIPTION: convert a flt to an integer, discarding the fractional part
556 */
f_ftoi(flt * a)557 static Int f_ftoi(flt *a)
558 {
559 Uint i;
560
561 if (a->exp < BIAS) {
562 return 0;
563 }
564 if (a->exp > BIAS + 30 &&
565 (a->sign == 0 || a->exp != BIAS + 31 || a->high != 0x4000 ||
566 a->low != 0)) {
567 f_erange();
568 }
569
570 i = (((Uint) a->high << 17) | (a->low >> 14)) >> (BIAS + 31 - a->exp);
571
572 return (a->sign == 0) ? i : -i;
573 }
574
575 /*
576 * NAME: f_ftoxf()
577 * DESCRIPTION: convert flt to xfloat
578 */
f_ftoxf(flt * a,xfloat * f)579 static void f_ftoxf(flt *a, xfloat *f)
580 {
581 unsigned short exp;
582 unsigned short high;
583 Uint low;
584
585 exp = a->exp;
586 if (exp == 0) {
587 /* zero */
588 f->high = 0;
589 f->low = 0;
590 return;
591 }
592 high = a->high;
593 low = a->low;
594
595 /* mantissa */
596 if ((Int) (low += 0x100) < 0) {
597 low = 0;
598 if ((short) ++high < 0) {
599 high >>= 1;
600 exp++;
601 }
602 }
603
604 /* exponent */
605 if (exp > BIAS + 1023) {
606 f_erange();
607 }
608 if (exp < BIAS - 1022) {
609 /* underflow */
610 f->high = 0;
611 f->low = 0;
612 return;
613 }
614
615 f->high = a->sign | ((exp - BIAS + 1023) << 4) | ((high >> 10) & 0x000f);
616 f->low = ((Uint) high << 22) | (low >> 9);
617 }
618
619 /*
620 * NAME: f_xftof()
621 * DESCRIPTION: convert xfloat to flt
622 */
f_xftof(xfloat * f,flt * a)623 static void f_xftof(xfloat *f, flt *a)
624 {
625 unsigned short exp;
626
627 a->sign = f->high & 0x8000;
628 exp = (f->high >> 4) & 0x07ff;
629 if (exp == 0) {
630 /* zero */
631 a->exp = 0;
632 return;
633 }
634 a->exp = exp + BIAS - 1023;
635 a->high = 0x4000 | ((f->high & 0x0f) << 10) | (f->low >> 22);
636 a->low = (f->low << 9) & 0x7fffffffL;
637 }
638
639
640 static flt tens[] = {
641 { 0, 0x8002, 0x5000, 0x00000000L }, /* 10 ** 1 */
642 { 0, 0x8005, 0x6400, 0x00000000L }, /* 10 ** 2 */
643 { 0, 0x800C, 0x4E20, 0x00000000L }, /* 10 ** 4 */
644 { 0, 0x8019, 0x5F5E, 0x08000000L }, /* 10 ** 8 */
645 { 0, 0x8034, 0x470D, 0x726FC100L }, /* 10 ** 16 */
646 { 0, 0x8069, 0x4EE2, 0x6B6A0ADCL }, /* 10 ** 32 */
647 { 0, 0x80D3, 0x613C, 0x07D27FF4L }, /* 10 ** 64 */
648 { 0, 0x81A8, 0x49DD, 0x11F2603CL }, /* 10 ** 128 */
649 { 0, 0x8351, 0x553F, 0x3AFEE780L }, /* 10 ** 256 */
650 { 0, 0x86A3, 0x718C, 0x682BA984L }, /* 10 ** 512 */
651 };
652
653 static flt tenths[] = {
654 { 0, 0x7FFB, 0x6666, 0x33333334L }, /* 10 ** -1 */
655 { 0, 0x7FF8, 0x51EB, 0x428F5C28L }, /* 10 ** -2 */
656 { 0, 0x7FF1, 0x68DB, 0x45D63888L }, /* 10 ** -4 */
657 { 0, 0x7FE4, 0x55E6, 0x1DC46118L }, /* 10 ** -8 */
658 { 0, 0x7FC9, 0x734A, 0x652FB114L }, /* 10 ** -16 */
659 { 0, 0x7F94, 0x67D8, 0x47AB5150L }, /* 10 ** -32 */
660 { 0, 0x7F2A, 0x543F, 0x7A89E950L }, /* 10 ** -64 */
661 { 0, 0x7E55, 0x6EE8, 0x119F1930L }, /* 10 ** -128 */
662 { 0, 0x7CAC, 0x6018, 0x50C958E0L }, /* 10 ** -256 */
663 { 0, 0x795A, 0x4824, 0x7B8CB6C8L }, /* 10 ** -512 */
664 };
665
666 /*
667 * NAME: float->atof()
668 * DESCRIPTION: Convert a string to a float. The string must be in the
669 * proper format. Return TRUE if the operation was successful,
670 * FALSE otherwise.
671 */
flt_atof(char ** s,xfloat * f)672 bool flt_atof(char **s, xfloat *f)
673 {
674 flt a = { 0 };
675 flt b, c, *t;
676 unsigned short e, h;
677 char *p;
678
679 p = *s;
680
681 /* sign */
682 if (*p == '-') {
683 a.sign = b.sign = 0x8000;
684 p++;
685 } else {
686 a.sign = b.sign = 0;
687 }
688
689 a.exp = 0;
690 b.low = 0;
691
692 /* digits before . */
693 while (isdigit(*p)) {
694 f_mult(&a, &tens[0]);
695 h = (*p++ - '0') << 12;
696 if (h != 0) {
697 e = BIAS + 3;
698 while ((short) h >= 0) {
699 h <<= 1;
700 --e;
701 }
702 b.exp = e;
703 b.high = h >> 1;
704 f_add(&a, &b);
705 }
706 if (a.exp > 0xffff - 10) {
707 return FALSE;
708 }
709 }
710
711 /* digits after . */
712 if (*p == '.') {
713 c = tenths[0];
714 while (isdigit(*++p)) {
715 if (c.exp > 10) {
716 h = (*p - '0') << 12;
717 if (h != 0) {
718 e = BIAS + 3;
719 while ((short) h >= 0) {
720 h <<= 1;
721 --e;
722 }
723 b.exp = e;
724 b.high = h >> 1;
725 b.low = 0;
726 f_mult(&b, &c);
727 f_add(&a, &b);
728 }
729 f_mult(&c, &tenths[0]);
730 }
731 }
732 }
733
734 /* exponent */
735 if (*p == 'e' || *p == 'E') {
736 /* sign of exponent */
737 if (*++p == '-') {
738 t = tenths;
739 p++;
740 } else {
741 t = tens;
742 if (*p == '+') {
743 p++;
744 }
745 }
746
747 /* get exponent */
748 e = 0;
749 do {
750 e *= 10;
751 e += *p++ - '0';
752 if (e >= 1024) {
753 return FALSE;
754 }
755 } while (isdigit(*p));
756
757 /* adjust number */
758 while (e != 0) {
759 if ((e & 1) != 0) {
760 f_mult(&a, t);
761 if (a.exp < 0x1000 || a.exp > 0xf000) {
762 break;
763 }
764 }
765 e >>= 1;
766 t++;
767 }
768 }
769 if (a.exp >= BIAS + 1023 &&
770 (a.exp > BIAS + 1023 || (a.high == 0x7fff && a.low >= 0x7fffff00L))) {
771 return FALSE;
772 }
773
774 f_ftoxf(&a, f);
775 *s = p;
776 return TRUE;
777 }
778
779 /*
780 * NAME: float->ftoa()
781 * DESCRIPTION: convert a float to a string
782 */
flt_ftoa(xfloat * f,char * buffer)783 void flt_ftoa(xfloat *f, char *buffer)
784 {
785 unsigned short i;
786 short e;
787 Uint n;
788 char *p;
789 flt *t, *t2;
790 char digits[10];
791 flt a;
792
793 f_xftof(f, &a);
794 if (a.exp == 0) {
795 strcpy(buffer, "0");
796 return;
797 }
798
799 if (a.sign != 0) {
800 *buffer++ = '-';
801 a.sign = 0;
802 }
803
804 /* reduce the float to range 1 .. 9.999999999, and extract exponent */
805 e = 0;
806 if (a.exp >= BIAS) {
807 /* >= 1 */
808 for (i = 10, t = &tens[9], t2 = &tenths[9]; i > 0; --i, --t, --t2) {
809 e <<= 1;
810 if (f_cmp(&a, t) >= 0) {
811 e |= 1;
812 f_mult(&a, t2);
813 }
814 }
815 } else {
816 /* < 1 */
817 for (i = 10, t = &tenths[9], t2 = &tens[9]; i > 0; --i, --t, --t2) {
818 e <<= 1;
819 if (f_cmp(&a, t) <= 0) {
820 e |= 1;
821 f_mult(&a, t2);
822 }
823 }
824 if (a.exp < BIAS) {
825 /* still < 1 */
826 f_mult(&a, &tens[0]);
827 e++;
828 }
829 e = -e;
830 }
831 f_mult(&a, &tens[3]);
832 f_37bits(&a);
833
834 /*
835 * obtain digits
836 */
837 f_add(&a, &half);
838 i = a.exp - BIAS + 1 - 15;
839 n = ((Uint) a.high << i) | (a.low >> (31 - i));
840 if (n == 1000000000L) {
841 p = digits + 8;
842 p[0] = '1';
843 p[1] = '\0';
844 i = 1;
845 e++;
846 } else {
847 while (n != 0 && n % 10 == 0) {
848 n /= 10;
849 }
850 p = digits + 9;
851 *p = '\0';
852 i = 0;
853 do {
854 i++;
855 *--p = '0' + n % 10;
856 n /= 10;
857 } while (n != 0);
858 }
859
860 if (e >= 9 || (e < -3 && i - e > 9)) {
861 buffer[0] = *p;
862 if (i != 1) {
863 buffer[1] = '.';
864 memcpy(buffer + 2, p + 1, i - 1);
865 i++;
866 }
867 buffer[i++] = 'e';
868 if (e >= 0) {
869 buffer[i] = '+';
870 } else {
871 buffer[i] = '-';
872 e = -e;
873 }
874 p = digits + 9;
875 do {
876 *--p = '0' + e % 10;
877 e /= 10;
878 } while (e != 0);
879 strcpy(buffer + i + 1, p);
880 } else if (e < 0) {
881 e = 1 - e;
882 memcpy(buffer, "0.0000000", e);
883 strcpy(buffer + e, p);
884 } else {
885 while (e >= 0) {
886 *buffer++ = (*p == '\0') ? '0' : *p++;
887 --e;
888 }
889 if (*p != '\0') {
890 *buffer = '.';
891 strcpy(buffer + 1, p);
892 } else {
893 *buffer = '\0';
894 }
895 }
896 }
897
898 /*
899 * NAME: float->itof()
900 * DESCRIPTION: convert an integer to a float
901 */
flt_itof(Int i,xfloat * f)902 void flt_itof(Int i, xfloat *f)
903 {
904 flt a;
905
906 f_itof(i, &a);
907 f_ftoxf(&a, f);
908 }
909
910 /*
911 * NAME: float->ftoi()
912 * DESCRIPTION: convert a float to an integer
913 */
flt_ftoi(xfloat * f)914 Int flt_ftoi(xfloat *f)
915 {
916 flt a;
917
918 f_xftof(f, &a);
919 f_round(&a);
920 return f_ftoi(&a);
921 }
922
923 /*
924 * NAME: float->add()
925 * DESCRIPTION: add two floats
926 */
flt_add(xfloat * f1,xfloat * f2)927 void flt_add(xfloat *f1, xfloat *f2)
928 {
929 flt a, b;
930
931 f_xftof(f2, &b);
932 f_xftof(f1, &a);
933 f_add(&a, &b);
934 f_ftoxf(&a, f1);
935 }
936
937 /*
938 * NAME: float->sub()
939 * DESCRIPTION: subtract a float from a float
940 */
flt_sub(xfloat * f1,xfloat * f2)941 void flt_sub(xfloat *f1, xfloat *f2)
942 {
943 flt a, b;
944
945 f_xftof(f2, &b);
946 f_xftof(f1, &a);
947 f_sub(&a, &b);
948 f_ftoxf(&a, f1);
949 }
950
951 /*
952 * NAME: float->mult()
953 * DESCRIPTION: multiply two floats
954 */
flt_mult(xfloat * f1,xfloat * f2)955 void flt_mult(xfloat *f1, xfloat *f2)
956 {
957 flt a, b;
958
959 f_xftof(f1, &a);
960 f_xftof(f2, &b);
961 f_mult(&a, &b);
962 f_ftoxf(&a, f1);
963 }
964
965 /*
966 * NAME: float->div()
967 * DESCRIPTION: divide a float by a float
968 */
flt_div(xfloat * f1,xfloat * f2)969 void flt_div(xfloat *f1, xfloat *f2)
970 {
971 flt a, b;
972
973 f_xftof(f2, &b);
974 f_xftof(f1, &a);
975 f_div(&a, &b);
976 f_ftoxf(&a, f1);
977 }
978
979 /*
980 * NAME: float->cmp()
981 * DESCRIPTION: compare two xfloats
982 */
flt_cmp(xfloat * f1,xfloat * f2)983 int flt_cmp(xfloat *f1, xfloat *f2)
984 {
985 if ((short) (f1->high ^ f2->high) < 0) {
986 return ((short) f1->high < 0) ? -1 : 1;
987 }
988
989 if (f1->high == f2->high && f1->low == f2->low) {
990 return 0;
991 }
992 if (f1->high <= f2->high && (f1->high < f2->high || f1->low < f2->low)) {
993 return ((short) f1->high < 0) ? 1 : -1;
994 }
995 return ((short) f1->high < 0) ? -1 : 1;
996 }
997
998 /*
999 * NAME: float->floor()
1000 * DESCRIPTION: round a float downwards
1001 */
flt_floor(xfloat * f)1002 void flt_floor(xfloat *f)
1003 {
1004 flt a, b;
1005
1006 f_xftof(f, &a);
1007 b = a;
1008 f_trunc(&b);
1009 if (b.sign != 0 && f_cmp(&a, &b) != 0) {
1010 f_sub(&b, &one);
1011 }
1012 f_ftoxf(&b, f);
1013 }
1014
1015 /*
1016 * NAME: float->ceil()
1017 * DESCRIPTION: round a float upwards
1018 */
flt_ceil(xfloat * f)1019 void flt_ceil(xfloat *f)
1020 {
1021 flt a, b;
1022
1023 f_xftof(f, &a);
1024 b = a;
1025 f_trunc(&b);
1026 if (b.sign == 0 && f_cmp(&a, &b) != 0) {
1027 f_add(&b, &one);
1028 }
1029 f_ftoxf(&b, f);
1030 }
1031
1032 /*
1033 * NAME: float->fmod()
1034 * DESCRIPTION: perform fmod
1035 */
flt_fmod(xfloat * f1,xfloat * f2)1036 void flt_fmod(xfloat *f1, xfloat *f2)
1037 {
1038 flt a, b, c;
1039 unsigned short sign;
1040
1041 f_xftof(f2, &b);
1042 if (b.exp == 0) {
1043 f_edom();
1044 }
1045 f_xftof(f1, &a);
1046 if (a.exp == 0) {
1047 return;
1048 }
1049
1050 sign = a.sign;
1051 a.sign = b.sign = 0;
1052 c = b;
1053 while (f_cmp(&a, &b) >= 0) {
1054 c.exp = a.exp;
1055 if (f_cmp(&a, &c) < 0) {
1056 c.exp--;
1057 }
1058 f_sub(&a, &c);
1059 }
1060
1061 a.sign = sign;
1062 f_ftoxf(&a, f1);
1063 }
1064
1065 /*
1066 * NAME: float->frexp()
1067 * DESCRIPTION: split a float into a fraction and an exponent
1068 */
flt_frexp(xfloat * f)1069 Int flt_frexp(xfloat *f)
1070 {
1071 short e;
1072
1073 if (f->high == 0) {
1074 return 0;
1075 }
1076 e = ((f->high & 0x7ff0) >> 4) - 1022;
1077 f->high = (f->high & 0x800f) | (1022 << 4);
1078 return e;
1079 }
1080
1081 /*
1082 * NAME: float->ldexp()
1083 * DESCRIPTION: make a float from a fraction and an exponent
1084 */
flt_ldexp(xfloat * f,Int exp)1085 void flt_ldexp(xfloat *f, Int exp)
1086 {
1087 if (f->high == 0) {
1088 return;
1089 }
1090 exp += (f->high & 0x7ff0) >> 4;
1091 if (exp <= 0) {
1092 f->high = 0;
1093 f->low = 0;
1094 return;
1095 }
1096 if (exp > 1023 + 1023) {
1097 f_erange();
1098 }
1099 f->high = (f->high & 0x800f) | (exp << 4);
1100 }
1101
1102 /*
1103 * NAME: float->modf()
1104 * DESCRIPTION: split float into fraction and integer part
1105 */
flt_modf(xfloat * f1,xfloat * f2)1106 void flt_modf(xfloat *f1, xfloat *f2)
1107 {
1108 flt a, b;
1109
1110 f_xftof(f1, &a);
1111 if (a.exp < BIAS) {
1112 b.exp = 0;
1113 } else {
1114 b = a;
1115 f_trunc(&b);
1116 f_sub(&a, &b);
1117 }
1118 f_ftoxf(&a, f1);
1119 f_ftoxf(&b, f2);
1120 }
1121
1122
1123 /*
1124 * The algorithms for much of the following are taken from the Cephes Math
1125 * Library 2.1, by Stephen L. Moshier.
1126 */
1127
1128 /*
1129 * NAME: f_poly()
1130 * DESCRIPTION: evaluate polynomial
1131 */
f_poly(flt * x,flt * coef,int n)1132 static void f_poly(flt *x, flt *coef, int n)
1133 {
1134 flt result;
1135
1136 result = *coef++;
1137 do {
1138 f_mult(&result, x);
1139 f_add(&result, coef++);
1140 } while (--n != 0);
1141
1142 *x = result;
1143 }
1144
1145 /*
1146 * NAME: f_poly1()
1147 * DESCRIPTION: evaluate polynomial with coefficient of x ** (n + 1) == 1.0.
1148 */
f_poly1(flt * x,flt * coef,int n)1149 static void f_poly1(flt *x, flt *coef, int n)
1150 {
1151 flt result;
1152
1153 result = *x;
1154 f_add(&result, coef++);
1155 do {
1156 f_mult(&result, x);
1157 f_add(&result, coef++);
1158 } while (--n != 0);
1159
1160 *x = result;
1161 }
1162
1163 /*
1164 * NAME: f_exp()
1165 * DESCRIPTION: internal version of exp(f)
1166 */
f_exp(flt * a)1167 static void f_exp(flt *a)
1168 {
1169 static flt p[] = {
1170 { 0x0000, 0x7ff2, 0x4228, 0x01073370L },
1171 { 0x0000, 0x7ff9, 0x7c1b, 0x4362a050L },
1172 { 0x0000, 0x7fff, 0x4000, 0x00000000L }
1173 };
1174 static flt q[] = {
1175 { 0x0000, 0x7fec, 0x64bd, 0x3130af58L },
1176 { 0x0000, 0x7ff6, 0x52b9, 0x2c76e408L },
1177 { 0x0000, 0x7ffc, 0x745c, 0x1b8352c0L },
1178 { 0x0000, 0x8000, 0x4000, 0x00000000L }
1179 };
1180 static flt log2e = { 0x0000, 0x7fff, 0x5c55, 0x0eca5704L };
1181 static flt c1 = { 0x0000, 0x7ffe, 0x58c0, 0x00000000L };
1182 static flt c2 = { 0x0000, 0x7ff2, 0x6f40, 0x20b8c218L };
1183 flt b, c;
1184 short n;
1185
1186 b = *a;
1187 f_mult(&b, &log2e);
1188 f_round(&b);
1189 n = f_ftoi(&b);
1190 c = b;
1191 f_mult(&c, &c1);
1192 f_sub(a, &c);
1193 f_mult(&b, &c2);
1194 f_add(a, &b);
1195
1196 b = *a;
1197 f_mult(&b, a);
1198 c = b;
1199 f_poly(&c, p, 2);
1200 f_mult(a, &c);
1201 f_poly(&b, q, 3);
1202 f_sub(&b, a);
1203 f_div(a, &b);
1204
1205 if (a->exp != 0) {
1206 a->exp++;
1207 }
1208 f_add(a, &one);
1209 if (a->exp != 0) {
1210 a->exp += n;
1211 }
1212 }
1213
1214 /*
1215 * NAME: float->exp()
1216 * DESCRIPTION: exp(f)
1217 */
flt_exp(xfloat * f)1218 void flt_exp(xfloat *f)
1219 {
1220 flt a;
1221
1222 f_xftof(f, &a);
1223 if (f_cmp(&a, &maxlog) > 0) {
1224 /* overflow */
1225 f_erange();
1226 }
1227 if (f_cmp(&a, &minlog) < 0) {
1228 /* underflow */
1229 a.exp = 0;
1230 } else {
1231 f_exp(&a);
1232 }
1233
1234 f_ftoxf(&a, f);
1235 }
1236
1237 static flt logp[] = {
1238 { 0x0000, 0x7ff0, 0x6026, 0x4ed4bf30L },
1239 { 0x0000, 0x7ffd, 0x7f9f, 0x5db2f2b4L },
1240 { 0x0000, 0x8001, 0x6902, 0x458cd8e8L },
1241 { 0x0000, 0x8003, 0x7726, 0x52fc7a84L },
1242 { 0x0000, 0x8004, 0x7939, 0x5ac9d7b8L },
1243 { 0x0000, 0x8004, 0x7178, 0x244a33a8L },
1244 { 0x0000, 0x8003, 0x4f8e, 0x4b136264L }
1245 };
1246 static flt logq[] = {
1247 { 0x0000, 0x8002, 0x7840, 0x2c1bf7a0L },
1248 { 0x0000, 0x8005, 0x52bd, 0x5a8f5cf4L },
1249 { 0x0000, 0x8006, 0x6e55, 0x0548968cL },
1250 { 0x0000, 0x8007, 0x4cd0, 0x22530620L },
1251 { 0x0000, 0x8006, 0x6b7a, 0x28551a68L },
1252 { 0x0000, 0x8004, 0x7755, 0x709d1394L }
1253 };
1254
1255 /*
1256 * NAME: float->log()
1257 * DESCRIPTION: log(f)
1258 */
flt_log(xfloat * f)1259 void flt_log(xfloat *f)
1260 {
1261 static flt r[] = {
1262 { 0x8000, 0x7ffe, 0x6510, 0x7bb8d81cL },
1263 { 0x0000, 0x8003, 0x418b, 0x78e604f8L },
1264 { 0x8000, 0x8005, 0x4024, 0x0c224454L }
1265 };
1266 static flt s[] = {
1267 { 0x8000, 0x8004, 0x4758, 0x1a87d8dcL },
1268 { 0x0000, 0x8007, 0x4e06, 0x002255c8L },
1269 { 0x8000, 0x8008, 0x6036, 0x12336680L }
1270 };
1271 static flt c1 = { 0x0000, 0x7ff2, 0x6f40, 0x20b8c218L };
1272 static flt c2 = { 0x0000, 0x7ffe, 0x58c0, 0x00000000L };
1273 flt a, b, c, d;
1274 short n;
1275
1276 f_xftof(f, &a);
1277 if (a.sign != 0 || a.exp == 0) {
1278 /* <= 0.0 */
1279 f_edom();
1280 }
1281
1282 n = a.exp - BIAS + 1;
1283 a.exp = BIAS - 1;
1284
1285 if (n > 2 || n < -2) {
1286 if (f_cmp(&a, &sqrth) < 0) {
1287 --n;
1288 f_sub(&a, &half);
1289 b = a;
1290 } else {
1291 b = a;
1292 f_sub(&a, &half);
1293 f_sub(&a, &half);
1294 }
1295 if (b.exp != 0) {
1296 --b.exp;
1297 }
1298 f_add(&b, &half);
1299
1300 f_div(&a, &b);
1301 b = a;
1302 f_mult(&b, &b);
1303 c = b;
1304 f_poly(&b, r, 2);
1305 f_mult(&b, &c);
1306 f_poly1(&c, s, 2);
1307 f_div(&b, &c);
1308 f_mult(&b, &a);
1309 f_add(&a, &b);
1310 } else {
1311 if (f_cmp(&a, &sqrth) < 0) {
1312 --n;
1313 a.exp++;
1314 }
1315 f_sub(&a, &one);
1316
1317 b = a;
1318 f_mult(&b, &a);
1319 c = a;
1320 f_poly(&c, logp, 6);
1321 f_mult(&c, &b);
1322 d = a;
1323 f_poly1(&d, logq, 5);
1324 f_div(&c, &d);
1325 f_mult(&c, &a);
1326 if (b.exp != 0) {
1327 --b.exp;
1328 }
1329 f_sub(&c, &b);
1330 f_add(&a, &c);
1331 }
1332
1333 if (n != 0) {
1334 f_itof((Int) n, &b);
1335 c = b;
1336 f_mult(&c, &c1);
1337 f_sub(&a, &c);
1338 f_mult(&b, &c2);
1339 f_add(&a, &b);
1340 }
1341
1342 f_ftoxf(&a, f);
1343 }
1344
1345 /*
1346 * NAME: float->log10()
1347 * DESCRIPTION: log10(f)
1348 */
flt_log10(xfloat * f)1349 void flt_log10(xfloat *f)
1350 {
1351 static flt l102a = { 0x0000, 0x7ffd, 0x4d00, 0x00000000L };
1352 static flt l102b = { 0x0000, 0x7ff3, 0x4135, 0x04fbcff8L };
1353 static flt l10ea = { 0x0000, 0x7ffd, 0x6f00, 0x00000000L };
1354 static flt l10eb = { 0x0000, 0x7ff4, 0x5bd8, 0x549b9438L };
1355 flt a, b, c, d;
1356 short n;
1357
1358 f_xftof(f, &a);
1359 if (a.sign != 0 || a.exp == 0) {
1360 /* <= 0.0 */
1361 f_edom();
1362 }
1363
1364 n = a.exp - BIAS + 1;
1365 a.exp = BIAS - 1;
1366
1367 if (f_cmp(&a, &sqrth) < 0) {
1368 --n;
1369 a.exp++;
1370 }
1371 f_sub(&a, &one);
1372
1373 b = a;
1374 f_mult(&b, &a);
1375 c = a;
1376 f_poly(&c, logp, 6);
1377 f_mult(&c, &b);
1378 d = a;
1379 f_poly1(&d, logq, 5);
1380 f_div(&c, &d);
1381 f_mult(&c, &a);
1382 if (b.exp != 0) {
1383 --b.exp;
1384 }
1385 f_sub(&c, &b);
1386
1387 b = a;
1388 f_add(&b, &c);
1389 f_mult(&b, &l10eb);
1390 f_mult(&a, &l10ea);
1391 f_add(&a, &b);
1392 f_mult(&c, &l10ea);
1393 f_add(&a, &c);
1394 f_itof((Int) n, &b);
1395 c = b;
1396 f_mult(&b, &l102b);
1397 f_add(&a, &b);
1398 f_mult(&c, &l102a);
1399 f_add(&a, &c);
1400
1401 f_ftoxf(&a, f);
1402 }
1403
1404 /*
1405 * NAME: f_powi()
1406 * DESCRIPTION: take a number to an integer power
1407 */
f_powi(flt * a,int n)1408 static void f_powi(flt *a, int n)
1409 {
1410 flt b;
1411 unsigned short sign;
1412 bool neg;
1413
1414 if (n == 0) {
1415 /* pow(x, 0.0) == 1.0 */
1416 *a = one;
1417 return;
1418 }
1419
1420 if (a->exp == 0) {
1421 if (n < 0) {
1422 /* negative power of 0.0 */
1423 f_edom();
1424 }
1425 /* pow(0.0, y) == 0.0 */
1426 return;
1427 }
1428
1429 sign = a->sign;
1430 a->sign = 0;
1431
1432 if (n < 0) {
1433 neg = TRUE;
1434 n = -n;
1435 } else {
1436 neg = FALSE;
1437 }
1438
1439 if (n & 1) {
1440 b = *a;
1441 } else {
1442 b = one;
1443 sign = 0;
1444 }
1445 while ((n >>= 1) != 0) {
1446 f_mult(a, a);
1447 if (a->exp > BIAS + 1023) {
1448 f_erange();
1449 }
1450 if (n & 1) {
1451 f_mult(&b, a);
1452 }
1453 }
1454 /* range of b is checked when converting back to xfloat */
1455
1456 b.sign = sign;
1457 if (neg) {
1458 *a = one;
1459 f_div(a, &b);
1460 } else {
1461 *a = b;
1462 }
1463 }
1464
1465 /*
1466 * NAME: float->pow()
1467 * DESCRIPTION: pow(f1, f2)
1468 */
flt_pow(xfloat * f1,xfloat * f2)1469 void flt_pow(xfloat *f1, xfloat *f2)
1470 {
1471 static flt p[] = {
1472 { 0x0000, 0x7ffd, 0x7f6e, 0x32feb6b8L },
1473 { 0x0000, 0x8000, 0x7777, 0x5fd53dc0L },
1474 { 0x0000, 0x8001, 0x7b32, 0x7afef1d8L },
1475 { 0x0000, 0x8001, 0x4aaa, 0x076cb938L }
1476 };
1477 static flt q[] = {
1478 { 0x0000, 0x8002, 0x4aaa, 0x69364124L },
1479 { 0x0000, 0x8003, 0x6fff, 0x7e838394L },
1480 { 0x0000, 0x8004, 0x4332, 0x78362ec8L },
1481 { 0x0000, 0x8002, 0x6fff, 0x0b2315d4L }
1482 };
1483 static flt aloga[] = {
1484 { 0x0000, 0x7fff, 0x4000, 0x00000000L },
1485 { 0x0000, 0x7ffe, 0x7a92, 0x5f454920L },
1486 { 0x0000, 0x7ffe, 0x7560, 0x31b9f748L },
1487 { 0x0000, 0x7ffe, 0x7066, 0x37bb0aa4L },
1488 { 0x0000, 0x7ffe, 0x6ba2, 0x3f32b5a8L },
1489 { 0x0000, 0x7ffe, 0x6712, 0x230547e0L },
1490 { 0x0000, 0x7ffe, 0x62b3, 0x4a845540L },
1491 { 0x0000, 0x7ffe, 0x5e84, 0x28e7d604L },
1492 { 0x0000, 0x7ffe, 0x5a82, 0x3cccfe78L },
1493 { 0x0000, 0x7ffe, 0x56ac, 0x0fba90a8L },
1494 { 0x0000, 0x7ffe, 0x52ff, 0x35aa6c54L },
1495 { 0x0000, 0x7ffe, 0x4f7a, 0x4c982468L },
1496 { 0x0000, 0x7ffe, 0x4c1b, 0x7c146370L },
1497 { 0x0000, 0x7ffe, 0x48e1, 0x74dceac4L },
1498 { 0x0000, 0x7ffe, 0x45ca, 0x7078faa4L },
1499 { 0x0000, 0x7ffe, 0x42d5, 0x30d9f314L },
1500 { 0x0000, 0x7ffe, 0x4000, 0x00000000L }
1501 };
1502 static flt alogb[] = {
1503 { 0x0000, 0x0000, 0x0000, 0x00000000L },
1504 { 0x0000, 0x7fc7, 0x4bb4, 0x05aeb670L },
1505 { 0x0000, 0x7fc8, 0x5e87, 0x1a68bb98L },
1506 { 0x0000, 0x7fc8, 0x5ba7, 0x62ad0c98L },
1507 { 0x8000, 0x7fc8, 0x6f74, 0x682764c8L },
1508 { 0x0000, 0x7fc6, 0x750e, 0x2f5fd884L },
1509 { 0x0000, 0x7fc7, 0x5bd1, 0x55a46304L },
1510 { 0x8000, 0x7fc7, 0x4641, 0x0373af14L },
1511 { 0x0000, 0x0000, 0x0000, 0x00000000L }
1512 };
1513 static flt r[] = {
1514 { 0x0000, 0x7fee, 0x7d8c, 0x0fafe528L },
1515 { 0x0000, 0x7ff2, 0x50be, 0x7cc1f924L },
1516 { 0x0000, 0x7ff5, 0x5761, 0x7d9095e0L },
1517 { 0x0000, 0x7ff8, 0x4eca, 0x56dde268L },
1518 { 0x0000, 0x7ffa, 0x71ac, 0x11ae0834L },
1519 { 0x0000, 0x7ffc, 0x7afe, 0x7bff058cL },
1520 { 0x0000, 0x7ffe, 0x58b9, 0x05fdf474L }
1521 };
1522 static flt log2ea = { 0x0000, 0x7ffd, 0x7154, 0x3b295c18L };
1523 static flt sixteenth = { 0x0000, 0x7ffb, 0x4000, 0x00000000L };
1524 flt a, b, c, d, e;
1525 int n, i;
1526 unsigned short sign;
1527
1528 f_xftof(f1, &a);
1529 f_xftof(f2, &b);
1530
1531 c = b;
1532 f_trunc(&c);
1533 if (f_cmp(&b, &c) == 0 && b.exp < 0x800e) {
1534 /* integer power < 32768 */
1535 f_powi(&a, (int) f_ftoi(&c));
1536 f_ftoxf(&a, f1);
1537 return;
1538 }
1539
1540 sign = a.sign;
1541 if (sign != 0) {
1542 if (f_cmp(&b, &c) != 0) {
1543 /* non-integer power of negative number */
1544 f_edom();
1545 }
1546 a.sign = 0;
1547 --c.exp;
1548 f_trunc(&c);
1549 if (f_cmp(&b, &c) == 0) {
1550 /* even power of negative number */
1551 sign = 0;
1552 }
1553 }
1554 if (a.exp == 0) {
1555 if (b.sign != 0) {
1556 /* negative power of 0.0 */
1557 f_edom();
1558 }
1559 /* pow(0.0, y) == 0.0 */
1560 return;
1561 }
1562
1563 n = a.exp - BIAS + 1;
1564 a.exp = BIAS - 1;
1565
1566 if (f_cmp(&a, &aloga[1]) >= 0) {
1567 i = 0;
1568 } else {
1569 i = 1;
1570 if (f_cmp(&a, &aloga[9]) <= 0) {
1571 i = 9;
1572 }
1573 if (f_cmp(&a, &aloga[i + 4]) <= 0) {
1574 i += 4;
1575 }
1576 if (f_cmp(&a, &aloga[i + 2]) <= 0) {
1577 i += 2;
1578 }
1579 i++;
1580 }
1581 f_sub(&a, &aloga[i]);
1582 f_sub(&a, &alogb[i >> 1]);
1583 f_div(&a, &aloga[i]);
1584
1585 c = a;
1586 f_mult(&c, &a);
1587 d = a;
1588 f_poly(&d, p, 3);
1589 f_mult(&d, &c);
1590 e = a;
1591 f_poly1(&e, q, 3);
1592 f_div(&d, &e);
1593 f_mult(&d, &a);
1594 if (c.exp != 0) {
1595 --c.exp;
1596 }
1597 f_sub(&d, &c);
1598
1599 c = d;
1600 f_mult(&d, &log2ea);
1601 f_add(&c, &d);
1602 d = a;
1603 f_mult(&d, &log2ea);
1604 f_add(&c, &d);
1605 f_add(&c, &a);
1606 f_mult(&c, &b);
1607
1608 f_itof((Int) -i, &d);
1609 if (d.exp != 0) {
1610 d.exp -= 4;
1611 }
1612 f_itof((Int) n, &e);
1613 f_add(&d, &e);
1614
1615 e = b;
1616 e.exp += 4;
1617 f_trunc(&e);
1618 if (e.exp != 0) {
1619 e.exp -= 4;
1620 }
1621 f_sub(&b, &e);
1622
1623 f_mult(&b, &d);
1624 f_add(&c, &b);
1625 b = c;
1626 if (b.exp != 0) {
1627 b.exp += 4;
1628 f_trunc(&b);
1629 if (b.exp != 0) {
1630 b.exp -= 4;
1631 }
1632 }
1633 f_sub(&c, &b);
1634
1635 f_mult(&d, &e);
1636 f_add(&b, &d);
1637 d = b;
1638 if (d.exp != 0) {
1639 d.exp += 4;
1640 f_trunc(&d);
1641 if (d.exp != 0) {
1642 d.exp -= 4;
1643 }
1644 }
1645 f_sub(&b, &d);
1646
1647 f_add(&b, &c);
1648 c = b;
1649 if (c.exp != 0) {
1650 c.exp += 4;
1651 f_trunc(&c);
1652 if (c.exp != 0) {
1653 c.exp -= 4;
1654 }
1655 }
1656 f_add(&d, &c);
1657 if (d.exp != 0) {
1658 d.exp += 4;
1659 }
1660
1661 if (d.exp >= 0x800d) {
1662 /* exponent >= 16384 */
1663 if (d.sign == 0) {
1664 /* overflow */
1665 f_erange();
1666 }
1667 /* underflow */
1668 a.exp = 0;
1669 f_ftoxf(&a, f1);
1670 return;
1671 }
1672 n = f_ftoi(&d);
1673 f_sub(&b, &c);
1674 if (b.sign == 0 && b.exp != 0) {
1675 n++;
1676 f_sub(&b, &sixteenth);
1677 }
1678
1679 a = b;
1680 f_poly(&a, r, 6);
1681 f_mult(&a, &b);
1682
1683 i = n / 16 + ((n < 0) ? 0 : 1);
1684 n = i * 16 - n;
1685 f_mult(&a, &aloga[n]);
1686 f_add(&a, &aloga[n]);
1687 if (a.exp != 0) {
1688 a.exp += i;
1689 }
1690 a.sign = sign;
1691
1692 f_ftoxf(&a, f1);
1693 }
1694
1695 /*
1696 * NAME: f_sqrt()
1697 * DESCRIPTION: internal version of sqrt(f)
1698 */
f_sqrt(flt * a)1699 static void f_sqrt(flt *a)
1700 {
1701 static flt c1 = { 0x0000, 0x7ffe, 0x4b8a, 0x371e5fa0L };
1702 static flt c2 = { 0x0000, 0x7ffd, 0x6ad4, 0x55de691cL };
1703 static flt sqrt2 = { 0x0000, 0x7fff, 0x5a82, 0x3cccfe78L };
1704 flt b, c;
1705 int n;
1706
1707 if (a->exp == 0) {
1708 return;
1709 }
1710
1711 b = *a;
1712 n = a->exp - BIAS + 1;
1713 a->exp = BIAS - 1;
1714 f_mult(a, &c1);
1715 f_add(a, &c2);
1716 if (n & 1) {
1717 f_mult(a, &sqrt2);
1718 }
1719 a->exp += n >> 1;
1720
1721 c = b;
1722 f_div(&c, a);
1723 f_add(a, &c);
1724 --a->exp;
1725 c = b;
1726 f_div(&c, a);
1727 f_add(a, &c);
1728 --a->exp;
1729 f_div(&b, a);
1730 f_add(a, &b);
1731 --a->exp;
1732 }
1733
1734 /*
1735 * NAME: float->sqrt()
1736 * DESCRIPTION: sqrt(f)
1737 */
flt_sqrt(xfloat * f)1738 void flt_sqrt(xfloat *f)
1739 {
1740 flt a;
1741
1742 f_xftof(f, &a);
1743 if (a.sign != 0) {
1744 f_edom();
1745 }
1746 f_sqrt(&a);
1747 f_ftoxf(&a, f);
1748 }
1749
1750 static flt sincof[] = {
1751 { 0x0000, 0x7fde, 0x5763, 0x7a3fa338L },
1752 { 0x8000, 0x7fe5, 0x6b97, 0x4b525240L },
1753 { 0x0000, 0x7fec, 0x5c77, 0x46acfa90L },
1754 { 0x8000, 0x7ff2, 0x6806, 0x40337fc0L },
1755 { 0x0000, 0x7ff8, 0x4444, 0x222221f0L },
1756 { 0x8000, 0x7ffc, 0x5555, 0x2aaaaaacL }
1757 };
1758 static flt coscof[] = {
1759 { 0x8000, 0x7fda, 0x63e9, 0x13410c34L },
1760 { 0x0000, 0x7fe2, 0x47ba, 0x3af69c80L },
1761 { 0x8000, 0x7fe9, 0x49f9, 0x1efd5898L },
1762 { 0x0000, 0x7fef, 0x6806, 0x40339088L },
1763 { 0x8000, 0x7ff5, 0x5b05, 0x582d82a0L },
1764 { 0x0000, 0x7ffa, 0x5555, 0x2aaaaaacL }
1765 };
1766 static flt sc1 = { 0x0000, 0x7ffe, 0x6487, 0x76800000L };
1767 static flt sc2 = { 0x0000, 0x7fe6, 0x5110, 0x5a000000L };
1768 static flt sc3 = { 0x0000, 0x7fce, 0x611a, 0x313198a4L };
1769
1770 /*
1771 * NAME: float->cos()
1772 * DESCRIPTION: cos(f)
1773 */
flt_cos(xfloat * f)1774 void flt_cos(xfloat *f)
1775 {
1776 flt a, b, c;
1777 int n;
1778 unsigned short sign;
1779
1780 f_xftof(f, &a);
1781 if (a.exp >= 0x801d) {
1782 f_edom();
1783 }
1784
1785 a.sign = sign = 0;
1786 b = a;
1787 f_div(&b, &pio4);
1788 f_trunc(&b);
1789 n = f_ftoi(&b);
1790 if (n & 1) {
1791 n++;
1792 f_add(&b, &one);
1793 }
1794 n &= 7;
1795 if (n > 3) {
1796 sign = 0x8000;
1797 n -= 4;
1798 }
1799 if (n > 1) {
1800 sign ^= 0x8000;
1801 }
1802
1803 c = b;
1804 f_mult(&c, &sc1);
1805 f_sub(&a, &c);
1806 c = b;
1807 f_mult(&c, &sc2);
1808 f_sub(&a, &c);
1809 f_mult(&b, &sc3);
1810 f_sub(&a, &b);
1811
1812 b = a;
1813 f_mult(&b, &a);
1814 if (n == 1 || n == 2) {
1815 c = b;
1816 f_mult(&b, &a);
1817 f_poly(&c, sincof, 5);
1818 } else {
1819 a = one;
1820 c = b;
1821 if (c.exp != 0) {
1822 --c.exp;
1823 }
1824 f_sub(&a, &c);
1825 c = b;
1826 f_mult(&b, &b);
1827 f_poly(&c, coscof, 5);
1828 }
1829 f_mult(&b, &c);
1830 f_add(&a, &b);
1831 a.sign ^= sign;
1832
1833 f_ftoxf(&a, f);
1834 }
1835
1836 /*
1837 * NAME: float->sin()
1838 * DESCRIPTION: sin(f)
1839 */
flt_sin(xfloat * f)1840 void flt_sin(xfloat *f)
1841 {
1842 flt a, b, c;
1843 int n;
1844 unsigned short sign;
1845
1846 f_xftof(f, &a);
1847 if (a.exp >= 0x801d) {
1848 f_edom();
1849 }
1850
1851 sign = a.sign;
1852 a.sign = 0;
1853 b = a;
1854 f_div(&b, &pio4);
1855 f_trunc(&b);
1856 n = f_ftoi(&b);
1857 if (n & 1) {
1858 n++;
1859 f_add(&b, &one);
1860 }
1861 n &= 7;
1862 if (n > 3) {
1863 sign ^= 0x8000;
1864 n -= 4;
1865 }
1866
1867 c = b;
1868 f_mult(&c, &sc1);
1869 f_sub(&a, &c);
1870 c = b;
1871 f_mult(&c, &sc2);
1872 f_sub(&a, &c);
1873 f_mult(&b, &sc3);
1874 f_sub(&a, &b);
1875
1876 b = a;
1877 f_mult(&b, &a);
1878 if (n == 1 || n == 2) {
1879 a = one;
1880 c = b;
1881 if (c.exp != 0) {
1882 --c.exp;
1883 }
1884 f_sub(&a, &c);
1885 c = b;
1886 f_mult(&b, &b);
1887 f_poly(&c, coscof, 5);
1888 } else {
1889 c = b;
1890 f_mult(&b, &a);
1891 f_poly(&c, sincof, 5);
1892 }
1893 f_mult(&b, &c);
1894 f_add(&a, &b);
1895 a.sign ^= sign;
1896
1897 f_ftoxf(&a, f);
1898 }
1899
1900 /*
1901 * NAME: float->tan()
1902 * DESCRIPTION: float(f)
1903 */
flt_tan(xfloat * f)1904 void flt_tan(xfloat *f)
1905 {
1906 static flt p[] = {
1907 { 0x8000, 0x800c, 0x664b, 0x31a49e80L },
1908 { 0x0000, 0x8013, 0x4667, 0x594bf93cL },
1909 { 0x8000, 0x8017, 0x447f, 0x55a65324L }
1910 };
1911 static flt q[] = {
1912 { 0x0000, 0x800c, 0x6ae2, 0x4bdd66ccL },
1913 { 0x8000, 0x8013, 0x509e, 0x78b05578L },
1914 { 0x0000, 0x8017, 0x5f66, 0x1f85d5b0L },
1915 { 0x8000, 0x8018, 0x66bf, 0x40797cb4L }
1916 };
1917 static flt p1 = { 0x0000, 0x7ffe, 0x6487, 0x76800000L };
1918 static flt p2 = { 0x0000, 0x7fe6, 0x5110, 0x5a000000L };
1919 static flt p3 = { 0x0000, 0x7fce, 0x611a, 0x313198a4L };
1920 flt a, b, c;
1921 int n;
1922 unsigned short sign;
1923
1924 f_xftof(f, &a);
1925 if (a.exp >= 0x801d) {
1926 f_edom();
1927 }
1928
1929 sign = a.sign;
1930 a.sign = 0;
1931 b = a;
1932 f_div(&b, &pio4);
1933 f_trunc(&b);
1934 n = f_ftoi(&b);
1935 if (n & 1) {
1936 n++;
1937 f_add(&b, &one);
1938 }
1939
1940 c = b;
1941 f_mult(&c, &p1);
1942 f_sub(&a, &c);
1943 c = b;
1944 f_mult(&c, &p2);
1945 f_sub(&a, &c);
1946 f_mult(&b, &p3);
1947 f_sub(&a, &b);
1948
1949 b = a;
1950 f_mult(&b, &a);
1951 if (b.exp > 0x7fd0) { /* ~1e-14 */
1952 c = b;
1953 f_poly(&b, p, 2);
1954 f_mult(&b, &c);
1955 f_poly1(&c, q, 3);
1956 f_div(&b, &c);
1957 f_mult(&b, &a);
1958 f_add(&a, &b);
1959 }
1960
1961 if (n & 2) {
1962 b = one;
1963 f_div(&b, &a);
1964 a = b;
1965 a.sign ^= 0x8000;
1966 }
1967 a.sign ^= sign;
1968
1969 f_ftoxf(&a, f);
1970 }
1971
1972 static flt ascp[] = {
1973 { 0x8000, 0x7ffe, 0x5931, 0x3dd1792cL },
1974 { 0x0000, 0x8002, 0x5147, 0x2a1c6244L },
1975 { 0x8000, 0x8004, 0x4f77, 0x6c7ab96cL },
1976 { 0x0000, 0x8004, 0x7295, 0x0d081500L },
1977 { 0x8000, 0x8003, 0x6da8, 0x634b0bb0L }
1978 };
1979 static flt ascq[] = {
1980 { 0x8000, 0x8003, 0x5f58, 0x7442cc70L },
1981 { 0x0000, 0x8006, 0x4b8c, 0x15abd1acL },
1982 { 0x8000, 0x8007, 0x5f95, 0x630cc2e0L },
1983 { 0x0000, 0x8007, 0x6871, 0x0dbab9b8L },
1984 { 0x8000, 0x8006, 0x523e, 0x4a7848c4L }
1985 };
1986
1987 /*
1988 * NAME: float->acos()
1989 * DESCRIPTION: acos(f)
1990 */
flt_acos(xfloat * f)1991 void flt_acos(xfloat *f)
1992 {
1993 flt a, b, c;
1994 unsigned short sign;
1995 bool flag;
1996
1997 f_xftof(f, &a);
1998 sign = a.sign;
1999 a.sign = 0;
2000 if (f_cmp(&a, &one) > 0) {
2001 f_edom();
2002 }
2003
2004 if (f_cmp(&a, &half) > 0) {
2005 b = half;
2006 f_sub(&b, &a);
2007 f_add(&b, &half);
2008 if (b.exp != 0) {
2009 --b.exp;
2010 }
2011 a = b;
2012 f_sqrt(&a);
2013 flag = TRUE;
2014 } else {
2015 b = a;
2016 f_mult(&b, &a);
2017 flag = FALSE;
2018 }
2019
2020 if (a.exp >= 0x7fe7) { /* ~1e-7 */
2021 c = b;
2022 f_poly(&c, ascp, 4);
2023 f_mult(&c, &b);
2024 f_poly1(&b, ascq, 4);
2025 f_div(&c, &b);
2026 f_mult(&c, &a);
2027 f_add(&a, &c);
2028 }
2029
2030 if (flag) {
2031 if (a.exp != 0) {
2032 a.exp++;
2033 }
2034 if (sign != 0) {
2035 b = pi;
2036 f_sub(&b, &a);
2037 a = b;
2038 }
2039 } else {
2040 if (sign != 0) {
2041 f_add(&a, &pio2);
2042 } else {
2043 b = pio2;
2044 f_sub(&b, &a);
2045 a = b;
2046 }
2047 }
2048
2049 f_ftoxf(&a, f);
2050 }
2051
2052 /*
2053 * NAME: float->asin()
2054 * DESCRIPTION: asin(f)
2055 */
flt_asin(xfloat * f)2056 void flt_asin(xfloat *f)
2057 {
2058 flt a, b, c;
2059 unsigned short sign;
2060 bool flag;
2061
2062 f_xftof(f, &a);
2063 sign = a.sign;
2064 a.sign = 0;
2065 if (f_cmp(&a, &one) > 0) {
2066 f_edom();
2067 }
2068
2069 if (f_cmp(&a, &half) > 0) {
2070 b = half;
2071 f_sub(&b, &a);
2072 f_add(&b, &half);
2073 if (b.exp != 0) {
2074 --b.exp;
2075 }
2076 a = b;
2077 f_sqrt(&a);
2078 flag = TRUE;
2079 } else {
2080 b = a;
2081 f_mult(&b, &a);
2082 flag = FALSE;
2083 }
2084
2085 if (a.exp >= 0x7fe7) { /* ~1e-7 */
2086 c = b;
2087 f_poly(&c, ascp, 4);
2088 f_mult(&c, &b);
2089 f_poly1(&b, ascq, 4);
2090 f_div(&c, &b);
2091 f_mult(&c, &a);
2092 f_add(&a, &c);
2093 }
2094
2095 if (flag) {
2096 if (a.exp != 0) {
2097 a.exp++;
2098 }
2099 b = pio2;
2100 f_sub(&b, &a);
2101 a = b;
2102 }
2103 a.sign ^= sign;
2104
2105 f_ftoxf(&a, f);
2106 }
2107
2108 static flt atp[] = {
2109 { 0x8000, 0x7ffe, 0x6ba5, 0x2175f64cL },
2110 { 0x8000, 0x8002, 0x46b5, 0x3c27114cL },
2111 { 0x8000, 0x8003, 0x5763, 0x7b6b8ba4L },
2112 { 0x8000, 0x8002, 0x76a5, 0x2457275cL }
2113 };
2114 static flt atq[] = {
2115 { 0x0000, 0x8002, 0x7bfa, 0x59b1a2acL },
2116 { 0x0000, 0x8004, 0x7d94, 0x68676274L },
2117 { 0x0000, 0x8005, 0x5c3c, 0x7b2444acL },
2118 { 0x0000, 0x8004, 0x58fb, 0x7b415d88L }
2119 };
2120 static flt t3p8 = { 0x0000, 0x8000, 0x4d41, 0x1e667f3cL };
2121 static flt tp8 = { 0x0000, 0x7ffd, 0x6a09, 0x7333f9e0L };
2122
2123 /*
2124 * NAME: float->atan()
2125 * DESCRIPTION: atan(f)
2126 */
flt_atan(xfloat * f)2127 void flt_atan(xfloat *f)
2128 {
2129 flt a, b, c, d, e;
2130 unsigned short sign;
2131
2132 f_xftof(f, &a);
2133 sign = a.sign;
2134 a.sign = 0;
2135
2136 if (f_cmp(&a, &t3p8) > 0) {
2137 b = pio2;
2138 c = one;
2139 f_div(&c, &a);
2140 a = c;
2141 a.sign = 0x8000;
2142 } else if (f_cmp(&a, &tp8) > 0) {
2143 b = pio4;
2144 c = a;
2145 f_sub(&a, &one);
2146 f_add(&c, &one);
2147 f_div(&a, &c);
2148 } else {
2149 b.exp = 0;
2150 }
2151
2152 c = a;
2153 f_mult(&c, &a);
2154 d = e = c;
2155 f_poly(&c, atp, 3);
2156 f_poly1(&d, atq, 3);
2157 f_div(&c, &d);
2158 f_mult(&c, &e);
2159 f_mult(&c, &a);
2160 f_add(&c, &b);
2161 f_add(&a, &c);
2162 a.sign ^= sign;
2163
2164 f_ftoxf(&a, f);
2165 }
2166
2167 /*
2168 * NAME: float->atan2()
2169 * DESCRIPTION: atan2(f)
2170 */
flt_atan2(xfloat * f1,xfloat * f2)2171 void flt_atan2(xfloat *f1, xfloat *f2)
2172 {
2173 flt a, b, c, d, e;
2174 unsigned short asign, bsign;
2175
2176 f_xftof(f1, &a);
2177 f_xftof(f2, &b);
2178
2179 if (b.exp == 0) {
2180 if (a.exp == 0) {
2181 /* atan2(0.0, 0.0); */
2182 return;
2183 }
2184 a.exp = pio2.exp;
2185 a.high = pio2.high;
2186 a.low = pio2.low;
2187 f_ftoxf(&a, f1);
2188 return;
2189 }
2190 if (a.exp == 0) {
2191 if (b.sign != 0) {
2192 a = pi;
2193 }
2194 f_ftoxf(&a, f1);
2195 return;
2196 }
2197
2198 asign = a.sign;
2199 bsign = b.sign;
2200 f_div(&a, &b);
2201 a.sign = 0;
2202
2203 if (f_cmp(&a, &t3p8) > 0) {
2204 b = pio2;
2205 c = one;
2206 f_div(&c, &a);
2207 a = c;
2208 a.sign = 0x8000;
2209 } else if (f_cmp(&a, &tp8) > 0) {
2210 b = pio4;
2211 c = a;
2212 f_sub(&a, &one);
2213 f_add(&c, &one);
2214 f_div(&a, &c);
2215 } else {
2216 b.exp = 0;
2217 }
2218
2219 c = a;
2220 f_mult(&c, &a);
2221 d = e = c;
2222 f_poly(&c, atp, 3);
2223 f_poly1(&d, atq, 3);
2224 f_div(&c, &d);
2225 f_mult(&c, &e);
2226 f_mult(&c, &a);
2227 f_add(&c, &b);
2228 f_add(&a, &c);
2229 a.sign ^= asign ^ bsign;
2230
2231 if (bsign != 0) {
2232 if (asign == 0) {
2233 f_add(&a, &pi);
2234 } else {
2235 f_sub(&a, &pi);
2236 }
2237 }
2238
2239 f_ftoxf(&a, f1);
2240 }
2241
2242 /*
2243 * NAME: float->cosh()
2244 * DESCRIPTION: cosh(f)
2245 */
flt_cosh(xfloat * f)2246 void flt_cosh(xfloat *f)
2247 {
2248 flt a, b;
2249
2250 f_xftof(f, &a);
2251 a.sign = 0;
2252 if (f_cmp(&a, &maxlog) > 0) {
2253 f_erange();
2254 }
2255
2256 f_exp(&a);
2257 b = one;
2258 f_div(&b, &a);
2259 f_add(&a, &b);
2260 --a.exp;
2261
2262 f_ftoxf(&a, f);
2263 }
2264
2265 /*
2266 * NAME: float->sinh()
2267 * DESCRIPTION: sinh(f)
2268 */
flt_sinh(xfloat * f)2269 void flt_sinh(xfloat *f)
2270 {
2271 static flt p[] = {
2272 { 0x8000, 0x7ffe, 0x650d, 0x3fd17678L },
2273 { 0x8000, 0x8006, 0x51dc, 0x74731fe8L },
2274 { 0x8000, 0x800c, 0x5a52, 0x718e3ac4L },
2275 { 0x8000, 0x8011, 0x55e0, 0x57b7ed58L }
2276 };
2277 static flt q[] = {
2278 { 0x8000, 0x8007, 0x456d, 0x412dd2c8L },
2279 { 0x0000, 0x800e, 0x469e, 0x74fdae44L },
2280 { 0x8000, 0x8014, 0x4068, 0x41c9f200L }
2281 };
2282 flt a, b, c, d;
2283 unsigned short sign;
2284
2285 f_xftof(f, &a);
2286 if (f_cmp(&a, &maxlog) > 0 || f_cmp(&a, &minlog) < 0) {
2287 f_erange();
2288 }
2289
2290 sign = a.sign;
2291 a.sign = 0;
2292
2293 if (f_cmp(&a, &one) > 0) {
2294 f_exp(&a);
2295 b = half;
2296 f_div(&b, &a);
2297 --a.exp;
2298 f_sub(&a, &b);
2299 a.sign ^= sign;
2300 } else {
2301 b = a;
2302 f_mult(&b, &a);
2303 c = d = b;
2304 f_poly(&c, p, 3);
2305 f_poly1(&d, q, 2);
2306 f_div(&c, &d);
2307 f_mult(&b, &a);
2308 f_mult(&b, &c);
2309 f_add(&a, &b);
2310 }
2311
2312 f_ftoxf(&a, f);
2313 }
2314
2315 /*
2316 * NAME: float->tanh()
2317 * DESCRIPTION: tanh(f)
2318 */
flt_tanh(xfloat * f)2319 void flt_tanh(xfloat *f)
2320 {
2321 static flt p[] = {
2322 { 0x8000, 0x7ffe, 0x7b71, 0x3755fae0L },
2323 { 0x8000, 0x8005, 0x6349, 0x541c4cd0L },
2324 { 0x8000, 0x8009, 0x64eb, 0x0060b00cL }
2325 };
2326 static flt q[] = {
2327 { 0x0000, 0x8005, 0x70cf, 0x6514b038L },
2328 { 0x0000, 0x800a, 0x45db, 0x741caa6cL },
2329 { 0x0000, 0x800b, 0x4bb0, 0x20488408L }
2330 };
2331 static flt mlog2 = { 0x0000, 0x8007, 0x58b9, 0x05fdf474L };
2332 static flt d625 = { 0x0000, 0x7ffe, 0x5000, 0x00000000L };
2333 static flt two = { 0x0000, 0x8000, 0x4000, 0x00000000L };
2334 flt a, b, c, d;
2335 unsigned short sign;
2336
2337 f_xftof(f, &a);
2338 sign = a.sign;
2339 a.sign = 0;
2340
2341 if (f_cmp(&a, &mlog2) > 0) {
2342 a.exp = one.exp;
2343 a.high = one.high;
2344 a.low = one.low;
2345 } else if (f_cmp(&a, &d625) >= 0) {
2346 a.exp++;
2347 f_exp(&a);
2348 f_add(&a, &one);
2349 b = two;
2350 f_div(&b, &a);
2351 a = one;
2352 f_sub(&a, &b);
2353 } else if (a.exp != 0) {
2354 b = a;
2355 f_mult(&b, &a);
2356 c = d = b;
2357 f_poly(&c, p, 2);
2358 f_poly1(&d, q, 2);
2359 f_div(&c, &d);
2360 f_mult(&b, &c);
2361 f_mult(&b, &a);
2362 f_add(&a, &b);
2363 }
2364 a.sign = sign;
2365
2366 f_ftoxf(&a, f);
2367 }
2368