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