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