+Fix32

8 files changed, 2031 insertions, 7 deletions

diff --git a/Client/ThirdParty/fix32/fix32.cpp b/Client/ThirdParty/fix32/fix32.cpp
deleted file mode 100644
index e69de29..0000000
--- a/Client/ThirdParty/fix32/fix32.cpp
+++ /dev/null
diff --git a/Client/ThirdParty/fix32/fix32.h b/Client/ThirdParty/fix32/fix32.h
deleted file mode 100644
index fb3d65c..0000000
--- a/Client/ThirdParty/fix32/fix32.h
+++ /dev/null
@@ -1,7 +0,0 @@
-#pragma once

-

-#include <stdint.h>

-

-// -2,147,483,648 to 2,147,483,647

-typedef int64_t fix32_t;

-

diff --git a/Client/ThirdParty/fix32/fix32.hpp b/Client/ThirdParty/fix32/fix32.hpp
new file mode 100644
index 0000000..7999ad4
--- /dev/null
+++ b/Client/ThirdParty/fix32/fix32.hpp
@@ -0,0 +1,161 @@
+#pragma once

+

+#include <stdint.h>

+extern "C" {

+#include "math-sll/math-sll.h"

+}

+

+typedef sll fixed32_t;

+

+// Q32.32

+class Fix32

+{

+public:

+    sll value;

+

+    inline Fix32() { value = 0; }

+    inline Fix32(const Fix32 &inValue) { value = inValue.value; }

+    inline Fix32(const float inValue) { value = dbl2sll(inValue); }

+    inline Fix32(const double inValue) { value = dbl2sll(inValue); }

+    inline Fix32(const int32_t inValue) { value = int2sll(inValue); }

+    inline Fix32(const sll inValue) { value = inValue; }

+

+    inline operator sll() const { return value; }

+    inline operator double()  const { return sll2dbl(value); }

+    inline operator float()   const { return (float)sll2dbl(value); }

+    inline operator int32_t() const { return (int32_t)sll2int(value); }

+

+    inline Fix32 & operator=(const Fix32 &rhs) { value = rhs.value;             return *this; }

+    inline Fix32 & operator=(const sll rhs) { value = rhs;                   return *this; }

+    inline Fix32 & operator=(const double rhs) { value = dbl2sll(rhs);   return *this; }

+    inline Fix32 & operator=(const float rhs) { value = (float)dbl2sll(rhs); return *this; }

+    inline Fix32 & operator=(const int32_t rhs) { value = int2sll(rhs);   return *this; }

+

+    inline Fix32 & operator+=(const Fix32 &rhs) { value += rhs.value;             return *this; }

+    inline Fix32 & operator+=(const sll rhs) { value += rhs;                   return *this; }

+    inline Fix32 & operator+=(const double rhs) { value += dbl2sll(rhs);   return *this; }

+    inline Fix32 & operator+=(const float rhs) { value += (float)dbl2sll(rhs); return *this; }

+    inline Fix32 & operator+=(const int32_t rhs) { value += int2sll(rhs);   return *this; }

+

+    inline Fix32 & operator-=(const Fix32 &rhs) { value -= rhs.value; return *this; }

+    inline Fix32 & operator-=(const sll rhs) { value -= rhs; return *this; }

+    inline Fix32 & operator-=(const double rhs) { value -= dbl2sll(rhs); return *this; }

+    inline Fix32 & operator-=(const float rhs) { value -= (float)dbl2sll(rhs); return *this; }

+    inline Fix32 & operator-=(const int32_t rhs) { value -= int2sll(rhs); return *this; }

+

+    inline Fix32 & operator*=(const Fix32 &rhs) { value = sllmul(value, rhs.value); return *this; }

+    inline Fix32 & operator*=(const sll rhs) { value = sllmul(value, rhs); return *this; }

+    inline Fix32 & operator*=(const double rhs) { value = sllmul(value, dbl2sll(rhs)); return *this; }

+    inline Fix32 & operator*=(const float rhs) { value = sllmul(value, (float)dbl2sll(rhs)); return *this; }

+    inline Fix32 & operator*=(const int32_t rhs) { value *= rhs; return *this; }

+

+    inline Fix32 & operator/=(const Fix32 &rhs) { value = slldiv(value, rhs.value); return *this; }

+    inline Fix32 & operator/=(const sll rhs) { value = slldiv(value, rhs); return *this; }

+    inline Fix32 & operator/=(const double rhs) { value = slldiv(value, dbl2sll(rhs)); return *this; }

+    inline Fix32 & operator/=(const float rhs) { value = slldiv(value, (float)dbl2sll(rhs)); return *this; }

+    inline Fix32 & operator/=(const int32_t rhs) { value /= rhs; return *this; }

+

+    inline const Fix32 operator+(const Fix32 &other) const { Fix32 ret = *this; ret += other; return ret; }

+    inline const Fix32 operator+(const sll other) const { Fix32 ret = *this; ret += other; return ret; }

+    inline const Fix32 operator+(const double other) const { Fix32 ret = *this; ret += other; return ret; }

+    inline const Fix32 operator+(const float other) const { Fix32 ret = *this; ret += other; return ret; }

+    inline const Fix32 operator+(const int32_t other) const { Fix32 ret = *this; ret += other; return ret; }

+

+#ifndef FIXMATH_NO_OVERFLOW

+    inline const Fix32 sadd(const Fix32 &other)  const { Fix32 ret = slladd(value, other.value);             return ret; }

+    inline const Fix32 sadd(const sll other) const { Fix32 ret = slladd(value, other);                   return ret; }

+    inline const Fix32 sadd(const double other)  const { Fix32 ret = slladd(value, dbl2sll(other));   return ret; }

+    inline const Fix32 sadd(const float other)   const { Fix32 ret = slladd(value, (float)dbl2sll(other)); return ret; }

+    inline const Fix32 sadd(const int32_t other) const { Fix32 ret = slladd(value, int2sll(other));   return ret; }

+#endif

+

+    inline const Fix32 operator-(const Fix32 &other) const { Fix32 ret = *this; ret -= other; return ret; }

+    inline const Fix32 operator-(const sll other) const { Fix32 ret = *this; ret -= other; return ret; }

+    inline const Fix32 operator-(const double other) const { Fix32 ret = *this; ret -= other; return ret; }

+    inline const Fix32 operator-(const float other) const { Fix32 ret = *this; ret -= other; return ret; }

+    inline const Fix32 operator-(const int32_t other) const { Fix32 ret = *this; ret -= other; return ret; }

+

+    inline const Fix32 operator-() const { Fix32 ret = sllsub(0, value); return ret; }

+

+#ifndef FIXMATH_NO_OVERFLOW

+    inline const Fix32 ssub(const Fix32 &other)  const { Fix32 ret = slladd(value, -other.value);             return ret; }

+    inline const Fix32 ssub(const sll other) const { Fix32 ret = slladd(value, -other);                   return ret; }

+    inline const Fix32 ssub(const double other)  const { Fix32 ret = slladd(value, -dbl2sll(other));   return ret; }

+    inline const Fix32 ssub(const float other)   const { Fix32 ret = slladd(value, -(float)dbl2sll(other)); return ret; }

+    inline const Fix32 ssub(const int32_t other) const { Fix32 ret = slladd(value, -int2sll(other));   return ret; }

+#endif

+

+    inline const Fix32 operator*(const Fix32 &other) const { Fix32 ret = *this; ret *= other; return ret; }

+    inline const Fix32 operator*(const sll other) const { Fix32 ret = *this; ret *= other; return ret; }

+    inline const Fix32 operator*(const double other) const { Fix32 ret = *this; ret *= other; return ret; }

+    inline const Fix32 operator*(const float other) const { Fix32 ret = *this; ret *= other; return ret; }

+    inline const Fix32 operator*(const int32_t other) const { Fix32 ret = *this; ret *= other; return ret; }

+

+#ifndef FIXMATH_NO_OVERFLOW

+    inline const Fix32 smul(const Fix32 &other)  const { Fix32 ret = sllmul(value, other.value);             return ret; }

+    inline const Fix32 smul(const sll other) const { Fix32 ret = sllmul(value, other);                   return ret; }

+    inline const Fix32 smul(const double other)  const { Fix32 ret = sllmul(value, dbl2sll(other));   return ret; }

+    inline const Fix32 smul(const float other)   const { Fix32 ret = sllmul(value, (float)dbl2sll(other)); return ret; }

+    inline const Fix32 smul(const int32_t other) const { Fix32 ret = sllmul(value, int2sll(other));   return ret; }

+#endif

+

+    inline const Fix32 operator/(const Fix32 &other) const { Fix32 ret = *this; ret /= other; return ret; }

+    inline const Fix32 operator/(const sll other) const { Fix32 ret = *this; ret /= other; return ret; }

+    inline const Fix32 operator/(const double other) const { Fix32 ret = *this; ret /= other; return ret; }

+    inline const Fix32 operator/(const float other) const { Fix32 ret = *this; ret /= other; return ret; }

+    inline const Fix32 operator/(const int32_t other) const { Fix32 ret = *this; ret /= other; return ret; }

+

+#ifndef FIXMATH_NO_OVERFLOW

+    inline const Fix32 sdiv(const Fix32 &other)  const { Fix32 ret = slldiv(value, other.value);             return ret; }

+    inline const Fix32 sdiv(const sll other) const { Fix32 ret = slldiv(value, other);                   return ret; }

+    inline const Fix32 sdiv(const double other)  const { Fix32 ret = slldiv(value, dbl2sll(other));   return ret; }

+    inline const Fix32 sdiv(const float other)   const { Fix32 ret = slldiv(value, (float)dbl2sll(other)); return ret; }

+    inline const Fix32 sdiv(const int32_t other) const { Fix32 ret = slldiv(value, int2sll(other));   return ret; }

+#endif

+

+    inline int operator==(const Fix32 &other)  const { return (value == other.value); }

+    inline int operator==(const sll other) const { return (value == other); }

+    inline int operator==(const double other)  const { return (value == dbl2sll(other)); }

+    inline int operator==(const float other)   const { return (value == (float)dbl2sll(other)); }

+    inline int operator==(const int32_t other) const { return (value == int2sll(other)); }

+

+    inline int operator!=(const Fix32 &other)  const { return (value != other.value); }

+    inline int operator!=(const sll other) const { return (value != other); }

+    inline int operator!=(const double other)  const { return (value != dbl2sll(other)); }

+    inline int operator!=(const float other)   const { return (value != (float)dbl2sll(other)); }

+    inline int operator!=(const int32_t other) const { return (value != int2sll(other)); }

+

+    inline int operator<=(const Fix32 &other)  const { return (value <= other.value); }

+    inline int operator<=(const sll other) const { return (value <= other); }

+    inline int operator<=(const double other)  const { return (value <= dbl2sll(other)); }

+    inline int operator<=(const float other)   const { return (value <= (float)dbl2sll(other)); }

+    inline int operator<=(const int32_t other) const { return (value <= int2sll(other)); }

+

+    inline int operator>=(const Fix32 &other)  const { return (value >= other.value); }

+    inline int operator>=(const sll other) const { return (value >= other); }

+    inline int operator>=(const double other)  const { return (value >= dbl2sll(other)); }

+    inline int operator>=(const float other)   const { return (value >= (float)dbl2sll(other)); }

+    inline int operator>=(const int32_t other) const { return (value >= int2sll(other)); }

+

+    inline int operator< (const Fix32 &other)  const { return (value < other.value); }

+    inline int operator< (const sll other) const { return (value < other); }

+    inline int operator< (const double other)  const { return (value < dbl2sll(other)); }

+    inline int operator< (const float other)   const { return (value < (float)dbl2sll(other)); }

+    inline int operator< (const int32_t other) const { return (value < int2sll(other)); }

+

+    inline int operator> (const Fix32 &other)  const { return (value > other.value); }

+    inline int operator> (const sll other) const { return (value > other); }

+    inline int operator> (const double other)  const { return (value > dbl2sll(other)); }

+    inline int operator> (const float other)   const { return (value > (float)dbl2sll(other)); }

+    inline int operator> (const int32_t other) const { return (value > int2sll(other)); }

+

+    inline Fix32  sin() const { return Fix32(sllsin(value)); }

+    inline Fix32  cos() const { return Fix32(sllcos(value)); }

+    inline Fix32  tan() const { return Fix32(slltan(value)); }

+    inline Fix32 asin() const { return Fix32(sllasin(value)); }

+    inline Fix32 acos() const { return Fix32(sllacos(value)); }

+    inline Fix32 atan() const { return Fix32(sllatan(value)); }

+    //inline Fix32 atan2(const Fix32 &inY) const { return Fix32(fix16_atan2(value, inY.value)); }

+    inline Fix32 sqrt() const { return Fix32(sllsqrt(value)); }

+

+};
\ No newline at end of file
diff --git a/Client/ThirdParty/math-sll/LICENSE b/Client/ThirdParty/math-sll/LICENSE
new file mode 100644
index 0000000..4e683cf
--- /dev/null
+++ b/Client/ThirdParty/math-sll/LICENSE
@@ -0,0 +1,23 @@
+
+Licensed under the terms of the MIT license:
+
+Copyright (c) 2000,2002,2006,2012,2016 Andrew E. Mileski <andrewm@isoar.ca>
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to
+deal in the Software without restriction, including without limitation the
+rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
+sell copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The copyright notice, and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+DEALINGS IN THE SOFTWARE.
+
diff --git a/Client/ThirdParty/math-sll/Makefile b/Client/ThirdParty/math-sll/Makefile
new file mode 100644
index 0000000..b7b512f
--- /dev/null
+++ b/Client/ThirdParty/math-sll/Makefile
@@ -0,0 +1,69 @@
+#
+# Licensed under the terms of the MIT license:
+#
+# Copyright (c) 2000,2002,2006,2012,2016 Andrew E. Mileski <andrewm@isoar.ca>
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to
+# deal in the Software without restriction, including without limitation the
+# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
+# sell copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The copyright notice, and this permission notice shall be included in all
+# copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+# DEALINGS IN THE SOFTWARE.
+#
+
+#
+# Installation directories
+#
+
+PREFIX	:= /usr/local
+INCDIR	:= $(PREFIX)/include
+LIBDIR	:= $(PREFIX)/lib
+
+#
+# Executables
+#
+
+AR	:= ar
+CC	:= gcc
+CFLAGS	:= -O2 -W -Wall
+INSTALL := install
+RANLIB	:= ranlib
+RM	:= rm -f
+STRIP	:= strip --strip-unneeded
+
+#
+# Recipes
+#
+
+OBJS	:= math-sll.o
+LIBS	:= math-sll.a
+
+.PHONY: all clean install
+
+all: $(LIBS)
+
+clean:
+	$(RM) $(LIBS) $(OBJS)
+
+install: $(LIBS) math-sll.h
+	$(INSTALL) -m a=rx,u+w math-sll.a $(LIBDIR)
+	$(INSTALL) -m a=r,u+w math-sll.h $(INCDIR)
+
+math-sll.o: math-sll.c math-sll.h
+
+math-sll.a: math-sll.o
+	$(STRIP) $<
+	$(AR) rcs $@ $<
+	$(RANLIB) $@
+
diff --git a/Client/ThirdParty/math-sll/README b/Client/ThirdParty/math-sll/README
new file mode 100644
index 0000000..e53ca36
--- /dev/null
+++ b/Client/ThirdParty/math-sll/README
@@ -0,0 +1,51 @@
+
+math-sll
+
+	A fixed point (32.32 bit) math library.
+
+	See math-sll.h for details.
+
+License
+
+	Licensed under the terms of the MIT license as of revision v1.19
+
+	See the LICENSE file for details.
+
+	Earlier revisions were licensed under the terms of GPL or LGPL.
+
+Installation
+
+	The source can be added as-is into your project, or a library can optionally
+	be built and installed, then statically-linked into your project.
+
+	To build:
+
+		make clean
+		make all
+
+	Optionally:
+
+		make install
+
+	The default installation path prefix is /usr/local
+
+	See the Makefile for details.
+
+Repository
+
+	A math-sll project is being maintained at:
+
+	http://github.com/aemileski/math-sll
+
+	You can also obtain a copy of the latest source code directly via GIT:
+
+	git clone git://github.com/aemileski/math-sll.git
+
+Contact
+
+	Feel free to contact me with questions, bug reports, patches, feature
+	requests, etc., or just to tell me about how you are using math-sll in
+	your project!
+
+	Andrew E. Mileski <andrewm@isoar.ca>
+
diff --git a/Client/ThirdParty/math-sll/math-sll.c b/Client/ThirdParty/math-sll/math-sll.c
new file mode 100644
index 0000000..ecc4238
--- /dev/null
+++ b/Client/ThirdParty/math-sll/math-sll.c
@@ -0,0 +1,957 @@
+/*
+ * Revision v1.24
+ *
+ * Credits
+ *
+ *	Maintained, conceived, written, and fiddled with by:
+ *
+ *		Andrew E. Mileski <andrewm@isoar.ca>
+ *	
+ *	Other source code contributors:
+ *
+ *		Kevin Rockel
+ *		Kevin Michael Woley
+ *		Mark Anthony Lisher
+ *		Nicolas Pitre
+ *		Anonymous
+ *
+ * License
+ *
+ *	Licensed under the terms of the MIT license:
+ *
+ * Copyright (c) 2000,2002,2006,2012,2016 Andrew E. Mileski <andrewm@isoar.ca>
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy
+ * of this software and associated documentation files (the "Software"), to
+ * deal in the Software without restriction, including without limitation the
+ * rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
+ * sell copies of the Software, and to permit persons to whom the Software is
+ * furnished to do so, subject to the following conditions:
+ *
+ * The copyright notice, and this permission notice shall be included in all
+ * copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+ * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ */
+
+/* See header for full details */
+#include "math-sll.h"
+
+/*
+ * Local prototypes
+ */
+
+static sll _sllcos(sll x);
+static sll _sllsin(sll x);
+
+static sll _sllexp(sll x);
+
+/*
+ * Unpack IEEE 754 floating point double format into fixed point sll format
+ *
+ * Description
+ *
+ *	IEEE 754 specifies the binary64 type ("double" in C) as having:
+ *
+ *	1 bit sign
+ *	11 bit exponent
+ *	53 bit significand
+ *
+ *	The first bit of the significand is an implied 1 which is not stored.
+ *	The decimal would be to the right of that implied 1, or to the left of
+ *	the stored significand.
+ *
+ *	The exponent is unsigned, and biased with an offset of 1023.
+ *
+ *	The IEEE 754 standard does not specify endianess, but the endian used is
+ *	traditionally the same endian that the processor uses.
+ */
+
+sll dbl2sll(double dbl)
+{
+	union {
+		double d;
+		unsigned u[2];
+		ull ull;
+		sll sll;
+	} in, retval;
+	register unsigned exp;
+
+	/* Move into memory as args might be passed in regs */
+	in.d = dbl;
+
+#if defined(BROKEN_IEEE754_DOUBLE)
+
+	exp = in.u[0];
+	in.u[0] = in.u[1];
+	in.u[1] = exp;
+
+#endif /* defined(BROKEN_IEEE754_DOUBLE) */
+
+	/* Leading 1 is assumed by IEEE */
+	retval.u[1] = 0x40000000;
+
+	/* Unpack the mantissa into the unsigned long */
+	retval.u[1] |= (in.u[1] << 10) & 0x3ffffc00;
+	retval.u[1] |= (in.u[0] >> 22) & 0x000003ff;
+	retval.u[0] = in.u[0] << 10;
+
+	/* Extract the exponent and align the decimals */
+	exp = (in.u[1] >> 20) & 0x7ff;
+	if (exp)
+		/* IEEE 754 decimal begins at right of bit position 30 */
+		retval.ull >>= (1023 + 30) - exp;
+	else
+		return 0L;
+
+	/* Negate if negative flag set */
+	if (in.u[1] & 0x80000000)
+		retval.sll = _sllneg(retval.sll);
+
+	return retval.sll;
+}
+
+/*
+ * Pack fixed point sll format into IEEE 754 floating point double format
+ *
+ * Description
+ *
+ *	IEEE 754 specifies the binary64 type ("double" in C) as having:
+ *
+ *	1 bit sign
+ *	11 bit exponent
+ *	53 bit significand
+ *
+ *	The first bit of the significand is an implied 1 which is not stored.
+ *	The decimal would be to the right of that implied 1, or to the left of
+ *	the stored significand.
+ *
+ *	The exponent is unsigned, and biased with an offset of 1023.
+ *
+ *	The IEEE 754 standard does not specify endianess, but the endian used is
+ *	traditionally the same endian that the processor uses.
+ */
+
+double sll2dbl(sll s)
+{
+	union {
+		double d;
+		unsigned u[2];
+		ull ull;
+		sll sll;
+	} in, retval;
+	register unsigned exp;
+	register unsigned flag;
+
+	if (s == 0)
+		return 0.0;
+
+	/* Move into memory as args might be passed in regs */
+	in.sll = s;
+
+	/* Handle the negative flag */
+	if (in.sll < 1) {
+		flag = 0x80000000;
+		in.ull = _sllneg(in.sll);
+	} else
+		flag = 0x00000000;
+
+	/*
+	 * Normalize
+	 *
+	 * IEEE 754 decimal-point begins at right of bit position 30
+	 */
+	for (exp = (1023 + 30); in.ull && (in.u[1] & 0x80000000) == 0; exp--) {
+		in.ull <<= 1;
+	}
+	in.ull <<= 1;
+	exp++;
+	in.ull >>= 12;
+	retval.ull = in.ull;
+	retval.u[1] |= flag | (exp << 20);
+
+#if defined(BROKEN_IEEE754_DOUBLE)
+
+	exp = retval.u[0];
+	retval.u[0] = retval.u[1];
+	retval.u[1] = exp;
+
+#endif /* defined(BROKEN_IEEE754_DOUBLE)  */
+
+	return retval.d;
+}
+
+/*
+ * Multiply two sll values
+ *
+ * Description
+ *
+ *	When multiplying two 64 bit sll numbers, the result is 128 bits, but there
+ *	is only room for a 64 bit result with sll!
+ *
+ *	The 128 bit result has 64 bits on either side of the decimal, so 32 bits
+ *	of overflow to the left of the decimal, and 32 bits of underflow to the
+ *	right of the decmial.
+ *
+ *	32.32 * 32.32 = 64.64 = overflow(32) + 32.32 + underflow(32)
+ *
+ *	However, a "long long" multiply has 64 bits of overflow to the left of the
+ *	decimal, resulting in the entire integer part being lost!
+ *
+ *	32.32 * 32.32 = 64.64 = overflow(64) + .64
+ *
+ *	Hence a custom multiply routine is required, to preserve the parts
+ *	of the result that sll needs.
+ *
+ *	Consider two sll numbers, x and y:
+ *
+ *	Let x = x_hi * 2^0 + x_lo * 2^(-32)
+ *	Let y = y_hi * 2^0 + y_lo * 2^(-32)
+ *
+ * 	Where:
+ *
+ *	*_hi is the signed 32 bit integer part to the left of the decimal
+ *	*_lo is the unsigned 32 bit fractional part to the right of the decimal
+ *
+ *	x * y = (x_hi * 2^0 + x_lo * 2^(-32))
+ *	      * (y_hi * 2^0 + y_lo * 2^(-32))
+ *
+ *	Expanding the terms, we get:
+ *
+ *	= x_hi * y_hi * 2^0 + x_hi * y_lo * 2^(-32)
+ *	+ x_lo * y_hi * 2^(-32) + x_lo * y_lo * 2^(-64)
+ *
+ *	Grouping by powers of 2, we get:
+ *
+ *	(x_hi * y_hi) * 2^0
+ *	We only need the low 32 bits of this term, as the rest is overflow
+ *
+ *	(x_hi * y_lo + x_lo * y_hi) * 2^-32
+ *	We need all bits of this term
+ *
+ *	x_lo * y_lo * 2^-64
+ *	We only need the high 32 bits of this term, as the rest is underflow
+ */
+
+sll sllmul(sll x, sll y)
+{
+	register unsigned int x_lo;
+	register signed int x_hi;
+
+	register unsigned int y_lo;
+	register signed int y_hi;
+
+	x_hi = (signed int) ((ull) x >> 32);	// Discard lower 32 bits
+	x_lo = (unsigned int) x;		// Discard upper 32 bits
+
+	y_hi = (signed int) ((ull) y >> 32);	// Discard lower 32 bits
+	y_lo = (unsigned int) y;		// Discard upper 32 bits
+
+	return (sll) (
+		  ((ull) (x_hi * y_hi) << 32)
+		+ ((ull) x_hi * y_lo + x_lo * (ull) y_hi)
+		+ (((ull) x_lo * y_lo) >> 32)
+	);
+}
+
+/*
+ * Calculate cos x where -pi/4 <= x <= pi/4
+ *
+ * Description
+ *
+ *	cos x = 1 - x^2 / 2! + x^4 / 4! - ... + x^(2N) / (2N)!
+ *	Note that (pi/4)^12 / 12! < 2^-32 which is the smallest possible number.
+ *
+ *	cos x = t0 + t1 + t2 + t3 + t4 + t5 + t6
+ *
+ *	Consider only the factorials:
+ *	f0 =  0! =  1
+ *	f1 =  2! =  2 *  1 * f0 =   2 * f0
+ *	f2 =  4! =  4 *  3 * f1 =  12 * f1
+ *	f3 =  6! =  6 *  5 * f2 =  30 * f2
+ *	f4 =  8! =  8 *  7 * f3 =  56 * f3
+ *	f5 = 10! = 10 *  9 * f4 =  90 * f4
+ *	f6 = 12! = 12 * 11 * f5 = 132 * f5
+ *
+ *	Now consider each term of the series:
+ *	t0 = 1
+ *	t1 = -t0 * x^2 / f1 = -t0 * x^2 * CONST_1_2
+ *	t2 = -t1 * x^2 / f2 = -t1 * x^2 * CONST_1_12
+ *	t3 = -t2 * x^2 / f3 = -t2 * x^2 * CONST_1_30
+ *	t4 = -t3 * x^2 / f4 = -t3 * x^2 * CONST_1_56
+ *	t5 = -t4 * x^2 / f5 = -t4 * x^2 * CONST_1_90
+ *	t6 = -t5 * x^2 / f6 = -t5 * x^2 * CONST_1_132
+ */
+
+sll _sllcos(sll x)
+{
+	sll retval;
+	sll x2;
+
+	x2 = sllmul(x, x);
+
+	retval = _sllsub(CONST_1, sllmul(x2, CONST_1_132));
+	retval = _sllsub(CONST_1, sllmul(sllmul(x2, retval), CONST_1_90));
+	retval = _sllsub(CONST_1, sllmul(sllmul(x2, retval), CONST_1_56));
+	retval = _sllsub(CONST_1, sllmul(sllmul(x2, retval), CONST_1_30));
+	retval = _sllsub(CONST_1, sllmul(sllmul(x2, retval), CONST_1_12));
+	retval = _sllsub(CONST_1, slldiv2(sllmul(x2, retval)));
+
+	return retval;
+}
+
+/*
+ * Calculate sin x where -pi/4 <= x <= pi/4
+ *
+ * Description
+ *
+ *	sin x = x - x^3 / 3! + x^5 / 5! - ... + x^(2N+1) / (2N+1)!
+ *	Note that (pi/4)^13 / 13! < 2^-32 which is the smallest possible number.
+ *
+ *	sin x = t0 + t1 + t2 + t3 + t4 + t5 + t6
+ *
+ *	Consider only the factorials:
+ *	f0 =  0! =  1
+ *	f1 =  3! =  3 *  2 * f0 =   6 * f0
+ *	f2 =  5! =  5 *  4 * f1 =  20 * f1
+ *	f3 =  7! =  7 *  6 * f2 =  42 * f2
+ *	f4 =  9! =  9 *  8 * f3 =  72 * f3
+ *	f5 = 11! = 11 * 10 * f4 = 110 * f4
+ *	f6 = 13! = 13 * 12 * f5 = 156 * f5
+ *
+ *	Now consider each term of the series:
+ *	t0 = 1
+ *	t1 = -t0 * x^2 /   6 = -t0 * x^2 * CONST_1_6
+ *	t2 = -t1 * x^2 /  20 = -t1 * x^2 * CONST_1_20
+ *	t3 = -t2 * x^2 /  42 = -t2 * x^2 * CONST_1_42
+ *	t4 = -t3 * x^2 /  72 = -t3 * x^2 * CONST_1_72
+ *	t5 = -t4 * x^2 / 110 = -t4 * x^2 * CONST_1_110
+ *	t6 = -t5 * x^2 / 156 = -t5 * x^2 * CONST_1_156
+ */
+
+sll _sllsin(sll x)
+{
+	sll retval;
+	sll x2;
+
+	x2 = sllmul(x, x);
+
+	retval = _sllsub(x, sllmul(x2, CONST_1_156));
+	retval = _sllsub(x, sllmul(sllmul(x2, retval), CONST_1_110));
+	retval = _sllsub(x, sllmul(sllmul(x2, retval), CONST_1_72));
+	retval = _sllsub(x, sllmul(sllmul(x2, retval), CONST_1_42));
+	retval = _sllsub(x, sllmul(sllmul(x2, retval), CONST_1_20));
+	retval = _sllsub(x, sllmul(sllmul(x2, retval), CONST_1_6));
+
+	return retval;
+}
+
+/*
+ * Calculate cos x for any value of x, by quadrant
+ */
+
+sll sllcos(sll x)
+{
+	int i;
+	sll retval;
+
+	/* Calculate cos (x - i * pi/2), where -pi/4 <= x - i * pi/2 <= pi/4  */
+	i = _sll2int(_slladd(sllmul(x, CONST_2_PI), CONST_1_2));
+	x = _sllsub(x, sllmul(_int2sll(i), CONST_PI_2));
+
+	/* Locate the quadrant */
+	switch (i & 3) {
+		default:
+		case 0:
+			retval = _sllcos(x);
+			break;
+		case 1:
+			retval = sllneg(_sllsin(x));
+			break;
+		case 2:
+			retval = sllneg(_sllcos(x));
+			break;
+		case 3:
+			retval = _sllsin(x);
+			break;
+	}
+
+	return retval;
+}
+
+/*
+ * Calculate sin x for any value of x, by quadrant
+ */
+
+sll sllsin(sll x)
+{
+	int i;
+	sll retval;
+
+	/* Calculate sin (x - n * pi/2), where -pi/4 <= x - i * pi/2 <= pi/4 */
+	i = _sll2int(_slladd(sllmul(x, CONST_2_PI), CONST_1_2));
+	x = _sllsub(x, sllmul(_int2sll(i), CONST_PI_2));
+
+	/* Locate the quadrant */
+	switch (i & 3) {
+		default:
+		case 0:
+			retval = _sllsin(x);
+			break;
+		case 1:
+			retval = _sllcos(x);
+			break;
+		case 2:
+			retval = sllneg(_sllsin(x));
+			break;
+		case 3:
+			retval = sllneg(_sllcos(x));
+			break;
+	}
+
+	return retval;
+}
+
+/*
+ * Calculate tan x for any value of x, by quadrant
+ */
+
+sll slltan(sll x)
+{
+	int i;
+	sll retval;
+
+	i = _sll2int(_slladd(sllmul(x, CONST_2_PI), CONST_1_2));
+	x = _sllsub(x, sllmul(_int2sll(i), CONST_PI_2));
+
+	/* Locate the quadrant */
+	switch (i & 3) {
+		default:
+		case 0:
+		case 2:
+			retval = slldiv(_sllsin(x), _sllcos(x));
+			break;
+		case 1:
+		case 3:
+			retval = _sllneg(slldiv(_sllcos(x), _sllsin(x)));
+			break;
+	}
+
+	return retval;
+}
+
+/*
+ *
+ * Calculate asin x, where |x| <= 1
+ *
+ * Description
+ *
+ *	asin x = SUM[n=0,) C(2 * n, n) * x^(2 * n + 1) / (4^n * (2 * n + 1)), |x| <= 1
+ *
+ *	where C(n, r) = nCr = n! / (r! * (n - r)!)
+ *
+ *	Using a two term approximation:
+ * [1]	a = x + x^3 / 6
+ *
+ *	Results in:
+ *	asin x = a + D
+ *	where D is the difference from the exact result
+ *
+ *	Letting D = asin d results in:
+ * [2]	asin x = a + asin d
+ *
+ *	Re-arranging:
+ *	asin x - a = asin d
+ *
+ *	Applying sin to both sides:
+ *	sin (asin x - a) = sin asin d
+ *	sin (asin x - a) = d
+ *	d = sin (asin x - a)
+ *
+ *	Applying the standard identity:
+ *	sin (u - v) = sin u * cos v - cos u * sin v
+ *
+ *	Results in:
+ *	d = sin asin x * cos a - cos asin x * sin a
+ *	d = x * cos a - cos asin x * sin a
+ *
+ *	Applying the standard identity:
+ *	cos asin u = (1 - u^2)^(1 / 2)
+ *
+ *	Results in:
+ * [3]	d = x * cos a - (1 - x^2)^(1 / 2) * sin a
+ *
+ *	Putting the pieces together:
+ * [1]	a = x + x^3 / 6
+ * [3]	d = x * cos a - (1 - x^2)^(1 / 2) * sin a
+ * [2]	asin x = a + asin d
+ *
+ *	The worst case is x = 1.0 which converges after 2 iterations.
+ */
+
+sll sllasin(sll x)
+{
+	int left_side;
+	sll a;
+	sll retval;
+
+	/* asin -x = -asin x */
+	if ((left_side = x < 0))
+		x = _sllneg(x);
+
+	/* Out-of-range */
+	if (x > CONST_1)
+		return 0;
+
+	/* Initial approximate */
+	a = sllmul(x, _slladd(CONST_1, sllmul(x, sllmul(x, CONST_1_6))));
+	retval = a;
+
+	/* First iteration */
+	x = _sllsub(sllmul(x, sllcos(a)), sllmul(sllsqrt(_sllsub(CONST_1, sllmul(x, x))), sllsin(a)));
+	a = sllmul(x, _slladd(CONST_1, sllmul(x, sllmul(x, CONST_1_6))));
+	retval = _slladd(retval, a);
+
+	/* Second iteration */
+	x = _sllsub(sllmul(x, sllcos(a)), sllmul(sllsqrt(_sllsub(CONST_1, sllmul(x, x))), sllsin(a)));
+	a = sllmul(x, _slladd(CONST_1, sllmul(x, sllmul(x, CONST_1_6))));
+	retval = _slladd(retval, a);
+
+	/* Negate result if necessary */
+	return (left_side ? _sllneg(retval): retval);
+}
+
+/*
+ * Calculate atan x
+ *
+ * Description
+ *
+ *	atan x = SUM[n=0,) (-1)^n * x^(2  * n + 1) / (2 * n + 1), |x| <= 1
+ *
+ *	Using a two term approximation:
+ * [1]	a = x - x^3 / 3
+ *
+ *	Results in:
+ * 	atan x = a + D
+ *	where D is the difference from the exact result
+ *
+ *	Letting D = atan d results in:
+ * [2]	atan x = a + atan d
+ *
+ *	Re-arranging:
+ *	atan x - a = atan d
+ *
+ *	Applying tan to both sides:
+ *	tan (atan x - a) = tan atan d
+ *	tan (atan x - a) = d
+ * 	d = tan (atan x - a)
+ *
+ *	Applying the standard identity:
+ *	tan (u - v) = (tan u - tan v) / (1 + tan u * tan v)
+ *
+ *	Results in:
+ *	d = tan (atan x - a) = (tan atan x - tan a) / (1 + tan atan x * tan a)
+ * 	d = tan (atan x - a) = (x - tan a) / (1 + x * tan a)
+ *
+ *	Let:
+ * [3]	t = tan a
+ *
+ *	Results in:
+ * [4]	d = (x - t) / (1 + x * t)
+ *
+ *	So putting the pieces together:
+ * [1]	a = x - x^3 / 3
+ * [3]	t = tan a
+ * [4]	d = (x - t) / (1 + x * t)
+ * [2]	atan x = a + atan d
+ * 	atan x = a + atan ((x - t) / (1 + x * t))
+ *
+ *	The worst case is x = 1.0 which converges after 2 iterations.
+ */
+
+sll sllatan(sll x)
+{
+	int side;
+	sll a;
+	sll t;
+	sll retval;
+
+
+	if (x < CONST_1) {
+
+		/* Left:  if (x < -1) then atan x = pi / 2 + atan 1 / x */
+		side = -1;
+		x = sllinv(x);
+
+	} else if (x > CONST_1) {
+
+		/* Right:  if (x > 1) then atan x = pi / 2 - atan 1 / x */
+		side = 1;
+		x = sllinv(x);
+
+	} else {
+		/* Middle:  -1 <= x <= 1 */
+		side = 0;
+	}
+
+	/* Initial approximate */
+	a = sllmul(x, _sllsub(CONST_1, sllmul(x, sllmul(x, CONST_1_3))));
+	retval = a;
+
+	/* First iteration */
+	t = _slldiv(_sllsin(a), _sllcos(a));
+	x = _slldiv(_sllsub(x, t), _slladd(CONST_1, sllmul(x, t)));
+	a = sllmul(x, _sllsub(CONST_1, sllmul(x, sllmul(x, CONST_1_3))));
+	retval = _slladd(retval, a);
+
+	/* Second iteration */
+	t = _slldiv(_sllsin(a), _sllcos(a));
+	x = _slldiv(_sllsub(x, t), _slladd(CONST_1, sllmul(x, t)));
+	a = sllmul(x, _sllsub(CONST_1, sllmul(x, sllmul(x, CONST_1_3))));
+	retval =  _slladd(retval, a);
+
+	if (side == -1) {
+
+		/* Left:  if (x < -1) then atan x = pi / 2 + atan 1 / x */
+		retval = _slladd(CONST_PI_2, retval);
+
+	} else if (side == 1) {
+
+		/* Right:  if (x > 1) then atan x = pi / 2 - atan 1 / x */
+		retval = _sllsub(CONST_PI_2, retval);
+	}
+
+	return retval;
+}
+
+/*
+ * Calculate e^x where -0.5 <= x <= 0.5
+ *
+ * Description:
+ *	e^x = x^0 / 0! + x^1 / 1! + ... + x^N / N!
+ *	Note that 0.5^11 / 11! < 2^-32 which is the smallest possible number.
+ */
+
+sll _sllexp(sll x)
+{
+	sll retval;
+
+	retval = CONST_1;
+
+	retval = _slladd(CONST_1, sllmul(retval, sllmul(x, CONST_1_11)));
+	retval = _slladd(CONST_1, sllmul(retval, sllmul(x, CONST_1_10)));
+	retval = _slladd(CONST_1, sllmul(retval, sllmul(x, CONST_1_9)));
+	retval = _slladd(CONST_1, sllmul(retval, slldiv2n(x, 3)));
+	retval = _slladd(CONST_1, sllmul(retval, sllmul(x, CONST_1_7)));
+	retval = _slladd(CONST_1, sllmul(retval, sllmul(x, CONST_1_6)));
+	retval = _slladd(CONST_1, sllmul(retval, sllmul(x, CONST_1_5)));
+	retval = _slladd(CONST_1, sllmul(retval, slldiv4(x)));
+	retval = _slladd(CONST_1, sllmul(retval, sllmul(x, CONST_1_3)));
+	retval = _slladd(CONST_1, sllmul(retval, slldiv2(x)));
+	retval = _slladd(CONST_1, sllmul(retval, x));
+
+	return retval;
+}
+
+/*
+ * Calculate e^x for any value of x
+ */
+
+sll sllexp(sll x)
+{
+	int i;
+	sll e;
+	sll retval;
+
+	e = CONST_E;
+
+	/* -0.5 <= x <= 0.5  */
+	i = _sll2int(_slladd(x, CONST_1_2));
+	retval = _sllexp(_sllsub(x, _int2sll(i)));
+
+	/* i >= 0 */
+	if (i < 0) {
+		i = -i;
+		e = CONST_1_E;
+	}
+
+	/* Scale the result */
+	for (; i; i >>= 1) {
+		if (i & 1)
+			retval = sllmul(retval, e);
+		e = sllmul(e, e);
+	}
+
+	return retval;
+}
+
+/*
+ * Calculate natural logarithm using Netwton-Raphson method
+ */
+
+sll slllog(sll x)
+{
+	sll x1;
+	sll ln;
+
+	ln = 0;
+
+	/* Scale: e^(-1/2) <= x <= e^(1/2) */
+	while (x < CONST_1_SQRTE) {
+		ln = _sllsub(ln, CONST_1);
+		x = sllmul(x, CONST_E);
+	}
+	while (x > CONST_SQRTE) {
+		ln = _slladd(ln, CONST_1);
+		x = sllmul(x, CONST_1_E);
+	}
+
+	/* First iteration */
+	x1 = sllmul(_sllsub(x, CONST_1), slldiv2(_sllsub(x, CONST_3)));
+	ln = _sllsub(ln, x1);
+	x = sllmul(x, _sllexp(x1));
+
+	/* Second iteration */
+	x1 = sllmul(_sllsub(x, CONST_1), slldiv2(_sllsub(x, CONST_3)));
+	ln = _sllsub(ln, x1);
+	x = sllmul(x, _sllexp(x1));
+
+	/* Third iteration */
+	x1 = sllmul(_sllsub(x, CONST_1), slldiv2(_sllsub(x, CONST_3)));
+	ln = _sllsub(ln, x1);
+
+	return ln;
+}
+
+/*
+ * Calculate the inverse for non-zero values
+ */
+
+sll sllinv(sll x)
+{
+	int sgn;
+	sll u;
+	ull s;
+
+	/* Use positive numbers, or the approximation won't work */
+	if (x < CONST_0) {
+		x = _sllneg(x);
+		sgn = 1;
+	} else {
+		sgn = 0;
+	}
+
+	/* Starting-point (gets shifted right to become positive) */
+	s = -1;
+
+	/* An approximation - must be larger than the actual value */
+	for (u = x; u; u = ((ull) u) >> 1)
+		s >>= 1;
+
+	/* Newton's Method */
+	u = sllmul(s, _sllsub(CONST_2, sllmul(x, s)));
+	u = sllmul(u, _sllsub(CONST_2, sllmul(x, u)));
+	u = sllmul(u, _sllsub(CONST_2, sllmul(x, u)));
+	u = sllmul(u, _sllsub(CONST_2, sllmul(x, u)));
+	u = sllmul(u, _sllsub(CONST_2, sllmul(x, u)));
+	u = sllmul(u, _sllsub(CONST_2, sllmul(x, u)));
+
+	return ((sgn) ? _sllneg(u): u);
+}
+
+/*
+ * Calculate x^y
+ *
+ * Description
+ *
+ *	The standard identity:
+ *	ln x^y = y * log x
+ *
+ *	Raising e to the power of either sides:
+ *	e^(ln x^y) = e^(y * log x)
+ *
+ *	Which simplifies to:
+ *	x^y = e^(y * ln x)
+ */
+
+sll sllpow(sll x, sll y)
+{
+	if (y == CONST_0)
+		return CONST_1;
+	if (y == CONST_1)
+		return x;
+	if (y == CONST_2)
+		return sllmul(x, x);
+	
+	return sllexp(sllmul(y, slllog(x)));
+}
+
+/*
+ * Calculate the square-root
+ *
+ * Description
+ *
+ *	Consider a parabola centered on the y-axis:
+ * 	y = a * x^2 + b
+ *
+ *	Has zeros (y = 0) located at:
+ *	a * x^2 + b = 0
+ *	a * x^2 = -b
+ *	x^2 = -b / a
+ *	x = +- (-b / a)^(1 / 2)
+ *
+ *	Letting a = 1 and b = -X results in:
+ *	y = x^2 - X
+ *	x = +- X^(1 / 2)
+ *
+ *	Which is a convenient result, since we want to find the square root of X,
+ *	and we can * use Newton's Method to find the zeros of any f(x):
+ *	xn = x - f(x) / f'(x)
+ *
+ *	Applying Newton's Method to our parabola:
+ *	f(x) = x^2 - X
+ *	xn = x - (x^2 - X) / (2 * x)
+ *	xn = x - (x - X / x) / 2
+ *
+ *	To make this converge quickly, we scale X so that:
+ *	X = 4^N * z
+ *
+ *	Taking the roots of both sides
+ *	X^(1 / 2) = (4^n * z)^(1 / 2)
+ *	X^(1 / 2) = 2^n * z^(1 / 2)
+ *
+ *	Letting N = 2^n results in:
+ *	x^(1 / 2) = N * z^(1 / 2)
+ *
+ *	We want this to converge to the positive root, so we must start at a point
+ *	0 < start <= x^(1 / 2)
+ * 	or
+ *	x^(1/2) <= start <= infinity
+ *
+ *	Since:
+ *	(1/2)^(1/2) = 0.707
+ *	2^(1/2) = 1.414
+ *
+ *	A good choice is 1 which lies in the middle, and takes 4 iterations to
+ *	converge from either extreme.
+ */
+
+sll sllsqrt(sll x)
+{
+	sll n;
+	sll xn;
+       
+	/* Quick solutions for the simple cases */
+	if (x <= CONST_0 || x == CONST_1)
+		return x;
+
+	/* Start with a scaling factor of 1 */
+	n = CONST_1;
+
+	/* Scale x so that 0.5 <= x < 2 */
+	while (x >= CONST_2) {
+		x = slldiv4(x);
+		n = sllmul2(n);
+	}
+	while (x < CONST_1_2) {
+		x = sllmul4(x);
+		n = slldiv2(n);
+	}
+
+	/* Simple solution if x = 4^n */
+	if (x == CONST_1)
+		return n;
+
+	/* The starting point */
+	xn = CONST_1;
+
+	/* Four iterations will be enough */
+	xn = _sllsub(xn, slldiv2(_sllsub(xn, slldiv(x, xn))));
+	xn = _sllsub(xn, slldiv2(_sllsub(xn, slldiv(x, xn))));
+	xn = _sllsub(xn, slldiv2(_sllsub(xn, slldiv(x, xn))));
+	xn = _sllsub(xn, slldiv2(_sllsub(xn, slldiv(x, xn))));
+
+	/* Scale the result */
+	return sllmul(n, xn);
+}
+
+sll slld2dsqrt(sll num)
+{
+	if (num <= CONST_0 || num == CONST_1)
+	{
+		return num;
+	}
+	sll result = 0;
+	sll bit;
+	
+	// Many numbers will be less than 15, so
+	// this gives a good balance between time spent
+	// in if vs. time spent in the while loop
+	// when searching for the starting value.
+	if (num & 0xFFFFFF0000000000LL)
+		bit = 0x4000000000000000LL; //(sll)1 << 62;
+	else
+		bit = 0x0000004000000000LL; //(sll)1 << 38;
+	
+	while (bit > num) bit >>= 2;
+	
+	while (bit)
+	{
+		if (num >= result + bit)
+		{
+			num -= result + bit;
+			result = (result >> 1) + bit;
+		}
+		else
+		{
+			result = (result >> 1);
+		}
+		bit >>= 2;
+	}
+		
+	// Then process it again to get the lowest 8 bits.
+	// 尾数定点数大于等于1
+	if (num > CONST_P9999)
+	{
+		// The remainder 'num' is too large to be shifted left
+		// by 32, so we have to add 0.5 to result manually and
+		// adjust 'num' accordingly. a是原本的数
+		// why 0.5 is enough? because:  result^2 <= a < (result+1)^2 
+		// num = a - (result + 0.5)^2   
+		//	 = num + result^2 - (result + 0.5)^2
+		//	 = num - result - 0.25
+		num -= result;
+		// 注意，严格来说，这里要左移16位才缩小到结果，但是为了后续小数不用再放大，这里提前做了放大
+		num = (num << 32) - CONST_1_4;
+		result = (result << 32) + CONST_1_2;
+	}
+	else
+	{
+		num <<= 32;
+		result <<= 32;
+	}
+	
+	bit = 0x0000000040000000LL; //(sll)1 << 30;
+
+	while (bit)
+	{
+		if (num >= result + bit)
+		{
+			num -= result + bit;
+			result = (result >> 1) + bit;
+		}
+		else
+		{
+			result = (result >> 1);
+		}
+		bit >>= 2;
+	}
+
+	return result;
+}
\ No newline at end of file
diff --git a/Client/ThirdParty/math-sll/math-sll.h b/Client/ThirdParty/math-sll/math-sll.h
new file mode 100644
index 0000000..b10dbe7
--- /dev/null
+++ b/Client/ThirdParty/math-sll/math-sll.h
@@ -0,0 +1,770 @@
+#if !defined(MATH_SLL_H)
+#  define MATH_SLL_H
+
+/*
+ * Revision v1.24
+ *
+ *	A fixed point (32.32 bit) math library.
+ *
+ * Description
+ *
+ *	Floating point packs the most accuracy in the available bits, but it
+ *	often provides more accuracy than is required.  It is time consuming to
+ *	carry the extra precision around, particularly on platforms that don't
+ *	have a dedicated floating point processor.
+ *
+ *	This library is a compromise.  All math is done using the 64 bit signed
+ *	"long long" format (sll), and is not intended to be portable, just as
+ *	simple and as fast as possible.
+ *
+ *	As some processors lack division instructions but have multiplication
+ *	instructions, multiplication is favored over division.  This can be a
+ *	penalty when used on a processor with a division instruction, so it is
+ *	recommended to modify the division functions and macros in that case.
+ *
+ *	On procesors without multiplication instructions, other algorithms, for
+ *	example CORDIC, are probably faster.
+ *
+ *	Since "long long" is a elementary type, it can be passed around without
+ *	resorting to the use of pointers.  Since the format used is fixed point,
+ *	there is never a need to do time consuming checks and adjustments to
+ *	maintain normalized numbers, as is the case in floating point.
+ *
+ *	Simply put, this library is limited to handling numbers with a whole
+ *	part of up to 2^31 - 1 = 2.147483647e9 in magnitude, and fractional
+ *	parts down to 2^-32 = 2.3283064365e-10 in magnitude.  This yields a
+ *	decent range and accuracy for many applications.
+ *
+ * IMPORTANT
+ *
+ *	No checking for arguments out of range (error).
+ *	No checking for divide by zero (error).
+ *	No checking for overflow (error).
+ *	No checking for underflow (warning).
+ *	Chops, doesn't round.
+ *
+ * Functions
+ *
+ *	sll dbl2sll(double d)			double to sll
+ *	double sll2dbl(sll s)			sll to double
+ *
+ *	sll int2sll(int i)			integer to sll
+ *	int sll2int(sll s)			sll to integer
+ *
+ *	sll sllint(sll s)			Set fractional-part to 0
+ *	sll sllfrac(sll s)			Set integer-part to 0
+ *
+ *	sll slladd(sll x, sll y)		x + y
+ *	sll sllneg(sll x)			-x
+ *	sll sllsub(sll x, sll y)		x - y
+ *
+ *	sll sllmul(sll x, sll y)		x * y
+ *	sll sllmul2(sll x)			x * 2
+ *	sll sllmul2n(sll x, int n)		x * 2^n, 0 <= n <= 31
+ *	sll sllmul4(sll x)			x * 4
+ *
+ *	sll slldiv(sll x, sll y)		x / y
+ *	sll slldiv2(sll x)			x / 2
+ *	sll slldiv2n(sll x, int n)		x / 2^n, 0 <= n <= 31
+ *	sll slldiv4(sll x)			x / 4
+ *
+ *	sll sllcos(sll x)			cos x
+ *	sll sllsin(sll x)			sin x
+ *	sll slltan(sll x)			tan x
+ *
+ *	sll sllsec(sll x)			sec x = 1 / cos x
+ *	sll sllcsc(sll x)			csc x = 1 / sin x
+ *	sll sllcot(sll x)			cot x = 1 / tan x = cos x / sin x
+ *
+ *	sll sllacos(sll x)			acos x
+ *	sll sllasin(sll x)			asin x
+ *	sll sllatan(sll x)			atan x
+ *
+ *	sll sllcosh(sll x)			cosh x
+ *	sll sllsinh(sll x)			sinh x
+ *	sll slltanh(sll x)			tanh x
+ *
+ *	sll sllsech(sll x)			sech x
+ *	sll sllcsch(sll x)			cosh x
+ *	sll sllcoth(sll x)			coth x
+ *
+ *	sll sllexp(sll x)			e^x
+ *	sll slllog(sll x)			ln x
+ *
+ *	sll sllinv(sll v)			1 / x
+ *	sll sllpow(sll x, sll y)		x^y
+ *	sll sllsqrt(sll x)			x^(1 / 2)
+ *
+ *	sll sllfloor(sll x)			floor x
+ *	sll sllceil(sll x)			ceiling x
+ *
+ * Macros
+ *
+ *	Use of the following macros is optional, but may be beneficial with
+ *	some compilers.  Using the non-macro versions is strongly recommended.
+ *
+ *	WARNING:  macros do not type-check their arguments!
+ *
+ *	_int2sll(X)				See function int2sll()
+ *	_sll2int(X)				See function sll2int()
+ *
+ *	_sllint(X)				See function sllint()
+ *	_sllfrac(X)				See function sllfrac()
+ *
+ *	_slladd(X,Y)				See function slladd()
+ *	_sllneg(X)				See function sllneg()
+ *	_sllsub(X,Y)				See function sllsub()
+ *
+ *	_sllmul2(X)				See function sllmul2()
+ *	_sllmul2n(X)				See function sllmul2n()
+ *	_sllmul4(X)				See function sllmul4()
+ *
+ *	_slldiv(X,Y)				See function slldiv()
+ *	_slldiv2(X,Y)				See function slldiv2()
+ *	_slldiv2n(X,Y)				See function slldiv2n()
+ *	_slldiv4(X,Y)				See function slldiv4()
+ *
+ * Credits
+ *
+ *	Maintained, conceived, written, and fiddled with by:
+ *
+ *		Andrew E. Mileski <andrewm@isoar.ca>
+ *
+ *	Other source code contributors:
+ *
+ *		Kevin Rockel
+ *		Kevin Michael Woley
+ *		Mark Anthony Lisher
+ *		Nicolas Pitre
+ *		Anonymous
+ *
+ * License
+ *
+ *	Licensed under the terms of the MIT license:
+ *
+ * Copyright (c) 2000,2002,2006,2012,2016 Andrew E. Mileski <andrewm@isoar.ca>
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy
+ * of this software and associated documentation files (the "Software"), to
+ * deal in the Software without restriction, including without limitation the
+ * rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
+ * sell copies of the Software, and to permit persons to whom the Software is
+ * furnished to do so, subject to the following conditions:
+ *
+ * The copyright notice, and this permission notice shall be included in all
+ * copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+ * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ */
+
+/* DEC SA-110 "StrongARM" (armv4l) architecture has a big-endian double */
+#if defined(__arm__)
+#  if (!defined(__BYTE_ORDER__) || (__BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__))
+#	define BROKEN_IEEE754_DOUBLE
+#  endif
+#endif
+
+#ifdef _MSC_VER
+# define __inline__ __inline
+# ifndef _MSC_STDINT_H_
+typedef signed __int64       int64_t;
+typedef unsigned __int64     uint64_t;
+#endif
+#define __extension__ 
+#else 
+#include <stdint.h>
+#if defined(__GNUC__)
+	
+#else
+	#define __inline__ inline
+#endif
+#endif
+/*
+ * Data types
+ */
+
+__extension__ typedef int64_t sll;
+__extension__ typedef uint64_t ull;
+
+/*
+ * Function prototypes
+ */
+
+sll dbl2sll(double d);
+double sll2dbl(sll s);
+
+static __inline__ sll int2sll(int i);
+static __inline__ int sll2int(sll s);
+
+static __inline__ sll sllint(sll s);
+static __inline__ sll sllfrac(sll s);
+
+static __inline__ sll slladd(sll x, sll y);
+static __inline__ sll sllneg(sll s);
+static __inline__ sll sllsub(sll x, sll y);
+sll sllmul(sll x, sll y);
+static __inline__ sll sllmul2(sll x);
+static __inline__ sll sllmul4(sll x);
+static __inline__ sll sllmul2n(sll x, int n);
+
+static __inline__ sll slldiv(sll x, sll y);
+static __inline__ sll slldiv2(sll x);
+static __inline__ sll slldiv4(sll x);
+static __inline__ sll slldiv2n(sll x, int n);
+
+sll sllcos(sll x);
+sll sllsin(sll x);
+sll slltan(sll x);
+
+static __inline__ sll sllacos(sll x);
+sll sllasin(sll x);
+sll sllatan(sll x);
+
+static __inline__ sll sllsec(sll x);
+static __inline__ sll sllcsc(sll x);
+static __inline__ sll sllcot(sll x);
+
+static __inline__ sll sllcosh(sll x);
+static __inline__ sll sllsinh(sll x);
+static __inline__ sll slltanh(sll x);
+
+static __inline__ sll sllsech(sll x);
+static __inline__ sll sllcsch(sll x);
+static __inline__ sll sllcoth(sll x);
+
+sll sllexp(sll x);
+sll slllog(sll x);
+
+sll sllpow(sll x, sll y);
+sll sllinv(sll v);
+sll sllsqrt(sll x);
+sll slld2dsqrt(sll x);
+
+static __inline__ sll sllfloor(sll x);
+static __inline__ sll sllceil(sll x);
+
+/*
+ * Macros
+ *
+ * WARNING - Macros don't type-check!
+ */
+
+#define _int2sll(X)	(((sll) (X)) << 32)
+#define _sll2int(X)	((int) ((X) >> 32))
+
+#define _sllint(X)	((X) & 0xffffffff00000000LL)
+#define _sllfrac(X)	((X) & 0x00000000ffffffffLL)
+
+#define _slladd(X,Y)	((X) + (Y))
+#define _sllneg(X)	(-(X))
+#define _sllsub(X,Y)	((X) - (Y))
+
+#define _sllmul2(X)	((X) << 1)
+#define _sllmul4(X)	((X) << 2)
+#define _sllmul2n(X,N)	((X) << (N))
+
+#define _slldiv(X,Y)	sllmul((X), sllinv(Y))
+#define _slldiv2(X)	((X) >> 1)
+#define _slldiv4(X)	((X) >> 2)
+#define _slldiv2n(X,N)	((X) >> (N))
+
+/*
+ * Constants (converted from double)
+*/
+
+#define CONST_0		0x0000000000000000LL	// 0.0
+#define CONST_1		0x0000000100000000LL	// 1.0
+#define CONST_neg1	0xffffffff00000000LL	// -1.0
+#define CONST_2		0x0000000200000000LL 	// 2.0
+#define CONST_3		0x0000000300000000LL	// 3.0
+#define CONST_4		0x0000000400000000LL	// 4.0
+#define CONST_10	0x0000000a00000000LL	// 10.0
+#define CONST_10	0x0000000a00000000LL	// 10.0
+#define CONST_1_2	0x0000000080000000LL	// 1.0 / 2.0
+#define CONST_1_3	0x0000000055555555LL	// 1.0 / 3.0
+#define CONST_1_4	0x0000000040000000LL	// 1.0 / 4.0
+#define CONST_1_5	0x0000000033333333LL	// 1.0 / 5.0
+#define CONST_1_6	0x000000002aaaaaaaLL	// 1.0 / 6.0
+#define CONST_1_7	0x0000000024924924LL	// 1.0 / 7.0
+#define CONST_1_8	0x0000000020000000LL	// 1.0 / 8.0
+#define CONST_1_9	0x000000001c71c71cLL	// 1.0 / 9.0
+#define CONST_1_10	0x0000000019999999LL	// 1.0 / 10.0
+#define CONST_1_11	0x000000001745d174LL	// 1.0 / 11.0
+#define CONST_1_12	0x0000000015555555LL	// 1.0 / 12.0
+#define CONST_1_20	0x000000000cccccccLL	// 1.0 / 20.0
+#define CONST_1_30	0x0000000008888888LL	// 1.0 / 30.0
+#define CONST_1_42	0x0000000006186186LL	// 1.0 / 42.0
+#define CONST_1_56	0x0000000004924924LL	// 1.0 / 56.0
+#define CONST_1_72	0x00000000038e38e3LL	// 1.0 / 72.0
+#define CONST_1_90	0x0000000002d82d82LL	// 1.0 / 90.0
+#define CONST_1_110	0x000000000253c825LL	// 1.0 / 110.0
+#define CONST_1_132	0x0000000001f07c1fLL	// 1.0 / 132.0
+#define CONST_1_156	0x0000000001a41a41LL	// 1.0 / 156.0
+#define CONST_P9999 0x00000000FFFFFFFFLL    // 0.999999999999
+
+#define CONST_E		0x00000002b7e15162LL	// E
+#define CONST_1_E	0x000000005e2d58d8LL	// 1 / E
+#define CONST_SQRTE	0x00000001a61298e1LL	// sqrt(E)
+#define CONST_1_SQRTE	0x000000009b4597e3LL	// 1 / sqrt(E)
+#define CONST_LOG2_E	0x0000000171547652LL	// ln(E)
+#define CONST_LOG10_E	0x000000006f2dec54LL	// log(E)
+#define CONST_LN2	0x00000000b17217f7LL	// ln(2)
+#define CONST_LN10	0x000000024d763776LL	// ln(10)
+
+#define CONST_PI	0x00000003243f6a88LL	// PI
+#define CONST_2PI	0x00000006487ED510LL	// PI
+#define CONST_PI_2	0x00000001921fb544LL	// PI / 2
+#define CONST_PI_4	0x00000000c90fdaa2LL	// PI / 4
+#define CONST_1_PI	0x00000000517cc1b7LL	// 1 / PI
+#define CONST_2_PI	0x00000000a2f9836eLL	// 2 / PI
+#define CONST_180_PI  	0x000000394BB834C7LL	// 180 / PI 246083499207.51537232162612011973
+#define CONST_PI_180	0x000000000477d1a8LL	// PI / 180 74961320.580677883103327681382757
+#define CONST_2_SQRTPI	0x0000000120dd7504LL	// 2 / sqrt(PI)
+#define CONST_SQRT2	0x000000016a09e667LL	// sqrt(2)
+#define CONST_1_SQRT2	0x00000000b504f333LL	// 1 / sqrt(2)
+
+#define CONST_FACT_0	0x0000000100000000LL	// 0!
+#define CONST_FACT_1	0x0000000100000000LL	// 1!
+#define CONST_FACT_2	0x0000000200000000LL	// 2!
+#define CONST_FACT_3	0x0000000600000000LL	// 3!
+#define CONST_FACT_4	0x0000001800000000LL	// 4!
+#define CONST_FACT_5	0x0000007800000000LL	// 5!
+#define CONST_FACT_6	0x000002d000000000LL	// 6!
+#define CONST_FACT_7	0x000013b000000000LL	// 7!
+#define CONST_FACT_8	0x00009d8000000000LL	// 8!
+#define CONST_FACT_9	0x0005898000000000LL	// 9!
+#define CONST_FACT_10	0x00375f0000000000LL	// 10!
+#define CONST_FACT_11	0x0261150000000000LL	// 11!
+#define CONST_FACT_12	0x1c8cfc0000000000LL	// 12!
+
+#define CONST_MAX		0x7FFFFFFFFFFFFFFFLL  // 最大
+#define CONST_MIN		0x8000000000000000LL	// 最小
+
+/*
+ * Convert integer to sll
+ */
+
+static __inline__ sll int2sll(int i)
+{
+	return _int2sll(i);
+}
+
+/*
+ * Convert sll to integer (truncates)
+ */
+
+static __inline__ int sll2int(sll s)
+{
+	return _sll2int(s);
+}
+
+/*
+ * Integer-part of sll (fractional-part set to 0)
+ */
+
+static __inline__ sll sllint(sll s)
+{
+	return _sllint(s);
+}
+
+/*
+ * Fractional-part of sll (integer-part set to 0)
+ */
+
+static __inline__ sll sllfrac(sll s)
+{
+	return _sllfrac(s);
+}
+
+/*
+ * Addition
+ */
+
+static __inline__ sll slladd(sll x, sll y)
+{
+	return _slladd(x, y);
+}
+
+/*
+ * Negate
+ */
+
+static __inline__ sll sllneg(sll s)
+{
+	return _sllneg(s);
+}
+
+/*
+ * Subtraction
+ */
+
+static __inline__ sll sllsub(sll x, sll y)
+{
+	return _sllsub(x, y);
+}
+
+/*
+ * Multiply two sll values
+ *
+ * Description
+ *
+ *	Let a = A * 2^32 + a_h * 2^0 + a_l * 2^(-32)
+ *	Let b = B * 2^32 + b_h * 2^0 + b_l * 2^(-32)
+ *
+ * 	Where:
+ *
+ *	*_h is the integer part
+ *	*_l the fractional part
+ *	A and B are the sign (0 for positive, -1 for negative).
+ *
+ *	a * b = (A * 2^32 + a_h * 2^0 + a_l * 2^-32)
+ *		* (B * 2^32 + b_h * 2^0 + b_l * 2^-32)
+ *
+ *	Expanding the terms, we get:
+ *
+ *	= A * B * 2^64 + A * b_h * 2^32 + A * b_l * 2^0
+ *	+ a_h * B * 2^32 + a_h * b_h * 2^0 + a_h * b_l * 2^-32
+ *	+ a_l * B * 2^0 + a_l * b_h * 2^-32 + a_l * b_l * 2^-64
+ *
+ *	Grouping by powers of 2, we get:
+ *
+ *	= A * B * 2^64
+ *	Meaningless overflow from sign extension - ignore
+ *
+ *	+ (A * b_h + a_h * B) * 2^32
+ *	Overflow which we can't handle - ignore
+ *
+ *	+ (A * b_l + a_h * b_h + a_l * B) * 2^0
+ *	We only need the low 32 bits of this term, as the rest is overflow
+ *
+ *	+ (a_h * b_l + a_l * b_h) * 2^-32
+ *	We need all 64 bits of this term
+ *
+ *	+  a_l * b_l * 2^-64
+ *	We only need the high 32 bits of this term, as the rest is underflow
+ *
+ *	Note that:
+ *	a > 0 && b > 0: A =  0, B =  0 and the third term is a_h * b_h
+ *	a < 0 && b > 0: A = -1, B =  0 and the third term is a_h * b_h - b_l
+ *	a > 0 && b < 0: A =  0, B = -1 and the third term is a_h * b_h - a_l
+ *	a < 0 && b < 0: A = -1, B = -1 and the third term is a_h * b_h - a_l - b_l
+ */
+
+/*
+ * Multiplication by 2
+ */
+
+static __inline__ sll sllmul2(sll x)
+{
+	return _sllmul2(x);
+}
+
+/*
+ * Multiplication by 4
+ */
+
+static __inline__ sll sllmul4(sll x)
+{
+	return _sllmul4(x);
+}
+
+/*
+ * Multiplication by power of 2
+ */
+
+static __inline__ sll sllmul2n(sll x, int n)
+{
+	return _sllmul2n(x, n);
+}
+
+/*
+ * Division
+ */
+
+static __inline__ sll slldiv(sll x, sll y)
+{
+	return _slldiv(x, y);
+}
+
+/*
+ * Division by 2
+ */
+
+static __inline__ sll slldiv2(sll x)
+{
+	return _slldiv2(x);
+}
+
+/*
+ * Division by 4
+ */
+
+static __inline__ sll slldiv4(sll x)
+{
+	return _slldiv4(x);
+}
+
+/*
+ * Division by power of 2
+ */
+
+static __inline__ sll slldiv2n(sll x, int n)
+{
+	return _slldiv2n(x, n);
+}
+
+/*
+ *
+ * Calculate acos x, where |x| <= 1
+ *
+ * Description
+ *
+ *	acos x = pi / 2 - asin x
+ *	acos x = pi / 2 - SUM[n=0,) C(2 * n, n) * x^(2 * n + 1) / (4^n * (2 * n + 1)), |x| <= 1
+ *
+ *	where C(n, r) = nCr = n! / (r! * (n - r)!)
+ */
+
+static __inline__ sll sllacos(sll x)
+{
+	return _sllsub(CONST_PI_2, sllasin(x));
+}
+
+/*
+ * Trigonometric secant
+ *
+ * Description
+ *
+ *	sec x = 1 / cos x
+ *
+ * An alternate algorithm, like a power series, would be more accurate.
+ */
+
+static __inline__ sll sllsec(sll x)
+{
+	return sllinv(sllcos(x));
+}
+
+/*
+ * Trigonometric cosecant
+ *
+ * Description
+ *
+ *	csc x = 1 / sin x
+ *
+ * An alternate algorithm, like a power series, would be more accurate.
+ */
+
+static __inline__ sll sllcsc(sll x)
+{
+	return sllinv(sllsin(x));
+}
+
+/*
+ * Trigonometric cotangent
+ *
+ * Description
+ *
+ *	cot x = 1 / tan x
+ *
+ *	cot x = cos x / sin x
+ *
+ * An alternate algorithm, like a power series, would be more accurate.
+ */
+
+static __inline__ sll sllcot(sll x)
+{
+	return _slldiv(sllcos(x), sllsin(x));
+}
+
+/*
+ * Hyperbolic cosine
+ *
+ * Description
+ *
+ *	cosh x = (e^x + e^(-x)) / 2
+ *
+ *	cosh x = 1 + x^2 / 2! + ... + x^(2 * N) / (2 * N)!
+ *
+ * An alternate algorithm, like a power series, would be more accurate.
+ */
+
+static __inline__ sll sllcosh(sll x)
+{
+	return _slldiv2(_slladd(sllexp(x), sllexp(_sllneg(x))));
+}
+
+/*
+ * Hyperbolic sine
+ *
+ * Description
+ *
+ *	sinh x = (e^x - e^(-x)) / 2
+ *
+ *	sinh x = 1 + x^3 / 3! + ... + x^(2 * N + 1) / (2 * N + 1)!
+ *
+ * An alternate algorithm, like a power series, would be more accurate.
+ */
+
+static __inline__ sll sllsinh(sll x)
+{
+	return _slldiv2(_sllsub(sllexp(x), sllexp(_sllneg(x))));
+}
+
+/*
+ * Hyperbolic tangent
+ *
+ * Description
+ *
+ *	tanh x = sinh x / cosh x
+ *
+ *	tanh x = (e^x - e^(-x)) / (e^x + e^(-x))
+ *
+ *	tanh x = (e^(2 * x) - 1) / (e^(2 * x) + 1)
+ *
+ * An alternate algorithm, like a power series, would be more accurate.
+ */
+
+static __inline__ sll slltanh(sll x)
+{
+	register sll e2x;
+
+	e2x = sllexp(_sllmul2(x));
+
+	return _slldiv(_sllsub(e2x, CONST_1), _slladd(e2x, CONST_1));
+}
+
+/*
+ * Hyperbolic secant
+ *
+ * Description
+ *
+ *	sech x = 1 / cosh x
+ *
+ *	sech x = 2 / (e^x + e^(-x))
+ *
+ *	sech x = 2 * e^x / (e^(2 * x) + 1)
+ *
+ * An alternate algorithm, like a power series, would be more accurate.
+ */
+
+static __inline__ sll sllsech(sll x)
+{
+	return _slldiv(_sllmul2(sllexp(x)), _slladd(sllexp(_sllmul2(x)), CONST_1));
+}
+
+/*
+ * Hyperbolic cosecant
+ *
+ * Description
+ *
+ *	csch x = = 1 / sinh x
+ *
+ *	csch x = 2 / (e^x - e^(-x))
+ *
+ *	csch x = 2 * e^x / (e^(2 * x) - 1)
+ *
+ * An alternate algorithm, like a power series, would be more accurate.
+ */
+
+static __inline__ sll sllcsch(sll x)
+{
+	return _slldiv(_sllmul2(sllexp(x)), _sllsub(sllexp(_sllmul2(x)), CONST_1));
+}
+
+/*
+ * Hyperbolic cotangent
+ *
+ * Description
+ *
+ *	coth x =  1 / tanh x
+ *
+ *	coth x = cosh x / sinh x
+ *
+ *	coth x = (e^x + e^(-x)) / (e^x - e^(-x))
+ *
+ *	coth x = (e^(2 * x) + 1) / (e^(2 * x) - 1)
+ *
+ * An alternate algorithm, like a power series, would be more accurate.
+ */
+
+static __inline__ sll sllcoth(sll x)
+{
+	register sll e2x;
+
+	e2x = sllexp(sllmul2(x));
+
+	return _slldiv(_slladd(e2x, CONST_1), _sllsub(e2x, CONST_1));
+}
+
+/*
+ * Floor
+ *
+ * Description
+ *
+ *	floor x = largest integer not larger than x
+ */
+
+static __inline__ sll sllfloor(sll x)
+{
+	register sll retval;
+
+	retval = _sllint(x);
+
+	return ((retval > x) ? _sllsub(retval, CONST_1): retval);
+}
+
+/*
+ * Ceiling
+ *
+ * Description
+ *
+ *	ceil x = smallest integer not smaller than x
+ */
+
+static __inline__ sll sllceil(sll x)
+{
+	register sll retval;
+
+	retval = _sllint(x);
+
+	return ((retval < x) ? _slladd(retval, CONST_1): retval);
+}
+
+#define sllabs(x) ((x) < 0 ? -(x) : (x))
+#define __VECTOR2_META__ "__VECTOR2_META__"
+#define __VECTOR3_META__ "__VECTOR3_META__"
+#define __ROT2_META__ "__ROT2_META__"
+#define __ROT4_META__ "__ROT4_META__"
+#define __METATABLE_NAME "__FIX_METATABLE__"
+
+static int _mul[] = {1, 10, 100, 1000, 10000, 100000, 1000000};
+
+typedef struct Vector2
+{
+	sll x;
+	sll y;
+}Vector2;
+
+typedef struct Vector3
+{
+	sll x;
+	sll y;
+	sll z;
+}Vector3;
+
+typedef struct Vector4
+{
+	sll x;
+	sll y;
+	sll z;
+	sll w;
+	
+}Vector4;
+#endif