/*-------------------------------------------------------------------------- ** ** Copyright (c) 2003, Steve Peltz and Tom Hunter (see license.txt) ** ** Name: float.c ** ** Description: ** Perform CDC 6600 floating point unit functions. ** ** Note: ** Lot of this source file has been "borrowed" from a floating point ** implementation done by Steve Peltz for his Zephyr emulator. ** **-------------------------------------------------------------------------- */ /* ** ------------- ** Include Files ** ------------- */ #include #include #include #include "const.h" #include "types.h" #include "proto.h" /* ** ----------------- ** Private Constants ** ----------------- */ #define ID ((CpWord)01777) #define OR ((CpWord)03777) #define IR(x) (((x) & ID) != ID) #define OVFL(s) ((OR ^ (s >> 48)) << 48) #define IND (ID << 48) /* ** ----------------------- ** Private Macro Functions ** ----------------------- */ /* ** The following code fragments are used throughout: ** ** sign = SignX(v, 60); extract sign bit ** v ^= sign; make absolute ** exponent = (v >> 48); extract exponent ** ** exponent -= 02000; un-bias exponent ** exponent -= (exponent >> 11); make two's complement ** ** exponent += 02000 + (exponent >> 11); make one's comp ** ** (if exponent is still less than zero (two's complement) then ** it is an underflow condition) */ #define SignX(v, bit) (((v) & ((CpWord)1 << ((bit) - 1))) == 0 ? 0 : Mask60) /* ** ----------------------------------------- ** Private Typedef and Structure Definitions ** ----------------------------------------- */ /* ** --------------------------- ** Private Function Prototypes ** --------------------------- */ /* ** ---------------- ** Public Variables ** ---------------- */ /* ** ----------------- ** Private Variables ** ----------------- */ /* **-------------------------------------------------------------------------- ** ** Public Functions ** **-------------------------------------------------------------------------- */ /* ** Visual C++ 6.0 Enterprise Edition appears to have a bug in its optimizer ** which results in incorrect maths. The line affected by this bug is marked ** with a comment "line triggers Visual C++ OPTIMIZER BUG". To avoid the problem ** optimization is turned off for the whole file. */ #if defined(_WIN32) #pragma optimize("", off) #endif /*-------------------------------------------------------------------------- ** Purpose: Floating point addition. ** Extract both signs, make absolute, extract exponents, ** check for special cases, shift smaller exponent to the right, ** do a 96 bit signed one's complement add, shift right by one ** if overflow. ** Return upper 48 bits and adjusted exponent for ** ** Parameters: Name Description. ** v1 First operand ** v2 Second operand ** doRound TRUE if rounding required, FALSE otherwise. ** doDouble TRUE if double precision required, FALSE otherwise. ** ** Returns: Upper 48 bits and adjusted exponent for single precision. ** Lower 48 bits and adjusted exponent for double precision. ** **------------------------------------------------------------------------*/ CpWord floatAdd(CpWord v1, CpWord v2, bool doRound, bool doDouble) { CpWord sign1; CpWord sign2; int exponent1; int exponent2; int shift; CpWord upper; /* upper half of 98 bit adder register - 50 bits */ CpWord lower; /* bottom half is 48 bits */ CpWord round = 0; CpWord result; sign1 = SignX(v1, 60); sign2 = SignX(v2, 60); v1 ^= sign1; v2 ^= sign2; exponent1 = (int)(v1 >> 48); exponent2 = (int)(v2 >> 48); if (!IR(exponent1)) { if ((exponent1 == ID) || (exponent2 == ID) || ((exponent2 == OR) && (sign1 != sign2))) { return IND; } return OVFL(sign1); } if (!IR(exponent2)) { if (exponent2 == ID) { return IND; } return OVFL(sign2); } exponent1 -= 02000; exponent2 -= 02000; exponent1 -= (exponent1 >> 11); exponent2 -= (exponent2 >> 11); /* ** pre-set round bit - second rounding bit is inserted if ** both values are normalized or the signs are different. ** bit 47 is the rounding bit. calculate this part here ** even though it won't be needed if shift is 48 or more, ** so the variables can be over-used. ** ** if both variables are normalized (both have bit 47 set), then ** the round bit will be set. if the signs are different, all 60 ** bits will be set, so the value must be masked later on if it is ** used. */ if (doRound) { round = (v1 & v2) | (sign1 ^ sign2); } sign1 >>= 10; sign2 >>= 10; /* ** move value with larger exponent to 98 bit "register", and ** leave smaller value in v1, larger exponent in exponent1, ** and calculate shift count. note that the lower 48 bits of ** the first value are extended with the sign. ** ** value with larger exponent always has a rounding bit ** inserted after the least significant bit (top bit of ** lower part of register). */ if (exponent1 > exponent2) { upper = (v1 & Mask48) ^ sign1; v1 = v2 & Mask48; lower = sign1 >> 2; if (doRound) { lower = Sign48 ^ lower; /* rounding bit */ } sign1 = sign2; shift = exponent1 - exponent2; } else { upper = (v2 & Mask48) ^ sign2; v1 &= Mask48; lower = sign2 >> 2; if (doRound) { lower = Sign48 ^ lower; /* rounding bit */ } shift = exponent2 - exponent1; exponent1 = exponent2; } /* ** sign1 is 50 bit sign extension of second (smaller exponent) value; ** make sign2 be the 48 bit version for the lower half. */ sign2 = sign1 >> 2; /* ** three possible choices - if shift count is less than 48 bits, ** add to both upper and lower (or just upper if shift count is ** zero, but that gets taken care of automatically); if shift count ** is 48 through 95, add to just the lower half; else just add in ** sign bits (result is shifted off the end of the register). ** ** if shift is less than 48, insert shifted rounding bit. */ if (shift < 48) { upper += (v1 >> shift) ^ sign1; if (doRound) { round = (round & Sign48) >> shift; lower += (((v1 << (48 - shift)) & Mask48) | round) ^ sign2; } else { lower += ((v1 << (48 - shift)) & Mask48) ^ sign2; } } else if (shift < 96) { upper += sign1; lower += (v1 >> (shift - 48)) ^ sign2; } else { upper += sign1; lower += sign2; } /* ** carry-out from lower to upper, mask off overflow, add one if ** adding one would cause an end-around carry, then carry-out ** again from lower to upper. this is the same algorithm as the ** 18 or 60 bit one's complement add uses, and adjusts for -0 ** the same way the cyber does (cyber actually inverts the second ** operand and subtracts, rather than adds) */ upper += lower >> 48; lower &= Mask48; lower += (upper + ((lower + 1) >> 48)) >> 50; // line triggers Visual C++ OPTIMIZER BUG upper += lower >> 48; upper &= Mask50; lower &= Mask48; /* ** get sign of result, and make result absolute */ sign1 = SignX(upper, 50); upper ^= (sign1 >> 10); lower ^= (sign1 >> 12); if (doDouble) { /* ** post-normalize - shift bottom bit of upper half to upper ** bit of bottom half (instead of shifting the upper half, as ** the other add routines do) */ if (upper >> 48) { lower = ((upper & 1) << 47) | (lower >> 1); exponent1 += 1; } /* ** check if exponent will be out of range after subtracting 48 (if ** resulting exponent is -1777 (octal), the result is an underflow, ** but the significant digits are still returned; if it is less than ** that, a positive zero is returned). -1777 is 0000 in one's complement ** biased form. */ if (exponent1 < -01717) { return 0; } exponent1 -= 48; exponent1 += 02000 + (exponent1 >> 11); result = ((((CpWord) exponent1) << 48) | (lower & Mask48)) ^ sign1; } else { /* ** post-normalize if necessary */ if (upper >> 48) { upper >>= 1; exponent1 += 1; } exponent1 += 02000 + (exponent1 >> 11); result = ((((CpWord) exponent1) << 48) | upper) ^ sign1; } return(result); } /*-------------------------------------------------------------------------- ** Purpose: Floating multiply ** do four 24 bit multiplies, combine, offset exponent and return ** upper 48 bits of result. does a one bit post-normalize if ** necessary if both values were normalized to begin with. ** ** Double precision multiply. returns the bottom half of the 96 ** bit product. also used for integer multiply by checking for ** both exponents zero and one or both values not normalized. ** ** Rounding multiply is almost identical to floating multiply, ** except that a single rounding bit is added to the result in ** bit 46 of the 96 bit product. ** ** Parameters: Name Description. ** v1 First operand ** v2 Second operand ** doRound TRUE if rounding required, FALSE otherwise. ** doDouble TRUE if double precision required, FALSE otherwise. ** ** Returns: Upper 48 bits and adjusted exponent for single precision. ** Lower 48 bits and adjusted exponent for double precision. ** **------------------------------------------------------------------------*/ CpWord floatMultiply(CpWord v1, CpWord v2, bool doRound, bool doDouble) { CpWord sign1; CpWord sign2; int exponent1; int exponent2; int norm; /* flag for post-normalize */ CpWord upper; /* upper 48 bits of product */ CpWord middle; /* middle cross-product */ CpWord lower; /* lower 48 bits of product */ sign1 = SignX(v1, 60); sign2 = SignX(v2, 60); v1 ^= sign1; v2 ^= sign2; /* ** get sign of result */ sign1 ^= sign2; exponent1 = (int)(v1 >> 48); exponent2 = (int)(v2 >> 48); if (!IR(exponent1)) { if ((exponent1 == ID) || (exponent2 == ID) || !exponent2) { return IND; } return OVFL(sign1); } if (!IR(exponent2)) { if ((exponent2 == ID) || !exponent1) { return IND; } return OVFL(sign1); } v1 &= Mask48; v2 &= Mask48; /* ** get the post-normalize flag */ norm = (int)((v1 & v2) >> 47); /* ** form middle cross-product, upper and lower product, and add them ** all together, with a carry from lower to upper. */ middle = (v1 & Mask24) * (v2 >> 24); if (doRound) { /* ** rounding bit (46) is bit 22 in the middle cross-product. */ middle += ((CpWord)1 << 22); } middle += (v1 >> 24) * (v2 & Mask24); lower = (v1 & Mask24) * (v2 & Mask24); lower += (middle & Mask24) << 24; upper = (v1 >> 24) * (v2 >> 24); upper += (middle >> 24) + (lower >> 48); /* ** do an integer multiply if one or both values are not normalized ** and both exponents are zero (this is really only specified for ** the double precision multiply, but the same check is done for ** floating and rounding as well - this is necessary to make -0 ** results come out correctly. */ if (doDouble) { if (!(norm | exponent1 | exponent2)) { return (lower & Mask48) ^ sign1; } if (!(exponent1 && exponent2)) { return 0; } upper = (v1 >> 24) * (v2 >> 24); upper += (middle >> 24) + (lower >> 48); exponent1 -= 02000; exponent2 -= 02000; exponent1 -= (exponent1 >> 11); exponent2 -= (exponent2 >> 11); exponent1 += exponent2; if (norm && !(upper >> 47)) { lower <<= 1; exponent1 -= 1; } /* ** since the bottom half is returned, exponent doesn't need to be ** offset by 48. */ if (exponent1 > 01777) { return OVFL(sign1); } exponent1 += 02000 + (exponent1 >> 11); if (exponent1 < 0) { return 0; } return ((((CpWord) exponent1) << 48) | (lower & Mask48)) ^ sign1; } if (!(norm | exponent1 | exponent2)) { return upper ^ sign1; } /* ** if not an integer multiply and one or both exponents are zero ** (underflow), return positive zero. */ if (!(exponent1 && exponent2)) return 0; exponent1 -= 02000; exponent2 -= 02000; exponent1 -= (exponent1 >> 11); exponent2 -= (exponent2 >> 11); exponent1 += exponent2; /* add exponents together for multiply */ /* ** post normalize if necessary */ if (norm && !(upper >> 47)) { upper = (upper << 1) | ((lower >> 47) & 1); exponent1 -= 1; } /* ** offset exponent by 48, since we ended up with a 96 bit product ** and are only returning the upper 48 bits. check for overflow ** values (biased values less than 0000 or greater than 3777). */ if (exponent1 > 01717) { return OVFL(sign1); } exponent1 += 48; exponent1 += 02000 + (exponent1 >> 11); if (exponent1 < 0) { return 0; } return ((((CpWord) exponent1) << 48) | upper) ^ sign1; } /*-------------------------------------------------------------------------- ** Purpose: Floating divide implemented using shift and subtract. ** ** Rounding divide is identical to floating divide, except ** that as the dividend gets shifted in, 1/3 is shifted in ** (1/3 is alternating bits: 25252525... octal). ** ** Parameters: Name Description. ** v1 First operand ** v2 Second operand ** doRound TRUE if rounding required, FALSE otherwise. ** ** Returns: Upper 48 bits and adjusted exponent for single precision. ** Lower 48 bits and adjusted exponent for double precision. ** **------------------------------------------------------------------------*/ CpWord floatDivide(CpWord v1, CpWord v2, bool doRound) { CpWord sign1; CpWord sign2; int round = 0; int exponent1; int exponent2; sign1 = SignX(v1, 60); sign2 = SignX(v2, 60); v1 ^= sign1; v2 ^= sign2; exponent1 = (int)(v1 >> 48); exponent2 = (int)(v2 >> 48); sign1 ^= sign2; /* ** indefinite divided by anything is indefinite ** anything divided by indefinite is indefinite ** infinite divided by infinite is indefinite ** infinite divided by anything else is infinite */ if (!IR(exponent1)) { if ((exponent1 == ID) || (exponent2 == ID) || (exponent2 == OR)) { return IND; } return OVFL(sign1); } if (!IR(exponent2)) { if (exponent2 == ID) { return IND; } return 0; } /* ** exponent = 0 is taken to mean value = 0 ** if non-zero divided by zero, return overflow ** if zero divided by non-zero, return positive zero ** if zero divided by zero, return positive indefinite */ if (!(exponent1 && exponent2)) { if (exponent1) { return OVFL(sign1); } if (exponent2) { return 0; } return IND; } v1 &= Mask48; v2 &= Mask48; /* ** if divisor is less than half of dividend, return indefinite - divisor ** should be normalized, but it isn't checked for explicitly. */ if (v1 >= (v2 << 1)) return IND; exponent1 -= 02000; exponent2 -= 02000; exponent1 -= (exponent1 >> 11); exponent2 -= (exponent2 >> 11); /* ** divide exponents by subtracting */ exponent1 -= exponent2; /* ** pre-normalize if necessary. this is guaranteed to make v1 >= v2 ** due to earlier check. */ if (v1 < v2) { v1 <<= 1; exponent1 -= 1; if (doRound) { round = 1; /* round bit (of zero) got shifted in */ } } /* ** figure out final exponent and check for overflow before ** actually doing the divides */ if (exponent1 > 02056) /* 1777 + 0057 (octal) */ { return OVFL(sign1); } exponent1 -= 47; exponent1 += 02000 + (exponent1 >> 11); if (exponent1 < 0) { return 0; } sign2 = 0; /* used to accumulate the result */ /* ** main divide loop - shift and subtract for 48 bits */ for (exponent2 = 47; exponent2 >= 0; exponent2--) { sign2 <<= 1; if (v1 >= v2) { v1 -= v2; sign2 += 1; } if (doRound) { v1 = (v1 << 1) | round; /* shift in rounding bit */ round = 1 - round; /* toggle round back and forth */ } else { v1 <<= 1; } } return ((((CpWord) exponent1) << 48) | sign2) ^ sign1; } /* **-------------------------------------------------------------------------- ** ** Private Functions ** **-------------------------------------------------------------------------- */ /*--------------------------- End Of File ------------------------------*/