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