| /* |
| ** emfloat.c |
| ** Source for emulated floating-point routines. |
| ** BYTEmark (tm) |
| ** BYTE's Native Mode Benchmarks |
| ** Rick Grehan, BYTE Magazine. |
| ** |
| ** Created: |
| ** Last update: 3/95 |
| ** |
| ** DISCLAIMER |
| ** The source, executable, and documentation files that comprise |
| ** the BYTEmark benchmarks are made available on an "as is" basis. |
| ** This means that we at BYTE Magazine have made every reasonable |
| ** effort to verify that the there are no errors in the source and |
| ** executable code. We cannot, however, guarantee that the programs |
| ** are error-free. Consequently, McGraw-HIll and BYTE Magazine make |
| ** no claims in regard to the fitness of the source code, executable |
| ** code, and documentation of the BYTEmark. |
| ** Furthermore, BYTE Magazine, McGraw-Hill, and all employees |
| ** of McGraw-Hill cannot be held responsible for any damages resulting |
| ** from the use of this code or the results obtained from using |
| ** this code. |
| */ |
| |
| |
| #include <stdio.h> |
| #include <string.h> |
| #include "nmglobal.h" |
| #include "emfloat.h" |
| |
| /* |
| ** Floating-point emulator. |
| ** These routines are only "sort of" IEEE-compliant. All work is |
| ** done using an internal representation. Also, the routines do |
| ** not check for many of the exceptions that might occur. |
| ** Still, the external formats produced are IEEE-compatible, |
| ** with the restriction that they presume a low-endian machine |
| ** (though the endianism will not effect the performance). |
| ** |
| ** Some code here was based on work done by Steve Snelgrove of |
| ** Orem, UT. Other code comes from routines presented in |
| ** the long-ago book: "Microprocessor Programming for |
| ** Computer Hobbyists" by Neill Graham. |
| */ |
| |
| /************************** |
| ** SetupCPUEmFloatArrays ** |
| *************************** |
| ** Set up the arrays that will be used in the emulated |
| ** floating-point tests. |
| ** This is done by loading abase and bbase elements with |
| ** random numbers. We use our long-to-floating point |
| ** routine to set them up. |
| ** NOTE: We really don't need the pointer to cbase...cbase |
| ** is overwritten in the benchmark. |
| */ |
| void SetupCPUEmFloatArrays(InternalFPF *abase, |
| InternalFPF *bbase, |
| InternalFPF *cbase, |
| ulong arraysize) |
| { |
| ulong i; |
| InternalFPF locFPF1,locFPF2; |
| /* |
| ** Reset random number generator so things repeat. Inserted by Uwe F. Mayer. |
| */ |
| extern int32 randnum(int32 lngval); |
| randnum((int32)13); |
| |
| for(i=0;i<arraysize;i++) |
| {/* LongToInternalFPF(randwc(50000L),&locFPF1); */ |
| Int32ToInternalFPF(randwc((int32)50000),&locFPF1); |
| /* LongToInternalFPF(randwc(50000L)+1L,&locFPF2); */ |
| Int32ToInternalFPF(randwc((int32)50000)+(int32)1,&locFPF2); |
| DivideInternalFPF(&locFPF1,&locFPF2,abase+i); |
| /* LongToInternalFPF(randwc(50000L)+1L,&locFPF2); */ |
| Int32ToInternalFPF(randwc((int32)50000)+(int32)1,&locFPF2); |
| DivideInternalFPF(&locFPF1,&locFPF2,bbase+i); |
| } |
| return; |
| } |
| |
| /*********************** |
| ** DoEmFloatIteration ** |
| ************************ |
| ** Perform an iteration of the emulated floating-point |
| ** benchmark. Note that "an iteration" can involve multiple |
| ** loops through the benchmark. |
| */ |
| ulong DoEmFloatIteration(InternalFPF *abase, |
| InternalFPF *bbase, |
| InternalFPF *cbase, |
| ulong arraysize, ulong loops) |
| { |
| ulong elapsed; /* For the stopwatch */ |
| static uchar jtable[16] = {0,0,0,0,1,1,1,1,2,2,2,2,2,3,3,3}; |
| ulong i; |
| #ifdef DEBUG |
| int number_of_loops; |
| #endif |
| /* |
| ** Begin timing |
| */ |
| elapsed=StartStopwatch(); |
| #ifdef DEBUG |
| number_of_loops=loops-1; /* the index of the first loop we run */ |
| #endif |
| |
| /* |
| ** Each pass through the array performs operations in |
| ** the followingratios: |
| ** 4 adds, 4 subtracts, 5 multiplies, 3 divides |
| ** (adds and subtracts being nearly the same operation) |
| */ |
| while(loops--) |
| { |
| for(i=0;i<arraysize;i++) |
| switch(jtable[i % 16]) |
| { |
| case 0: /* Add */ |
| AddSubInternalFPF(0,abase+i, |
| bbase+i, |
| cbase+i); |
| break; |
| case 1: /* Subtract */ |
| AddSubInternalFPF(1,abase+i, |
| bbase+i, |
| cbase+i); |
| break; |
| case 2: /* Multiply */ |
| MultiplyInternalFPF(abase+i, |
| bbase+i, |
| cbase+i); |
| break; |
| case 3: /* Divide */ |
| DivideInternalFPF(abase+i, |
| bbase+i, |
| cbase+i); |
| break; |
| } |
| #ifdef DEBUG |
| { |
| ulong j[8]; /* we test 8 entries */ |
| int k; |
| ulong i; |
| char buffer[1024]; |
| if (number_of_loops==loops) /* the first loop */ |
| { |
| j[0]=(ulong)2; |
| j[1]=(ulong)6; |
| j[2]=(ulong)10; |
| j[3]=(ulong)14; |
| j[4]=(ulong)(arraysize-14); |
| j[5]=(ulong)(arraysize-10); |
| j[6]=(ulong)(arraysize-6); |
| j[7]=(ulong)(arraysize-2); |
| for(k=0;k<8;k++){ |
| i=j[k]; |
| InternalFPFToString(buffer,abase+i); |
| printf("%6ld: (%s) ",i,buffer); |
| switch(jtable[i % 16]) |
| { |
| case 0: strcpy(buffer,"+"); break; |
| case 1: strcpy(buffer,"-"); break; |
| case 2: strcpy(buffer,"*"); break; |
| case 3: strcpy(buffer,"/"); break; |
| } |
| printf("%s ",buffer); |
| InternalFPFToString(buffer,bbase+i); |
| printf("(%s) = ",buffer); |
| InternalFPFToString(buffer,cbase+i); |
| printf("%s\n",buffer); |
| } |
| } |
| } |
| #endif |
| } |
| return(StopStopwatch(elapsed)); |
| } |
| |
| /*********************** |
| ** SetInternalFPFZero ** |
| ************************ |
| ** Set an internal floating-point-format number to zero. |
| ** sign determines the sign of the zero. |
| */ |
| static void SetInternalFPFZero(InternalFPF *dest, |
| uchar sign) |
| { |
| int i; /* Index */ |
| |
| dest->type=IFPF_IS_ZERO; |
| dest->sign=sign; |
| dest->exp=MIN_EXP; |
| for(i=0;i<INTERNAL_FPF_PRECISION;i++) |
| dest->mantissa[i]=0; |
| return; |
| } |
| |
| /*************************** |
| ** SetInternalFPFInfinity ** |
| **************************** |
| ** Set an internal floating-point-format number to infinity. |
| ** This can happen if the exponent exceeds MAX_EXP. |
| ** As above, sign picks the sign of infinity. |
| */ |
| static void SetInternalFPFInfinity(InternalFPF *dest, |
| uchar sign) |
| { |
| int i; /* Index */ |
| |
| dest->type=IFPF_IS_INFINITY; |
| dest->sign=sign; |
| dest->exp=MIN_EXP; |
| for(i=0;i<INTERNAL_FPF_PRECISION;i++) |
| dest->mantissa[i]=0; |
| return; |
| } |
| |
| /********************** |
| ** SetInternalFPFNaN ** |
| *********************** |
| ** Set an internal floating-point-format number to Nan |
| ** (not a number). Note that we "emulate" an 80x87 as far |
| ** as the mantissa bits go. |
| */ |
| static void SetInternalFPFNaN(InternalFPF *dest) |
| { |
| int i; /* Index */ |
| |
| dest->type=IFPF_IS_NAN; |
| dest->exp=MAX_EXP; |
| dest->sign=1; |
| dest->mantissa[0]=0x4000; |
| for(i=1;i<INTERNAL_FPF_PRECISION;i++) |
| dest->mantissa[i]=0; |
| |
| return; |
| } |
| |
| /******************* |
| ** IsMantissaZero ** |
| ******************** |
| ** Pass this routine a pointer to an internal floating point format |
| ** number's mantissa. It checks for an all-zero mantissa. |
| ** Returns 0 if it is NOT all zeros, !=0 otherwise. |
| */ |
| static int IsMantissaZero(u16 *mant) |
| { |
| int i; /* Index */ |
| int n; /* Return value */ |
| |
| n=0; |
| for(i=0;i<INTERNAL_FPF_PRECISION;i++) |
| n|=mant[i]; |
| |
| return(!n); |
| } |
| |
| /************** |
| ** Add16Bits ** |
| *************** |
| ** Add b, c, and carry. Retult in a. New carry in carry. |
| */ |
| static void Add16Bits(u16 *carry, |
| u16 *a, |
| u16 b, |
| u16 c) |
| { |
| u32 accum; /* Accumulator */ |
| |
| /* |
| ** Do the work in the 32-bit accumulator so we can return |
| ** the carry. |
| */ |
| accum=(u32)b; |
| accum+=(u32)c; |
| accum+=(u32)*carry; |
| *carry=(u16)((accum & 0x00010000) ? 1 : 0); /* New carry */ |
| *a=(u16)(accum & 0xFFFF); /* Result is lo 16 bits */ |
| return; |
| } |
| |
| /************** |
| ** Sub16Bits ** |
| *************** |
| ** Additive inverse of above. |
| */ |
| static void Sub16Bits(u16 *borrow, |
| u16 *a, |
| u16 b, |
| u16 c) |
| { |
| u32 accum; /* Accumulator */ |
| |
| accum=(u32)b; |
| accum-=(u32)c; |
| accum-=(u32)*borrow; |
| *borrow=(u32)((accum & 0x00010000) ? 1 : 0); /* New borrow */ |
| *a=(u16)(accum & 0xFFFF); |
| return; |
| } |
| |
| /******************* |
| ** ShiftMantLeft1 ** |
| ******************** |
| ** Shift a vector of 16-bit numbers left 1 bit. Also provides |
| ** a carry bit, which is shifted in at the beginning, and |
| ** shifted out at the end. |
| */ |
| static void ShiftMantLeft1(u16 *carry, |
| u16 *mantissa) |
| { |
| int i; /* Index */ |
| int new_carry; |
| u16 accum; /* Temporary holding placed */ |
| |
| for(i=INTERNAL_FPF_PRECISION-1;i>=0;i--) |
| { accum=mantissa[i]; |
| new_carry=accum & 0x8000; /* Get new carry */ |
| accum=accum<<1; /* Do the shift */ |
| if(*carry) |
| accum|=1; /* Insert previous carry */ |
| *carry=new_carry; |
| mantissa[i]=accum; /* Return shifted value */ |
| } |
| return; |
| } |
| |
| /******************** |
| ** ShiftMantRight1 ** |
| ********************* |
| ** Shift a mantissa right by 1 bit. Provides carry, as |
| ** above |
| */ |
| static void ShiftMantRight1(u16 *carry, |
| u16 *mantissa) |
| { |
| int i; /* Index */ |
| int new_carry; |
| u16 accum; |
| |
| for(i=0;i<INTERNAL_FPF_PRECISION;i++) |
| { accum=mantissa[i]; |
| new_carry=accum & 1; /* Get new carry */ |
| accum=accum>>1; |
| if(*carry) |
| accum|=0x8000; |
| *carry=new_carry; |
| mantissa[i]=accum; |
| } |
| return; |
| } |
| |
| |
| /***************************** |
| ** StickyShiftMantRight ** |
| ****************************** |
| ** This is a shift right of the mantissa with a "sticky bit". |
| ** I.E., if a carry of 1 is shifted out of the least significant |
| ** bit, the least significant bit is set to 1. |
| */ |
| static void StickyShiftRightMant(InternalFPF *ptr, |
| int amount) |
| { |
| int i; /* Index */ |
| u16 carry; /* Self-explanatory */ |
| u16 *mantissa; |
| |
| mantissa=ptr->mantissa; |
| |
| if(ptr->type!=IFPF_IS_ZERO) /* Don't bother shifting a zero */ |
| { |
| /* |
| ** If the amount of shifting will shift everyting |
| ** out of existence, then just clear the whole mantissa |
| ** and set the lowmost bit to 1. |
| */ |
| if(amount>=INTERNAL_FPF_PRECISION * 16) |
| { |
| for(i=0;i<INTERNAL_FPF_PRECISION-1;i++) |
| mantissa[i]=0; |
| mantissa[INTERNAL_FPF_PRECISION-1]=1; |
| } |
| else |
| for(i=0;i<amount;i++) |
| { |
| carry=0; |
| ShiftMantRight1(&carry,mantissa); |
| if(carry) |
| mantissa[INTERNAL_FPF_PRECISION-1] |= 1; |
| } |
| } |
| return; |
| } |
| |
| |
| /************************************************** |
| ** POST ARITHMETIC PROCESSING ** |
| ** (NORMALIZE, ROUND, OVERFLOW, AND UNDERFLOW) ** |
| **************************************************/ |
| |
| /************** |
| ** normalize ** |
| *************** |
| ** Normalize an internal-representation number. Normalization |
| ** discards empty most-significant bits. |
| */ |
| static void normalize(InternalFPF *ptr) |
| { |
| u16 carry; |
| |
| /* |
| ** As long as there's a highmost 0 bit, shift the significand |
| ** left 1 bit. Each time you do this, though, you've |
| ** gotta decrement the exponent. |
| */ |
| while ((ptr->mantissa[0] & 0x8000) == 0) |
| { |
| carry = 0; |
| ShiftMantLeft1(&carry, ptr->mantissa); |
| ptr->exp--; |
| } |
| return; |
| } |
| |
| /**************** |
| ** denormalize ** |
| ***************** |
| ** Denormalize an internal-representation number. This means |
| ** shifting it right until its exponent is equivalent to |
| ** minimum_exponent. (You have to do this often in order |
| ** to perform additions and subtractions). |
| */ |
| static void denormalize(InternalFPF *ptr, |
| int minimum_exponent) |
| { |
| long exponent_difference; |
| |
| if (IsMantissaZero(ptr->mantissa)) |
| { |
| printf("Error: zero significand in denormalize\n"); |
| } |
| |
| exponent_difference = ptr->exp-minimum_exponent; |
| if (exponent_difference < 0) |
| { |
| /* |
| ** The number is subnormal |
| */ |
| exponent_difference = -exponent_difference; |
| if (exponent_difference >= (INTERNAL_FPF_PRECISION * 16)) |
| { |
| /* Underflow */ |
| SetInternalFPFZero(ptr, ptr->sign); |
| } |
| else |
| { |
| ptr->exp+=exponent_difference; |
| StickyShiftRightMant(ptr, exponent_difference); |
| } |
| } |
| return; |
| } |
| |
| |
| /********************* |
| ** RoundInternalFPF ** |
| ********************** |
| ** Round an internal-representation number. |
| ** The kind of rounding we do here is simplest...referred to as |
| ** "chop". "Extraneous" rightmost bits are simply hacked off. |
| */ |
| void RoundInternalFPF(InternalFPF *ptr) |
| { |
| /* int i; */ |
| |
| if (ptr->type == IFPF_IS_NORMAL || |
| ptr->type == IFPF_IS_SUBNORMAL) |
| { |
| denormalize(ptr, MIN_EXP); |
| if (ptr->type != IFPF_IS_ZERO) |
| { |
| |
| /* clear the extraneous bits */ |
| ptr->mantissa[3] &= 0xfff8; |
| /* for (i=4; i<INTERNAL_FPF_PRECISION; i++) |
| { |
| ptr->mantissa[i] = 0; |
| } |
| */ |
| /* |
| ** Check for overflow |
| */ |
| /* Does not do anything as ptr->exp is a short and MAX_EXP=37268 |
| if (ptr->exp > MAX_EXP) |
| { |
| SetInternalFPFInfinity(ptr, ptr->sign); |
| } |
| */ |
| } |
| } |
| return; |
| } |
| |
| /******************************************************* |
| ** ARITHMETIC OPERATIONS ON INTERNAL REPRESENTATION ** |
| *******************************************************/ |
| |
| /*************** |
| ** choose_nan ** |
| **************** |
| ** Called by routines that are forced to perform math on |
| ** a pair of NaN's. This routine "selects" which NaN is |
| ** to be returned. |
| */ |
| static void choose_nan(InternalFPF *x, |
| InternalFPF *y, |
| InternalFPF *z, |
| int intel_flag) |
| { |
| int i; |
| |
| /* |
| ** Compare the two mantissas, |
| ** return the larger. Note that we will be emulating |
| ** an 80387 in this operation. |
| */ |
| for (i=0; i<INTERNAL_FPF_PRECISION; i++) |
| { |
| if (x->mantissa[i] > y->mantissa[i]) |
| { |
| memmove((void *)x,(void *)z,sizeof(InternalFPF)); |
| return; |
| } |
| if (x->mantissa[i] < y->mantissa[i]) |
| { |
| memmove((void *)y,(void *)z,sizeof(InternalFPF)); |
| return; |
| } |
| } |
| |
| /* |
| ** They are equal |
| */ |
| if (!intel_flag) |
| /* if the operation is addition */ |
| memmove((void *)x,(void *)z,sizeof(InternalFPF)); |
| else |
| /* if the operation is multiplication */ |
| memmove((void *)y,(void *)z,sizeof(InternalFPF)); |
| return; |
| } |
| |
| |
| /********************** |
| ** AddSubInternalFPF ** |
| *********************** |
| ** Adding or subtracting internal-representation numbers. |
| ** Internal-representation numbers pointed to by x and y are |
| ** added/subtracted and the result returned in z. |
| */ |
| static void AddSubInternalFPF(uchar operation, |
| InternalFPF *x, |
| InternalFPF *y, |
| InternalFPF *z) |
| { |
| int exponent_difference; |
| u16 borrow; |
| u16 carry; |
| int i; |
| InternalFPF locx,locy; /* Needed since we alter them */ |
| |
| /* |
| ** Following big switch statement handles the |
| ** various combinations of operand types. |
| */ |
| switch ((x->type * IFPF_TYPE_COUNT) + y->type) |
| { |
| case ZERO_ZERO: |
| memmove((void *)x,(void *)z,sizeof(InternalFPF)); |
| if (x->sign ^ y->sign ^ operation) |
| { |
| z->sign = 0; /* positive */ |
| } |
| break; |
| |
| case NAN_ZERO: |
| case NAN_SUBNORMAL: |
| case NAN_NORMAL: |
| case NAN_INFINITY: |
| case SUBNORMAL_ZERO: |
| case NORMAL_ZERO: |
| case INFINITY_ZERO: |
| case INFINITY_SUBNORMAL: |
| case INFINITY_NORMAL: |
| memmove((void *)x,(void *)z,sizeof(InternalFPF)); |
| break; |
| |
| |
| case ZERO_NAN: |
| case SUBNORMAL_NAN: |
| case NORMAL_NAN: |
| case INFINITY_NAN: |
| memmove((void *)y,(void *)z,sizeof(InternalFPF)); |
| break; |
| |
| case ZERO_SUBNORMAL: |
| case ZERO_NORMAL: |
| case ZERO_INFINITY: |
| case SUBNORMAL_INFINITY: |
| case NORMAL_INFINITY: |
| memmove((void *)y,(void *)z,sizeof(InternalFPF)); |
| z->sign ^= operation; |
| break; |
| |
| case SUBNORMAL_SUBNORMAL: |
| case SUBNORMAL_NORMAL: |
| case NORMAL_SUBNORMAL: |
| case NORMAL_NORMAL: |
| /* |
| ** Copy x and y to locals, since we may have |
| ** to alter them. |
| */ |
| memmove((void *)&locx,(void *)x,sizeof(InternalFPF)); |
| memmove((void *)&locy,(void *)y,sizeof(InternalFPF)); |
| |
| /* compute sum/difference */ |
| exponent_difference = locx.exp-locy.exp; |
| if (exponent_difference == 0) |
| { |
| /* |
| ** locx.exp == locy.exp |
| ** so, no shifting required |
| */ |
| if (locx.type == IFPF_IS_SUBNORMAL || |
| locy.type == IFPF_IS_SUBNORMAL) |
| z->type = IFPF_IS_SUBNORMAL; |
| else |
| z->type = IFPF_IS_NORMAL; |
| |
| /* |
| ** Assume that locx.mantissa > locy.mantissa |
| */ |
| z->sign = locx.sign; |
| z->exp= locx.exp; |
| } |
| else |
| if (exponent_difference > 0) |
| { |
| /* |
| ** locx.exp > locy.exp |
| */ |
| StickyShiftRightMant(&locy, |
| exponent_difference); |
| z->type = locx.type; |
| z->sign = locx.sign; |
| z->exp = locx.exp; |
| } |
| else /* if (exponent_difference < 0) */ |
| { |
| /* |
| ** locx.exp < locy.exp |
| */ |
| StickyShiftRightMant(&locx, |
| -exponent_difference); |
| z->type = locy.type; |
| z->sign = locy.sign ^ operation; |
| z->exp = locy.exp; |
| } |
| |
| if (locx.sign ^ locy.sign ^ operation) |
| { |
| /* |
| ** Signs are different, subtract mantissas |
| */ |
| borrow = 0; |
| for (i=(INTERNAL_FPF_PRECISION-1); i>=0; i--) |
| Sub16Bits(&borrow, |
| &z->mantissa[i], |
| locx.mantissa[i], |
| locy.mantissa[i]); |
| |
| if (borrow) |
| { |
| /* The y->mantissa was larger than the |
| ** x->mantissa leaving a negative |
| ** result. Change the result back to |
| ** an unsigned number and flip the |
| ** sign flag. |
| */ |
| z->sign = locy.sign ^ operation; |
| borrow = 0; |
| for (i=(INTERNAL_FPF_PRECISION-1); i>=0; i--) |
| { |
| Sub16Bits(&borrow, |
| &z->mantissa[i], |
| 0, |
| z->mantissa[i]); |
| } |
| } |
| else |
| { |
| /* The assumption made above |
| ** (i.e. x->mantissa >= y->mantissa) |
| ** was correct. Therefore, do nothing. |
| ** z->sign = x->sign; |
| */ |
| } |
| |
| if (IsMantissaZero(z->mantissa)) |
| { |
| z->type = IFPF_IS_ZERO; |
| z->sign = 0; /* positive */ |
| } |
| else |
| if (locx.type == IFPF_IS_NORMAL || |
| locy.type == IFPF_IS_NORMAL) |
| { |
| normalize(z); |
| } |
| } |
| else |
| { |
| /* signs are the same, add mantissas */ |
| carry = 0; |
| for (i=(INTERNAL_FPF_PRECISION-1); i>=0; i--) |
| { |
| Add16Bits(&carry, |
| &z->mantissa[i], |
| locx.mantissa[i], |
| locy.mantissa[i]); |
| } |
| |
| if (carry) |
| { |
| z->exp++; |
| carry=0; |
| ShiftMantRight1(&carry,z->mantissa); |
| z->mantissa[0] |= 0x8000; |
| z->type = IFPF_IS_NORMAL; |
| } |
| else |
| if (z->mantissa[0] & 0x8000) |
| z->type = IFPF_IS_NORMAL; |
| } |
| break; |
| |
| case INFINITY_INFINITY: |
| SetInternalFPFNaN(z); |
| break; |
| |
| case NAN_NAN: |
| choose_nan(x, y, z, 1); |
| break; |
| } |
| |
| /* |
| ** All the math is done; time to round. |
| */ |
| RoundInternalFPF(z); |
| return; |
| } |
| |
| |
| /************************ |
| ** MultiplyInternalFPF ** |
| ************************* |
| ** Two internal-representation numbers x and y are multiplied; the |
| ** result is returned in z. |
| */ |
| static void MultiplyInternalFPF(InternalFPF *x, |
| InternalFPF *y, |
| InternalFPF *z) |
| { |
| int i; |
| int j; |
| u16 carry; |
| u16 extra_bits[INTERNAL_FPF_PRECISION]; |
| InternalFPF locy; /* Needed since this will be altered */ |
| /* |
| ** As in the preceding function, this large switch |
| ** statement selects among the many combinations |
| ** of operands. |
| */ |
| switch ((x->type * IFPF_TYPE_COUNT) + y->type) |
| { |
| case INFINITY_SUBNORMAL: |
| case INFINITY_NORMAL: |
| case INFINITY_INFINITY: |
| case ZERO_ZERO: |
| case ZERO_SUBNORMAL: |
| case ZERO_NORMAL: |
| memmove((void *)x,(void *)z,sizeof(InternalFPF)); |
| z->sign ^= y->sign; |
| break; |
| |
| case SUBNORMAL_INFINITY: |
| case NORMAL_INFINITY: |
| case SUBNORMAL_ZERO: |
| case NORMAL_ZERO: |
| memmove((void *)y,(void *)z,sizeof(InternalFPF)); |
| z->sign ^= x->sign; |
| break; |
| |
| case ZERO_INFINITY: |
| case INFINITY_ZERO: |
| SetInternalFPFNaN(z); |
| break; |
| |
| case NAN_ZERO: |
| case NAN_SUBNORMAL: |
| case NAN_NORMAL: |
| case NAN_INFINITY: |
| memmove((void *)x,(void *)z,sizeof(InternalFPF)); |
| break; |
| |
| case ZERO_NAN: |
| case SUBNORMAL_NAN: |
| case NORMAL_NAN: |
| case INFINITY_NAN: |
| memmove((void *)y,(void *)z,sizeof(InternalFPF)); |
| break; |
| |
| |
| case SUBNORMAL_SUBNORMAL: |
| case SUBNORMAL_NORMAL: |
| case NORMAL_SUBNORMAL: |
| case NORMAL_NORMAL: |
| /* |
| ** Make a local copy of the y number, since we will be |
| ** altering it in the process of multiplying. |
| */ |
| memmove((void *)&locy,(void *)y,sizeof(InternalFPF)); |
| |
| /* |
| ** Check for unnormal zero arguments |
| */ |
| if (IsMantissaZero(x->mantissa) || IsMantissaZero(y->mantissa)) |
| SetInternalFPFInfinity(z, 0); |
| |
| /* |
| ** Initialize the result |
| */ |
| if (x->type == IFPF_IS_SUBNORMAL || |
| y->type == IFPF_IS_SUBNORMAL) |
| z->type = IFPF_IS_SUBNORMAL; |
| else |
| z->type = IFPF_IS_NORMAL; |
| |
| z->sign = x->sign ^ y->sign; |
| z->exp = x->exp + y->exp ; |
| for (i=0; i<INTERNAL_FPF_PRECISION; i++) |
| { |
| z->mantissa[i] = 0; |
| extra_bits[i] = 0; |
| } |
| |
| for (i=0; i<(INTERNAL_FPF_PRECISION*16); i++) |
| { |
| /* |
| ** Get rightmost bit of the multiplier |
| */ |
| carry = 0; |
| ShiftMantRight1(&carry, locy.mantissa); |
| if (carry) |
| { |
| /* |
| ** Add the multiplicand to the product |
| */ |
| carry = 0; |
| for (j=(INTERNAL_FPF_PRECISION-1); j>=0; j--) |
| Add16Bits(&carry, |
| &z->mantissa[j], |
| z->mantissa[j], |
| x->mantissa[j]); |
| } |
| else |
| { |
| carry = 0; |
| } |
| |
| /* |
| ** Shift the product right. Overflow bits get |
| ** shifted into extra_bits. We'll use it later |
| ** to help with the "sticky" bit. |
| */ |
| ShiftMantRight1(&carry, z->mantissa); |
| ShiftMantRight1(&carry, extra_bits); |
| } |
| |
| /* |
| ** Normalize |
| ** Note that we use a "special" normalization routine |
| ** because we need to use the extra bits. (These are |
| ** bits that may have been shifted off the bottom that |
| ** we want to reclaim...if we can. |
| */ |
| while ((z->mantissa[0] & 0x8000) == 0) |
| { |
| carry = 0; |
| ShiftMantLeft1(&carry, extra_bits); |
| ShiftMantLeft1(&carry, z->mantissa); |
| z->exp--; |
| } |
| |
| /* |
| ** Set the sticky bit if any bits set in extra bits. |
| */ |
| if (IsMantissaZero(extra_bits)) |
| { |
| z->mantissa[INTERNAL_FPF_PRECISION-1] |= 1; |
| } |
| break; |
| |
| case NAN_NAN: |
| choose_nan(x, y, z, 0); |
| break; |
| } |
| |
| /* |
| ** All math done...do rounding. |
| */ |
| RoundInternalFPF(z); |
| return; |
| } |
| |
| |
| /********************** |
| ** DivideInternalFPF ** |
| *********************** |
| ** Divide internal FPF number x by y. Return result in z. |
| */ |
| static void DivideInternalFPF(InternalFPF *x, |
| InternalFPF *y, |
| InternalFPF *z) |
| { |
| int i; |
| int j; |
| u16 carry; |
| u16 extra_bits[INTERNAL_FPF_PRECISION]; |
| InternalFPF locx; /* Local for x number */ |
| |
| /* |
| ** As with preceding function, the following switch |
| ** statement selects among the various possible |
| ** operands. |
| */ |
| switch ((x->type * IFPF_TYPE_COUNT) + y->type) |
| { |
| case ZERO_ZERO: |
| case INFINITY_INFINITY: |
| SetInternalFPFNaN(z); |
| break; |
| |
| case ZERO_SUBNORMAL: |
| case ZERO_NORMAL: |
| if (IsMantissaZero(y->mantissa)) |
| { |
| SetInternalFPFNaN(z); |
| break; |
| } |
| |
| case ZERO_INFINITY: |
| case SUBNORMAL_INFINITY: |
| case NORMAL_INFINITY: |
| SetInternalFPFZero(z, x->sign ^ y->sign); |
| break; |
| |
| case SUBNORMAL_ZERO: |
| case NORMAL_ZERO: |
| if (IsMantissaZero(x->mantissa)) |
| { |
| SetInternalFPFNaN(z); |
| break; |
| } |
| |
| case INFINITY_ZERO: |
| case INFINITY_SUBNORMAL: |
| case INFINITY_NORMAL: |
| SetInternalFPFInfinity(z, 0); |
| z->sign = x->sign ^ y->sign; |
| break; |
| |
| case NAN_ZERO: |
| case NAN_SUBNORMAL: |
| case NAN_NORMAL: |
| case NAN_INFINITY: |
| memmove((void *)x,(void *)z,sizeof(InternalFPF)); |
| break; |
| |
| case ZERO_NAN: |
| case SUBNORMAL_NAN: |
| case NORMAL_NAN: |
| case INFINITY_NAN: |
| memmove((void *)y,(void *)z,sizeof(InternalFPF)); |
| break; |
| |
| case SUBNORMAL_SUBNORMAL: |
| case NORMAL_SUBNORMAL: |
| case SUBNORMAL_NORMAL: |
| case NORMAL_NORMAL: |
| /* |
| ** Make local copy of x number, since we'll be |
| ** altering it in the process of dividing. |
| */ |
| memmove((void *)&locx,(void *)x,sizeof(InternalFPF)); |
| |
| /* |
| ** Check for unnormal zero arguments |
| */ |
| if (IsMantissaZero(locx.mantissa)) |
| { |
| if (IsMantissaZero(y->mantissa)) |
| SetInternalFPFNaN(z); |
| else |
| SetInternalFPFZero(z, 0); |
| break; |
| } |
| if (IsMantissaZero(y->mantissa)) |
| { |
| SetInternalFPFInfinity(z, 0); |
| break; |
| } |
| |
| /* |
| ** Initialize the result |
| */ |
| z->type = x->type; |
| z->sign = x->sign ^ y->sign; |
| z->exp = x->exp - y->exp + |
| ((INTERNAL_FPF_PRECISION * 16 * 2)); |
| for (i=0; i<INTERNAL_FPF_PRECISION; i++) |
| { |
| z->mantissa[i] = 0; |
| extra_bits[i] = 0; |
| } |
| |
| while ((z->mantissa[0] & 0x8000) == 0) |
| { |
| carry = 0; |
| ShiftMantLeft1(&carry, locx.mantissa); |
| ShiftMantLeft1(&carry, extra_bits); |
| |
| /* |
| ** Time to subtract yet? |
| */ |
| if (carry == 0) |
| for (j=0; j<INTERNAL_FPF_PRECISION; j++) |
| { |
| if (y->mantissa[j] > extra_bits[j]) |
| { |
| carry = 0; |
| goto no_subtract; |
| } |
| if (y->mantissa[j] < extra_bits[j]) |
| break; |
| } |
| /* |
| ** Divisor (y) <= dividend (x), subtract |
| */ |
| carry = 0; |
| for (j=(INTERNAL_FPF_PRECISION-1); j>=0; j--) |
| Sub16Bits(&carry, |
| &extra_bits[j], |
| extra_bits[j], |
| y->mantissa[j]); |
| carry = 1; /* 1 shifted into quotient */ |
| no_subtract: |
| ShiftMantLeft1(&carry, z->mantissa); |
| z->exp--; |
| } |
| break; |
| |
| case NAN_NAN: |
| choose_nan(x, y, z, 0); |
| break; |
| } |
| |
| /* |
| ** Math complete...do rounding |
| */ |
| RoundInternalFPF(z); |
| } |
| |
| /********************** |
| ** LongToInternalFPF ** |
| ** Int32ToInternalFPF ** |
| *********************** |
| ** Convert a signed (long) 32-bit integer into an internal FPF number. |
| */ |
| /* static void LongToInternalFPF(long mylong, */ |
| static void Int32ToInternalFPF(int32 mylong, |
| InternalFPF *dest) |
| { |
| int i; /* Index */ |
| u16 myword; /* Used to hold converted stuff */ |
| /* |
| ** Save the sign and get the absolute value. This will help us |
| ** with 64-bit machines, since we use only the lower 32 |
| ** bits just in case. (No longer necessary after we use int32.) |
| */ |
| /* if(mylong<0L) */ |
| if(mylong<(int32)0) |
| { dest->sign=1; |
| mylong=(int32)0-mylong; |
| } |
| else |
| dest->sign=0; |
| /* |
| ** Prepare the destination floating point number |
| */ |
| dest->type=IFPF_IS_NORMAL; |
| for(i=0;i<INTERNAL_FPF_PRECISION;i++) |
| dest->mantissa[i]=0; |
| |
| /* |
| ** See if we've got a zero. If so, make the resultant FP |
| ** number a true zero and go home. |
| */ |
| if(mylong==0) |
| { dest->type=IFPF_IS_ZERO; |
| dest->exp=0; |
| return; |
| } |
| |
| /* |
| ** Not a true zero. Set the exponent to 32 (internal FPFs have |
| ** no bias) and load the low and high words into their proper |
| ** locations in the mantissa. Then normalize. The action of |
| ** normalizing slides the mantissa bits into place and sets |
| ** up the exponent properly. |
| */ |
| dest->exp=32; |
| myword=(u16)((mylong >> 16) & 0xFFFFL); |
| dest->mantissa[0]=myword; |
| myword=(u16)(mylong & 0xFFFFL); |
| dest->mantissa[1]=myword; |
| normalize(dest); |
| return; |
| } |
| |
| #ifdef DEBUG |
| /************************ |
| ** InternalFPFToString ** |
| ************************* |
| ** FOR DEBUG PURPOSES |
| ** This routine converts an internal floating point representation |
| ** number to a string. Used in debugging the package. |
| ** Returns length of converted number. |
| ** NOTE: dest must point to a buffer big enough to hold the |
| ** result. Also, this routine does append a null (an effect |
| ** of using the sprintf() function). It also returns |
| ** a length count. |
| ** NOTE: This routine returns 5 significant digits. Thats |
| ** about all I feel safe with, given the method of |
| ** conversion. It should be more than enough for programmers |
| ** to determine whether the package is properly ported. |
| */ |
| static int InternalFPFToString(char *dest, |
| InternalFPF *src) |
| { |
| InternalFPF locFPFNum; /* Local for src (will be altered) */ |
| InternalFPF IFPF10; /* Floating-point 10 */ |
| InternalFPF IFPFComp; /* For doing comparisons */ |
| int msign; /* Holding for mantissa sign */ |
| int expcount; /* Exponent counter */ |
| int ccount; /* Character counter */ |
| int i,j,k; /* Index */ |
| u16 carryaccum; /* Carry accumulator */ |
| u16 mycarry; /* Local for carry */ |
| |
| /* |
| ** Check first for the simple things...Nan, Infinity, Zero. |
| ** If found, copy the proper string in and go home. |
| */ |
| switch(src->type) |
| { |
| case IFPF_IS_NAN: |
| memcpy(dest,"NaN",3); |
| return(3); |
| |
| case IFPF_IS_INFINITY: |
| if(src->sign==0) |
| memcpy(dest,"+Inf",4); |
| else |
| memcpy(dest,"-Inf",4); |
| return(4); |
| |
| case IFPF_IS_ZERO: |
| if(src->sign==0) |
| memcpy(dest,"+0",2); |
| else |
| memcpy(dest,"-0",2); |
| return(2); |
| } |
| |
| /* |
| ** Move the internal number into our local holding area, since |
| ** we'll be altering it to print it out. |
| */ |
| memcpy((void *)&locFPFNum,(void *)src,sizeof(InternalFPF)); |
| |
| /* |
| ** Set up a floating-point 10...which we'll use a lot in a minute. |
| */ |
| /* LongToInternalFPF(10L,&IFPF10); */ |
| Int32ToInternalFPF((int32)10,&IFPF10); |
| |
| /* |
| ** Save the mantissa sign and make it positive. |
| */ |
| msign=src->sign; |
| |
| /* src->sign=0 */ /* bug, fixed Nov. 13, 1997 */ |
| (&locFPFNum)->sign=0; |
| |
| expcount=0; /* Init exponent counter */ |
| |
| /* |
| ** See if the number is less than 10. If so, multiply |
| ** the number repeatedly by 10 until it's not. For each |
| ** multiplication, decrement a counter so we can keep track |
| ** of the exponent. |
| */ |
| |
| while(1) |
| { AddSubInternalFPF(1,&locFPFNum,&IFPF10,&IFPFComp); |
| if(IFPFComp.sign==0) break; |
| MultiplyInternalFPF(&locFPFNum,&IFPF10,&IFPFComp); |
| expcount--; |
| memcpy((void *)&locFPFNum,(void *)&IFPFComp,sizeof(InternalFPF)); |
| } |
| /* |
| ** Do the reverse of the above. As long as the number is |
| ** greater than or equal to 10, divide it by 10. Increment the |
| ** exponent counter for each multiplication. |
| */ |
| |
| while(1) |
| { |
| AddSubInternalFPF(1,&locFPFNum,&IFPF10,&IFPFComp); |
| if(IFPFComp.sign!=0) break; |
| DivideInternalFPF(&locFPFNum,&IFPF10,&IFPFComp); |
| expcount++; |
| memcpy((void *)&locFPFNum,(void *)&IFPFComp,sizeof(InternalFPF)); |
| } |
| |
| /* |
| ** About time to start storing things. First, store the |
| ** mantissa sign. |
| */ |
| ccount=1; /* Init character counter */ |
| if(msign==0) |
| *dest++='+'; |
| else |
| *dest++='-'; |
| |
| /* |
| ** At this point we know that the number is in the range |
| ** 10 > n >=1. We need to "strip digits" out of the |
| ** mantissa. We do this by treating the mantissa as |
| ** an integer and multiplying by 10. (Not a floating-point |
| ** 10, but an integer 10. Since this is debug code and we |
| ** could care less about speed, we'll do it the stupid |
| ** way and simply add the number to itself 10 times. |
| ** Anything that makes it to the left of the implied binary point |
| ** gets stripped off and emitted. We'll do this for |
| ** 5 significant digits (which should be enough to |
| ** verify things). |
| */ |
| /* |
| ** Re-position radix point |
| */ |
| carryaccum=0; |
| while(locFPFNum.exp>0) |
| { |
| mycarry=0; |
| ShiftMantLeft1(&mycarry,locFPFNum.mantissa); |
| carryaccum=(carryaccum<<1); |
| if(mycarry) carryaccum++; |
| locFPFNum.exp--; |
| } |
| |
| while(locFPFNum.exp<0) |
| { |
| mycarry=0; |
| ShiftMantRight1(&mycarry,locFPFNum.mantissa); |
| locFPFNum.exp++; |
| } |
| |
| for(i=0;i<6;i++) |
| if(i==1) |
| { /* Emit decimal point */ |
| *dest++='.'; |
| ccount++; |
| } |
| else |
| { /* Emit a digit */ |
| *dest++=('0'+carryaccum); |
| ccount++; |
| |
| carryaccum=0; |
| memcpy((void *)&IFPF10, |
| (void *)&locFPFNum, |
| sizeof(InternalFPF)); |
| |
| /* Do multiply via repeated adds */ |
| for(j=0;j<9;j++) |
| { |
| mycarry=0; |
| for(k=(INTERNAL_FPF_PRECISION-1);k>=0;k--) |
| Add16Bits(&mycarry,&(IFPFComp.mantissa[k]), |
| locFPFNum.mantissa[k], |
| IFPF10.mantissa[k]); |
| carryaccum+=mycarry ? 1 : 0; |
| memcpy((void *)&locFPFNum, |
| (void *)&IFPFComp, |
| sizeof(InternalFPF)); |
| } |
| } |
| |
| /* |
| ** Now move the 'E', the exponent sign, and the exponent |
| ** into the string. |
| */ |
| *dest++='E'; |
| |
| /* sprint is supposed to return an integer, but it caused problems on SunOS |
| * with the native cc. Hence we force it. |
| * Uwe F. Mayer |
| */ |
| ccount+=(int)sprintf(dest,"%4d",expcount); |
| |
| /* |
| ** All done, go home. |
| */ |
| return(ccount); |
| |
| } |
| |
| #endif |
