/*-- Float.c --*/

/*
Based on the SoftFloat IEC/IEEE Floating-Point Arithmetic
Package, Release 2b.
SoftFloat written by John R. Hauser.
Porting to C4Script by Nicolas Hake.

SoftFloat license and disclaimer follows.

============================================================================

This C source fragment is part of the SoftFloat IEC/IEEE Floating-point
Arithmetic Package, Release 2b.

Written by John R. Hauser.  This work was made possible in part by the
International Computer Science Institute, located at Suite 600, 1947 Center
Street, Berkeley, California 94704.  Funding was partially provided by the
National Science Foundation under grant MIP-9311980.  The original version
of this code was written as part of a project to build a fixed-point vector
processor in collaboration with the University of California at Berkeley,
overseen by Profs. Nelson Morgan and John Wawrzynek.  More information
is available through the Web page `http://www.cs.berkeley.edu/~jhauser/
arithmetic/SoftFloat.html'.

THIS SOFTWARE IS DISTRIBUTED AS IS, FOR FREE.  Although reasonable effort has
been made to avoid it, THIS SOFTWARE MAY CONTAIN FAULTS THAT WILL AT TIMES
RESULT IN INCORRECT BEHAVIOR.  USE OF THIS SOFTWARE IS RESTRICTED TO PERSONS
AND ORGANIZATIONS WHO CAN AND WILL TAKE FULL RESPONSIBILITY FOR ALL LOSSES,
COSTS, OR OTHER PROBLEMS THEY INCUR DUE TO THE SOFTWARE, AND WHO FURTHERMORE
EFFECTIVELY INDEMNIFY JOHN HAUSER AND THE INTERNATIONAL COMPUTER SCIENCE
INSTITUTE (possibly via similar legal warning) AGAINST ALL LOSSES, COSTS, OR
OTHER PROBLEMS INCURRED BY THEIR CUSTOMERS AND CLIENTS DUE TO THE SOFTWARE.

Derivative works are acceptable, even for commercial purposes, so long as
(1) the source code for the derivative work includes prominent notice that
the work is derivative, and (2) the source code includes prominent notice with
these four paragraphs for those parts of this code that are retained.

=============================================================================

End of SoftFloat license and disclaimer.
*/

#strict

// Rounding modes.
static const fRoundToNearestEven = 0;
static const fRoundDown = 1;
static const fRoundUp = 2;
static const fRoundToZero = 3;

// Exception codes.
static const fFlagInvalid = 1;
static const fFlagDivZero = 4;
static const fFlagOverflow = 8;
static const fFlagUnderflow = 16;
static const fFlagInexact = 32;

// 1. Infty is returned on x / +-0, if x != 0.
// Infty may also be returned on overflowing or underflowing arithmetic operations.

// 2. The following operations return a QNaN:
// a) any operation involving a QNaN
// b) +-0 / +-0
// c) +-Infty / +-Infty
// d) +Infty + (-Infty)
// e) +-0 * +-Infty
// f) sqrt(x), if x < 0.

// Representation of FP numbers:
// seee eeee emmm mmmm mmmm mmmm mmmm mmmm
// where s is the sign bit (set means "negative") (bit 31 = 1 bit)
// e are the exponent bits (bits 30 through 23 = 8 bits)
// m are the mantissa bits (bits 22 through 0 = 23 bits)

/** Converts an integer to a float.
 * Flags:
 *	This function does not set any flags.
 */
global func IntToFloat(int number) {
	// Quick check: If number is 0, no calculation is needed.
	if(number == 0) return 0;
	if(number == 0x80000000) return __fJoinSEM(1, 0x9E, 0);
	
	return __fNormRoundPack(number < 0, 0x9C, Abs(number));
}

/** Converts a float to an integer.
 * Flags:
 *	fFlagInvalid	Set if the number is too large to be represented as 32bit integer.
 *	fFlagInexact	Set if the result was rounded.
 */
global func FloatToInt(int number) {
	var s, e, m, z, mExtra;
	__fSplitSEM(number, s, e, m);
	var shift = e - 0x96;
	
	if(shift >= 0) {
		if(e >= 0x9E) {
			if(number != 0xCF000000) {
				__fFlags |= fFlagInvalid;
				if(!s || ((e == 0xFF) && m)) {
					return 0x7FFFFFFF;
				}
			}
			return 0x80000000;
		}
		z = (m |0x800000) << shift;
		if(s) z = -z;
	} else {
		if(e < 0x7E) {
			mExtra = e | m;
		} else {
			m |= 0x800000;
			mExtra = m << (shift & 31);
			z = __fShr(m, -shift);
		}
		
		if(mExtra) __fFlags |= fFlagInexact;
		
		if(__fRound == fRoundToNearestEven) {
			if(mExtra < 0) {
				++z;
				if((mExtra << 1) == 0) z &= ~1;
			}
			if(s) z = -z;
		} else {
			mExtra = !!mExtra;
			if(s) {
				z += (__fRound == fRoundDown) && mExtra;
				z = -z;
			} else {
				z += (__fRound == fRoundUp) && mExtra;
			}
		}
	}
	return z;
}

/** Returns the absolute value of a floating-point number.
 * Flags:
 *	This function does not set any flags.
 */
global func fAbs(int number) {
	return number & ~__fSignMask;
}

/** Returns the negated value of a floating-point number.
 * Flags:
 *	This function does not set any flags.
 */
global func fNeg(int number) {
	// Toggle sign bit
	return number ^ __fSignMask;
}

/** Returns true if the argument is a NaN.
 * Flags:
 *	This function does not set any flags.
 */
global func fIsNaN(int number) {
	// NaNs have all exponent bits and at least one mantissa bit set.
	if((number & __fExpMask) != __fExpMask) return false;
	return !!(number & __fMantMask);
}

/** Returns true if the argument is a signalling NaN.
 * Flags:
 *	This function does not set any flags.
 */
global func fIsSNaN(int number) {
	// SNaNs have all exponent bits set, the MSB of the mantissa cleared, and at least one other mantissa bit set.
	return ((number & 0x7FC00000 == 0x7F800000) && (number & __fMantMask));
}

/** Adds two floating-point numbers.
 * Flags:
 *	fFlagInvalid	Either one of the arguments was a signalling NaN, or the arguments were infinity values with different signs.
 *	fFlagInexact	The result was rounded.
 *	fFlagOverflow	The result was too big to be stored in a float.
 *	fFlagUnderflow	The result was too small to be stored in a float.
 */
global func fAdd(int a, int b) {
	var as = !!(a & __fSignMask);
	var bs = !!(b & __fSignMask);
	if(as == bs) {
		return __fAddSgn(a, b, as);
	} else {
		return __fSubSgn(a, b, as);
	}
}

/** Subtracts two floating-point numbers.
 * Flags:
 *	fFlagInvalid	Either one of the arguments was a signalling NaN, or both of the arguments were infinity values.
 *	fFlagInexact	The result was rounded.
 *	fFlagOverflow	The result was too big to be stored in a float.
 *	fFlagUnderflow	The result was too small to be stored in a float.
 */
global func fSub(int a, int b) {
	var as = !!(a & __fSignMask);
	var bs = !!(b & __fSignMask);
	if(as == bs) {
		return __fSubSgn(a, b, as);
	} else {
		return __fAddSgn(a, b, as);
	}
}

/** Multiplies two floating-point numbers.
 * Flags:
 *	fFlagInvalid	Either one of the arguments was a signalling NaN, or the arguments were Infinity and zero.
 *	fFlagInexact	The result was rounded.
 *	fFlagOverflow	The result was too big to be stored in a float.
 *	fFlagUnderflow	The result was too small to be stored in a float.
 */
global func fMul(int a, int b) {
	var as, ae, am;
	var bs, be, bm;
	var zs, ze, zm;
	
	__fSplitSEM(a, as, ae, am);
	__fSplitSEM(b, bs, be, bm);
	zs = as ^ bs;
	
	// Infinity/NaN
	if(ae == 0xFF) { 
		if(am || ((be == 0xFF) && bm))
			return __fPropNaN(a, b);
		else if(!be && !bs) {
			__fFlags |= fFlagInvalid;
			return __fQNaN;
		}
		return __fJoinSEM(zs, 0xFF, 0);
	}
	if(be == 0xFF) { 
		if(bm)
			return __fPropNaN(a, b);
		else if(!ae && !as) {
			__fFlags |= fFlagInvalid;
			return __fQNaN;
		}
		return __fJoinSEM(zs, 0xFF, 0);
	}
	
	if(ae == 0) {
		if(am == 0) return __fJoinSEM(zs, 0, 0);
		// Normalize
		__fSplitSEM(__fNormalizeSubnormal(a), as, ae, am);
	}
	if(be == 0) {
		if(bm == 0) return __fJoinSEM(zs, 0, 0);
		__fSplitSEM(__fNormalizeSubnormal(b), bs, be, bm);
	}
	
	ze = ae + be - 0x7F;
	am = (am | 0x00800000) << 7;
	bm = (bm | 0x00800000) << 8;
	
	var z0, z1;
	__fMul64(am, bm, z0, z1);
	z0 |= !!z1;
	if((z0 << 1) >= 0) {
		z0 = z0 << 1;
		--ze;
	}
	
	return __fRoundPack(zs, ze, z0);
}

/** Divides two floating-point numbers.
 * Flags:
 *	fFlagInvalid	Either one of the arguments was a signalling NaN, or both of the arguments were Infinity, or both of the arguments were zero.
 *	fFlagInexact	The result was rounded.
 *	fFlagOverflow	The result was too big to be stored in a float.
 *	fFlagUnderflow	The result was too small to be stored in a float.
 *	fFlagDivZero	The divisor was zero.
 */
global func fDiv(int a, int b) {
	var as, ae, am;
	var bs, be, bm;
	__fSplitSEM(a, as, ae, am);
	__fSplitSEM(b, bs, be, bm);
	
	var zs, ze, zm;
	zs = as ^ bs;
	
	// Handle Infty and NaNs
	if(ae == 0xFF) {
		if(am) return __fPropNaN(a, b);
		else if(be == 0xFF) {
			if(bm) return __fPropNaN(a, b);
			__fFlags |= fFlagInvalid;
			return __fQNaN;
		}
		return __fJoinSEM(zs, 0xFF, 0);
	} else if(be == 0xFF) {
		if(bm) return __fPropNaN(a, b);
		return __fJoinSEM(zs, 0, 0);
	}
	
	if(be == 0) {
		if(bm == 0) {
			if(!ae && !am) {
				__fFlags |= fFlagInvalid;
				return __fQNaN;
			}
			__fFlags |= fFlagDivZero;
			return __fJoinSEM(zs, 0xFF, 0);
		}
		__fSplitSEM(__fNormalizeSubnormal(b), bs, be, bm);
	}
	
	if(ae == 0) {
		if(am == 0) return __fJoinSEM(zs, 0, 0);
		__fSplitSEM(__fNormalizeSubnormal(a), as, ae, am);
	}
	
//	Log("a: %d/%d/%d", as, ae, am);
//	Log("b: %d/%d/%d", bs, be, bm);
	
	ze = ae - be + 0x7D;
	am = (am | 0x800000) << 7;
	bm = (bm | 0x800000) << 8;
	
	if(__fUnsignedLE(bm, (am << 1))) {
		am = __fShr(am, 1);
		++ze;
	}
	
	zm = __fEstDiv64(am, 0, bm);
	if((zm & 0x3F) <= 2) {
		var t0, t1, r0, r1;
		__fMul64(bm, zm, t0, t1);
		__fSub64(am, 0, t0, t1, r0, r1);
		while(r0 < 0) {
			--zm;
			__fAdd64(r0, r1, 0, bm, r0, r1);
		}
		zm |= !!r1;
	}
	
//	Log("fDiv => %d/%d/%d", zs, ze, zm);
	return __fRoundPack(zs, ze, zm);
}

/** Returns the remainder of a division of two floating-point numbers.
 * Flags:
 *	fFlagInvalid	Either one of the arguments was a signalling NaN, or both of the arguments were Infinity, or the second argument was zero.
 *	fFlagInexact	The result was rounded.
 *	fFlagOverflow	The result was too big to be stored in a float.
 *	fFlagUnderflow	The result was too small to be stored in a float.
 *	fFlagDivZero	The divisor was zero.
 */
global func fRem(int a, int b) {
	var as, ae, am;
	var bs, be, bm;
	__fSplitSEM(a, as, ae, am);
	__fSplitSEM(b, bs, be, bm);
	
	// Handle Infty and NaNs
	if(ae == 0xFF) {
		if(am || ((be == 0xFF) && bm)) {
			return __fPropNaN(a, b);
		}
		__fFlags |= fFlagInvalid;
		return __fQNaN;
	}
	if(be == 0xFF) {
		if(bm) return __fPropNaN(a, b);
		return a;
	}
	if(be == 0) {
		if(bm == 0) {
			__fFlags |= fFlagInvalid;
			return __fQNaN;
		}
		__fSplitSEM(__fNormalizeSubnormal(b), bs, be, bm);
	}
	if(ae == 0) {
		if(am == 0) return a;
		__fSplitSEM(__fNormalizeSubnormal(a), as, ae, am);
	}
	
	var de = ae - be;
	am = (am | 0x800000) << 8;
	bm = (bm | 0x800000) << 8;
	if(de < 0) {
		if(de < -1) return a;
		am = __fShr(am, 1);
	}
	
	var q = __fUnsignedLE(bm, am);
	if(q) am -= bm;
	de -= 32;
	
	while(de > 0) {
		q = __fEstDiv64(am, 0, bm);
		if(q > 2) q -= 2;
		else q = 0;
		am = -__fUnsignedIntMul(__fShr(bm, 2), q);
		de -= 30;
	}
	de += 32;
	
	if(de > 0) {
		q = __fEstDiv64(am, 0, bm);
		if(q > 2) q -= 2;
		else q = 0;
		q = __fShr(q, 32 - de);
		bm = __fShr(bm, 2);
		am = __fUnsignedIntMul((__fShr(am, 1) << (de - 1)) - bm, q);
	} else {
		am = __fShr(am, 2);
		bm = __fShr(bm, 2);
	}
	
	var altAm;
	while(true) {
		altAm = am;
		++q;
		am -= bm;
		if(am < 0) break;
	}
	
	var mMean = am + altAm;
	if((mMean < 0) || ((mMean == 0) && (q & 1))) {
		am = altAm;
	}
	
	if(am < 0) {
		am = -am;
		as = !as;
	}
	return __fNormRoundPack(as, be, am);
}

/** Returns the square root of a floating-point number.
 * Flags:
 *	fFlagInvalid	Either the argument was a signalling NaN, or the argument was less than zero.
 *	fFlagInexact	The result was rounded.
 *	fFlagOverflow	The result was too big to be stored in a float.
 */
global func fSqrt(int number) {
	var as, ae, am;
	__fSplitSEM(number, as, ae, am);
	
	// Infty/NaNs
	if(ae == 0xFF) {
		if(am) return __fPropNaN(number, 0);
		if(!as) return number;
		__fFlags |= fFlagInvalid;
		return __fQNaN;
	}
	
	// sqrt(x < 0)
	if(as) {
		if(!(ae|am)) return number;
		__fFlags |= fFlagInvalid;
		return __fQNaN;
	}
	
	if(ae == 0) {
		if(am == 0) return 0;
		__fSplitSEM(__fNormalizeSubnormal(number), as, ae, am);
	}
	
	var ze = __fShr(ae - 0x7F, 1) + 0x7E;
	am = (am | 0x800000) << 8;
	var zm = __fEstSqrt(ae, am) + 2;
	if((zm & 0x7F) <= 5) {
		if(__fUnsignedLT(zm, 2)) {
			return __fRoundPack(0, ze, 0x7FFFFFFF);
		} else {
			if(ae & 1)
				am = __fShr(am, 1);
			var t0, t1, r0, r1;
			__fMul64(zm, zm, t0, t1);
			__fSub64(am, 0, t0, t1, r0, r1);
			while(r0 < 0) {
				--zm;
				__fShl64(0, zm, 1, t0, t1);
				t1 |= 1;
				__fAdd64(r0, r1, t0, t1, r0, r1);
			}
			if(r0||r1)
				zm |= 1;
		}
	}
	zm = __fShrJam(zm, 1);
	return __fRoundPack(0, ze, zm);
}

/** Rounds a float to an integral value.
 * Flags:
 *	fFlagInvalid	The argument was a signalling NaN.
 *	fFlagInexact	The result was rounded.
 */
global func fRoundToInt(int number) {
	var as, ae, am;
	__fSplitSEM(number, as, ae, am);
	
	if(ae >= 0x96) {
		if(ae == 0xFF && am) {
			return __fPropNaN(number, number);
		}
		return number;
	}
	if(ae <= 0x7E) {
		if(!(number << 1)) return number;
		__fFlags |= fFlagInexact;
		if(__fRound == fRoundToNearestEven) {
			if(ae == 0x7E && am)
				return __fJoinSEM(as, 0x7F, 0);
		} else if(__fRound == fRoundDown) {
			if(as) return 0xBF800000;
			else return 0;
		} else if(__fRound == fRoundUp) {
			if(as) return 0x80000000;
			else return 0x3F800000;
		}
		return __fJoinSEM(as, 0, 0);
	}
	
	var lastBM = 1 << (0x96 - ae);
	var roundBM = lastBM - 1;
	var z = number;
	
	if(__fRound == fRoundToNearestEven) {
		z += __fShr(lastBM, 1);
		if(!(z & roundBM)) z &= ~lastBM;
	} else if(__fRound != fRoundToZero) {
		if(as ^ (__fRound == fRoundUp)) {
			z += roundBM;
		}
	}
	
	z &= ~roundBM;
	if(z != number) __fFlags |= fFlagInexact;
	return z;
}

/** Returns true if both arguments are equal.
 * Note that this is different than a bitwise compare. Any NaN is not equal to any other value, even other NaNs with the same representation. Also, -0 and +0 are equal, although they differ in sign.
 * Flags:
 *	fFlagInvalid	One of the arguments was a signalling NaN.
 */
global func fEq(int a, int b) {
	if(((a & __fExpMask) == __fExpMask && (a & __fMantMask))
		|| ((b & __fExpMask) == __fExpMask && (b & __fMantMask)))
	{
		if(fIsSNaN(a) || fIsSNaN(b))
			__fFlags |= fFlagInvalid;
		return false;
	}
	
	return ((a == b) || ((a | b) << 1) == 0);
}

/** Returns true if a is less than or equal to b.
 * Flags:
 *	fFlagInvalid	One of the arguments was a NaN.
 */
global func fLe(int a, int b) {
	var as, ae, am;
	var bs, be, bm;
	__fSplitSEM(a, as, ae, am);
	__fSplitSEM(b, bs, be, bm);

	if((ae == 0xFF && !am) || (be == 0xFF && !bm)) {
		__fFlags |= fFlagInvalid;
		return false;
	}
	
	if(as != bs)
		return as || ((a | b) << 1) == 0;
	return ((a == b) || (as ^ __fUnsignedLT(a, b)));
}

/** Returns true if a is less than b.
 * Flags:
 *	fFlagInvalid	One of the arguments was a NaN.
 */
global func fLt(int a, int b) {
	var as, ae, am;
	var bs, be, bm;
	__fSplitSEM(a, as, ae, am);
	__fSplitSEM(b, bs, be, bm);

	if((ae == 0xFF && !am) || (be == 0xFF && !bm)) {
		__fFlags |= fFlagInvalid;
		return false;
	}
	
	if(as != bs)
		return as && ((a | b) << 1);
	return ((a != b) && (as ^ __fUnsignedLT(a, b)));
}

/** Returns the error flags of the FPU.
 */
global func fFlags() { return __fFlags; }

/** Clears the error flags of the FPU.
 */
global func fClearFlags() { __fFlags = 0; }

/** Returns the rounding mode of the FPU.
 */
global func fRoundingMode() { return __fRound; }

/** Sets the rounding mode of the FPU.
 */
global func fSetRoundingMode(int round) { __fRound = round; }

/* === END OF INTERFACE === */

global func __fUnsignedLE(int a, int b) {
	if(a == b) return true;
	return __fUnsignedLT(a, b);
}

global func __fUnsignedLT(int a, int b) {
	if((a & 0x80000000) && !(b & 0x80000000)) return false;
	if((b & 0x80000000) && !(a & 0x80000000)) return true;
	return (a & 0x7FFFFFFF) < (b & 0x7FFFFFFF);
}

global func __fAdd64(int a0, int a1, int b0, int b1, &out0, &out1) {
	var z = a1 + b1;
	out1 = z;
	out0 = a0 + b0 + (z < a1);
}

global func __fSub64(int a0, int a1, int b0, int b1, &out0, &out1) {
	out1 = a1 - b1;
	out0 = a0 - b0 - (a1 < b1);
}

global func __fEstDiv64(int a0, int a1, int b) {
	var b0, b1;
	var r0, r1, t0, t1;
	var z;
	
	if(__fUnsignedLE(b, a0)) {
//		Log("__fEstDiv64(%d, %d, %d) = %d", a0, a1, b, 0xFFFFFFFF);
		return 0xFFFFFFFF;
	}
	b0 = __fShr(b, 16);
	if(__fUnsignedLE(b0 << 16, a0))
		z = 0xFFFF0000;
	else
		z = __fUnsignedIntDiv(a0, b0) << 16;
	
	__fMul64(b, z, t0, t1);
	__fSub64(a0, a1, t0, t1, r0, r1);
	while(r0 < 0) {
		z -= 0x10000;
		b1 = b << 16;
		__fAdd64(r0, r1, b0, b1, r0, r1);
	}
	
	r0 = (r0 << 16) | __fShr(r1, 16);
	if(__fUnsignedLE(b0 << 16, r0))
		z |= 0xFFFF;
	else
		z |= __fUnsignedIntDiv(r0, b0);
	
//	Log("__fEstDiv64(%d, %d, %d) = %d", a0, a1, b, z);
	return z;
}

global func __fUnsignedIntDiv(int a, int b) {
	if(!((a|b)&0x80000000)) return a/b;
	if(__fUnsignedLT(a, b)) return 0;
	
	var rem, res;
	for(var ctr = 32; ctr; --ctr) {
		rem = rem << 1;
		res = res << 1;
		rem |= !!(a&0x80000000);
		a = a << 1;
		if(__fUnsignedLE(b, rem)) {
			rem -= b;
			res |= 1;
		}
	}
	return res;
}

global func __fShl64(int a0, int a1, int count, &out0, &out1) {
	out1 = a1 << count;
	if(count == 0) out0 = a0;
	else out0 = (a0 << count) | __fShr(a1, (-count) & 31);
}

global func __fEstSqrt(int ae, int am) {
	var oddAdjust = [0x0004, 0x0022, 0x005D, 0x00B1, 0x011D, 0x019F, 0x0236, 0x02E0, 0x039C, 0x0468, 0x0545, 0x0631, 0x072B, 0x0832, 0x0946, 0x0A67];
	var evenAdjust = [0x0A2D, 0x08AF, 0x075A, 0x0629, 0x051A, 0x0429, 0x0356, 0x029E, 0x0200, 0x0179, 0x0109, 0x00AF, 0x0068, 0x0034, 0x0012, 0x0002];
	
//	Log("__fES: in: %d", FloatToInt(__fJoinSEM(0, ae, am)));
	
	var idx = (am >> 27) & 15;
	var z;
	if(ae & 1) {
//		Log("Odd path");
		z = 0x4000 + __fShr(am, 17) - oddAdjust[idx];
		z = (__fUnsignedIntDiv(am, z) << 14) + (z << 15);
		am = __fShr(am, 1);
	} else {
//		Log("Even path");
		z = 0x8000 + __fShr(am, 17) - evenAdjust[idx];
		z += __fUnsignedIntDiv(am, z);
		if(__fUnsignedLE(0x20000, z))
			z = 0xFFFF8000;
		else
			z = z << 15;
		if(__fUnsignedLE(z, am))
			return __fShr(am, 1);
	}
	
//	Log("out: %d", __fShr(__fEstDiv64(am, 0, z), 1) + __fShr(z, 1));
	return __fShr(__fEstDiv64(am, 0, z), 1) + __fShr(z, 1);
}

global func __fPrintFlags() {
	var f;
	if(__fFlags & fFlagInvalid)
		f = "?";
	else
		f = ".";
	if(__fFlags & fFlagDivZero)
		f = Format("%s%s", f, "0");
	else
		f = Format("%s.", f);
	if(__fFlags & fFlagOverflow)
		f = Format("%s%s", f, ">");
	else
		f = Format("%s.", f);
	if(__fFlags & fFlagUnderflow)
		f = Format("%s%s", f, "<");
	else
		f = Format("%s.", f);
	if(__fFlags & fFlagInexact)
		f = Format("%s%s", f, "~");
	else
		f = Format("%s.", f);
	
	Log("%s", f);
}

global func __fShr(int number, int shift) {
	// Need to work around sign extension.
	number = (number >> 1) & ~0x80000000;
	number = number >> (shift - 1);
	return number;
}

global func __fShrJam(int number, int shift) {
	if(!shift) return number;
	if(shift < 32) {
		return ((__fShr(number, shift) | !!(number << ((-shift) & 31))));
	}
	return 0+!!number;
}

global func __fCountZeros(int number) {
	var zeros;
	if(!number) return 32;
	for(zeros = 0; !(number & 0x80000000); ++zeros)
		number = number << 1;
	return zeros;
}

global func __fNormalizeSubnormal(int number) {
	var s, e, m;
	__fSplitSEM(number, s, e, m);
	var shift = __fCountZeros(m) - 8;
	m = m << shift;
	e = 1 - shift;
	return __fJoinSEM(s, e, m);
}

global func __fSplitSEM(int in, &s, &e, &m) {
	s = !!(in & __fSignMask);
	e = (in & __fExpMask) >> 23;
	m = in & __fMantMask;
}

global func __fJoinSEM(int s, int e, int m) {
	return (s << 31) + (e << 23) + m;
}

global func __fRoundPack(int s, int e, int m) {
	var rInc, rBits;
	var tiny;
	
	rInc = 0x40;
	if(__fRound != fRoundToNearestEven) {
		if(__fRound == fRoundToZero)
			rInc = 0;
		else {
			rInc = 0x7F;
			if(s) {
				if(__fRound == fRoundUp) rInc = 0;
			} else
				if(__fRound == fRoundDown) rInc = 0;
		}
	}
	
	rBits = m & 0x7F;
	if(e >= 0xFD) {
		if((e > 0xFD) || ((e == 0xFD) && ((m + rInc) < 0))) {
			__fFlags |= fFlagOverflow | fFlagInexact;
			// Return +-(Infty - 1) when rounding towards zero, +-Infty otherwise
			return __fJoinSEM(s, 0xFF, 0) - !!rInc;
		}
		if(e < 0) {
			tiny = (e < -1) || (m + rInc < 0x80000000);
			m = __fShrJam(m, -e);
			e = 0;
			rBits = m & 0x7F;
			if(tiny && rBits) __fFlags |= fFlagUnderflow;
		}
	}
	
	if(rBits) __fFlags |= fFlagInexact;
	m = (m + rInc) >> 7;
	m &= ~(((rBits ^ 0x40) == 0) & (__fRound == fRoundToNearestEven));
	if(!m) e = 0;
	
	return __fJoinSEM(s, e, m);
}

global func __fNormRoundPack(int s, int e, int m) {
	var shift = __fCountZeros(m) - 1;
	return __fRoundPack(s, e - shift, m << shift);
}

global func __fPropNaN(int a, int b) {
	var aNaN, aSNaN, bNaN, bSNaN;
	aNaN = fIsNaN(a);
	aSNaN = fIsSNaN(a);
	bNaN = fIsNaN(b);
	bSNaN = fIsSNaN(b);
	a |= 0x00400000;
	b |= 0x00400000;
	
	// SNaNs always raise the invalid flag.
	if(aSNaN || bSNaN) __fFlags |= fFlagInvalid;
	if(aSNaN) {
		if(!bSNaN) {
			// If any operand is a NaN, this is returned. QNaNs override SNaNs.
			if(bNaN) return b;
			return a;
		}
	} else if(aNaN) {
		if(bSNaN || !bNaN) return a;
	} else {
		return b;
	}
	if((a << 1) < (b << 1)) return b;
	if((a << 1) < (b << 1)) return a;
	return Min(a, b);
}

global func __fAddSgn(int a, int b, int sgn) {
	// Find out which number has the larger exponent.
	var as, ae, am, bs, be, bm;
	__fSplitSEM(a, as, ae, am);
	__fSplitSEM(b, bs, be, bm);
	
	var de = ae - be;
	var ze, zm;
	am = am << 6;
	bm = bm << 6;
	
	if(de > 0) {
		if(ae == 0xFF) {
			if(am) return __fPropNaN(a, b);
			return a;
		}
		if(be == 0) {
			--de;
		} else {
			bm |= 0x20000000;
		}
		bm = __fShrJam(bm, de);
		ze = ae;
	} else if(de < 0) {
		if(be == 0xFF) {
			if(bm) return __fPropNaN(a, b);
			return __fJoinSEM(sgn, 0xFF, 0);
		}
		if(ae == 0) {
			++de;
		} else {
			am |= 0x20000000;
		}
		am = __fShrJam(am, -de);
		ze = be;
	} else {
		if(ae == 0xFF) {
			if(am || bm) return __fPropNaN(a, b);
			return a;
		}
		if(ae == 0) return __fJoinSEM(sgn, 0, __fShr(am + bm, 6));
		zm = 0x40000000 + am + bm;
		ze = ae;
		return __fRoundPack(sgn, ze, zm);
	}
	
	am |= 0x20000000;
	zm = (am + bm) << 1;
	--ze;
	if(zm < 0) {
		zm = am + bm;
		++ze;
	}
	return __fRoundPack(sgn, ze, zm);
}

global func __fSubSgn(int a, int b, int sgn) {
	var as, ae, am, bs, be, bm;
	__fSplitSEM(a, as, ae, am);
	__fSplitSEM(b, bs, be, bm);
	
	var de = ae - be;
	var ze, zm;
	am = am << 7;
	bm = bm << 7;
	
	if(de == 0) {
		if(ae == 0xFF) {
			if(am || bm) return __fPropNaN(a, b);
			__fFlags |= fFlagInvalid;
			return __fQNaN;
		}
		if(ae == 0) {
			ae = be = 1;
		}
		if(am == bm) return __fJoinSEM(__fRound == fRoundDown, 0, 0);
		
		else if(am > bm) {
			return __fNormRoundPack(sgn, ae - 1, am - bm);
		}
		else /* if(am < bm) */ {
			return __fNormRoundPack(!sgn, be - 1, bm - am);
		}
	} else if(de > 0) {
		if(ae == 0xFF) {
			if(am) return __fPropNaN(a, b);
			return a;
		} else if(be == 0) {
			--de;
		} else {
			bm |= 0x40000000;
		}
		bm = __fShrJam(bm, de);
		am |= 0x40000000;
		return __fNormRoundPack(sgn, ae - 1, am - bm);
	} else /* if(de < 0) */ {
		if(be == 0xFF) {
			if(bm) return __fPropNaN(a, b);
			return __fJoinSEM(!sgn, 0xFF, 0);
		}
		else if(ae == 0) {
			++de;
		} else {
			am |= 0x40000000;
		}
		am = __fShrJam(am, -de);
		bm |= 0x40000000;
		return __fNormRoundPack(!sgn, be - 1, bm - am);
	}
}

global func __fMul64(int a, int b, &out0, &out1) {
	var al, ah, bl, bh;
	var zMid0, zMid1;
	
	al = a & 0xFFFF;
	ah = __fShr(a, 16);
	bl = b & 0xFFFF;
	bh = __fShr(b, 16);
	out1 = al * bl;
	zMid0 = al * bh;
	zMid1 = ah * bl;
	out0 = ah * bh;
	zMid0 += zMid1;
	out0 += ((zMid0 < zMid1) << 16) + (__fShr(zMid0, 16));
	zMid0 = zMid0 << 16;
	out1 += zMid0;
	out0 += (out1 < zMid0);
}

global func __fUnsignedIntMul(int a, int b) {
	var o0, o1;
	__fMul64(a, b, o0, o1);
	return o1;
}

// Some special numbers.
static const __fSignMask = 0x80000000;	// Sign bit mask
static const __fExpMask = 0x7F800000;	// Exponent bits mask
static const __fMantMask = 0x7FFFFF;	// Mantissa bits mask
static const __fQNaN = 0x7FC00000;	// Quiet NaN

static __fRound;	// Defaults to fRoundToNearestEven
static __fFlags;