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 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 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 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 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 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 112 static Int vexxx_strlen ( const HChar* str ) 113 { 114 Int i = 0; 115 while (str[i] != 0) i++; 116 return i; 117 } 118 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 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 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 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 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 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 555 static int32 randwc(int32 num) 556 { 557 return(randnum((int32)0)%num); 558 } 559 560 static int32 randw[2] = { (int32)13 , (int32)117 }; 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 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 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 */ 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 */ 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 */ 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 */ 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 */ 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 */ 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 */ 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 */ 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 */ 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 */ 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 */ 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 */ 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 */ 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 */ 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 */ 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 */ 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, */ 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 */ 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 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 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 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 ////////////////// 1938 void entry ( HWord(*f)(HWord,HWord) ) 1939 { 1940 serviceFn = f; 1941 vexxx_printf("starting\n"); 1942 DoEmFloat(); 1943 (*serviceFn)(0,0); 1944 } 1945
