1 /*
2 ** emfloat.c
3 ** Source for emulated floating-point routines.
4 ** BYTEmark (tm)
5 ** BYTE's Native Mode Benchmarks
6 ** Rick Grehan, BYTE Magazine.
7 **
8 ** Created:
9 ** Last update: 3/95
10 **
11 ** DISCLAIMER
12 ** The source, executable, and documentation files that comprise
13 ** the BYTEmark benchmarks are made available on an "as is" basis.
14 ** This means that we at BYTE Magazine have made every reasonable
15 ** effort to verify that the there are no errors in the source and
16 ** executable code. We cannot, however, guarantee that the programs
17 ** are error-free. Consequently, McGraw-HIll and BYTE Magazine make
18 ** no claims in regard to the fitness of the source code, executable
19 ** code, and documentation of the BYTEmark.
20 ** Furthermore, BYTE Magazine, McGraw-Hill, and all employees
21 ** of McGraw-Hill cannot be held responsible for any damages resulting
22 ** from the use of this code or the results obtained from using
23 ** this code.
24 */
25
26 #include "../pub/libvex_basictypes.h"
27
28 static HWord (*serviceFn)(HWord,HWord) = 0;
29
30
31 /////////////////////////////////////////////////////////////////////
32 /////////////////////////////////////////////////////////////////////
33
my_strcpy(char * dest,const char * src)34 static char* my_strcpy ( char* dest, const char* src )
35 {
36 char* dest_orig = dest;
37 while (*src) *dest++ = *src++;
38 *dest = 0;
39 return dest_orig;
40 }
41
my_memcpy(void * dest,const void * src,int sz)42 static void* my_memcpy ( void *dest, const void *src, int sz )
43 {
44 const char *s = (const char *)src;
45 char *d = (char *)dest;
46
47 while (sz--)
48 *d++ = *s++;
49
50 return dest;
51 }
52
my_memmove(void * dst,const void * src,unsigned int len)53 static void* my_memmove( void *dst, const void *src, unsigned int len )
54 {
55 register char *d;
56 register char *s;
57 if ( dst > src ) {
58 d = (char *)dst + len - 1;
59 s = (char *)src + len - 1;
60 while ( len >= 4 ) {
61 *d-- = *s--;
62 *d-- = *s--;
63 *d-- = *s--;
64 *d-- = *s--;
65 len -= 4;
66 }
67 while ( len-- ) {
68 *d-- = *s--;
69 }
70 } else if ( dst < src ) {
71 d = (char *)dst;
72 s = (char *)src;
73 while ( len >= 4 ) {
74 *d++ = *s++;
75 *d++ = *s++;
76 *d++ = *s++;
77 *d++ = *s++;
78 len -= 4;
79 }
80 while ( len-- ) {
81 *d++ = *s++;
82 }
83 }
84 return dst;
85 }
86
87 /////////////////////////////////////////////////////////////////////
88
vexxx_log_bytes(char * p,int n)89 static void vexxx_log_bytes ( char* p, int n )
90 {
91 int i;
92 for (i = 0; i < n; i++)
93 (*serviceFn)( 1, (int)p[i] );
94 }
95
96 /*---------------------------------------------------------*/
97 /*--- vexxx_printf ---*/
98 /*---------------------------------------------------------*/
99
100 /* This should be the only <...> include in the entire VEXXX library.
101 New code for vexxx_util.c should go above this point. */
102 #include <stdarg.h>
103
vexxx_toupper(HChar c)104 static HChar vexxx_toupper ( HChar c )
105 {
106 if (c >= 'a' && c <= 'z')
107 return toHChar(c + ('A' - 'a'));
108 else
109 return c;
110 }
111
vexxx_strlen(const HChar * str)112 static Int vexxx_strlen ( const HChar* str )
113 {
114 Int i = 0;
115 while (str[i] != 0) i++;
116 return i;
117 }
118
vexxx_streq(const HChar * s1,const HChar * s2)119 Bool vexxx_streq ( const HChar* s1, const HChar* s2 )
120 {
121 while (True) {
122 if (*s1 == 0 && *s2 == 0)
123 return True;
124 if (*s1 != *s2)
125 return False;
126 s1++;
127 s2++;
128 }
129 }
130
131 /* Some flags. */
132 #define VG_MSG_SIGNED 1 /* The value is signed. */
133 #define VG_MSG_ZJUSTIFY 2 /* Must justify with '0'. */
134 #define VG_MSG_LJUSTIFY 4 /* Must justify on the left. */
135 #define VG_MSG_PAREN 8 /* Parenthesize if present (for %y) */
136 #define VG_MSG_COMMA 16 /* Add commas to numbers (for %d, %u) */
137
138 /* Copy a string into the buffer. */
139 static UInt
myvprintf_str(void (* send)(HChar),Int flags,Int width,HChar * str,Bool capitalise)140 myvprintf_str ( void(*send)(HChar), Int flags, Int width, HChar* str,
141 Bool capitalise )
142 {
143 # define MAYBE_TOUPPER(ch) toHChar(capitalise ? vexxx_toupper(ch) : (ch))
144 UInt ret = 0;
145 Int i, extra;
146 Int len = vexxx_strlen(str);
147
148 if (width == 0) {
149 ret += len;
150 for (i = 0; i < len; i++)
151 send(MAYBE_TOUPPER(str[i]));
152 return ret;
153 }
154
155 if (len > width) {
156 ret += width;
157 for (i = 0; i < width; i++)
158 send(MAYBE_TOUPPER(str[i]));
159 return ret;
160 }
161
162 extra = width - len;
163 if (flags & VG_MSG_LJUSTIFY) {
164 ret += extra;
165 for (i = 0; i < extra; i++)
166 send(' ');
167 }
168 ret += len;
169 for (i = 0; i < len; i++)
170 send(MAYBE_TOUPPER(str[i]));
171 if (!(flags & VG_MSG_LJUSTIFY)) {
172 ret += extra;
173 for (i = 0; i < extra; i++)
174 send(' ');
175 }
176
177 # undef MAYBE_TOUPPER
178
179 return ret;
180 }
181
182 /* Write P into the buffer according to these args:
183 * If SIGN is true, p is a signed.
184 * BASE is the base.
185 * If WITH_ZERO is true, '0' must be added.
186 * WIDTH is the width of the field.
187 */
188 static UInt
myvprintf_int64(void (* send)(HChar),Int flags,Int base,Int width,ULong pL)189 myvprintf_int64 ( void(*send)(HChar), Int flags, Int base, Int width, ULong pL)
190 {
191 HChar buf[40];
192 Int ind = 0;
193 Int i, nc = 0;
194 Bool neg = False;
195 HChar *digits = "0123456789ABCDEF";
196 UInt ret = 0;
197 UInt p = (UInt)pL;
198
199 if (base < 2 || base > 16)
200 return ret;
201
202 if ((flags & VG_MSG_SIGNED) && (Int)p < 0) {
203 p = - (Int)p;
204 neg = True;
205 }
206
207 if (p == 0)
208 buf[ind++] = '0';
209 else {
210 while (p > 0) {
211 if ((flags & VG_MSG_COMMA) && 10 == base &&
212 0 == (ind-nc) % 3 && 0 != ind)
213 {
214 buf[ind++] = ',';
215 nc++;
216 }
217 buf[ind++] = digits[p % base];
218 p /= base;
219 }
220 }
221
222 if (neg)
223 buf[ind++] = '-';
224
225 if (width > 0 && !(flags & VG_MSG_LJUSTIFY)) {
226 for(; ind < width; ind++) {
227 //vassert(ind < 39);
228 buf[ind] = toHChar((flags & VG_MSG_ZJUSTIFY) ? '0': ' ');
229 }
230 }
231
232 /* Reverse copy to buffer. */
233 ret += ind;
234 for (i = ind -1; i >= 0; i--) {
235 send(buf[i]);
236 }
237 if (width > 0 && (flags & VG_MSG_LJUSTIFY)) {
238 for(; ind < width; ind++) {
239 ret++;
240 send(' '); // Never pad with zeroes on RHS -- changes the value!
241 }
242 }
243 return ret;
244 }
245
246
247 /* A simple vprintf(). */
248 static
vprintf_wrk(void (* send)(HChar),const HChar * format,va_list vargs)249 UInt vprintf_wrk ( void(*send)(HChar), const HChar *format, va_list vargs )
250 {
251 UInt ret = 0;
252 int i;
253 int flags;
254 int width;
255 Bool is_long;
256
257 /* We assume that vargs has already been initialised by the
258 caller, using va_start, and that the caller will similarly
259 clean up with va_end.
260 */
261
262 for (i = 0; format[i] != 0; i++) {
263 if (format[i] != '%') {
264 send(format[i]);
265 ret++;
266 continue;
267 }
268 i++;
269 /* A '%' has been found. Ignore a trailing %. */
270 if (format[i] == 0)
271 break;
272 if (format[i] == '%') {
273 /* `%%' is replaced by `%'. */
274 send('%');
275 ret++;
276 continue;
277 }
278 flags = 0;
279 is_long = False;
280 width = 0; /* length of the field. */
281 if (format[i] == '(') {
282 flags |= VG_MSG_PAREN;
283 i++;
284 }
285 /* If ',' follows '%', commas will be inserted. */
286 if (format[i] == ',') {
287 flags |= VG_MSG_COMMA;
288 i++;
289 }
290 /* If '-' follows '%', justify on the left. */
291 if (format[i] == '-') {
292 flags |= VG_MSG_LJUSTIFY;
293 i++;
294 }
295 /* If '0' follows '%', pads will be inserted. */
296 if (format[i] == '0') {
297 flags |= VG_MSG_ZJUSTIFY;
298 i++;
299 }
300 /* Compute the field length. */
301 while (format[i] >= '0' && format[i] <= '9') {
302 width *= 10;
303 width += format[i++] - '0';
304 }
305 while (format[i] == 'l') {
306 i++;
307 is_long = True;
308 }
309
310 switch (format[i]) {
311 case 'd': /* %d */
312 flags |= VG_MSG_SIGNED;
313 if (is_long)
314 ret += myvprintf_int64(send, flags, 10, width,
315 (ULong)(va_arg (vargs, Long)));
316 else
317 ret += myvprintf_int64(send, flags, 10, width,
318 (ULong)(va_arg (vargs, Int)));
319 break;
320 case 'u': /* %u */
321 if (is_long)
322 ret += myvprintf_int64(send, flags, 10, width,
323 (ULong)(va_arg (vargs, ULong)));
324 else
325 ret += myvprintf_int64(send, flags, 10, width,
326 (ULong)(va_arg (vargs, UInt)));
327 break;
328 case 'p': /* %p */
329 ret += 2;
330 send('0');
331 send('x');
332 ret += myvprintf_int64(send, flags, 16, width,
333 (ULong)((HWord)va_arg (vargs, void *)));
334 break;
335 case 'x': /* %x */
336 if (is_long)
337 ret += myvprintf_int64(send, flags, 16, width,
338 (ULong)(va_arg (vargs, ULong)));
339 else
340 ret += myvprintf_int64(send, flags, 16, width,
341 (ULong)(va_arg (vargs, UInt)));
342 break;
343 case 'c': /* %c */
344 ret++;
345 send(toHChar(va_arg (vargs, int)));
346 break;
347 case 's': case 'S': { /* %s */
348 char *str = va_arg (vargs, char *);
349 if (str == (char*) 0) str = "(null)";
350 ret += myvprintf_str(send, flags, width, str,
351 toBool(format[i]=='S'));
352 break;
353 }
354 # if 0
355 case 'y': { /* %y - print symbol */
356 Addr a = va_arg(vargs, Addr);
357
358 HChar *name;
359 if (VG_(get_fnname_w_offset)(a, &name)) {
360 HChar buf[1 + VG_strlen(name) + 1 + 1];
361 if (flags & VG_MSG_PAREN) {
362 VG_(sprintf)(str, "(%s)", name):
363 } else {
364 VG_(sprintf)(str, "%s", name):
365 }
366 ret += myvprintf_str(send, flags, width, buf, 0);
367 }
368 break;
369 }
370 # endif
371 default:
372 break;
373 }
374 }
375 return ret;
376 }
377
378
379 /* A general replacement for printf(). Note that only low-level
380 debugging info should be sent via here. The official route is to
381 to use vg_message(). This interface is deprecated.
382 */
383 static HChar myprintf_buf[1000];
384 static Int n_myprintf_buf;
385
add_to_myprintf_buf(HChar c)386 static void add_to_myprintf_buf ( HChar c )
387 {
388 if (c == '\n' || n_myprintf_buf >= 1000-10 /*paranoia*/ ) {
389 (*vexxx_log_bytes)( myprintf_buf, vexxx_strlen(myprintf_buf) );
390 n_myprintf_buf = 0;
391 myprintf_buf[n_myprintf_buf] = 0;
392 }
393 myprintf_buf[n_myprintf_buf++] = c;
394 myprintf_buf[n_myprintf_buf] = 0;
395 }
396
vexxx_printf(const char * format,...)397 static UInt vexxx_printf ( const char *format, ... )
398 {
399 UInt ret;
400 va_list vargs;
401 va_start(vargs,format);
402
403 n_myprintf_buf = 0;
404 myprintf_buf[n_myprintf_buf] = 0;
405 ret = vprintf_wrk ( add_to_myprintf_buf, format, vargs );
406
407 if (n_myprintf_buf > 0) {
408 (*vexxx_log_bytes)( myprintf_buf, n_myprintf_buf );
409 }
410
411 va_end(vargs);
412
413 return ret;
414 }
415
416 /*---------------------------------------------------------------*/
417 /*--- end vexxx_util.c ---*/
418 /*---------------------------------------------------------------*/
419
420
421 /////////////////////////////////////////////////////////////////////
422 /////////////////////////////////////////////////////////////////////
423
424 //#include <stdio.h>
425 //#include <string.h>
426 //#include <malloc.h>
427
428 typedef unsigned char uchar;
429 typedef unsigned int uint;
430 typedef unsigned short ushort;
431 typedef unsigned long ulong;
432 typedef int int32; /* Signed 32 bit integer */
433
434 #define INTERNAL_FPF_PRECISION 4
435 #define CPUEMFLOATLOOPMAX 500000L
436 #define EMFARRAYSIZE 3000L
437
438 typedef struct {
439 int adjust; /* Set adjust code */
440 ulong request_secs; /* # of seconds requested */
441 ulong arraysize; /* Size of array */
442 ulong loops; /* Loops per iterations */
443 double emflops; /* Results */
444 } EmFloatStruct;
445
446
447
448 /* Is this a 64 bit architecture? If so, this will define LONG64 */
449 /* Uwe F. Mayer 15 November 1997 */
450 // #include "pointer.h"
451
452 #define u8 unsigned char
453 #define u16 unsigned short
454 #ifdef LONG64
455 #define u32 unsigned int
456 #else
457 #define u32 unsigned long
458 #endif
459 #define uchar unsigned char
460 #define ulong unsigned long
461
462 #define MAX_EXP 32767L
463 #define MIN_EXP (-32767L)
464
465 #define IFPF_IS_ZERO 0
466 #define IFPF_IS_SUBNORMAL 1
467 #define IFPF_IS_NORMAL 2
468 #define IFPF_IS_INFINITY 3
469 #define IFPF_IS_NAN 4
470 #define IFPF_TYPE_COUNT 5
471
472 #define ZERO_ZERO 0
473 #define ZERO_SUBNORMAL 1
474 #define ZERO_NORMAL 2
475 #define ZERO_INFINITY 3
476 #define ZERO_NAN 4
477
478 #define SUBNORMAL_ZERO 5
479 #define SUBNORMAL_SUBNORMAL 6
480 #define SUBNORMAL_NORMAL 7
481 #define SUBNORMAL_INFINITY 8
482 #define SUBNORMAL_NAN 9
483
484 #define NORMAL_ZERO 10
485 #define NORMAL_SUBNORMAL 11
486 #define NORMAL_NORMAL 12
487 #define NORMAL_INFINITY 13
488 #define NORMAL_NAN 14
489
490 #define INFINITY_ZERO 15
491 #define INFINITY_SUBNORMAL 16
492 #define INFINITY_NORMAL 17
493 #define INFINITY_INFINITY 18
494 #define INFINITY_NAN 19
495
496 #define NAN_ZERO 20
497 #define NAN_SUBNORMAL 21
498 #define NAN_NORMAL 22
499 #define NAN_INFINITY 23
500 #define NAN_NAN 24
501 #define OPERAND_ZERO 0
502 #define OPERAND_SUBNORMAL 1
503 #define OPERAND_NORMAL 2
504 #define OPERAND_INFINITY 3
505 #define OPERAND_NAN 4
506
507 typedef struct
508 {
509 u8 type; /* Indicates, NORMAL, SUBNORMAL, etc. */
510 u8 sign; /* Mantissa sign */
511 short exp; /* Signed exponent...no bias */
512 u16 mantissa[INTERNAL_FPF_PRECISION];
513 } InternalFPF;
514
515 static
516 void SetupCPUEmFloatArrays(InternalFPF *abase,
517 InternalFPF *bbase, InternalFPF *cbase, ulong arraysize);
518 static
519 ulong DoEmFloatIteration(InternalFPF *abase,
520 InternalFPF *bbase, InternalFPF *cbase,
521 ulong arraysize, ulong loops);
522
523 static void SetInternalFPFZero(InternalFPF *dest,
524 uchar sign);
525 static void SetInternalFPFInfinity(InternalFPF *dest,
526 uchar sign);
527 static void SetInternalFPFNaN(InternalFPF *dest);
528 static int IsMantissaZero(u16 *mant);
529 static void Add16Bits(u16 *carry,u16 *a,u16 b,u16 c);
530 static void Sub16Bits(u16 *borrow,u16 *a,u16 b,u16 c);
531 static void ShiftMantLeft1(u16 *carry,u16 *mantissa);
532 static void ShiftMantRight1(u16 *carry,u16 *mantissa);
533 static void StickyShiftRightMant(InternalFPF *ptr,int amount);
534 static void normalize(InternalFPF *ptr);
535 static void denormalize(InternalFPF *ptr,int minimum_exponent);
536 static void RoundInternalFPF(InternalFPF *ptr);
537 static void choose_nan(InternalFPF *x,InternalFPF *y,InternalFPF *z,
538 int intel_flag);
539 static void AddSubInternalFPF(uchar operation,InternalFPF *x,
540 InternalFPF *y,InternalFPF *z);
541 static void MultiplyInternalFPF(InternalFPF *x,InternalFPF *y,
542 InternalFPF *z);
543 static void DivideInternalFPF(InternalFPF *x,InternalFPF *y,
544 InternalFPF *z);
545
546 static void Int32ToInternalFPF(int32 mylong,
547 InternalFPF *dest);
548 static int InternalFPFToString(char *dest,
549 InternalFPF *src);
550
551 static int32 randnum(int32 lngval);
552
randwc(int32 num)553 static int32 randwc(int32 num)
554 {
555 return(randnum((int32)0)%num);
556 }
557
558 static int32 randw[2] = { (int32)13 , (int32)117 };
randnum(int32 lngval)559 static int32 randnum(int32 lngval)
560 {
561 register int32 interm;
562
563 if (lngval!=(int32)0)
564 { randw[0]=(int32)13; randw[1]=(int32)117; }
565
566 interm=(randw[0]*(int32)254754+randw[1]*(int32)529562)%(int32)999563;
567 randw[1]=randw[0];
568 randw[0]=interm;
569 return(interm);
570 }
571
572
573 static
SetupCPUEmFloatArrays(InternalFPF * abase,InternalFPF * bbase,InternalFPF * cbase,ulong arraysize)574 void SetupCPUEmFloatArrays(InternalFPF *abase,
575 InternalFPF *bbase,
576 InternalFPF *cbase,
577 ulong arraysize)
578 {
579 ulong i;
580 InternalFPF locFPF1,locFPF2;
581
582 randnum((int32)13);
583
584 for(i=0;i<arraysize;i++)
585 {/* LongToInternalFPF(randwc(50000L),&locFPF1); */
586 Int32ToInternalFPF(randwc((int32)50000),&locFPF1);
587 /* LongToInternalFPF(randwc(50000L)+1L,&locFPF2); */
588 Int32ToInternalFPF(randwc((int32)50000)+(int32)1,&locFPF2);
589 DivideInternalFPF(&locFPF1,&locFPF2,abase+i);
590 /* LongToInternalFPF(randwc(50000L)+1L,&locFPF2); */
591 Int32ToInternalFPF(randwc((int32)50000)+(int32)1,&locFPF2);
592 DivideInternalFPF(&locFPF1,&locFPF2,bbase+i);
593 }
594 return;
595 }
596
597
598 static char* str1 = "loops %d\n";
599 static
DoEmFloatIteration(InternalFPF * abase,InternalFPF * bbase,InternalFPF * cbase,ulong arraysize,ulong loops)600 ulong DoEmFloatIteration(InternalFPF *abase,
601 InternalFPF *bbase,
602 InternalFPF *cbase,
603 ulong arraysize, ulong loops)
604 {
605 static uchar jtable[16] = {0,0,0,0,1,1,1,1,2,2,2,2,2,3,3,3};
606 ulong i;
607 int number_of_loops;
608 loops = 100;
609 number_of_loops=loops-1; /* the index of the first loop we run */
610
611 vexxx_printf(str1, (int)loops);
612
613 /*
614 ** Each pass through the array performs operations in
615 ** the followingratios:
616 ** 4 adds, 4 subtracts, 5 multiplies, 3 divides
617 ** (adds and subtracts being nearly the same operation)
618 */
619
620 {
621 for(i=0;i<arraysize;i++)
622 switch(jtable[i % 16])
623 {
624 case 0: /* Add */
625 AddSubInternalFPF(0,abase+i,
626 bbase+i,
627 cbase+i);
628 break;
629 case 1: /* Subtract */
630 AddSubInternalFPF(1,abase+i,
631 bbase+i,
632 cbase+i);
633 break;
634 case 2: /* Multiply */
635 MultiplyInternalFPF(abase+i,
636 bbase+i,
637 cbase+i);
638 break;
639 case 3: /* Divide */
640 DivideInternalFPF(abase+i,
641 bbase+i,
642 cbase+i);
643 break;
644 }
645 {
646 ulong j[8]; /* we test 8 entries */
647 int k;
648 ulong i;
649 char buffer[1024];
650 if (100==loops) /* the first loop */
651 {
652 j[0]=(ulong)2;
653 j[1]=(ulong)6;
654 j[2]=(ulong)10;
655 j[3]=(ulong)14;
656 j[4]=(ulong)(arraysize-14);
657 j[5]=(ulong)(arraysize-10);
658 j[6]=(ulong)(arraysize-6);
659 j[7]=(ulong)(arraysize-2);
660 for(k=0;k<8;k++){
661 i=j[k];
662 InternalFPFToString(buffer,abase+i);
663 vexxx_printf("%6d: (%s) ",i,buffer);
664 switch(jtable[i % 16])
665 {
666 case 0: my_strcpy(buffer,"+"); break;
667 case 1: my_strcpy(buffer,"-"); break;
668 case 2: my_strcpy(buffer,"*"); break;
669 case 3: my_strcpy(buffer,"/"); break;
670 }
671 vexxx_printf("%s ",buffer);
672 InternalFPFToString(buffer,bbase+i);
673 vexxx_printf("(%s) = ",buffer);
674 InternalFPFToString(buffer,cbase+i);
675 vexxx_printf("%s\n",buffer);
676 }
677 return 0;
678 }
679 }
680 }
681 return 0;
682 }
683
684 /***********************
685 ** SetInternalFPFZero **
686 ************************
687 ** Set an internal floating-point-format number to zero.
688 ** sign determines the sign of the zero.
689 */
SetInternalFPFZero(InternalFPF * dest,uchar sign)690 static void SetInternalFPFZero(InternalFPF *dest,
691 uchar sign)
692 {
693 int i; /* Index */
694
695 dest->type=IFPF_IS_ZERO;
696 dest->sign=sign;
697 dest->exp=MIN_EXP;
698 for(i=0;i<INTERNAL_FPF_PRECISION;i++)
699 dest->mantissa[i]=0;
700 return;
701 }
702
703 /***************************
704 ** SetInternalFPFInfinity **
705 ****************************
706 ** Set an internal floating-point-format number to infinity.
707 ** This can happen if the exponent exceeds MAX_EXP.
708 ** As above, sign picks the sign of infinity.
709 */
SetInternalFPFInfinity(InternalFPF * dest,uchar sign)710 static void SetInternalFPFInfinity(InternalFPF *dest,
711 uchar sign)
712 {
713 int i; /* Index */
714
715 dest->type=IFPF_IS_INFINITY;
716 dest->sign=sign;
717 dest->exp=MIN_EXP;
718 for(i=0;i<INTERNAL_FPF_PRECISION;i++)
719 dest->mantissa[i]=0;
720 return;
721 }
722
723 /**********************
724 ** SetInternalFPFNaN **
725 ***********************
726 ** Set an internal floating-point-format number to Nan
727 ** (not a number). Note that we "emulate" an 80x87 as far
728 ** as the mantissa bits go.
729 */
SetInternalFPFNaN(InternalFPF * dest)730 static void SetInternalFPFNaN(InternalFPF *dest)
731 {
732 int i; /* Index */
733
734 dest->type=IFPF_IS_NAN;
735 dest->exp=MAX_EXP;
736 dest->sign=1;
737 dest->mantissa[0]=0x4000;
738 for(i=1;i<INTERNAL_FPF_PRECISION;i++)
739 dest->mantissa[i]=0;
740
741 return;
742 }
743
744 /*******************
745 ** IsMantissaZero **
746 ********************
747 ** Pass this routine a pointer to an internal floating point format
748 ** number's mantissa. It checks for an all-zero mantissa.
749 ** Returns 0 if it is NOT all zeros, !=0 otherwise.
750 */
IsMantissaZero(u16 * mant)751 static int IsMantissaZero(u16 *mant)
752 {
753 int i; /* Index */
754 int n; /* Return value */
755
756 n=0;
757 for(i=0;i<INTERNAL_FPF_PRECISION;i++)
758 n|=mant[i];
759
760 return(!n);
761 }
762
763 /**************
764 ** Add16Bits **
765 ***************
766 ** Add b, c, and carry. Retult in a. New carry in carry.
767 */
Add16Bits(u16 * carry,u16 * a,u16 b,u16 c)768 static void Add16Bits(u16 *carry,
769 u16 *a,
770 u16 b,
771 u16 c)
772 {
773 u32 accum; /* Accumulator */
774
775 /*
776 ** Do the work in the 32-bit accumulator so we can return
777 ** the carry.
778 */
779 accum=(u32)b;
780 accum+=(u32)c;
781 accum+=(u32)*carry;
782 *carry=(u16)((accum & 0x00010000) ? 1 : 0); /* New carry */
783 *a=(u16)(accum & 0xFFFF); /* Result is lo 16 bits */
784 return;
785 }
786
787 /**************
788 ** Sub16Bits **
789 ***************
790 ** Additive inverse of above.
791 */
Sub16Bits(u16 * borrow,u16 * a,u16 b,u16 c)792 static void Sub16Bits(u16 *borrow,
793 u16 *a,
794 u16 b,
795 u16 c)
796 {
797 u32 accum; /* Accumulator */
798
799 accum=(u32)b;
800 accum-=(u32)c;
801 accum-=(u32)*borrow;
802 *borrow=(u32)((accum & 0x00010000) ? 1 : 0); /* New borrow */
803 *a=(u16)(accum & 0xFFFF);
804 return;
805 }
806
807 /*******************
808 ** ShiftMantLeft1 **
809 ********************
810 ** Shift a vector of 16-bit numbers left 1 bit. Also provides
811 ** a carry bit, which is shifted in at the beginning, and
812 ** shifted out at the end.
813 */
ShiftMantLeft1(u16 * carry,u16 * mantissa)814 static void ShiftMantLeft1(u16 *carry,
815 u16 *mantissa)
816 {
817 int i; /* Index */
818 int new_carry;
819 u16 accum; /* Temporary holding placed */
820
821 for(i=INTERNAL_FPF_PRECISION-1;i>=0;i--)
822 { accum=mantissa[i];
823 new_carry=accum & 0x8000; /* Get new carry */
824 accum=accum<<1; /* Do the shift */
825 if(*carry)
826 accum|=1; /* Insert previous carry */
827 *carry=new_carry;
828 mantissa[i]=accum; /* Return shifted value */
829 }
830 return;
831 }
832
833 /********************
834 ** ShiftMantRight1 **
835 *********************
836 ** Shift a mantissa right by 1 bit. Provides carry, as
837 ** above
838 */
ShiftMantRight1(u16 * carry,u16 * mantissa)839 static void ShiftMantRight1(u16 *carry,
840 u16 *mantissa)
841 {
842 int i; /* Index */
843 int new_carry;
844 u16 accum;
845
846 for(i=0;i<INTERNAL_FPF_PRECISION;i++)
847 { accum=mantissa[i];
848 new_carry=accum & 1; /* Get new carry */
849 accum=accum>>1;
850 if(*carry)
851 accum|=0x8000;
852 *carry=new_carry;
853 mantissa[i]=accum;
854 }
855 return;
856 }
857
858
859 /*****************************
860 ** StickyShiftMantRight **
861 ******************************
862 ** This is a shift right of the mantissa with a "sticky bit".
863 ** I.E., if a carry of 1 is shifted out of the least significant
864 ** bit, the least significant bit is set to 1.
865 */
StickyShiftRightMant(InternalFPF * ptr,int amount)866 static void StickyShiftRightMant(InternalFPF *ptr,
867 int amount)
868 {
869 int i; /* Index */
870 u16 carry; /* Self-explanatory */
871 u16 *mantissa;
872
873 mantissa=ptr->mantissa;
874
875 if(ptr->type!=IFPF_IS_ZERO) /* Don't bother shifting a zero */
876 {
877 /*
878 ** If the amount of shifting will shift everyting
879 ** out of existence, then just clear the whole mantissa
880 ** and set the lowmost bit to 1.
881 */
882 if(amount>=INTERNAL_FPF_PRECISION * 16)
883 {
884 for(i=0;i<INTERNAL_FPF_PRECISION-1;i++)
885 mantissa[i]=0;
886 mantissa[INTERNAL_FPF_PRECISION-1]=1;
887 }
888 else
889 for(i=0;i<amount;i++)
890 {
891 carry=0;
892 ShiftMantRight1(&carry,mantissa);
893 if(carry)
894 mantissa[INTERNAL_FPF_PRECISION-1] |= 1;
895 }
896 }
897 return;
898 }
899
900
901 /**************************************************
902 ** POST ARITHMETIC PROCESSING **
903 ** (NORMALIZE, ROUND, OVERFLOW, AND UNDERFLOW) **
904 **************************************************/
905
906 /**************
907 ** normalize **
908 ***************
909 ** Normalize an internal-representation number. Normalization
910 ** discards empty most-significant bits.
911 */
normalize(InternalFPF * ptr)912 static void normalize(InternalFPF *ptr)
913 {
914 u16 carry;
915
916 /*
917 ** As long as there's a highmost 0 bit, shift the significand
918 ** left 1 bit. Each time you do this, though, you've
919 ** gotta decrement the exponent.
920 */
921 while ((ptr->mantissa[0] & 0x8000) == 0)
922 {
923 carry = 0;
924 ShiftMantLeft1(&carry, ptr->mantissa);
925 ptr->exp--;
926 }
927 return;
928 }
929
930 /****************
931 ** denormalize **
932 *****************
933 ** Denormalize an internal-representation number. This means
934 ** shifting it right until its exponent is equivalent to
935 ** minimum_exponent. (You have to do this often in order
936 ** to perform additions and subtractions).
937 */
denormalize(InternalFPF * ptr,int minimum_exponent)938 static void denormalize(InternalFPF *ptr,
939 int minimum_exponent)
940 {
941 long exponent_difference;
942
943 if (IsMantissaZero(ptr->mantissa))
944 {
945 vexxx_printf("Error: zero significand in denormalize\n");
946 }
947
948 exponent_difference = ptr->exp-minimum_exponent;
949 if (exponent_difference < 0)
950 {
951 /*
952 ** The number is subnormal
953 */
954 exponent_difference = -exponent_difference;
955 if (exponent_difference >= (INTERNAL_FPF_PRECISION * 16))
956 {
957 /* Underflow */
958 SetInternalFPFZero(ptr, ptr->sign);
959 }
960 else
961 {
962 ptr->exp+=exponent_difference;
963 StickyShiftRightMant(ptr, exponent_difference);
964 }
965 }
966 return;
967 }
968
969
970 /*********************
971 ** RoundInternalFPF **
972 **********************
973 ** Round an internal-representation number.
974 ** The kind of rounding we do here is simplest...referred to as
975 ** "chop". "Extraneous" rightmost bits are simply hacked off.
976 */
RoundInternalFPF(InternalFPF * ptr)977 void RoundInternalFPF(InternalFPF *ptr)
978 {
979 /* int i; */
980
981 if (ptr->type == IFPF_IS_NORMAL ||
982 ptr->type == IFPF_IS_SUBNORMAL)
983 {
984 denormalize(ptr, MIN_EXP);
985 if (ptr->type != IFPF_IS_ZERO)
986 {
987
988 /* clear the extraneous bits */
989 ptr->mantissa[3] &= 0xfff8;
990 /* for (i=4; i<INTERNAL_FPF_PRECISION; i++)
991 {
992 ptr->mantissa[i] = 0;
993 }
994 */
995 /*
996 ** Check for overflow
997 */
998 /* Does not do anything as ptr->exp is a short and MAX_EXP=37268
999 if (ptr->exp > MAX_EXP)
1000 {
1001 SetInternalFPFInfinity(ptr, ptr->sign);
1002 }
1003 */
1004 }
1005 }
1006 return;
1007 }
1008
1009 /*******************************************************
1010 ** ARITHMETIC OPERATIONS ON INTERNAL REPRESENTATION **
1011 *******************************************************/
1012
1013 /***************
1014 ** choose_nan **
1015 ****************
1016 ** Called by routines that are forced to perform math on
1017 ** a pair of NaN's. This routine "selects" which NaN is
1018 ** to be returned.
1019 */
choose_nan(InternalFPF * x,InternalFPF * y,InternalFPF * z,int intel_flag)1020 static void choose_nan(InternalFPF *x,
1021 InternalFPF *y,
1022 InternalFPF *z,
1023 int intel_flag)
1024 {
1025 int i;
1026
1027 /*
1028 ** Compare the two mantissas,
1029 ** return the larger. Note that we will be emulating
1030 ** an 80387 in this operation.
1031 */
1032 for (i=0; i<INTERNAL_FPF_PRECISION; i++)
1033 {
1034 if (x->mantissa[i] > y->mantissa[i])
1035 {
1036 my_memmove((void *)x,(void *)z,sizeof(InternalFPF));
1037 return;
1038 }
1039 if (x->mantissa[i] < y->mantissa[i])
1040 {
1041 my_memmove((void *)y,(void *)z,sizeof(InternalFPF));
1042 return;
1043 }
1044 }
1045
1046 /*
1047 ** They are equal
1048 */
1049 if (!intel_flag)
1050 /* if the operation is addition */
1051 my_memmove((void *)x,(void *)z,sizeof(InternalFPF));
1052 else
1053 /* if the operation is multiplication */
1054 my_memmove((void *)y,(void *)z,sizeof(InternalFPF));
1055 return;
1056 }
1057
1058
1059 /**********************
1060 ** AddSubInternalFPF **
1061 ***********************
1062 ** Adding or subtracting internal-representation numbers.
1063 ** Internal-representation numbers pointed to by x and y are
1064 ** added/subtracted and the result returned in z.
1065 */
AddSubInternalFPF(uchar operation,InternalFPF * x,InternalFPF * y,InternalFPF * z)1066 static void AddSubInternalFPF(uchar operation,
1067 InternalFPF *x,
1068 InternalFPF *y,
1069 InternalFPF *z)
1070 {
1071 int exponent_difference;
1072 u16 borrow;
1073 u16 carry;
1074 int i;
1075 InternalFPF locx,locy; /* Needed since we alter them */
1076
1077 /*
1078 ** Following big switch statement handles the
1079 ** various combinations of operand types.
1080 */
1081 switch ((x->type * IFPF_TYPE_COUNT) + y->type)
1082 {
1083 case ZERO_ZERO:
1084 my_memmove((void *)x,(void *)z,sizeof(InternalFPF));
1085 if (x->sign ^ y->sign ^ operation)
1086 {
1087 z->sign = 0; /* positive */
1088 }
1089 break;
1090
1091 case NAN_ZERO:
1092 case NAN_SUBNORMAL:
1093 case NAN_NORMAL:
1094 case NAN_INFINITY:
1095 case SUBNORMAL_ZERO:
1096 case NORMAL_ZERO:
1097 case INFINITY_ZERO:
1098 case INFINITY_SUBNORMAL:
1099 case INFINITY_NORMAL:
1100 my_memmove((void *)x,(void *)z,sizeof(InternalFPF));
1101 break;
1102
1103
1104 case ZERO_NAN:
1105 case SUBNORMAL_NAN:
1106 case NORMAL_NAN:
1107 case INFINITY_NAN:
1108 my_memmove((void *)y,(void *)z,sizeof(InternalFPF));
1109 break;
1110
1111 case ZERO_SUBNORMAL:
1112 case ZERO_NORMAL:
1113 case ZERO_INFINITY:
1114 case SUBNORMAL_INFINITY:
1115 case NORMAL_INFINITY:
1116 my_memmove((void *)y,(void *)z,sizeof(InternalFPF));
1117 z->sign ^= operation;
1118 break;
1119
1120 case SUBNORMAL_SUBNORMAL:
1121 case SUBNORMAL_NORMAL:
1122 case NORMAL_SUBNORMAL:
1123 case NORMAL_NORMAL:
1124 /*
1125 ** Copy x and y to locals, since we may have
1126 ** to alter them.
1127 */
1128 my_memmove((void *)&locx,(void *)x,sizeof(InternalFPF));
1129 my_memmove((void *)&locy,(void *)y,sizeof(InternalFPF));
1130
1131 /* compute sum/difference */
1132 exponent_difference = locx.exp-locy.exp;
1133 if (exponent_difference == 0)
1134 {
1135 /*
1136 ** locx.exp == locy.exp
1137 ** so, no shifting required
1138 */
1139 if (locx.type == IFPF_IS_SUBNORMAL ||
1140 locy.type == IFPF_IS_SUBNORMAL)
1141 z->type = IFPF_IS_SUBNORMAL;
1142 else
1143 z->type = IFPF_IS_NORMAL;
1144
1145 /*
1146 ** Assume that locx.mantissa > locy.mantissa
1147 */
1148 z->sign = locx.sign;
1149 z->exp= locx.exp;
1150 }
1151 else
1152 if (exponent_difference > 0)
1153 {
1154 /*
1155 ** locx.exp > locy.exp
1156 */
1157 StickyShiftRightMant(&locy,
1158 exponent_difference);
1159 z->type = locx.type;
1160 z->sign = locx.sign;
1161 z->exp = locx.exp;
1162 }
1163 else /* if (exponent_difference < 0) */
1164 {
1165 /*
1166 ** locx.exp < locy.exp
1167 */
1168 StickyShiftRightMant(&locx,
1169 -exponent_difference);
1170 z->type = locy.type;
1171 z->sign = locy.sign ^ operation;
1172 z->exp = locy.exp;
1173 }
1174
1175 if (locx.sign ^ locy.sign ^ operation)
1176 {
1177 /*
1178 ** Signs are different, subtract mantissas
1179 */
1180 borrow = 0;
1181 for (i=(INTERNAL_FPF_PRECISION-1); i>=0; i--)
1182 Sub16Bits(&borrow,
1183 &z->mantissa[i],
1184 locx.mantissa[i],
1185 locy.mantissa[i]);
1186
1187 if (borrow)
1188 {
1189 /* The y->mantissa was larger than the
1190 ** x->mantissa leaving a negative
1191 ** result. Change the result back to
1192 ** an unsigned number and flip the
1193 ** sign flag.
1194 */
1195 z->sign = locy.sign ^ operation;
1196 borrow = 0;
1197 for (i=(INTERNAL_FPF_PRECISION-1); i>=0; i--)
1198 {
1199 Sub16Bits(&borrow,
1200 &z->mantissa[i],
1201 0,
1202 z->mantissa[i]);
1203 }
1204 }
1205 else
1206 {
1207 /* The assumption made above
1208 ** (i.e. x->mantissa >= y->mantissa)
1209 ** was correct. Therefore, do nothing.
1210 ** z->sign = x->sign;
1211 */
1212 }
1213
1214 if (IsMantissaZero(z->mantissa))
1215 {
1216 z->type = IFPF_IS_ZERO;
1217 z->sign = 0; /* positive */
1218 }
1219 else
1220 if (locx.type == IFPF_IS_NORMAL ||
1221 locy.type == IFPF_IS_NORMAL)
1222 {
1223 normalize(z);
1224 }
1225 }
1226 else
1227 {
1228 /* signs are the same, add mantissas */
1229 carry = 0;
1230 for (i=(INTERNAL_FPF_PRECISION-1); i>=0; i--)
1231 {
1232 Add16Bits(&carry,
1233 &z->mantissa[i],
1234 locx.mantissa[i],
1235 locy.mantissa[i]);
1236 }
1237
1238 if (carry)
1239 {
1240 z->exp++;
1241 carry=0;
1242 ShiftMantRight1(&carry,z->mantissa);
1243 z->mantissa[0] |= 0x8000;
1244 z->type = IFPF_IS_NORMAL;
1245 }
1246 else
1247 if (z->mantissa[0] & 0x8000)
1248 z->type = IFPF_IS_NORMAL;
1249 }
1250 break;
1251
1252 case INFINITY_INFINITY:
1253 SetInternalFPFNaN(z);
1254 break;
1255
1256 case NAN_NAN:
1257 choose_nan(x, y, z, 1);
1258 break;
1259 }
1260
1261 /*
1262 ** All the math is done; time to round.
1263 */
1264 RoundInternalFPF(z);
1265 return;
1266 }
1267
1268
1269 /************************
1270 ** MultiplyInternalFPF **
1271 *************************
1272 ** Two internal-representation numbers x and y are multiplied; the
1273 ** result is returned in z.
1274 */
MultiplyInternalFPF(InternalFPF * x,InternalFPF * y,InternalFPF * z)1275 static void MultiplyInternalFPF(InternalFPF *x,
1276 InternalFPF *y,
1277 InternalFPF *z)
1278 {
1279 int i;
1280 int j;
1281 u16 carry;
1282 u16 extra_bits[INTERNAL_FPF_PRECISION];
1283 InternalFPF locy; /* Needed since this will be altered */
1284 /*
1285 ** As in the preceding function, this large switch
1286 ** statement selects among the many combinations
1287 ** of operands.
1288 */
1289 switch ((x->type * IFPF_TYPE_COUNT) + y->type)
1290 {
1291 case INFINITY_SUBNORMAL:
1292 case INFINITY_NORMAL:
1293 case INFINITY_INFINITY:
1294 case ZERO_ZERO:
1295 case ZERO_SUBNORMAL:
1296 case ZERO_NORMAL:
1297 my_memmove((void *)x,(void *)z,sizeof(InternalFPF));
1298 z->sign ^= y->sign;
1299 break;
1300
1301 case SUBNORMAL_INFINITY:
1302 case NORMAL_INFINITY:
1303 case SUBNORMAL_ZERO:
1304 case NORMAL_ZERO:
1305 my_memmove((void *)y,(void *)z,sizeof(InternalFPF));
1306 z->sign ^= x->sign;
1307 break;
1308
1309 case ZERO_INFINITY:
1310 case INFINITY_ZERO:
1311 SetInternalFPFNaN(z);
1312 break;
1313
1314 case NAN_ZERO:
1315 case NAN_SUBNORMAL:
1316 case NAN_NORMAL:
1317 case NAN_INFINITY:
1318 my_memmove((void *)x,(void *)z,sizeof(InternalFPF));
1319 break;
1320
1321 case ZERO_NAN:
1322 case SUBNORMAL_NAN:
1323 case NORMAL_NAN:
1324 case INFINITY_NAN:
1325 my_memmove((void *)y,(void *)z,sizeof(InternalFPF));
1326 break;
1327
1328
1329 case SUBNORMAL_SUBNORMAL:
1330 case SUBNORMAL_NORMAL:
1331 case NORMAL_SUBNORMAL:
1332 case NORMAL_NORMAL:
1333 /*
1334 ** Make a local copy of the y number, since we will be
1335 ** altering it in the process of multiplying.
1336 */
1337 my_memmove((void *)&locy,(void *)y,sizeof(InternalFPF));
1338
1339 /*
1340 ** Check for unnormal zero arguments
1341 */
1342 if (IsMantissaZero(x->mantissa) || IsMantissaZero(y->mantissa))
1343 SetInternalFPFInfinity(z, 0);
1344
1345 /*
1346 ** Initialize the result
1347 */
1348 if (x->type == IFPF_IS_SUBNORMAL ||
1349 y->type == IFPF_IS_SUBNORMAL)
1350 z->type = IFPF_IS_SUBNORMAL;
1351 else
1352 z->type = IFPF_IS_NORMAL;
1353
1354 z->sign = x->sign ^ y->sign;
1355 z->exp = x->exp + y->exp ;
1356 for (i=0; i<INTERNAL_FPF_PRECISION; i++)
1357 {
1358 z->mantissa[i] = 0;
1359 extra_bits[i] = 0;
1360 }
1361
1362 for (i=0; i<(INTERNAL_FPF_PRECISION*16); i++)
1363 {
1364 /*
1365 ** Get rightmost bit of the multiplier
1366 */
1367 carry = 0;
1368 ShiftMantRight1(&carry, locy.mantissa);
1369 if (carry)
1370 {
1371 /*
1372 ** Add the multiplicand to the product
1373 */
1374 carry = 0;
1375 for (j=(INTERNAL_FPF_PRECISION-1); j>=0; j--)
1376 Add16Bits(&carry,
1377 &z->mantissa[j],
1378 z->mantissa[j],
1379 x->mantissa[j]);
1380 }
1381 else
1382 {
1383 carry = 0;
1384 }
1385
1386 /*
1387 ** Shift the product right. Overflow bits get
1388 ** shifted into extra_bits. We'll use it later
1389 ** to help with the "sticky" bit.
1390 */
1391 ShiftMantRight1(&carry, z->mantissa);
1392 ShiftMantRight1(&carry, extra_bits);
1393 }
1394
1395 /*
1396 ** Normalize
1397 ** Note that we use a "special" normalization routine
1398 ** because we need to use the extra bits. (These are
1399 ** bits that may have been shifted off the bottom that
1400 ** we want to reclaim...if we can.
1401 */
1402 while ((z->mantissa[0] & 0x8000) == 0)
1403 {
1404 carry = 0;
1405 ShiftMantLeft1(&carry, extra_bits);
1406 ShiftMantLeft1(&carry, z->mantissa);
1407 z->exp--;
1408 }
1409
1410 /*
1411 ** Set the sticky bit if any bits set in extra bits.
1412 */
1413 if (IsMantissaZero(extra_bits))
1414 {
1415 z->mantissa[INTERNAL_FPF_PRECISION-1] |= 1;
1416 }
1417 break;
1418
1419 case NAN_NAN:
1420 choose_nan(x, y, z, 0);
1421 break;
1422 }
1423
1424 /*
1425 ** All math done...do rounding.
1426 */
1427 RoundInternalFPF(z);
1428 return;
1429 }
1430
1431
1432 /**********************
1433 ** DivideInternalFPF **
1434 ***********************
1435 ** Divide internal FPF number x by y. Return result in z.
1436 */
DivideInternalFPF(InternalFPF * x,InternalFPF * y,InternalFPF * z)1437 static void DivideInternalFPF(InternalFPF *x,
1438 InternalFPF *y,
1439 InternalFPF *z)
1440 {
1441 int i;
1442 int j;
1443 u16 carry;
1444 u16 extra_bits[INTERNAL_FPF_PRECISION];
1445 InternalFPF locx; /* Local for x number */
1446
1447 /*
1448 ** As with preceding function, the following switch
1449 ** statement selects among the various possible
1450 ** operands.
1451 */
1452 switch ((x->type * IFPF_TYPE_COUNT) + y->type)
1453 {
1454 case ZERO_ZERO:
1455 case INFINITY_INFINITY:
1456 SetInternalFPFNaN(z);
1457 break;
1458
1459 case ZERO_SUBNORMAL:
1460 case ZERO_NORMAL:
1461 if (IsMantissaZero(y->mantissa))
1462 {
1463 SetInternalFPFNaN(z);
1464 break;
1465 }
1466
1467 case ZERO_INFINITY:
1468 case SUBNORMAL_INFINITY:
1469 case NORMAL_INFINITY:
1470 SetInternalFPFZero(z, x->sign ^ y->sign);
1471 break;
1472
1473 case SUBNORMAL_ZERO:
1474 case NORMAL_ZERO:
1475 if (IsMantissaZero(x->mantissa))
1476 {
1477 SetInternalFPFNaN(z);
1478 break;
1479 }
1480
1481 case INFINITY_ZERO:
1482 case INFINITY_SUBNORMAL:
1483 case INFINITY_NORMAL:
1484 SetInternalFPFInfinity(z, 0);
1485 z->sign = x->sign ^ y->sign;
1486 break;
1487
1488 case NAN_ZERO:
1489 case NAN_SUBNORMAL:
1490 case NAN_NORMAL:
1491 case NAN_INFINITY:
1492 my_memmove((void *)x,(void *)z,sizeof(InternalFPF));
1493 break;
1494
1495 case ZERO_NAN:
1496 case SUBNORMAL_NAN:
1497 case NORMAL_NAN:
1498 case INFINITY_NAN:
1499 my_memmove((void *)y,(void *)z,sizeof(InternalFPF));
1500 break;
1501
1502 case SUBNORMAL_SUBNORMAL:
1503 case NORMAL_SUBNORMAL:
1504 case SUBNORMAL_NORMAL:
1505 case NORMAL_NORMAL:
1506 /*
1507 ** Make local copy of x number, since we'll be
1508 ** altering it in the process of dividing.
1509 */
1510 my_memmove((void *)&locx,(void *)x,sizeof(InternalFPF));
1511
1512 /*
1513 ** Check for unnormal zero arguments
1514 */
1515 if (IsMantissaZero(locx.mantissa))
1516 {
1517 if (IsMantissaZero(y->mantissa))
1518 SetInternalFPFNaN(z);
1519 else
1520 SetInternalFPFZero(z, 0);
1521 break;
1522 }
1523 if (IsMantissaZero(y->mantissa))
1524 {
1525 SetInternalFPFInfinity(z, 0);
1526 break;
1527 }
1528
1529 /*
1530 ** Initialize the result
1531 */
1532 z->type = x->type;
1533 z->sign = x->sign ^ y->sign;
1534 z->exp = x->exp - y->exp +
1535 ((INTERNAL_FPF_PRECISION * 16 * 2));
1536 for (i=0; i<INTERNAL_FPF_PRECISION; i++)
1537 {
1538 z->mantissa[i] = 0;
1539 extra_bits[i] = 0;
1540 }
1541
1542 while ((z->mantissa[0] & 0x8000) == 0)
1543 {
1544 carry = 0;
1545 ShiftMantLeft1(&carry, locx.mantissa);
1546 ShiftMantLeft1(&carry, extra_bits);
1547
1548 /*
1549 ** Time to subtract yet?
1550 */
1551 if (carry == 0)
1552 for (j=0; j<INTERNAL_FPF_PRECISION; j++)
1553 {
1554 if (y->mantissa[j] > extra_bits[j])
1555 {
1556 carry = 0;
1557 goto no_subtract;
1558 }
1559 if (y->mantissa[j] < extra_bits[j])
1560 break;
1561 }
1562 /*
1563 ** Divisor (y) <= dividend (x), subtract
1564 */
1565 carry = 0;
1566 for (j=(INTERNAL_FPF_PRECISION-1); j>=0; j--)
1567 Sub16Bits(&carry,
1568 &extra_bits[j],
1569 extra_bits[j],
1570 y->mantissa[j]);
1571 carry = 1; /* 1 shifted into quotient */
1572 no_subtract:
1573 ShiftMantLeft1(&carry, z->mantissa);
1574 z->exp--;
1575 }
1576 break;
1577
1578 case NAN_NAN:
1579 choose_nan(x, y, z, 0);
1580 break;
1581 }
1582
1583 /*
1584 ** Math complete...do rounding
1585 */
1586 RoundInternalFPF(z);
1587 }
1588
1589 /**********************
1590 ** LongToInternalFPF **
1591 ** Int32ToInternalFPF **
1592 ***********************
1593 ** Convert a signed (long) 32-bit integer into an internal FPF number.
1594 */
1595 /* static void LongToInternalFPF(long mylong, */
Int32ToInternalFPF(int32 mylong,InternalFPF * dest)1596 static void Int32ToInternalFPF(int32 mylong,
1597 InternalFPF *dest)
1598 {
1599 int i; /* Index */
1600 u16 myword; /* Used to hold converted stuff */
1601 /*
1602 ** Save the sign and get the absolute value. This will help us
1603 ** with 64-bit machines, since we use only the lower 32
1604 ** bits just in case. (No longer necessary after we use int32.)
1605 */
1606 /* if(mylong<0L) */
1607 if(mylong<(int32)0)
1608 { dest->sign=1;
1609 mylong=(int32)0-mylong;
1610 }
1611 else
1612 dest->sign=0;
1613 /*
1614 ** Prepare the destination floating point number
1615 */
1616 dest->type=IFPF_IS_NORMAL;
1617 for(i=0;i<INTERNAL_FPF_PRECISION;i++)
1618 dest->mantissa[i]=0;
1619
1620 /*
1621 ** See if we've got a zero. If so, make the resultant FP
1622 ** number a true zero and go home.
1623 */
1624 if(mylong==0)
1625 { dest->type=IFPF_IS_ZERO;
1626 dest->exp=0;
1627 return;
1628 }
1629
1630 /*
1631 ** Not a true zero. Set the exponent to 32 (internal FPFs have
1632 ** no bias) and load the low and high words into their proper
1633 ** locations in the mantissa. Then normalize. The action of
1634 ** normalizing slides the mantissa bits into place and sets
1635 ** up the exponent properly.
1636 */
1637 dest->exp=32;
1638 myword=(u16)((mylong >> 16) & 0xFFFFL);
1639 dest->mantissa[0]=myword;
1640 myword=(u16)(mylong & 0xFFFFL);
1641 dest->mantissa[1]=myword;
1642 normalize(dest);
1643 return;
1644 }
1645
1646 #if 1
1647 /************************
1648 ** InternalFPFToString **
1649 *************************
1650 ** FOR DEBUG PURPOSES
1651 ** This routine converts an internal floating point representation
1652 ** number to a string. Used in debugging the package.
1653 ** Returns length of converted number.
1654 ** NOTE: dest must point to a buffer big enough to hold the
1655 ** result. Also, this routine does append a null (an effect
1656 ** of using the sprintf() function). It also returns
1657 ** a length count.
1658 ** NOTE: This routine returns 5 significant digits. Thats
1659 ** about all I feel safe with, given the method of
1660 ** conversion. It should be more than enough for programmers
1661 ** to determine whether the package is properly ported.
1662 */
InternalFPFToString(char * dest,InternalFPF * src)1663 static int InternalFPFToString(char *dest,
1664 InternalFPF *src)
1665 {
1666 InternalFPF locFPFNum; /* Local for src (will be altered) */
1667 InternalFPF IFPF10; /* Floating-point 10 */
1668 InternalFPF IFPFComp; /* For doing comparisons */
1669 int msign; /* Holding for mantissa sign */
1670 int expcount; /* Exponent counter */
1671 int ccount; /* Character counter */
1672 int i,j,k; /* Index */
1673 u16 carryaccum; /* Carry accumulator */
1674 u16 mycarry; /* Local for carry */
1675
1676 /*
1677 ** Check first for the simple things...Nan, Infinity, Zero.
1678 ** If found, copy the proper string in and go home.
1679 */
1680 switch(src->type)
1681 {
1682 case IFPF_IS_NAN:
1683 my_memcpy(dest,"NaN",3);
1684 return(3);
1685
1686 case IFPF_IS_INFINITY:
1687 if(src->sign==0)
1688 my_memcpy(dest,"+Inf",4);
1689 else
1690 my_memcpy(dest,"-Inf",4);
1691 return(4);
1692
1693 case IFPF_IS_ZERO:
1694 if(src->sign==0)
1695 my_memcpy(dest,"+0",2);
1696 else
1697 my_memcpy(dest,"-0",2);
1698 return(2);
1699 }
1700
1701 /*
1702 ** Move the internal number into our local holding area, since
1703 ** we'll be altering it to print it out.
1704 */
1705 my_memcpy((void *)&locFPFNum,(void *)src,sizeof(InternalFPF));
1706
1707 /*
1708 ** Set up a floating-point 10...which we'll use a lot in a minute.
1709 */
1710 /* LongToInternalFPF(10L,&IFPF10); */
1711 Int32ToInternalFPF((int32)10,&IFPF10);
1712
1713 /*
1714 ** Save the mantissa sign and make it positive.
1715 */
1716 msign=src->sign;
1717
1718 /* src->sign=0 */ /* bug, fixed Nov. 13, 1997 */
1719 (&locFPFNum)->sign=0;
1720
1721 expcount=0; /* Init exponent counter */
1722
1723 /*
1724 ** See if the number is less than 10. If so, multiply
1725 ** the number repeatedly by 10 until it's not. For each
1726 ** multiplication, decrement a counter so we can keep track
1727 ** of the exponent.
1728 */
1729
1730 while(1)
1731 { AddSubInternalFPF(1,&locFPFNum,&IFPF10,&IFPFComp);
1732 if(IFPFComp.sign==0) break;
1733 MultiplyInternalFPF(&locFPFNum,&IFPF10,&IFPFComp);
1734 expcount--;
1735 my_memcpy((void *)&locFPFNum,(void *)&IFPFComp,sizeof(InternalFPF));
1736 }
1737 /*
1738 ** Do the reverse of the above. As long as the number is
1739 ** greater than or equal to 10, divide it by 10. Increment the
1740 ** exponent counter for each multiplication.
1741 */
1742
1743 while(1)
1744 {
1745 AddSubInternalFPF(1,&locFPFNum,&IFPF10,&IFPFComp);
1746 if(IFPFComp.sign!=0) break;
1747 DivideInternalFPF(&locFPFNum,&IFPF10,&IFPFComp);
1748 expcount++;
1749 my_memcpy((void *)&locFPFNum,(void *)&IFPFComp,sizeof(InternalFPF));
1750 }
1751
1752 /*
1753 ** About time to start storing things. First, store the
1754 ** mantissa sign.
1755 */
1756 ccount=1; /* Init character counter */
1757 if(msign==0)
1758 *dest++='+';
1759 else
1760 *dest++='-';
1761
1762 /*
1763 ** At this point we know that the number is in the range
1764 ** 10 > n >=1. We need to "strip digits" out of the
1765 ** mantissa. We do this by treating the mantissa as
1766 ** an integer and multiplying by 10. (Not a floating-point
1767 ** 10, but an integer 10. Since this is debug code and we
1768 ** could care less about speed, we'll do it the stupid
1769 ** way and simply add the number to itself 10 times.
1770 ** Anything that makes it to the left of the implied binary point
1771 ** gets stripped off and emitted. We'll do this for
1772 ** 5 significant digits (which should be enough to
1773 ** verify things).
1774 */
1775 /*
1776 ** Re-position radix point
1777 */
1778 carryaccum=0;
1779 while(locFPFNum.exp>0)
1780 {
1781 mycarry=0;
1782 ShiftMantLeft1(&mycarry,locFPFNum.mantissa);
1783 carryaccum=(carryaccum<<1);
1784 if(mycarry) carryaccum++;
1785 locFPFNum.exp--;
1786 }
1787
1788 while(locFPFNum.exp<0)
1789 {
1790 mycarry=0;
1791 ShiftMantRight1(&mycarry,locFPFNum.mantissa);
1792 locFPFNum.exp++;
1793 }
1794
1795 for(i=0;i<6;i++)
1796 if(i==1)
1797 { /* Emit decimal point */
1798 *dest++='.';
1799 ccount++;
1800 }
1801 else
1802 { /* Emit a digit */
1803 *dest++=('0'+carryaccum);
1804 ccount++;
1805
1806 carryaccum=0;
1807 my_memcpy((void *)&IFPF10,
1808 (void *)&locFPFNum,
1809 sizeof(InternalFPF));
1810
1811 /* Do multiply via repeated adds */
1812 for(j=0;j<9;j++)
1813 {
1814 mycarry=0;
1815 for(k=(INTERNAL_FPF_PRECISION-1);k>=0;k--)
1816 Add16Bits(&mycarry,&(IFPFComp.mantissa[k]),
1817 locFPFNum.mantissa[k],
1818 IFPF10.mantissa[k]);
1819 carryaccum+=mycarry ? 1 : 0;
1820 my_memcpy((void *)&locFPFNum,
1821 (void *)&IFPFComp,
1822 sizeof(InternalFPF));
1823 }
1824 }
1825
1826 /*
1827 ** Now move the 'E', the exponent sign, and the exponent
1828 ** into the string.
1829 */
1830 *dest++='E';
1831
1832 /* sprint is supposed to return an integer, but it caused problems on SunOS
1833 * with the native cc. Hence we force it.
1834 * Uwe F. Mayer
1835 */
1836 if (expcount < 0) {
1837 *dest++ = '-';
1838 expcount =- expcount;
1839 }
1840 else *dest++ = ' ';
1841
1842 *dest++ = (char)(expcount + '0');
1843 *dest++ = 0;
1844
1845 ccount += 3;
1846 /*
1847 ** All done, go home.
1848 */
1849 return(ccount);
1850
1851 }
1852
1853 #endif
1854
1855
1856
1857 ////////////////////////////////////////////////////////////////////////
1858 static
AllocateMemory(unsigned long n,int * p)1859 void* AllocateMemory ( unsigned long n, int* p )
1860 {
1861 *p = 0;
1862 void* r = (void*) (*serviceFn)(2,n);
1863 return r;
1864 }
1865 static
FreeMemory(void * p,int * zz)1866 void FreeMemory ( void* p, int* zz )
1867 {
1868 *zz = 0;
1869 // free(p);
1870 }
1871
1872
1873
1874 /**************
1875 ** DoEmFloat **
1876 ***************
1877 ** Perform the floating-point emulation routines portion of the
1878 ** CPU benchmark. Returns the operations per second.
1879 */
1880 static
DoEmFloat(void)1881 void DoEmFloat(void)
1882 {
1883 EmFloatStruct *locemfloatstruct; /* Local structure */
1884 InternalFPF *abase; /* Base of A array */
1885 InternalFPF *bbase; /* Base of B array */
1886 InternalFPF *cbase; /* Base of C array */
1887 ulong tickcount; /* # of ticks */
1888 char *errorcontext; /* Error context string pointer */
1889 int systemerror; /* For holding error code */
1890 ulong loops; /* # of loops */
1891
1892 /*
1893 ** Link to global structure
1894 */
1895 EmFloatStruct global_emfloatstruct;
1896 global_emfloatstruct.adjust = 0;
1897 global_emfloatstruct.request_secs = 0;
1898 global_emfloatstruct.arraysize = 100;
1899 global_emfloatstruct.loops = 1;
1900 global_emfloatstruct.emflops = 0.0;
1901 locemfloatstruct=&global_emfloatstruct;
1902
1903 /*
1904 ** Set the error context
1905 */
1906 errorcontext="CPU:Floating Emulation";
1907
1908
1909 abase=(InternalFPF *)AllocateMemory(locemfloatstruct->arraysize*sizeof(InternalFPF),
1910 &systemerror);
1911
1912 bbase=(InternalFPF *)AllocateMemory(locemfloatstruct->arraysize*sizeof(InternalFPF),
1913 &systemerror);
1914
1915 cbase=(InternalFPF *)AllocateMemory(locemfloatstruct->arraysize*sizeof(InternalFPF),
1916 &systemerror);
1917
1918 /*
1919 ** Set up the arrays
1920 */
1921 SetupCPUEmFloatArrays(abase,bbase,cbase,locemfloatstruct->arraysize);
1922
1923 loops=100;
1924 tickcount=DoEmFloatIteration(abase,bbase,cbase,
1925 locemfloatstruct->arraysize,
1926 loops);
1927
1928 FreeMemory((void *)abase,&systemerror);
1929 FreeMemory((void *)bbase,&systemerror);
1930 FreeMemory((void *)cbase,&systemerror);
1931
1932 return;
1933 }
1934
1935 //////////////////
entry(HWord (* f)(HWord,HWord))1936 void entry ( HWord(*f)(HWord,HWord) )
1937 {
1938 serviceFn = f;
1939 vexxx_printf("starting\n");
1940 DoEmFloat();
1941 (*serviceFn)(0,0);
1942 }
1943