Print done when done.

1 files changed, 402 insertions, 0 deletions

diff --git a/indra/llmath/llmath.h b/indra/llmath/llmath.h
new file mode 100644
index 0000000000..ad8ced9e1a
--- /dev/null
+++ b/indra/llmath/llmath.h
@@ -0,0 +1,402 @@
+/** 
+ * @file llmath.h
+ * @brief Useful math constants and macros.
+ *
+ * Copyright (c) 2000-$CurrentYear$, Linden Research, Inc.
+ * $License$
+ */
+
+#ifndef LLMATH_H
+#define LLMATH_H
+
+#include <cmath>
+#include <math.h>
+#include <stdlib.h>
+
+#include "lldefs.h"
+
+// work around for Windows & older gcc non-standard function names.
+#if LL_WINDOWS
+#define llisnan(val)	_isnan(val)
+#define llfinite(val)	_finite(val)
+#elif (LL_LINUX && __GNUC__ <= 2)
+#define llisnan(val)	isnan(val)
+#define llfinite(val)	isfinite(val)
+#else
+#define llisnan(val)	std::isnan(val)
+#define llfinite(val)	std::isfinite(val)
+#endif
+
+// Single Precision Floating Point Routines
+#ifndef fsqrtf
+#define fsqrtf(x)		((F32)sqrt((F64)(x)))
+#endif
+#ifndef sqrtf
+#define sqrtf(x)		((F32)sqrt((F64)(x)))
+#endif
+
+#ifndef cosf
+#define cosf(x)		((F32)cos((F64)(x)))
+#endif
+#ifndef sinf
+#define sinf(x)		((F32)sin((F64)(x)))
+#endif
+#ifndef tanf
+#define tanf(x)		((F32)tan((F64)(x)))
+#endif
+#ifndef acosf
+#define acosf(x)		((F32)acos((F64)(x)))
+#endif
+
+#ifndef powf
+#define powf(x,y) ((F32)pow((F64)(x),(F64)(y)))
+#endif
+
+const F32	GRAVITY			= -9.8f;
+
+// mathematical constants
+const F32	F_PI		= 3.1415926535897932384626433832795f;
+const F32	F_TWO_PI	= 6.283185307179586476925286766559f;
+const F32	F_PI_BY_TWO	= 1.5707963267948966192313216916398f;
+const F32	F_SQRT2		= 1.4142135623730950488016887242097f;
+const F32	F_SQRT3		= 1.73205080756888288657986402541f;
+const F32	OO_SQRT2	= 0.7071067811865475244008443621049f;
+const F32	DEG_TO_RAD	= 0.017453292519943295769236907684886f;
+const F32	RAD_TO_DEG	= 57.295779513082320876798154814105f;
+const F32	F_APPROXIMATELY_ZERO = 0.00001f;
+const F32	F_LN2		= 0.69314718056f;
+const F32	OO_LN2		= 1.4426950408889634073599246810019f;
+
+// BUG: Eliminate in favor of F_APPROXIMATELY_ZERO above?
+const F32 FP_MAG_THRESHOLD = 0.0000001f;
+
+// TODO: Replace with logic like is_approx_equal
+inline BOOL is_approx_zero( F32 f ) { return (-F_APPROXIMATELY_ZERO < f) && (f < F_APPROXIMATELY_ZERO); }
+
+inline BOOL is_approx_equal(F32 x, F32 y)
+{
+	const S32 COMPARE_MANTISSA_UP_TO_BIT = 0x02;
+	return (abs((S32) ((U32&)x - (U32&)y) ) < COMPARE_MANTISSA_UP_TO_BIT);
+}
+
+inline S32 llabs(const S32 a)
+{
+	return S32(labs(a));
+}
+
+inline F32 llabs(const F32 a)
+{
+	return F32(fabs(a));
+}
+
+inline F64 llabs(const F64 a)
+{
+	return F64(fabs(a));
+}
+
+inline S32 lltrunc( F32 f )
+{
+#if LL_WINDOWS && !defined( __INTEL_COMPILER )
+		// Avoids changing the floating point control word.
+		// Add or subtract 0.5 - epsilon and then round
+		const static U32 zpfp[] = { 0xBEFFFFFF, 0x3EFFFFFF };
+		S32 result;
+		__asm {
+			fld		f
+			mov		eax,	f
+			shr		eax,	29
+			and		eax,	4
+			fadd	dword ptr [zpfp + eax]
+			fistp	result
+		}
+		return result;
+#else
+		return (S32)f;
+#endif
+}
+
+inline S32 lltrunc( F64 f )
+{
+	return (S32)f;
+}
+
+inline S32 llfloor( F32 f )
+{
+#if LL_WINDOWS && !defined( __INTEL_COMPILER )
+		// Avoids changing the floating point control word.
+		// Accurate (unlike Stereopsis version) for all values between S32_MIN and S32_MAX and slightly faster than Stereopsis version.
+		// Add -(0.5 - epsilon) and then round
+		const U32 zpfp = 0xBEFFFFFF;
+		S32 result;
+		__asm { 
+			fld		f
+			fadd	dword ptr [zpfp]
+			fistp	result
+		}
+		return result;
+#else
+		return (S32)floor(f);
+#endif
+}
+
+
+inline S32 llceil( F32 f )
+{
+	// This could probably be optimized, but this works.
+	return (S32)ceil(f);
+}
+
+
+#ifndef BOGUS_ROUND
+// Use this round.  Does an arithmetic round (0.5 always rounds up)
+inline S32 llround(const F32 val)
+{
+	return llfloor(val + 0.5f);
+}
+
+#else // BOGUS_ROUND
+// Old llround implementation - does banker's round (toward nearest even in the case of a 0.5.
+// Not using this because we don't have a consistent implementation on both platforms, use
+// llfloor(val + 0.5f), which is consistent on all platforms.
+inline S32 llround(const F32 val)
+{
+	#if LL_WINDOWS
+		// Note: assumes that the floating point control word is set to rounding mode (the default)
+		S32 ret_val;
+		_asm fld	val
+		_asm fistp	ret_val;
+		return ret_val;
+	#elif LL_LINUX
+		// Note: assumes that the floating point control word is set
+		// to rounding mode (the default)
+		S32 ret_val;
+		__asm__ __volatile__( "flds %1    \n\t"
+							  "fistpl %0  \n\t"
+							  : "=m" (ret_val)
+							  : "m" (val) );
+		return ret_val;
+	#else
+		return llfloor(val + 0.5f);
+	#endif
+}
+
+// A fast arithmentic round on intel, from Laurent de Soras http://ldesoras.free.fr
+inline int round_int(double x)
+{
+	const float round_to_nearest = 0.5f;
+	int i;
+	__asm
+	{
+		fld x
+		fadd st, st (0)
+		fadd round_to_nearest
+		fistp i
+		sar i, 1
+	}
+	return (i);
+}
+#endif // BOGUS_ROUND
+
+inline F32 llround( F32 val, F32 nearest )
+{
+	return F32(floor(val * (1.0f / nearest) + 0.5f)) * nearest;
+}
+
+inline F64 llround( F64 val, F64 nearest )
+{
+	return F64(floor(val * (1.0 / nearest) + 0.5)) * nearest;
+}
+
+// these provide minimum peak error
+//
+// avg  error = -0.013049 
+// peak error = -31.4 dB
+// RMS  error = -28.1 dB
+
+const F32 FAST_MAG_ALPHA = 0.960433870103f;
+const F32 FAST_MAG_BETA = 0.397824734759f;
+
+// these provide minimum RMS error
+//
+// avg  error = 0.000003 
+// peak error = -32.6 dB
+// RMS  error = -25.7 dB
+//
+//const F32 FAST_MAG_ALPHA = 0.948059448969f;
+//const F32 FAST_MAG_BETA = 0.392699081699f;
+
+inline F32 fastMagnitude(F32 a, F32 b)
+{ 
+	a = (a > 0) ? a : -a;
+	b = (b > 0) ? b : -b;
+	return(FAST_MAG_ALPHA * llmax(a,b) + FAST_MAG_BETA * llmin(a,b));
+}
+
+
+
+////////////////////
+//
+// Fast F32/S32 conversions
+//
+// Culled from www.stereopsis.com/FPU.html
+
+const F64 LL_DOUBLE_TO_FIX_MAGIC	= 68719476736.0*1.5;     //2^36 * 1.5,  (52-_shiftamt=36) uses limited precisicion to floor
+const S32 LL_SHIFT_AMOUNT			= 16;                    //16.16 fixed point representation,
+
+// Endian dependent code
+#ifdef LL_LITTLE_ENDIAN
+	#define LL_EXP_INDEX				1
+	#define LL_MAN_INDEX				0
+#else
+	#define LL_EXP_INDEX				0
+	#define LL_MAN_INDEX				1
+#endif
+
+/* Deprecated: use llround(), lltrunc(), or llfloor() instead
+// ================================================================================================
+// Real2Int
+// ================================================================================================
+inline S32 F64toS32(F64 val)
+{
+	val		= val + LL_DOUBLE_TO_FIX_MAGIC;
+	return ((S32*)&val)[LL_MAN_INDEX] >> LL_SHIFT_AMOUNT; 
+}
+
+// ================================================================================================
+// Real2Int
+// ================================================================================================
+inline S32 F32toS32(F32 val)
+{
+	return F64toS32 ((F64)val);
+}
+*/
+
+////////////////////////////////////////////////
+//
+// Fast exp and log
+//
+
+// Implementation of fast exp() approximation (from a paper by Nicol N. Schraudolph
+// http://www.inf.ethz.ch/~schraudo/pubs/exp.pdf
+static union
+{
+	double d;
+	struct
+	{
+#ifdef LL_LITTLE_ENDIAN
+		S32 j, i;
+#else
+		S32 i, j;
+#endif
+	} n;
+} LLECO; // not sure what the name means
+
+#define LL_EXP_A (1048576 * OO_LN2) // use 1512775 for integer
+#define LL_EXP_C (60801)			// this value of C good for -4 < y < 4
+
+#define LL_FAST_EXP(y) (LLECO.n.i = llround(F32(LL_EXP_A*(y))) + (1072693248 - LL_EXP_C), LLECO.d)
+
+
+
+inline F32 llfastpow(const F32 x, const F32 y)
+{
+	return (F32)(LL_FAST_EXP(y * log(x)));
+}
+
+
+inline F32 snap_to_sig_figs(F32 foo, S32 sig_figs)
+{
+	// compute the power of ten
+	F32 bar = 1.f;
+	for (S32 i = 0; i < sig_figs; i++)
+	{
+		bar *= 10.f;
+	}
+
+	foo = (F32)llround(foo * bar);
+
+	// shift back
+	foo /= bar;
+	return foo;
+}
+
+inline F32 lerp(F32 a, F32 b, F32 u) 
+{
+	return a + ((b - a) * u);
+}
+
+inline F32 lerp2d(F32 x00, F32 x01, F32 x10, F32 x11, F32 u, F32 v)
+{
+	F32 a = x00 + (x01-x00)*u;
+	F32 b = x10 + (x11-x10)*u;
+	F32 r = a + (b-a)*v;
+	return r;
+}
+
+inline F32 ramp(F32 x, F32 a, F32 b)
+{
+	return (a == b) ? 0.0f : ((a - x) / (a - b));
+}
+
+inline F32 rescale(F32 x, F32 x1, F32 x2, F32 y1, F32 y2)
+{
+	return lerp(y1, y2, ramp(x, x1, x2));
+}
+
+inline F32 clamp_rescale(F32 x, F32 x1, F32 x2, F32 y1, F32 y2)
+{
+	if (y1 < y2)
+	{
+		return llclamp(rescale(x,x1,x2,y1,y2),y1,y2);
+	}
+	else
+	{
+		return llclamp(rescale(x,x1,x2,y1,y2),y2,y1);
+	}
+}
+
+
+inline F32 cubic_step( F32 x, F32 x0, F32 x1, F32 s0, F32 s1 )
+{
+	if (x <= x0)
+		return s0;
+
+	if (x >= x1)
+		return s1;
+
+	F32 f = (x - x0) / (x1 - x0);
+
+	return	s0 + (s1 - s0) * (f * f) * (3.0f - 2.0f * f);
+}
+
+inline F32 cubic_step( F32 x )
+{
+	x = llclampf(x);
+
+	return	(x * x) * (3.0f - 2.0f * x);
+}
+
+inline F32 quadratic_step( F32 x, F32 x0, F32 x1, F32 s0, F32 s1 )
+{
+	if (x <= x0)
+		return s0;
+
+	if (x >= x1)
+		return s1;
+
+	F32 f = (x - x0) / (x1 - x0);
+	F32 f_squared = f * f;
+
+	return	(s0 * (1.f - f_squared)) + ((s1 - s0) * f_squared);
+}
+
+inline F32 llsimple_angle(F32 angle)
+{
+	while(angle <= -F_PI)
+		angle += F_TWO_PI;
+	while(angle >  F_PI)
+		angle -= F_TWO_PI;
+	return angle;
+}
+
+#endif