Index: ps/trunk/source/maths/Fixed.cpp
===================================================================
--- ps/trunk/source/maths/Fixed.cpp (revision 7689)
+++ ps/trunk/source/maths/Fixed.cpp (revision 7690)
@@ -1,140 +1,183 @@
/* Copyright (C) 2010 Wildfire Games.
* This file is part of 0 A.D.
*
* 0 A.D. is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 2 of the License, or
* (at your option) any later version.
*
* 0 A.D. is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with 0 A.D. If not, see .
*/
#include "precompiled.h"
#include "Fixed.h"
#include "ps/CStr.h"
template<>
CFixed_15_16 CFixed_15_16::FromString(const CStr8& s)
{
// Parse a superset of the xsd:decimal syntax: [-+]?\d*(\.\d*)?
if (s.empty())
return CFixed_15_16::Zero();
bool neg = false;
CFixed_15_16 r;
const char* c = &s[0];
if (*c == '+')
{
++c;
}
else if (*c == '-')
{
++c;
neg = true;
}
while (true)
{
// Integer part:
if (*c >= '0' && *c <= '9')
{
r = r * 10; // TODO: handle overflow gracefully, maybe
r += CFixed_15_16::FromInt(*c - '0');
++c;
}
else if (*c == '.')
{
++c;
u32 frac = 0;
u32 div = 1;
// Fractional part
while (*c >= '0' && *c <= '9')
{
frac *= 10;
frac += (*c - '0');
div *= 10;
++c;
if (div >= 100000)
{
// any further digits will be too small to have any major effect
break;
}
}
// too many digits or invalid character or end of string - add the fractional part and stop
r += CFixed_15_16(((u64)frac << 16) / div);
break;
}
else
{
// invalid character or end of string
break;
}
}
return (neg ? -r : r);
}
template<>
CFixed_15_16 CFixed_15_16::FromString(const CStrW& s)
{
return FromString(CStr8(s));
}
// Based on http://www.dspguru.com/dsp/tricks/fixed-point-atan2-with-self-normalization
CFixed_15_16 atan2_approx(CFixed_15_16 y, CFixed_15_16 x)
{
CFixed_15_16 zero;
// Special case to avoid division-by-zero
if (x.IsZero() && y.IsZero())
return zero;
CFixed_15_16 c1;
c1.SetInternalValue(51472); // pi/4 << 16
CFixed_15_16 c2;
c2.SetInternalValue(154415); // 3*pi/4 << 16
CFixed_15_16 abs_y = y.Absolute();
CFixed_15_16 angle;
if (x >= zero)
{
CFixed_15_16 r = (x - abs_y) / (x + abs_y);
angle = c1 - c1.Multiply(r);
}
else
{
CFixed_15_16 r = (x + abs_y) / (abs_y - x);
angle = c2 - c1.Multiply(r);
}
if (y < zero)
return -angle;
else
return angle;
}
template<>
CFixed_15_16 CFixed_15_16::Pi()
{
return CFixed_15_16(205887); // = pi << 16
}
void sincos_approx(CFixed_15_16 a, CFixed_15_16& sin_out, CFixed_15_16& cos_out)
{
- // XXX: mustn't use floating-point here - need a fixed-point emulation
+ // Based on http://www.coranac.com/2009/07/sines/
- sin_out = CFixed_15_16::FromDouble(sin(a.ToDouble()));
- cos_out = CFixed_15_16::FromDouble(cos(a.ToDouble()));
+ // TODO: this could be made a bit more precise by being careful about scaling
+
+ typedef CFixed_15_16 fixed;
+
+ fixed c2_pi;
+ c2_pi.SetInternalValue(41721); // = 2/pi << 16
+
+ // Map radians onto the range [0, 4)
+ fixed z = a.Multiply(c2_pi) % fixed::FromInt(4);
+
+ // Map z onto the range [-1, +1] for sin, and the same with z+1 to compute cos
+ fixed sz, cz;
+ if (z >= fixed::FromInt(3)) // [3, 4)
+ {
+ sz = z - fixed::FromInt(4);
+ cz = z - fixed::FromInt(3);
+ }
+ else if (z >= fixed::FromInt(2)) // [2, 3)
+ {
+ sz = fixed::FromInt(2) - z;
+ cz = z - fixed::FromInt(3);
+ }
+ else if (z >= fixed::FromInt(1)) // [1, 2)
+ {
+ sz = fixed::FromInt(2) - z;
+ cz = fixed::FromInt(1) - z;
+ }
+ else // [0, 1)
+ {
+ sz = z;
+ cz = fixed::FromInt(1) - z;
+ }
+
+ // Third-order (max absolute error ~0.02)
+
+// sin_out = (sz / 2).Multiply(fixed::FromInt(3) - sz.Multiply(sz));
+// cos_out = (cz / 2).Multiply(fixed::FromInt(3) - cz.Multiply(cz));
+
+ // Fifth-order (max absolute error ~0.0005)
+
+ fixed sz2 = sz.Multiply(sz);
+ sin_out = sz.Multiply(fixed::Pi() - sz2.Multiply(fixed::Pi()*2 - fixed::FromInt(5) - sz2.Multiply(fixed::Pi() - fixed::FromInt(3)))) / 2;
+
+ fixed cz2 = cz.Multiply(cz);
+ cos_out = cz.Multiply(fixed::Pi() - cz2.Multiply(fixed::Pi()*2 - fixed::FromInt(5) - cz2.Multiply(fixed::Pi() - fixed::FromInt(3)))) / 2;
}
Index: ps/trunk/source/maths/Fixed.h
===================================================================
--- ps/trunk/source/maths/Fixed.h (revision 7689)
+++ ps/trunk/source/maths/Fixed.h (revision 7690)
@@ -1,313 +1,326 @@
/* Copyright (C) 2010 Wildfire Games.
* This file is part of 0 A.D.
*
* 0 A.D. is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 2 of the License, or
* (at your option) any later version.
*
* 0 A.D. is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with 0 A.D. If not, see .
*/
#ifndef INCLUDED_FIXED
#define INCLUDED_FIXED
#include "lib/types.h"
#include "maths/Sqrt.h"
class CStr8;
class CStrW;
#ifndef NDEBUG
#define USE_FIXED_OVERFLOW_CHECKS
#endif // NDEBUG
//define overflow macros
#ifndef USE_FIXED_OVERFLOW_CHECKS
#define CheckSignedSubtractionOverflow(type, left, right, overflowWarning, underflowWarning)
#define CheckSignedAdditionOverflow(type, left, right, overflowWarning, underflowWarning)
#define CheckCastOverflow(var, targetType, overflowWarning, underflowWarning)
#define CheckU32CastOverflow(var, targetType, overflowWarning)
#define CheckUnsignedAdditionOverflow(result, operand, overflowWarning)
#define CheckUnsignedSubtractionOverflow(result, operand, overflowWarning)
#define CheckNegationOverflow(var, type, overflowWarning)
#define CheckMultiplicationOverflow(type, left, right, overflowWarning, underflowWarning)
#define CheckDivisionOverflow(type, left, right, overflowWarning)
#else // USE_FIXED_OVERFLOW_CHECKS
#define CheckSignedSubtractionOverflow(type, left, right, overflowWarning, underflowWarning) \
if(left > 0 && right < 0 && left > std::numeric_limits::max() + right) \
debug_warn(overflowWarning); \
else if(left < 0 && right > 0 && left < std::numeric_limits::min() + right) \
debug_warn(underflowWarning);
#define CheckSignedAdditionOverflow(type, left, right, overflowWarning, underflowWarning) \
if(left > 0 && right > 0 && std::numeric_limits::max() - left < right) \
debug_warn(overflowWarning); \
else if(left < 0 && right < 0 && std::numeric_limits::min() - left > right) \
debug_warn(underflowWarning);
#define CheckCastOverflow(var, targetType, overflowWarning, underflowWarning) \
if(var > std::numeric_limits::max()) \
debug_warn(overflowWarning); \
else if(var < std::numeric_limits::min()) \
debug_warn(underflowWarning);
#define CheckU32CastOverflow(var, targetType, overflowWarning) \
if(var > (u32)std::numeric_limits::max()) \
debug_warn(overflowWarning);
#define CheckUnsignedAdditionOverflow(result, operand, overflowWarning) \
if(result < operand) \
debug_warn(overflowWarning);
#define CheckUnsignedSubtractionOverflow(result, left, overflowWarning) \
if(result > left) \
debug_warn(overflowWarning);
#define CheckNegationOverflow(var, type, overflowWarning) \
if(value == std::numeric_limits::min()) \
debug_warn(overflowWarning);
#define CheckMultiplicationOverflow(type, left, right, overflowWarning, underflowWarning) \
i64 res##left = (i64)left * (i64)right; \
CheckCastOverflow(res##left, type, overflowWarning, underflowWarning)
#define CheckDivisionOverflow(type, left, right, overflowWarning) \
if(right == -1) { CheckNegationOverflow(left, type, overflowWarning) }
#endif // USE_FIXED_OVERFLOW_CHECKS
template
inline T round_away_from_zero(float value)
{
return (T)(value >= 0 ? value + 0.5f : value - 0.5f);
}
/**
* A simple fixed-point number class.
*
* Use 'fixed' rather than using this class directly.
*/
template
class CFixed
{
private:
T value;
explicit CFixed(T v) : value(v) { }
public:
enum { fract_bits = fract_bits_ };
CFixed() : value(0) { }
static CFixed Zero() { return CFixed(0); }
static CFixed Pi();
T GetInternalValue() const { return value; }
void SetInternalValue(T n) { value = n; }
// Conversion to/from primitive types:
static CFixed FromInt(int n)
{
return CFixed(n << fract_bits);
}
static CFixed FromFloat(float n)
{
if (!isfinite(n))
return CFixed(0);
float scaled = n * fract_pow2;
return CFixed(round_away_from_zero(scaled));
}
static CFixed FromDouble(double n)
{
if (!isfinite(n))
return CFixed(0);
double scaled = n * fract_pow2;
return CFixed(round_away_from_zero(scaled));
}
static CFixed FromString(const CStr8& s);
static CFixed FromString(const CStrW& s);
float ToFloat() const
{
return value / (float)fract_pow2;
}
double ToDouble() const
{
return value / (double)fract_pow2;
}
int ToInt_RoundToZero() const
{
if (value > 0)
return value >> fract_bits;
else
return (value + fract_pow2 - 1) >> fract_bits;
}
/// Returns true if the number is precisely 0.
bool IsZero() const { return value == 0; }
/// Equality.
bool operator==(CFixed n) const { return (value == n.value); }
/// Inequality.
bool operator!=(CFixed n) const { return (value != n.value); }
/// Numeric comparison.
bool operator<=(CFixed n) const { return (value <= n.value); }
/// Numeric comparison.
bool operator<(CFixed n) const { return (value < n.value); }
/// Numeric comparison.
bool operator>=(CFixed n) const { return (value >= n.value); }
/// Numeric comparison.
bool operator>(CFixed n) const { return (value > n.value); }
// Basic arithmetic:
/// Add a CFixed. Might overflow.
CFixed operator+(CFixed n) const
{
CheckSignedAdditionOverflow(T, value, n.value, L"Overflow in CFixed::operator+(CFixed n)", L"Underflow in CFixed::operator+(CFixed n)")
return CFixed(value + n.value);
}
/// Subtract a CFixed. Might overflow.
CFixed operator-(CFixed n) const
{
CheckSignedSubtractionOverflow(T, value, n.value, L"Overflow in CFixed::operator-(CFixed n)", L"Underflow in CFixed::operator-(CFixed n)")
return CFixed(value - n.value);
}
/// Add a CFixed. Might overflow.
CFixed& operator+=(CFixed n) { *this = *this + n; return *this; }
/// Subtract a CFixed. Might overflow.
CFixed& operator-=(CFixed n) { *this = *this - n; return *this; }
/// Negate a CFixed.
CFixed operator-() const
{
CheckNegationOverflow(value, T, L"Overflow in CFixed::operator-()")
return CFixed(-value);
}
/// Divide by a CFixed. Must not have n.IsZero(). Might overflow.
CFixed operator/(CFixed n) const
{
i64 t = (i64)value << fract_bits;
i64 result = t / (i64)n.value;
CheckCastOverflow(result, T, L"Overflow in CFixed::operator/(CFixed n)", L"Underflow in CFixed::operator/(CFixed n)")
return CFixed((T)result);
}
/// Multiply by an integer. Might overflow.
CFixed operator*(int n) const
{
CheckMultiplicationOverflow(T, value, n, L"Overflow in CFixed::operator*(int n)", L"Underflow in CFixed::operator*(int n)")
return CFixed(value * n);
}
/// Divide by an integer. Must not have n == 0. Cannot overflow unless n == -1.
CFixed operator/(int n) const
{
CheckDivisionOverflow(T, value, n, L"Overflow in CFixed::operator/(int n)")
return CFixed(value / n);
}
+ /// Mod by a fixed. Must not have n == 0. Result has the same sign as n.
+ CFixed operator%(CFixed n) const
+ {
+ T t = value % n.value;
+ if (n.value > 0 && t < 0)
+ t += n.value;
+ else if (n.value < 0 && t > 0)
+ t += n.value;
+
+ return CFixed(t);
+ }
+
CFixed Absolute() const { return CFixed(abs(value)); }
/**
* Multiply by a CFixed. Likely to overflow if both numbers are large,
* so we use an ugly name instead of operator* to make it obvious.
*/
CFixed Multiply(CFixed n) const
{
i64 t = (i64)value * (i64)n.value;
t >>= fract_bits;
CheckCastOverflow(t, T, L"Overflow in CFixed::Multiply(CFixed n)", L"Underflow in CFixed::Multiply(CFixed n)")
return CFixed((T)t);
}
/**
* Compute this*m/d. Must not have d == 0. Won't overflow if the result can be represented as a CFixed.
*/
CFixed MulDiv(CFixed m, CFixed d) const
{
i64 t = ((i64)value * (i64)m.value) / (i64)d.value;
CheckCastOverflow(t, T, L"Overflow in CFixed::Multiply(CFixed n)", L"Underflow in CFixed::Multiply(CFixed n)")
return CFixed((T)t);
}
CFixed Sqrt() const
{
if (value <= 0)
return CFixed(0);
u32 s = isqrt64((u64)value << fract_bits);
return CFixed(s);
}
private:
// Prevent dangerous accidental implicit conversions of floats to ints in certain operations
CFixed operator*(float n) const;
CFixed operator/(float n) const;
};
/**
* A fixed-point number class with 1-bit sign, 15-bit integral part, 16-bit fractional part.
*/
typedef CFixed CFixed_15_16;
/**
* Default fixed-point type used by the engine.
*/
typedef CFixed_15_16 fixed;
namespace std
{
/**
* std::numeric_limits specialisation, currently just providing min and max
*/
template
class numeric_limits >
{
typedef CFixed fixed;
public:
static const bool is_specialized = true;
static fixed min() throw() { fixed f; f.SetInternalValue(std::numeric_limits::min()); return f; }
static fixed max() throw() { fixed f; f.SetInternalValue(std::numeric_limits::max()); return f; }
};
}
/**
* Inaccurate approximation of atan2 over fixed-point numbers.
* Maximum error is almost 0.08 radians (4.5 degrees).
*/
CFixed_15_16 atan2_approx(CFixed_15_16 y, CFixed_15_16 x);
/**
* Compute sin(a) and cos(a).
+ * Maximum error for -2pi < a < 2pi is almost 0.0005.
*/
void sincos_approx(CFixed_15_16 a, CFixed_15_16& sin_out, CFixed_15_16& cos_out);
#endif // INCLUDED_FIXED
Index: ps/trunk/source/maths/tests/test_Fixed.h
===================================================================
--- ps/trunk/source/maths/tests/test_Fixed.h (revision 7689)
+++ ps/trunk/source/maths/tests/test_Fixed.h (revision 7690)
@@ -1,191 +1,249 @@
/* Copyright (C) 2010 Wildfire Games.
* This file is part of 0 A.D.
*
* 0 A.D. is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 2 of the License, or
* (at your option) any later version.
*
* 0 A.D. is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with 0 A.D. If not, see .
*/
#include "lib/self_test.h"
#include "maths/Fixed.h"
#include "maths/MathUtil.h"
class TestFixed : public CxxTest::TestSuite
{
public:
void test_basic()
{
fixed a = fixed::FromInt(12345);
TS_ASSERT_EQUALS((a + a).ToDouble(), 24690);
TS_ASSERT_EQUALS((a - a).ToDouble(), 0);
TS_ASSERT_EQUALS((-a).ToDouble(), -12345);
TS_ASSERT_EQUALS((a / a).ToDouble(), 1);
TS_ASSERT_EQUALS((a / 823).ToDouble(), 15);
TS_ASSERT_EQUALS((a * 2).ToDouble(), 24690);
}
void test_FromInt()
{
fixed a = fixed::FromInt(123);
TS_ASSERT_EQUALS(a.ToFloat(), 123.0f);
TS_ASSERT_EQUALS(a.ToDouble(), 123.0);
TS_ASSERT_EQUALS(a.ToInt_RoundToZero(), 123);
}
void test_FromFloat()
{
fixed a = fixed::FromFloat(123.125f);
TS_ASSERT_EQUALS(a.ToFloat(), 123.125f);
TS_ASSERT_EQUALS(a.ToDouble(), 123.125);
fixed b = fixed::FromFloat(-123.125f);
TS_ASSERT_EQUALS(b.ToFloat(), -123.125f);
TS_ASSERT_EQUALS(b.ToDouble(), -123.125);
fixed c = fixed::FromFloat(INFINITY);
TS_ASSERT_EQUALS(c.GetInternalValue(), (i32)0);
fixed d = fixed::FromFloat(-INFINITY);
TS_ASSERT_EQUALS(d.GetInternalValue(), (i32)0);
fixed e = fixed::FromFloat(NAN);
TS_ASSERT_EQUALS(e.GetInternalValue(), (i32)0);
}
void test_FromDouble()
{
fixed a = fixed::FromDouble(123.125);
TS_ASSERT_EQUALS(a.ToFloat(), 123.125f);
TS_ASSERT_EQUALS(a.ToDouble(), 123.125);
fixed b = fixed::FromDouble(-123.125);
TS_ASSERT_EQUALS(b.ToFloat(), -123.125f);
TS_ASSERT_EQUALS(b.ToDouble(), -123.125);
fixed c = fixed::FromDouble(INFINITY);
TS_ASSERT_EQUALS(c.GetInternalValue(), (i32)0);
fixed d = fixed::FromDouble(-INFINITY);
TS_ASSERT_EQUALS(d.GetInternalValue(), (i32)0);
fixed e = fixed::FromDouble(NAN);
TS_ASSERT_EQUALS(e.GetInternalValue(), (i32)0);
}
void test_FromFloat_Rounding()
{
double eps = pow(2.0, -16.0);
TS_ASSERT_EQUALS(fixed::FromFloat(eps).ToDouble(), eps);
TS_ASSERT_EQUALS(fixed::FromFloat(eps * 0.5625).ToDouble(), eps);
TS_ASSERT_EQUALS(fixed::FromFloat(eps * 0.5).ToDouble(), eps);
TS_ASSERT_EQUALS(fixed::FromFloat(eps * 0.4375).ToDouble(), 0.0);
TS_ASSERT_EQUALS(fixed::FromFloat(-eps).ToDouble(), -eps);
TS_ASSERT_EQUALS(fixed::FromFloat(-eps * 0.5625).ToDouble(), -eps);
TS_ASSERT_EQUALS(fixed::FromFloat(-eps * 0.5).ToDouble(), -eps);
TS_ASSERT_EQUALS(fixed::FromFloat(-eps * 0.4375).ToDouble(), 0.0);
}
void test_FromString()
{
TS_ASSERT_EQUALS(fixed::FromString("").ToDouble(), 0.0);
TS_ASSERT_EQUALS(fixed::FromString("123").ToDouble(), 123.0);
TS_ASSERT_EQUALS(fixed::FromString("+123").ToDouble(), 123.0);
TS_ASSERT_EQUALS(fixed::FromString("-123").ToDouble(), -123.0);
TS_ASSERT_EQUALS(fixed::FromString("123.5").ToDouble(), 123.5);
TS_ASSERT_EQUALS(fixed::FromString("-123.5").ToDouble(), -123.5);
TS_ASSERT_EQUALS(fixed::FromString(".5").ToDouble(), 0.5);
TS_ASSERT_EQUALS(fixed::FromString("5.").ToDouble(), 5.0);
TS_ASSERT_EQUALS(fixed::FromString("254.391").GetInternalValue(), 16671768);
TS_ASSERT_EQUALS(fixed::FromString("123x456").ToDouble(), 123.0);
TS_ASSERT_EQUALS(fixed::FromString("1.00002").ToDouble(), 1.0 + 1.0/65536.0);
TS_ASSERT_EQUALS(fixed::FromString("1.0000200000000000000000").ToDouble(), 1.0 + 1.0/65536.0);
TS_ASSERT_EQUALS(fixed::FromString("1.000009").ToDouble(), 1.0);
TS_ASSERT_EQUALS(fixed::FromString("0.9999999999999999999999").ToDouble(), 1.0 - 1.0/65536.0);
// TODO: could test more large/small numbers, errors, etc
}
void test_RoundToZero()
{
TS_ASSERT_EQUALS(fixed::FromFloat(10.f).ToInt_RoundToZero(), 10);
TS_ASSERT_EQUALS(fixed::FromFloat(10.1f).ToInt_RoundToZero(), 10);
TS_ASSERT_EQUALS(fixed::FromFloat(10.5f).ToInt_RoundToZero(), 10);
TS_ASSERT_EQUALS(fixed::FromFloat(10.99f).ToInt_RoundToZero(), 10);
TS_ASSERT_EQUALS(fixed::FromFloat(0.1f).ToInt_RoundToZero(), 0);
TS_ASSERT_EQUALS(fixed::FromFloat(0.0f).ToInt_RoundToZero(), 0);
TS_ASSERT_EQUALS(fixed::FromFloat(-0.1f).ToInt_RoundToZero(), 0);
TS_ASSERT_EQUALS(fixed::FromFloat(-0.99f).ToInt_RoundToZero(), 0);
TS_ASSERT_EQUALS(fixed::FromFloat(-1.0f).ToInt_RoundToZero(), -1);
TS_ASSERT_EQUALS(fixed::FromFloat(-2.0f).ToInt_RoundToZero(), -2);
TS_ASSERT_EQUALS(fixed::FromFloat(-2.5f).ToInt_RoundToZero(), -2);
TS_ASSERT_EQUALS(fixed::FromFloat(-2.99f).ToInt_RoundToZero(), -2);
}
// TODO: test all the arithmetic operators
+ void test_Mod()
+ {
+ TS_ASSERT_EQUALS(fixed::FromInt(0) % fixed::FromInt(4), fixed::FromInt(0));
+ TS_ASSERT_EQUALS(fixed::FromInt(0) % fixed::FromInt(-4), fixed::FromInt(0));
+ TS_ASSERT_EQUALS(fixed::FromInt(5) % fixed::FromInt(4), fixed::FromInt(1));
+ TS_ASSERT_EQUALS(fixed::FromInt(5) % fixed::FromInt(-4), fixed::FromInt(-3));
+ TS_ASSERT_EQUALS(fixed::FromInt(-5) % fixed::FromInt(4), fixed::FromInt(3));
+ TS_ASSERT_EQUALS(fixed::FromInt(-5) % fixed::FromInt(-4), fixed::FromInt(-1));
+
+ TS_ASSERT_EQUALS((fixed::FromDouble(5.5) % fixed::FromInt(4)).ToDouble(), 1.5);
+
+ TS_ASSERT_EQUALS((fixed::FromDouble(1.75) % fixed::FromDouble(0.5)).ToDouble(), 0.25);
+ }
+
void test_Sqrt()
{
TS_ASSERT_EQUALS(fixed::FromDouble(1.0).Sqrt().ToDouble(), 1.0);
TS_ASSERT_DELTA(fixed::FromDouble(2.0).Sqrt().ToDouble(), sqrt(2.0), 1.0/65536.0);
TS_ASSERT_EQUALS(fixed::FromDouble(32400.0).Sqrt().ToDouble(), 180.0);
TS_ASSERT_EQUALS(fixed::FromDouble(0.0625).Sqrt().ToDouble(), 0.25);
TS_ASSERT_EQUALS(fixed::FromDouble(-1.0).Sqrt().ToDouble(), 0.0);
TS_ASSERT_EQUALS(fixed::FromDouble(0.0).Sqrt().ToDouble(), 0.0);
}
void test_Atan2()
{
// Special cases from atan2 man page:
TS_ASSERT_DELTA(atan2_approx(fixed::FromInt(0), fixed::FromInt(-1)).ToDouble(), M_PI, 0.01);
TS_ASSERT_EQUALS(atan2_approx(fixed::FromInt(0), fixed::FromInt(+1)).ToDouble(), 0);
TS_ASSERT_DELTA(atan2_approx(fixed::FromInt(-1), fixed::FromInt(0)).ToDouble(), -M_PI_2, 0.01);
TS_ASSERT_DELTA(atan2_approx(fixed::FromInt(+1), fixed::FromInt(0)).ToDouble(), +M_PI_2, 0.01);
TS_ASSERT_EQUALS(atan2_approx(fixed::FromInt(0), fixed::FromInt(0)).ToDouble(), 0);
// Test that it approximately matches libc's atan2
#define T(y, x) TS_ASSERT_DELTA(atan2_approx(fixed::FromDouble(y), fixed::FromDouble(x)).ToDouble(), atan2((double)(y), (double)(x)), 0.072)
// Quadrants
T(1, 1);
T(-1, 1);
T(1, -1);
T(-1, -1);
// Scale
T(10, 10);
T(100, 100);
T(1000, 1000);
T(10000, 10000);
// Some arbitrary cases
T(2, 1);
T(3, 1);
T(4, 1);
T(10, 1);
T(1000, 1);
T(1, 2);
T(1, 3);
T(1, 4);
T(1, 10);
T(1, 1000);
// The highest-error region
for (double x = 2.0; x < 4.0; x += 0.01)
{
T(1, x);
}
#undef T
}
+
+ void test_SinCos()
+ {
+ fixed s, c;
+
+ sincos_approx(fixed::FromInt(0), s, c);
+ TS_ASSERT_EQUALS(s, fixed::FromInt(0));
+ TS_ASSERT_EQUALS(c, fixed::FromInt(1));
+
+ sincos_approx(fixed::Pi(), s, c);
+ TS_ASSERT_DELTA(s.ToDouble(), 0.0, 0.00005);
+ TS_ASSERT_EQUALS(c, fixed::FromInt(-1));
+
+ sincos_approx(fixed::FromDouble(M_PI*2.0), s, c);
+ TS_ASSERT_DELTA(s.ToDouble(), 0.0, 0.0001);
+ TS_ASSERT_DELTA(c.ToDouble(), 1.0, 0.0001);
+
+ sincos_approx(fixed::FromDouble(M_PI*100.0), s, c);
+ TS_ASSERT_DELTA(s.ToDouble(), 0.0, 0.004);
+ TS_ASSERT_DELTA(c.ToDouble(), 1.0, 0.004);
+
+ sincos_approx(fixed::FromDouble(M_PI*-100.0), s, c);
+ TS_ASSERT_DELTA(s.ToDouble(), 0.0, 0.004);
+ TS_ASSERT_DELTA(c.ToDouble(), 1.0, 0.004);
+
+/*
+ for (double a = 0.0; a < 6.28; a += 0.1)
+ {
+ sincos_approx(fixed::FromDouble(a), s, c);
+ printf("%f: sin = %f / %f; cos = %f / %f\n", a, s.ToDouble(), sin(a), c.ToDouble(), cos(a));
+ }
+*/
+
+ double err = 0.0;
+ for (double a = -6.28; a < 6.28; a += 0.001)
+ {
+ sincos_approx(fixed::FromDouble(a), s, c);
+ err = std::max(err, fabs(sin(a) - s.ToDouble()));
+ err = std::max(err, fabs(cos(a) - c.ToDouble()));
+ }
+// printf("### Max error %f = %f/2^16\n", err, err*65536.0);
+
+ TS_ASSERT_LESS_THAN(err, 0.00046);
+ }
};