1 /* Copyright (C) 1989, 1996, 1998 Aladdin Enterprises. All rights reserved.
2
3 This software is provided AS-IS with no warranty, either express or
4 implied.
5
6 This software is distributed under license and may not be copied,
7 modified or distributed except as expressly authorized under the terms
8 of the license contained in the file LICENSE in this distribution.
9
10 For more information about licensing, please refer to
11 http://www.ghostscript.com/licensing/. For information on
12 commercial licensing, go to http://www.artifex.com/licensing/ or
13 contact Artifex Software, Inc., 101 Lucas Valley Road #110,
14 San Rafael, CA 94903, U.S.A., +1(415)492-9861.
15 */
16
17 /* $Id: gsfemu.c,v 1.4 2002/02/21 22:24:52 giles Exp $ */
18 /* Floating point emulator for gcc */
19
20 /* We actually only need arch.h + uint and ulong, but because signal.h */
21 /* may include <sys/types.h>, we have to include std.h to handle (avoid) */
22 /* redefinition of type names. */
23 #include "std.h"
24 #include <signal.h>
25
26 /*#define TEST */
27
28 /*
29 * This package is not a fully general IEEE floating point implementation.
30 * It implements only one rounding mode (round to nearest) and does not
31 * generate or properly handle NaNs.
32 *
33 * The names and interfaces of the routines in this package were designed to
34 * work with gcc. This explains some peculiarities such as the presence of
35 * single-precision arithmetic (forbidden by the ANSI standard) and the
36 * omission of unsigned-to-float conversions (which gcc implements with
37 * truly grotesque inline code).
38 *
39 * The following routines do not yet handle denormalized numbers:
40 * addd3 (operands or result)
41 * __muldf3 (operands or result)
42 * __divdf3 (operands or result)
43 * __truncdfsf2 (result)
44 * __extendsfdf2 (operand)
45 */
46
47 /*
48 * IEEE single-precision floats have the format:
49 * sign(1) exponent(8) mantissa(23)
50 * The exponent is biased by +127.
51 * The mantissa is a fraction with an implicit integer part of 1,
52 * unless the exponent is zero.
53 */
54 #define fx(ia) (((ia) >> 23) & 0xff)
55 #define fx_bias 127
56 /*
57 * IEEE double-precision floats have the format:
58 * sign(1) exponent(11) mantissa(52)
59 * The exponent is biased by +1023.
60 */
61 #define dx(ld) ((ld[msw] >> 20) & 0x7ff)
62 #define dx_bias 1023
63
64 #if arch_is_big_endian
65 # define msw 0
66 # define lsw 1
67 #else
68 # define msw 1
69 # define lsw 0
70 #endif
71 /* Double arguments/results */
72 #define la ((const long *)&a)
73 #define ula ((const ulong *)&a)
74 #define lb ((const long *)&b)
75 #define ulb ((const ulong *)&b)
76 #define dc (*(const double *)lc)
77 /* Float arguments/results */
78 #define ia (*(const long *)&a)
79 #define ua (*(const ulong *)&a)
80 #define ib (*(const long *)&b)
81 #define ub (*(const ulong *)&b)
82 #define fc (*(const float *)&lc)
83
84 /* Round a double-length fraction by adding 1 in the lowest bit and */
85 /* then shifting right by 1. */
86 #define roundr1(ms, uls)\
87 if ( uls == 0xffffffff ) ms++, uls = 0;\
88 else uls++;\
89 uls = (uls >> 1) + (ms << 31);\
90 ms >>= 1
91
92 /* Extend a non-zero float to a double. */
93 #define extend(lc, ia)\
94 ((lc)[msw] = ((ia) & 0x80000000) + (((ia) & 0x7fffffff) >> 3) + 0x38000000,\
95 (lc)[lsw] = (ia) << 29)
96
97 /* ---------------- Arithmetic ---------------- */
98
99 /* -------- Add/subtract/negate -------- */
100
101 double
__negdf2(double a)102 __negdf2(double a)
103 {
104 long lc[2];
105
106 lc[msw] = la[msw] ^ 0x80000000;
107 lc[lsw] = la[lsw];
108 return dc;
109 }
110 float
__negsf2(float a)111 __negsf2(float a)
112 {
113 long lc = ia ^ 0x80000000;
114
115 return fc;
116 }
117
118 double
__adddf3(double a,double b)119 __adddf3(double a, double b)
120 {
121 long lc[2];
122 int expt = dx(la);
123 int shift = expt - dx(lb);
124 long sign;
125 ulong msa, lsa;
126 ulong msb, lsb;
127
128 if (shift < 0) { /* Swap a and b so that expt(a) >= expt(b). */
129 double temp = a;
130
131 a = b, b = temp;
132 expt += (shift = -shift);
133 }
134 if (shift >= 54) /* also picks up most cases where b == 0 */
135 return a;
136 if (!(lb[msw] & 0x7fffffff))
137 return a;
138 sign = la[msw] & 0x80000000;
139 msa = (la[msw] & 0xfffff) + 0x100000, lsa = la[lsw];
140 msb = (lb[msw] & 0xfffff) + 0x100000, lsb = lb[lsw];
141 if ((la[msw] ^ lb[msw]) >= 0) { /* Adding numbers of the same sign. */
142 if (shift >= 32)
143 lsb = msb, msb = 0, shift -= 32;
144 if (shift) {
145 --shift;
146 lsb = (lsb >> shift) + (msb << (32 - shift));
147 msb >>= shift;
148 roundr1(msb, lsb);
149 }
150 if (lsb > (ulong) 0xffffffff - lsa)
151 msa++;
152 lsa += lsb;
153 msa += msb;
154 if (msa > 0x1fffff) {
155 roundr1(msa, lsa);
156 /* In principle, we have to worry about exponent */
157 /* overflow here, but we don't. */
158 ++expt;
159 }
160 } else { /* Adding numbers of different signs. */
161 if (shift > 53)
162 return a; /* b can't affect the result, even rounded */
163 if (shift == 0 && (msb > msa || (msb == msa && lsb >= lsa))) { /* This is the only case where the sign of the result */
164 /* differs from the sign of the first operand. */
165 sign ^= 0x80000000;
166 msa = msb - msa;
167 if (lsb < lsa)
168 msa--;
169 lsa = lsb - lsa;
170 } else {
171 if (shift >= 33) {
172 lsb = ((msb >> (shift - 32)) + 1) >> 1; /* round */
173 msb = 0;
174 } else if (shift) {
175 lsb = (lsb >> (shift - 1)) + (msb << (33 - shift));
176 msb >>= shift - 1;
177 roundr1(msb, lsb);
178 }
179 msa -= msb;
180 if (lsb > lsa)
181 msa--;
182 lsa -= lsb;
183 }
184 /* Now renormalize the result. */
185 /* For the moment, we do this the slow way. */
186 if (!(msa | lsa))
187 return 0;
188 while (msa < 0x100000) {
189 msa = (msa << 1) + (lsa >> 31);
190 lsa <<= 1;
191 expt -= 1;
192 }
193 if (expt <= 0) { /* Underflow. Return 0 rather than a denorm. */
194 lc[msw] = sign;
195 lc[lsw] = 0;
196 return dc;
197 }
198 }
199 lc[msw] = sign + ((ulong) expt << 20) + (msa & 0xfffff);
200 lc[lsw] = lsa;
201 return dc;
202 }
203 double
__subdf3(double a,double b)204 __subdf3(double a, double b)
205 {
206 long nb[2];
207
208 nb[msw] = lb[msw] ^ 0x80000000;
209 nb[lsw] = lb[lsw];
210 return a + *(const double *)nb;
211 }
212
213 float
__addsf3(float a,float b)214 __addsf3(float a, float b)
215 {
216 long lc;
217 int expt = fx(ia);
218 int shift = expt - fx(ib);
219 long sign;
220 ulong ma, mb;
221
222 if (shift < 0) { /* Swap a and b so that expt(a) >= expt(b). */
223 long temp = ia;
224
225 *(long *)&a = ib;
226 *(long *)&b = temp;
227 expt += (shift = -shift);
228 }
229 if (shift >= 25) /* also picks up most cases where b == 0 */
230 return a;
231 if (!(ib & 0x7fffffff))
232 return a;
233 sign = ia & 0x80000000;
234 ma = (ia & 0x7fffff) + 0x800000;
235 mb = (ib & 0x7fffff) + 0x800000;
236 if ((ia ^ ib) >= 0) { /* Adding numbers of the same sign. */
237 if (shift) {
238 --shift;
239 mb = ((mb >> shift) + 1) >> 1;
240 }
241 ma += mb;
242 if (ma > 0xffffff) {
243 ma = (ma + 1) >> 1;
244 /* In principle, we have to worry about exponent */
245 /* overflow here, but we don't. */
246 ++expt;
247 }
248 } else { /* Adding numbers of different signs. */
249 if (shift > 24)
250 return a; /* b can't affect the result, even rounded */
251 if (shift == 0 && mb >= ma) {
252 /* This is the only case where the sign of the result */
253 /* differs from the sign of the first operand. */
254 sign ^= 0x80000000;
255 ma = mb - ma;
256 } else {
257 if (shift) {
258 --shift;
259 mb = ((mb >> shift) + 1) >> 1;
260 }
261 ma -= mb;
262 }
263 /* Now renormalize the result. */
264 /* For the moment, we do this the slow way. */
265 if (!ma)
266 return 0;
267 while (ma < 0x800000) {
268 ma <<= 1;
269 expt -= 1;
270 }
271 if (expt <= 0) {
272 /* Underflow. Return 0 rather than a denorm. */
273 lc = sign;
274 return fc;
275 }
276 }
277 lc = sign + ((ulong)expt << 23) + (ma & 0x7fffff);
278 return fc;
279 }
280
281 float
__subsf3(float a,float b)282 __subsf3(float a, float b)
283 {
284 long lc = ib ^ 0x80000000;
285
286 return a + fc;
287 }
288
289 /* -------- Multiplication -------- */
290
291 double
__muldf3(double a,double b)292 __muldf3(double a, double b)
293 {
294 long lc[2];
295 ulong sign;
296 uint H, I, h, i;
297 ulong p0, p1, p2;
298 ulong expt;
299
300 if (!(la[msw] & 0x7fffffff) || !(lb[msw] & 0x7fffffff))
301 return 0;
302 /*
303 * We now have to multiply two 53-bit fractions to produce a 53-bit
304 * result. Since the idiots who specified the standard C libraries
305 * don't allow us to use the 32 x 32 => 64 multiply that every
306 * 32-bit CPU provides, we have to chop these 53-bit numbers up into
307 * 14-bit chunks so we can use 32 x 32 => 32 multiplies.
308 */
309 #define chop_ms(ulx, h, i)\
310 h = ((ulx[msw] >> 7) & 0x1fff) | 0x2000,\
311 i = ((ulx[msw] & 0x7f) << 7) | (ulx[lsw] >> 25)
312 #define chop_ls(ulx, j, k)\
313 j = (ulx[lsw] >> 11) & 0x3fff,\
314 k = (ulx[lsw] & 0x7ff) << 3
315 chop_ms(ula, H, I);
316 chop_ms(ulb, h, i);
317 #undef chop
318 #define prod(m, n) ((ulong)(m) * (n)) /* force long multiply */
319 p0 = prod(H, h);
320 p1 = prod(H, i) + prod(I, h);
321 /* If these doubles were produced from floats or ints, */
322 /* all the other products will be zero. It's worth a check. */
323 if ((ula[lsw] | ulb[lsw]) & 0x1ffffff) { /* No luck. */
324 /* We can almost always avoid computing p5 and p6, */
325 /* but right now it's not worth the trouble to check. */
326 uint J, K, j, k;
327
328 chop_ls(ula, J, K);
329 chop_ls(ulb, j, k);
330 {
331 ulong p6 = prod(K, k);
332 ulong p5 = prod(J, k) + prod(K, j) + (p6 >> 14);
333 ulong p4 = prod(I, k) + prod(J, j) + prod(K, i) + (p5 >> 14);
334 ulong p3 = prod(H, k) + prod(I, j) + prod(J, i) + prod(K, h) +
335 (p4 >> 14);
336
337 p2 = prod(H, j) + prod(I, i) + prod(J, h) + (p3 >> 14);
338 }
339 } else {
340 p2 = prod(I, i);
341 }
342 /* Now p0 through p2 hold the result. */
343 /****** ASSUME expt IS 32 BITS WIDE. ******/
344 expt = (la[msw] & 0x7ff00000) + (lb[msw] & 0x7ff00000) -
345 (dx_bias << 20);
346 /* Now expt is in the range [-1023..3071] << 20. */
347 /* Values outside the range [1..2046] are invalid. */
348 p1 += p2 >> 14;
349 p0 += p1 >> 14;
350 /* Now the 56-bit result consists of p0, the low 14 bits of p1, */
351 /* and the low 14 bits of p2. */
352 /* p0 can be anywhere between 2^26 and almost 2^28, so we might */
353 /* need to adjust the exponent by 1. */
354 if (p0 < 0x8000000) {
355 p0 = (p0 << 1) + ((p1 >> 13) & 1);
356 p1 = (p1 << 1) + ((p2 >> 13) & 1);
357 p2 <<= 1;
358 } else
359 expt += 0x100000;
360 /* The range of expt is now [-1023..3072] << 20. */
361 /* Round the mantissa to 53 bits. */
362 if (!((p2 += 4) & 0x3ffc) && !(++p1 & 0x3fff) && ++p0 >= 0x10000000) {
363 p0 >>= 1, p1 = 0x2000;
364 /* Check for exponent overflow, just in case. */
365 if ((ulong) expt < 0xc0000000)
366 expt += 0x100000;
367 }
368 sign = (la[msw] ^ lb[msw]) & 0x80000000;
369 if (expt == 0) { /* Underflow. Return 0 rather than a denorm. */
370 lc[msw] = sign;
371 lc[lsw] = 0;
372 } else if ((ulong) expt >= 0x7ff00000) { /* Overflow or underflow. */
373 if ((ulong) expt <= 0xc0000000) { /* Overflow. */
374 raise(SIGFPE);
375 lc[msw] = sign + 0x7ff00000;
376 lc[lsw] = 0;
377 } else { /* Underflow. See above. */
378 lc[msw] = sign;
379 lc[lsw] = 0;
380 }
381 } else {
382 lc[msw] = sign + expt + ((p0 >> 7) & 0xfffff);
383 lc[lsw] = (p0 << 25) | ((p1 & 0x3fff) << 11) | ((p2 >> 3) & 0x7ff);
384 }
385 return dc;
386 #undef prod
387 }
388 float
__mulsf3(float a,float b)389 __mulsf3(float a, float b)
390 {
391 uint au, al, bu, bl, cu, cl, sign;
392 long lc;
393 uint expt;
394
395 if (!(ia & 0x7fffffff) || !(ib & 0x7fffffff))
396 return 0;
397 au = ((ia >> 8) & 0x7fff) | 0x8000;
398 bu = ((ib >> 8) & 0x7fff) | 0x8000;
399 /* Since 0x8000 <= au,bu <= 0xffff, 0x40000000 <= cu <= 0xfffe0001. */
400 cu = au * bu;
401 if ((al = ia & 0xff) != 0) {
402 cl = bu * al;
403 } else
404 cl = 0;
405 if ((bl = ib & 0xff) != 0) {
406 cl += au * bl;
407 if (al)
408 cl += (al * bl) >> 8;
409 }
410 cu += cl >> 8;
411 sign = (ia ^ ib) & 0x80000000;
412 expt = (ia & 0x7f800000) + (ib & 0x7f800000) - (fx_bias << 23);
413 /* expt is now in the range [-127..383] << 23. */
414 /* Values outside [1..254] are invalid. */
415 if (cu <= 0x7fffffff)
416 cu <<= 1;
417 else
418 expt += 1 << 23;
419 cu = ((cu >> 7) + 1) >> 1;
420 if (expt < 1 << 23)
421 lc = sign; /* underflow */
422 else if (expt > (uint)(254 << 23)) {
423 if (expt <= 0xc0000000) { /* overflow */
424 raise(SIGFPE);
425 lc = sign + 0x7f800000;
426 } else { /* underflow */
427 lc = sign;
428 }
429 } else
430 lc = sign + expt + cu - 0x800000;
431 return fc;
432 }
433
434 /* -------- Division -------- */
435
436 double
__divdf3(double a,double b)437 __divdf3(double a, double b)
438 {
439 long lc[2];
440
441 /*
442 * Multi-precision division is really, really awful.
443 * We generate the result more or less brute force,
444 * 11 bits at a time.
445 */
446 ulong sign = (la[msw] ^ lb[msw]) & 0x80000000;
447 ulong msa = (la[msw] & 0xfffff) | 0x100000, lsa = la[lsw];
448 ulong msb = (lb[msw] & 0xfffff) | 0x100000, lsb = lb[lsw];
449 uint qn[5];
450 int i;
451 ulong msq, lsq;
452 int expt = dx(la) - dx(lb) + dx_bias;
453
454 if (!(lb[msw] & 0x7fffffff)) { /* Division by zero. */
455 raise(SIGFPE);
456 lc[lsw] = 0;
457 lc[msw] =
458 (la[msw] & 0x7fffffff ?
459 sign + 0x7ff00000 /* infinity */ :
460 0x7ff80000 /* NaN */ );
461 return dc;
462 }
463 if (!(la[msw] & 0x7fffffff))
464 return 0;
465 for (i = 0; i < 5; ++i) {
466 uint q;
467 ulong msp, lsp;
468
469 msa = (msa << 11) + (lsa >> 21);
470 lsa <<= 11;
471 q = msa / msb;
472 /* We know that 2^20 <= msb < 2^21; q could be too high by 1, */
473 /* but it can't be too low. */
474 /* Set p = d * q. */
475 msp = q * msb;
476 lsp = q * (lsb & 0x1fffff);
477 {
478 ulong midp = q * (lsb >> 21);
479
480 msp += (midp + (lsp >> 21)) >> 11;
481 lsp += midp << 21;
482 }
483 /* Set a = a - p, i.e., the tentative remainder. */
484 if (msp > msa || (lsp > lsa && msp == msa)) { /* q was too large. */
485 --q;
486 if (lsb > lsp)
487 msp--;
488 lsp -= lsb;
489 msp -= msb;
490 }
491 if (lsp > lsa)
492 msp--;
493 lsa -= lsp;
494 msa -= msp;
495 qn[i] = q;
496 }
497 /* Pack everything together again. */
498 msq = (qn[0] << 9) + (qn[1] >> 2);
499 lsq = (qn[1] << 30) + (qn[2] << 19) + (qn[3] << 8) + (qn[4] >> 3);
500 if (msq < 0x100000) {
501 msq = (msq << 1) + (lsq >> 31);
502 lsq <<= 1;
503 expt--;
504 }
505 /* We should round the quotient, but we don't. */
506 if (expt <= 0) { /* Underflow. Return 0 rather than a denorm. */
507 lc[msw] = sign;
508 lc[lsw] = 0;
509 } else if (expt >= 0x7ff) { /* Overflow. */
510 raise(SIGFPE);
511 lc[msw] = sign + 0x7ff00000;
512 lc[lsw] = 0;
513 } else {
514 lc[msw] = sign + (expt << 20) + (msq & 0xfffff);
515 lc[lsw] = lsq;
516 }
517 return dc;
518 }
519 float
__divsf3(float a,float b)520 __divsf3(float a, float b)
521 {
522 return (float)((double)a / (double)b);
523 }
524
525 /* ---------------- Comparisons ---------------- */
526
527 static int
compared2(const long * pa,const long * pb)528 compared2(const long *pa, const long *pb)
529 {
530 #define upa ((const ulong *)pa)
531 #define upb ((const ulong *)pb)
532 if (pa[msw] == pb[msw]) {
533 int result = (upa[lsw] < upb[lsw] ? -1 :
534 upa[lsw] > upb[lsw] ? 1 : 0);
535
536 return (pa[msw] < 0 ? -result : result);
537 }
538 if ((pa[msw] & pb[msw]) < 0)
539 return (pa[msw] < pb[msw] ? 1 : -1);
540 /* We have to treat -0 and +0 specially. */
541 else if (!((pa[msw] | pb[msw]) & 0x7fffffff) && !(pa[lsw] | pb[lsw]))
542 return 0;
543 else
544 return (pa[msw] > pb[msw] ? 1 : -1);
545 #undef upa
546 #undef upb
547 }
548
549 /* 0 means true */
550 int
__eqdf2(double a,double b)551 __eqdf2(double a, double b)
552 {
553 return compared2(la, lb);
554 }
555
556 /* !=0 means true */
557 int
__nedf2(double a,double b)558 __nedf2(double a, double b)
559 {
560 return compared2(la, lb);
561 }
562
563 /* >0 means true */
564 int
__gtdf2(double a,double b)565 __gtdf2(double a, double b)
566 {
567 return compared2(la, lb);
568 }
569
570 /* >=0 means true */
571 int
__gedf2(double a,double b)572 __gedf2(double a, double b)
573 {
574 return compared2(la, lb);
575 }
576
577 /* <0 means true */
578 int
__ltdf2(double a,double b)579 __ltdf2(double a, double b)
580 {
581 return compared2(la, lb);
582 }
583
584 /* <=0 means true */
585 int
__ledf2(double a,double b)586 __ledf2(double a, double b)
587 {
588 return compared2(la, lb);
589 }
590
591 static int
comparef2(long va,long vb)592 comparef2(long va, long vb)
593 { /* We have to treat -0 and +0 specially. */
594 if (va == vb)
595 return 0;
596 if ((va & vb) < 0)
597 return (va < vb ? 1 : -1);
598 else if (!((va | vb) & 0x7fffffff))
599 return 0;
600 else
601 return (va > vb ? 1 : -1);
602 }
603
604 /* See the __xxdf2 functions above for the return values. */
605 int
__eqsf2(float a,float b)606 __eqsf2(float a, float b)
607 {
608 return comparef2(ia, ib);
609 }
610 int
__nesf2(float a,float b)611 __nesf2(float a, float b)
612 {
613 return comparef2(ia, ib);
614 }
615 int
__gtsf2(float a,float b)616 __gtsf2(float a, float b)
617 {
618 return comparef2(ia, ib);
619 }
620 int
__gesf2(float a,float b)621 __gesf2(float a, float b)
622 {
623 return comparef2(ia, ib);
624 }
625 int
__ltsf2(float a,float b)626 __ltsf2(float a, float b)
627 {
628 return comparef2(ia, ib);
629 }
630 int
__lesf2(float a,float b)631 __lesf2(float a, float b)
632 {
633 return comparef2(ia, ib);
634 }
635
636 /* ---------------- Conversion ---------------- */
637
638 long
__fixdfsi(double a)639 __fixdfsi(double a)
640 { /* This routine does check for overflow. */
641 long i = (la[msw] & 0xfffff) + 0x100000;
642 int expt = dx(la) - dx_bias;
643
644 if (expt < 0)
645 return 0;
646 if (expt <= 20)
647 i >>= 20 - expt;
648 else if (expt >= 31 &&
649 (expt > 31 || i != 0x100000 || la[msw] >= 0 ||
650 ula[lsw] >= 1L << 21)
651 ) {
652 raise(SIGFPE);
653 i = (la[msw] < 0 ? 0x80000000 : 0x7fffffff);
654 } else
655 i = (i << (expt - 20)) + (ula[lsw] >> (52 - expt));
656 return (la[msw] < 0 ? -i : i);
657 }
658
659 long
__fixsfsi(float a)660 __fixsfsi(float a)
661 { /* This routine does check for overflow. */
662 long i = (ia & 0x7fffff) + 0x800000;
663 int expt = fx(ia) - fx_bias;
664
665 if (expt < 0)
666 return 0;
667 if (expt <= 23)
668 i >>= 23 - expt;
669 else if (expt >= 31 && (expt > 31 || i != 0x800000 || ia >= 0)) {
670 raise(SIGFPE);
671 i = (ia < 0 ? 0x80000000 : 0x7fffffff);
672 } else
673 i <<= expt - 23;
674 return (ia < 0 ? -i : i);
675 }
676
677 double
__floatsidf(long i)678 __floatsidf(long i)
679 {
680 long msc;
681 ulong v;
682 long lc[2];
683
684 if (i > 0)
685 msc = 0x41e00000 - 0x100000, v = i;
686 else if (i < 0)
687 msc = 0xc1e00000 - 0x100000, v = -i;
688 else /* i == 0 */
689 return 0;
690 while (v < 0x01000000)
691 v <<= 8, msc -= 0x00800000;
692 if (v < 0x10000000)
693 v <<= 4, msc -= 0x00400000;
694 while (v < 0x80000000)
695 v <<= 1, msc -= 0x00100000;
696 lc[msw] = msc + (v >> 11);
697 lc[lsw] = v << 21;
698 return dc;
699 }
700
701 float
__floatsisf(long i)702 __floatsisf(long i)
703 {
704 long lc;
705
706 if (i == 0)
707 lc = 0;
708 else {
709 ulong v;
710
711 if (i < 0)
712 lc = 0xcf000000, v = -i;
713 else
714 lc = 0x4f000000, v = i;
715 while (v < 0x01000000)
716 v <<= 8, lc -= 0x04000000;
717 while (v < 0x80000000)
718 v <<= 1, lc -= 0x00800000;
719 v = ((v >> 7) + 1) >> 1;
720 if (v > 0xffffff)
721 v >>= 1, lc += 0x00800000;
722 lc += v & 0x7fffff;
723 }
724 return fc;
725 }
726
727 float
__truncdfsf2(double a)728 __truncdfsf2(double a)
729 { /* This routine does check for overflow, but it converts */
730 /* underflows to zero rather than to a denormalized number. */
731 long lc;
732
733 if ((la[msw] & 0x7ff00000) < 0x38100000)
734 lc = la[msw] & 0x80000000;
735 else if ((la[msw] & 0x7ff00000) >= 0x47f00000) {
736 raise(SIGFPE);
737 lc = (la[msw] & 0x80000000) + 0x7f800000; /* infinity */
738 } else {
739 lc = (la[msw] & 0xc0000000) +
740 ((la[msw] & 0x07ffffff) << 3) +
741 (ula[lsw] >> 29);
742 /* Round the mantissa. Simply adding 1 works even if the */
743 /* mantissa overflows, because it increments the exponent and */
744 /* sets the mantissa to zero! */
745 if (ula[lsw] & 0x10000000)
746 ++lc;
747 }
748 return fc;
749 }
750
751 double
__extendsfdf2(float a)752 __extendsfdf2(float a)
753 {
754 long lc[2];
755
756 if (!(ia & 0x7fffffff))
757 lc[msw] = ia, lc[lsw] = 0;
758 else
759 extend(lc, ia);
760 return dc;
761 }
762
763 /* ---------------- Test program ---------------- */
764
765 #ifdef TEST
766
767 #include <stdio.h>
768 #include <stdlib.h>
769
770 int
test(double v1)771 test(double v1)
772 {
773 double v3 = v1 * 3;
774 double vh = v1 / 2;
775 double vd = v3 - vh;
776 double vdn = v1 - v3;
777
778 printf("%g=1 %g=3 %g=0.5 %g=2.5 %g=-2\n", v1, v3, vh, vd, vdn);
779 return 0;
780 }
781
782 float
randf(void)783 randf(void)
784 {
785 int v = rand();
786
787 v = (v << 16) ^ rand();
788 if (!(v & 0x7f800000))
789 return 0;
790 if ((v & 0x7f800000) == 0x7f800000)
791 return randf();
792 return *(float *)&v;
793 }
794
795 int
main(int argc,char * argv[])796 main(int argc, char *argv[])
797 {
798 int i;
799
800 test(1.0);
801 for (i = 0; i < 10; ++i) {
802 float a = randf(), b = randf(), r;
803 int c;
804
805 switch ((rand() >> 12) & 3) {
806 case 0:
807 r = a + b;
808 c = '+';
809 break;
810 case 1:
811 r = a - b;
812 c = '-';
813 break;
814 case 2:
815 r = a * b;
816 c = '*';
817 break;
818 case 3:
819 if (b == 0)
820 continue;
821 r = a / b;
822 c = '/';
823 break;
824 }
825 printf("0x%08x %c 0x%08x = 0x%08x\n",
826 *(int *)&a, c, *(int *)&b, *(int *)&r);
827 }
828 }
829
830 /*
831 Results from compiling with hardware floating point:
832
833 1=1 3=3 0.5=0.5 2.5=2.5 -2=-2
834 0x6f409924 - 0x01faa67a = 0x6f409924
835 0x00000000 + 0x75418107 = 0x75418107
836 0xe32fab71 - 0xc88f7816 = 0xe32fab71
837 0x94809067 * 0x84ddaeee = 0x00000000
838 0x2b0a5b7d + 0x38f70f50 = 0x38f70f50
839 0xa5efcef3 - 0xc5dc1a2c = 0x45dc1a2c
840 0xe7124521 * 0x3f4206d2 = 0xe6ddb891
841 0x8d175f17 * 0x2ed2c1c0 = 0x80007c9f
842 0x419e7a6d / 0xf2f95a35 = 0x8e22b404
843 0xe0b2f48f * 0xc72793fc = 0x686a49f8
844
845 */
846
847
848 #endif
849