celestar@9876: /* $Id */
celestar@9876: 
celestar@9876: /**
celestar@9876:  * @file
celestar@9876:  * Implements a Fixed-point data type.
celestar@9876:  * This files defines a new data type
celestar@9876:  * "FixedT<Tstorage, Tdec_bits>", which can be used to compute fractional
celestar@9876:  * numbers without having to worry about network stability (because
celestar@9876:  * there is no real rounding taking place.
celestar@9876:  * Use "Fixed" for the defaults values (48 bits integral part, 16 bits
celestar@9876:  * fractional part) or "FixedT<typename, int>" for flexible sizes. "int"
celestar@9876:  * number of bits are taken from a variable of "typename" for the fractional
celestar@9876:  * representation, the rest stays available for the integral part.
celestar@9876:  * @warning More than 16 bits fractional should not be used when working
celestar@9876:  * with large integers.
celestar@9876:  * @warning More than 31 fractional bits are not supported, and will
celestar@9876:  * trigger a compiler warning with appropriate settings.
celestar@9876:  * @warning There are no warnings for buffer overflows. Those are very likely
celestar@9876:  * to happen with divisions, as those first shift the numerator by the number
celestar@9876:  * of fractional bits. Use multiplication with the reciprocal whereever
celestar@9876:  * it is possible.
celestar@9876:  * @warning Please use signed variables for storage at the moment ONLY
celestar@9876:  *
celestar@9876:  * Variables can be initialized:
celestar@9876:  * 	- with integers: "FixedT<int, 9> a = -5"
celestar@9876:  * 	- with fractionals: "FixedT<int64, 11> b(7, 2);" where the first number is
celestar@9876:  * the numerator and the second is the denominator.
celestar@9876:  * 	- not at all: "Fixed c;", which will automatically set c to zero.
celestar@9876:  *
celestar@9876:  * @example fixed_ex.cpp How to initialize Fixed-point data types
celestar@9876:  */
celestar@9876: 
celestar@9876: /* Included to overload the stream operator << */
celestar@9876: #include <iostream>
celestar@9876: /* Needed for some ottd-specific data types */
celestar@9876: #include "stdafx.h"
celestar@9876: 
celestar@9876: /**
celestar@9876:  * Base Class for storing fixed-point data types.
celestar@9876:  * Fixed-point data types are stored in a single integer-type variable
celestar@9876:  * (Tstorage) of which a certain number of bits (Tdec_bits) is used to
celestar@9876:  * represent the fractional part of the number.
celestar@9876:  */
celestar@9876: template <typename Tstorage, int Tdec_bits>
celestar@9876: class FixedRawT {
celestar@9876: public:
celestar@9876: 	/** The number of bits that represent the fractional */
celestar@9876: 	static const int dec_bits = Tdec_bits;
celestar@9876: 
celestar@9876: 	/** The storage of the number itself */
celestar@9876: 	Tstorage m_data;
celestar@9876: 
celestar@9876: 	/**
celestar@9876: 	 * Basic constructor for integer arguments.
celestar@9876: 	 * @param value The raw value to which we want to initialize
celestar@9876: 	 */
celestar@9876: 	template <typename T>
celestar@9876: 	FixedRawT(T value) : m_data(value) {}
celestar@9876: 
celestar@9876: 	/**
celestar@9876: 	 * Make sure that a Fixed variable is always inizialized to zero by default.
celestar@9876: 	 * We need this ctor so that we can declare variable without
celestar@9876: 	 * the need to initialize it explicitly
celestar@9876: 	 */
celestar@9876: 	FixedRawT() : m_data (0) {}
celestar@9876: 
celestar@9876: 	/**
celestar@9876: 	 * Returns the biggest integral number we can represent
celestar@9876: 	 */
celestar@9876: 	int64 FIXED_MAX() const { return m_maximum; }
celestar@9876: 
celestar@9876: 	/**
celestar@9876: 	 * Returns the biggest integral number we can represent
celestar@9876: 	 */
celestar@9876: 	int64 FIXED_MIN() const { return m_minimum; }
celestar@9876: 
celestar@9876: private:
celestar@9876: 	/** The largest number we can hold */
celestar@9876: 	static const int64 m_maximum = +(1ULL << ((sizeof(Tstorage) * 8) - Tdec_bits - 1)) - 1;
celestar@9876: 
celestar@9876: 	/** The smallest number we can hold */
celestar@9876: 	static const int64 m_minimum = -(1ULL << ((sizeof(Tstorage) * 8) - Tdec_bits - 1)) - 0;
celestar@9876: };
celestar@9876: 
celestar@9876: /* forward-declare some structs */
celestar@9876: template <typename Tstorage, int Tdec_bits> class FixedT;
celestar@9876: template <typename T> struct FixedHelperT;
celestar@9876: 
celestar@9876: /**
celestar@9876:  * Specialization of FixedHelperT. Used to decimal-align two variables of Fixed type.
celestar@9876:  */
celestar@9876: template <typename Tstorage, int Tdec_bits> struct FixedHelperT<FixedRawT<Tstorage, Tdec_bits> >
celestar@9876: {
celestar@9876: 	/** The number of bits used for the fraction */
celestar@9876: 	static const int dec_bits = Tdec_bits;
celestar@9876: 
celestar@9876: 	/**
celestar@9876: 	 * Gives the raw data of a FixedRaw
celestar@9876: 	 * @param t The number to return
celestar@9876: 	 */
celestar@9876: 	static int64 Raw(const FixedRawT<Tstorage, Tdec_bits>& t) {return (int64)t.m_data;}
celestar@9876: 
celestar@9876: 	/**
celestar@9876: 	 * Returns the data from FixedRawT aligned so that it is aligned to a given number
celestar@9876: 	 * @param t    The number to be aligned
celestar@9876: 	 * @param bits The number of bits in the fractional part to be aligned to
celestar@9876: 	 */
celestar@9876: 	static int64 Raw(const FixedRawT<Tstorage, Tdec_bits>& t, int bits)
celestar@9876: 	{
celestar@9876: 		return (((int64)t.m_data) << (bits > Tdec_bits ? bits - Tdec_bits : 0)) >> (bits < Tdec_bits ? Tdec_bits - bits : 0);
celestar@9876: 	}
celestar@9876: };
celestar@9876: 
celestar@9876: /**
celestar@9876:  * Specialization of FixedHelperT. Used to decimal-align two variables of Fixed type.
celestar@9876:  */
celestar@9876: template <typename Tstorage, int Tdec_bits> struct FixedHelperT<FixedT<Tstorage, Tdec_bits> > : public FixedHelperT<FixedRawT<Tstorage, Tdec_bits> > {};
celestar@9876: 
celestar@9876: /**
celestar@9876:  * General implementation of FixedHelperT. This tempate makes sure that
celestar@9876:  * some number supplied is properly aligned at the decimal
celestar@9876:  */
celestar@9876: template <typename T> struct FixedHelperT
celestar@9876: {
celestar@9876: 	/** This version of FixedHelperT is only used for full integers, so assume the number of dec_bits to be zero */
celestar@9876: 	static const int dec_bits = 0;
celestar@9876: 	/** Converts a full integer to one that has a given number of fractional bits, assumes
celestar@9876: 	 * the number of fractional bits to be zero.
celestar@9876: 	 * @param t the number to be converted
celestar@9876: 	 */
celestar@9876: 	static int64 Raw(const T& t) {return (int64)t;}
celestar@9876: 
celestar@9876: 	/**
celestar@9876: 	 * Converts a full integer to one that has a given number of fractional bits
celestar@9876: 	 * @param bits the number of bits for the fraction
celestar@9876: 	 * @param t    the number to be converted
celestar@9876: 	 */
celestar@9876: 	static int64 Raw(const T& t, int bits) {return ((int64)t) << bits;}
celestar@9876: };
celestar@9876: 
celestar@9876: 
celestar@9876: /**
celestar@9876:  * A class that defines a fixed-point data type, which a variable length and precision.
celestar@9876:  * The data type that is defined is a fixed-point variable that has and number of Tdec_bits
celestar@9876:  * bits to represent the fractional part, the remaining bits of Tstorage are then used for
celestar@9876:  * the integer part of the number. This means, we always need to make a trade-off between
celestar@9876:  * the precision we want (higher number in Tdec_bits) and the range (smaller number of
celestar@9876:  * Tdec_bits). This class here only defines the constructors and the operators, everything
celestar@9876:  * else is happning in helper classes.
celestar@9876:  * @note for usage of fprintf and similar, explicit casts are needed (double or int64).
celestar@9876:  */
celestar@9876: template <typename Tstorage, int Tdec_bits> class FixedT : public FixedRawT<Tstorage, Tdec_bits> {
celestar@9876: /* private block up here, because we need the typedef later on */
celestar@9876: private:
celestar@9876: 	/** We shortcut the underlying data type to "Raw", to save typing */
celestar@9876: 	typedef FixedRawT<Tstorage, Tdec_bits> Raw;
celestar@9876: 
celestar@9876: public:
celestar@9876: 	/**
celestar@9876: 	 * Ctor for assignment with other, non floating, variable types
celestar@9876: 	 * @param value The Value we should assign to the Fixed
celestar@9876: 	 */
celestar@9876: 	template <typename T> FixedT(T value) : Raw(FixedHelperT<T>::Raw(value, Tdec_bits)) {}
celestar@9876: 
celestar@9876: 	/** As a constructor without initializing, just use the one that Raw uses. Sets a variable to zero */
celestar@9876: 	FixedT() : Raw() {}
celestar@9876: 
celestar@9876: 	/**
celestar@9876: 	 * Ctor for use with a fraction, useful for initing a variable to a non-integer number at declaration.
celestar@9876: 	 * @param numerator   The Nominator of the fraction
celestar@9876: 	 * @param denominator The Denominator of the fraction
celestar@9876: 	 */
celestar@9876: 	FixedT(int numerator, int denominator) {Raw::m_data = ((int64)numerator << Raw::dec_bits) / (int64)denominator;}
celestar@9876: 
celestar@9876: 	/**
celestar@9879: 	 * Cast to a double (for outputs)
celestar@9879: 	 */
celestar@9879: 	operator double() const
celestar@9879: 	{
celestar@9879: 		return (double)Raw::m_data / (1 << Raw::dec_bits);
celestar@9879: 	}
celestar@9879: 
celestar@9879: 	/**
celestar@9876: 	 * Equality operator
celestar@9876: 	 * @param comparator The non-floating point variable we want to compare against
celestar@9876: 	 */
celestar@9876: 	template <typename T>
celestar@9876: 	bool operator ==(const T &comparator) const { return Raw::m_data == (Tstorage)FixedHelperT<T>::Raw(comparator, Tdec_bits); }
celestar@9876: 
celestar@9876: 	/**
celestar@9876: 	 * Inequality operator
celestar@9876: 	 * @param comparator The non-floating point variable we want to compare against
celestar@9876: 	 */
celestar@9876: 	template <typename T>
celestar@9876: 	bool operator !=(const T &comparator) const { return Raw::m_data != (Tstorage)FixedHelperT<T>::Raw(comparator, Tdec_bits); }
celestar@9876: 
celestar@9876: 	/**
celestar@9876: 	 * Greater or equal operator
celestar@9876: 	 * @param comparator The non-floating point variable we want to compare against
celestar@9876: 	 */
celestar@9876: 	template <typename T>
celestar@9876: 	bool operator >=(const T &comparator) const { return Raw::m_data >= (Tstorage)FixedHelperT<T>::Raw(comparator, Tdec_bits); }
celestar@9876: 
celestar@9876: 	/**
celestar@9876: 	 * Less or equal operator
celestar@9876: 	 * @param comparator The non-floating point variable we want to compare against
celestar@9876: 	 */
celestar@9876: 	template <typename T>
celestar@9876: 	bool operator <=(const T &comparator) const { return Raw::m_data <= (Tstorage)FixedHelperT<T>::Raw(comparator, Tdec_bits); }
celestar@9876: 
celestar@9876: 	/**
celestar@9876: 	 * Greater than operator
celestar@9876: 	 * @param comparator The non-floating point variable we want to compare against
celestar@9876: 	 */
celestar@9876: 	template <typename T>
celestar@9876: 	bool operator >(const T &comparator) const { return Raw::m_data > (Tstorage)FixedHelperT<T>::Raw(comparator, Tdec_bits); }
celestar@9876: 
celestar@9876: 	/**
celestar@9876: 	 * Less than operator
celestar@9876: 	 * @param comparator The non-floating point variable we want to compare against
celestar@9876: 	 */
celestar@9876: 	template <typename T>
celestar@9876: 	bool operator <(const T &comparator) const { return Raw::m_data < (Tstorage)FixedHelperT<T>::Raw(comparator, Tdec_bits); }
celestar@9876: 
celestar@9876: 	/**
celestar@9876: 	 * Addition for Fixed-point data types
celestar@9876: 	 * @param value The non-floating point number to add
celestar@9876: 	 */
celestar@9876: 	template <typename T> FixedT operator +(const T &value) const
celestar@9876: 	{
celestar@9876: 		return Raw(Raw::m_data + (Tstorage)FixedHelperT<T>::Raw(value, Tdec_bits));
celestar@9876: 	}
celestar@9876: 
celestar@9876: 	/**
celestar@9876: 	 * Subtraction for Fixed-point data types
celestar@9876: 	 * @param value The non-floating point number to subtract
celestar@9876: 	 */
celestar@9876: 	template <typename T> FixedT operator -(const T &value) const
celestar@9876: 	{
celestar@9876: 		return Raw(Raw::m_data - (Tstorage)FixedHelperT<T>::Raw(value, Tdec_bits));
celestar@9876: 	}
celestar@9876: 
celestar@9876: 	/**
celestar@9876: 	 * A simple multiplication for non-floating data types
celestar@9876: 	 * @param value The non floating-point factor
celestar@9876: 	 */
celestar@9876: 	template <typename T> FixedT operator *(const T &value) const
celestar@9876: 	{
celestar@9876: 		return Raw( (((int64)Raw::m_data) * FixedHelperT<T>::Raw(value)) >> FixedHelperT<T>::dec_bits);
celestar@9876: 	}
celestar@9876: 
celestar@9876: 	/**
celestar@9876: 	 * A simple division for non-floating data types
celestar@9876: 	 * @param value The non floating-point divisor
celestar@9876: 	 */
celestar@9876: 	template <typename T> FixedT operator /(const T &value) const
celestar@9876: 	{
celestar@9876: 		return Raw( (((int64)Raw::m_data) << FixedHelperT<T>::dec_bits) / FixedHelperT<T>::Raw(value));
celestar@9876: 	}
celestar@9876: 
celestar@9876: 	/**
celestar@9876: 	 * Addition-assignment for Fixed-point data types
celestar@9876: 	 * @param value The non-floating point number to add
celestar@9876: 	 */
celestar@9876: 	template <typename T> FixedT& operator +=(const T &value)
celestar@9876: 	{
celestar@9876: 		Raw::m_data = Raw::m_data + (Tstorage)FixedHelperT<T>::Raw(value, Tdec_bits);
celestar@9876: 		return *this;
celestar@9876: 	}
celestar@9876: 
celestar@9876: 	/**
celestar@9876: 	 * Subtract-assignment for Fixed-point data types
celestar@9876: 	 * @param value The non-floating point number to subtract
celestar@9876: 	 */
celestar@9876: 	template <typename T> FixedT& operator -=(const T &value)
celestar@9876: 	{
celestar@9876: 		Raw::m_data = Raw::m_data - (Tstorage)FixedHelperT<T>::Raw(value, Tdec_bits);
celestar@9876: 		return *this;
celestar@9876: 	}
celestar@9876: 
celestar@9876: 	/**
celestar@9876: 	 * Multiply-assignment for Fixed-point data types
celestar@9876: 	 * @param value The non-floating point number to multiply
celestar@9876: 	 */
celestar@9876: 	template <typename T> FixedT& operator *=(const T &value)
celestar@9876: 	{
celestar@9882: 		Raw::m_data = ((int64)Raw::m_data * (Tstorage)FixedHelperT<T>::Raw(value)) >> FixedHelperT<T>::dec_bits;
celestar@9876: 		return *this;
celestar@9876: 	}
celestar@9876: 
celestar@9876: 	/**
celestar@9876: 	 * Divide-assignment for Fixed-point data types
celestar@9876: 	 * @param value The non-floating point number to use as divisor
celestar@9876: 	 */
celestar@9876: 	template <typename T> FixedT& operator /=(const T &value)
celestar@9876: 	{
celestar@9876: 		Raw::m_data = (Raw::m_data << FixedHelperT<T>::dec_bits) / (Tstorage)FixedHelperT<T>::Raw(value);
celestar@9876: 		return *this;
celestar@9876: 	}
celestar@9876: 
celestar@9876: 	/**
celestar@9877: 	 * Unary minus operator
celestar@9877: 	 */
celestar@9877: 	FixedT operator -() const { return Raw(-Raw::m_data); }
celestar@9877: 
celestar@9877: 	/**
celestar@9876: 	 * Stream operator, used for floating point output
celestar@9876: 	 * @param os    The stream we are going to write to
celestar@9876: 	 * @param value The Fixed-point variable we want to write in the stream
celestar@9876: 	 */
celestar@9876: 	friend std::ostream& operator << (std::ostream &os, const FixedT &value) { os << (double)value.m_data / (1ULL << Raw::dec_bits); return os; }
celestar@9887: };
celestar@9876: 
celestar@9887: /** The value of \f$\pi\f$ with ample precision for our computations */
celestar@9887: static const FixedT<int64, 16> PI(3141592, 1000000);
celestar@9887: /** The number of elements used in Taylor approximations */
celestar@9887: static const int PRECISION = 5;
celestar@9876: 
celestar@9880: /**
celestar@9880:  * Computes a integral, positive power of a FixedT. Uses an optimized algorithm
celestar@9880:  * to keep computational requirement at bay, by decomposing the power into powers
celestar@9880:  * of two itself.
celestar@9880:  *
celestar@9880:  * @param arg   The number to compute the power for
celestar@9880:  * @param pow   the power
celestar@9887:  * @return      arg^pow
celestar@9880:  * @todo        Add a nice LaTeX forumla to the documentation
celestar@9880:  * @todo        Implement negative powers
celestar@9880:  */
celestar@9880: template <typename Tstorage, int Tdec_bits>
celestar@9880: FixedT<Tstorage, Tdec_bits> pow(const FixedT<Tstorage, Tdec_bits> &arg, int pow)
celestar@9880: {
celestar@9880: 	FixedT<int64, Tdec_bits> temp = 1;
celestar@9880: 	FixedT<int64, Tdec_bits> mul = arg;
celestar@9880: 	if (pow == 0) return (FixedT<Tstorage, Tdec_bits>)1;
celestar@9880: 	for (int i = 1; pow > 0; pow >>= 1, i <<= 1, mul *= mul) if (pow & 1) temp *= mul;
celestar@9880: 
celestar@9880: 	return temp;
celestar@9880: }
celestar@9880: 
celestar@9887: /**
celestar@9887:  * Computes the n-th root of a number using an iterative scheme
celestar@9887:  * @param arg  the number to compute the root from
celestar@9887:  * @param root The "n" in the n-th root
celestar@9887:  * @return     the root
celestar@9887:  * @todo       Abort not after a fixed number of iterations but dynamically by step size
celestar@9887:  * @todo       Make some better first-guess.
celestar@9887:  */
celestar@9887: template <typename Tstorage, int Tdec_bits>
celestar@9887: FixedT<Tstorage, Tdec_bits> nroot(const FixedT<Tstorage, Tdec_bits> &arg, int root)
celestar@9887: {
celestar@9887: 	FixedT<int64, 20> one = 1;
celestar@9887: 	FixedT<int64, 20> guess(6, 5);
celestar@9887: 	for (int i = 0; i < 150; i++) guess = (one / root) * ( guess * (root - 1) + arg / pow(guess, root - 1));
celestar@9887: 
celestar@9887: 	return guess;
celestar@9887: }
celestar@9876: 
celestar@9878: /**
celestar@9878:  * Computes a single element of the Taylor series for the approximation of the cosine.
celestar@9878:  * A single element of the Taylor series computes as:
celestar@9878:  * \f[
celestar@9878:  *    \frac{(-1)^nx^{2n}}{(2n)!}
celestar@9878:  * \f]
celestar@9878:  * In order to prevent overflows, underflows and other nasty stuff, we do not compute
celestar@9878:  * the numerator and denominator separately in the formula but build a product of
celestar@9878:  * fractions:
celestar@9878:  * \f[
celestar@9878:  *    \frac{x}{2n} * \frac{x}{2n-1} * \frac{x}{2n-2} ...
celestar@9878:  * \f]
celestar@9878:  * @param arg The value of "x" in above formula
celestar@9878:  * @param pow The value of "2n" in above formula
celestar@9878:  * @return    The value of the element of the Taylor series
celestar@9878:  */
celestar@9878: template <typename Tstorage, int Tdec_bits>
celestar@9878: FixedT<Tstorage, Tdec_bits> trigbase(const FixedT<Tstorage, Tdec_bits> &arg, int pow)
celestar@9878: {
celestar@9878: 	/* Shortcut for the very first element */
celestar@9878: 	if (pow == 0) return (FixedT<Tstorage, Tdec_bits>)1;
celestar@9878: 
celestar@9878: 	/* find out the (-1)^n part, element should be negative if n is non-even
celestar@9878: 	 * but "pow" already holds 2n, account for that
celestar@9878: 	 */
celestar@9878: 	bool neg = ( (pow / 2) % 2 > 0);
celestar@9878: 
celestar@9878: 	FixedT<Tstorage, Tdec_bits> element = arg / pow;
celestar@9878: 	for (--pow; pow > 0; pow --) element *= arg / pow;
celestar@9878: 
celestar@9878: 	return (neg) ? -element : element;
celestar@9878: }
celestar@9878: 
celestar@9878: /**
celestar@9878:  * Computes the cosine of an argument using a Taylor series.
celestar@9878:  * @param arg The number to compute the cosine from, assumes arg to be in radians.
celestar@9887:  * @return    The cosine
celestar@9878:  * @todo      Optimize the two while loops.
celestar@9878:  * @note      It is possible to factor out the series to optimize even further, but is it
celestar@9878:  *            worth the effort? It could reduce the opcount a little.
celestar@9878:  * @see       trigbase, PRECISION
celestar@9878:  */
celestar@9878: template <typename Tstorage, int Tdec_bits>
celestar@9878: FixedT<Tstorage, Tdec_bits> cos(const FixedT<Tstorage, Tdec_bits> &arg)
celestar@9878: {
celestar@9878: 	FixedT<Tstorage, Tdec_bits> l_arg = arg;
celestar@9878: 	FixedT<Tstorage, Tdec_bits> cosine = 0;
celestar@9878: 
celestar@9878: 	/* We've got to adjust the argument around zero, otherwise we will need
celestar@9878: 	 * (literally) hundreds of elements
celestar@9878: 	 */
celestar@9878: 	while (l_arg < -PI) l_arg += PI * 2;
celestar@9878: 	while (l_arg >  PI) l_arg -= PI * 2;
celestar@9878: 
celestar@9878: 	/* sum up the Taylor elements */
celestar@9878: 	for(int i = 0; i < PRECISION; i++) cosine += trigbase(l_arg, 2 * i);
celestar@9878: 
celestar@9878: 	return cosine;
celestar@9878: }
celestar@9878: 
celestar@9878: /**
celestar@9878:  * Computes the sine of an argument using a Taylor series, uses the cosine
celestar@9878:  * computation in fact, as \f$sin(x) = cos(x - \frac{\pi}{2}))\f$
celestar@9878:  * @param arg The number to compute the cosine from, assumes arg to be in radians.
celestar@9887:  * @return    The sine
celestar@9878:  */
celestar@9878: template <typename Tstorage, int Tdec_bits>
celestar@9878: FixedT<Tstorage, Tdec_bits> sin(const FixedT<Tstorage, Tdec_bits> &arg)
celestar@9878: {
celestar@9878: 	return cos(arg - PI / 2);
celestar@9878: }
celestar@9878: