1

/* $Id */

     2

     3

/**

     4

 * @file

     5

 * Implements a Fixed-point data type.

     6

 * This files defines a new data type

     7

 * "FixedT<Tstorage, Tdec_bits>", which can be used to compute fractional

     8

 * numbers without having to worry about network stability (because

     9

 * there is no real rounding taking place.

    10

 * Use "Fixed" for the defaults values (48 bits integral part, 16 bits

    11

 * fractional part) or "FixedT<typename, int>" for flexible sizes. "int"

    12

 * number of bits are taken from a variable of "typename" for the fractional

    13

 * representation, the rest stays available for the integral part.

    14

 * @warning More than 16 bits fractional should not be used when working

    15

 * with large integers.

    16

 * @warning More than 31 fractional bits are not supported, and will

    17

 * trigger a compiler warning with appropriate settings.

    18

 * @warning There are no warnings for buffer overflows. Those are very likely

    19

 * to happen with divisions, as those first shift the numerator by the number

    20

 * of fractional bits. Use multiplication with the reciprocal whereever

    21

 * it is possible.

    22

 * @warning Please use signed variables for storage at the moment ONLY

    23

 *

    24

 * Variables can be initialized:

    25

 * 	- with integers: "FixedT<int, 9> a = -5"

    26

 * 	- with fractionals: "FixedT<int64, 11> b(7, 2);" where the first number is

    27

 * the numerator and the second is the denominator.

    28

 * 	- not at all: "Fixed c;", which will automatically set c to zero.

    29

 *

    30

 * @example fixed_ex.cpp How to initialize Fixed-point data types

    31

 */

    32

    33

/* Included to overload the stream operator << */

    34

#include <iostream>

    35

/* Needed for some ottd-specific data types */

    36

#include "stdafx.h"

    37

    38

/**

    39

 * Base Class for storing fixed-point data types.

    40

 * Fixed-point data types are stored in a single integer-type variable

    41

 * (Tstorage) of which a certain number of bits (Tdec_bits) is used to

    42

 * represent the fractional part of the number.

    43

 */

    44

template <typename Tstorage, int Tdec_bits>

    45

class FixedRawT {

    46

public:

    47

	/** The number of bits that represent the fractional */

    48

	static const int dec_bits = Tdec_bits;

    49

    50

	/** The storage of the number itself */

    51

	Tstorage m_data;

    52

    53

	/**

    54

	 * Basic constructor for integer arguments.

    55

	 * @param value The raw value to which we want to initialize

    56

	 */

    57

	template <typename T>

    58

	FixedRawT(T value) : m_data(value) {}

    59

    60

	/**

    61

	 * Make sure that a Fixed variable is always inizialized to zero by default.

    62

	 * We need this ctor so that we can declare variable without

    63

	 * the need to initialize it explicitly

    64

	 */

    65

	FixedRawT() : m_data (0) {}

    66

    67

	/**

    68

	 * Returns the biggest integral number we can represent

    69

	 */

    70

	int64 FIXED_MAX() const { return m_maximum; }

    71

    72

	/**

    73

	 * Returns the biggest integral number we can represent

    74

	 */

    75

	int64 FIXED_MIN() const { return m_minimum; }

    76

    77

private:

    78

	/** The largest number we can hold */

    79

	static const int64 m_maximum = +(1ULL << ((sizeof(Tstorage) * 8) - Tdec_bits - 1)) - 1;

    80

    81

	/** The smallest number we can hold */

    82

	static const int64 m_minimum = -(1ULL << ((sizeof(Tstorage) * 8) - Tdec_bits - 1)) - 0;

    83

};

    84

    85

/* forward-declare some structs */

    86

template <typename Tstorage, int Tdec_bits> class FixedT;

    87

template <typename T> struct FixedHelperT;

    88

    89

/**

    90

 * Specialization of FixedHelperT. Used to decimal-align two variables of Fixed type.

    91

 */

    92

template <typename Tstorage, int Tdec_bits> struct FixedHelperT<FixedRawT<Tstorage, Tdec_bits> >

    93

{

    94

	/** The number of bits used for the fraction */

    95

	static const int dec_bits = Tdec_bits;

    96

    97

	/**

    98

	 * Gives the raw data of a FixedRaw

    99

	 * @param t The number to return

   100

	 */

   101

	static int64 Raw(const FixedRawT<Tstorage, Tdec_bits>& t) {return (int64)t.m_data;}

   102

   103

	/**

   104

	 * Returns the data from FixedRawT aligned so that it is aligned to a given number

   105

	 * @param t    The number to be aligned

   106

	 * @param bits The number of bits in the fractional part to be aligned to

   107

	 */

   108

	static int64 Raw(const FixedRawT<Tstorage, Tdec_bits>& t, int bits)

   109

	{

   110

		return (((int64)t.m_data) << (bits > Tdec_bits ? bits - Tdec_bits : 0)) >> (bits < Tdec_bits ? Tdec_bits - bits : 0);

   111

	}

   112

};

   113

   114

/**

   115

 * Specialization of FixedHelperT. Used to decimal-align two variables of Fixed type.

   116

 */

   117

template <typename Tstorage, int Tdec_bits> struct FixedHelperT<FixedT<Tstorage, Tdec_bits> > : public FixedHelperT<FixedRawT<Tstorage, Tdec_bits> > {};

   118

   119

/**

   120

 * General implementation of FixedHelperT. This tempate makes sure that

   121

 * some number supplied is properly aligned at the decimal

   122

 */

   123

template <typename T> struct FixedHelperT

   124

{

   125

	/** This version of FixedHelperT is only used for full integers, so assume the number of dec_bits to be zero */

   126

	static const int dec_bits = 0;

   127

	/** Converts a full integer to one that has a given number of fractional bits, assumes

   128

	 * the number of fractional bits to be zero.

   129

	 * @param t the number to be converted

   130

	 */

   131

	static int64 Raw(const T& t) {return (int64)t;}

   132

   133

	/**

   134

	 * Converts a full integer to one that has a given number of fractional bits

   135

	 * @param bits the number of bits for the fraction

   136

	 * @param t    the number to be converted

   137

	 */

   138

	static int64 Raw(const T& t, int bits) {return ((int64)t) << bits;}

   139

};

   140

   141

   142

/**

   143

 * A class that defines a fixed-point data type, which a variable length and precision.

   144

 * The data type that is defined is a fixed-point variable that has and number of Tdec_bits

   145

 * bits to represent the fractional part, the remaining bits of Tstorage are then used for

   146

 * the integer part of the number. This means, we always need to make a trade-off between

   147

 * the precision we want (higher number in Tdec_bits) and the range (smaller number of

   148

 * Tdec_bits). This class here only defines the constructors and the operators, everything

   149

 * else is happning in helper classes.

   150

 * @note for usage of fprintf and similar, explicit casts are needed (double or int64).

   151

 */

   152

template <typename Tstorage, int Tdec_bits> class FixedT : public FixedRawT<Tstorage, Tdec_bits> {

   153

/* private block up here, because we need the typedef later on */

   154

private:

   155

	/** We shortcut the underlying data type to "Raw", to save typing */

   156

	typedef FixedRawT<Tstorage, Tdec_bits> Raw;

   157

   158

public:

   159

	/**

   160

	 * Ctor for assignment with other, non floating, variable types

   161

	 * @param value The Value we should assign to the Fixed

   162

	 */

   163

	template <typename T> FixedT(T value) : Raw(FixedHelperT<T>::Raw(value, Tdec_bits)) {}

   164

   165

	/** As a constructor without initializing, just use the one that Raw uses. Sets a variable to zero */

   166

	FixedT() : Raw() {}

   167

   168

	/**

   169

	 * Ctor for use with a fraction, useful for initing a variable to a non-integer number at declaration.

   170

	 * @param numerator   The Nominator of the fraction

   171

	 * @param denominator The Denominator of the fraction

   172

	 */

   173

	FixedT(int numerator, int denominator) {Raw::m_data = ((int64)numerator << Raw::dec_bits) / (int64)denominator;}

   174

   175

	/**

   176

	 * Cast to a double (for outputs)

   177

	 */

   178

	operator double() const

   179

	{

   180

		return (double)Raw::m_data / (1 << Raw::dec_bits);

   181

	}

   182

   183

	/**

   184

	 * Equality operator

   185

	 * @param comparator The non-floating point variable we want to compare against

   186

	 */

   187

	template <typename T>

   188

	bool operator ==(const T &comparator) const { return Raw::m_data == (Tstorage)FixedHelperT<T>::Raw(comparator, Tdec_bits); }

   189

   190

	/**

   191

	 * Inequality operator

   192

	 * @param comparator The non-floating point variable we want to compare against

   193

	 */

   194

	template <typename T>

   195

	bool operator !=(const T &comparator) const { return Raw::m_data != (Tstorage)FixedHelperT<T>::Raw(comparator, Tdec_bits); }

   196

   197

	/**

   198

	 * Greater or equal operator

   199

	 * @param comparator The non-floating point variable we want to compare against

   200

	 */

   201

	template <typename T>

   202

	bool operator >=(const T &comparator) const { return Raw::m_data >= (Tstorage)FixedHelperT<T>::Raw(comparator, Tdec_bits); }

   203

   204

	/**

   205

	 * Less or equal operator

   206

	 * @param comparator The non-floating point variable we want to compare against

   207

	 */

   208

	template <typename T>

   209

	bool operator <=(const T &comparator) const { return Raw::m_data <= (Tstorage)FixedHelperT<T>::Raw(comparator, Tdec_bits); }

   210

   211

	/**

   212

	 * Greater than operator

   213

	 * @param comparator The non-floating point variable we want to compare against

   214

	 */

   215

	template <typename T>

   216

	bool operator >(const T &comparator) const { return Raw::m_data > (Tstorage)FixedHelperT<T>::Raw(comparator, Tdec_bits); }

   217

   218

	/**

   219

	 * Less than operator

   220

	 * @param comparator The non-floating point variable we want to compare against

   221

	 */

   222

	template <typename T>

   223

	bool operator <(const T &comparator) const { return Raw::m_data < (Tstorage)FixedHelperT<T>::Raw(comparator, Tdec_bits); }

   224

   225

	/**

   226

	 * Addition for Fixed-point data types

   227

	 * @param value The non-floating point number to add

   228

	 */

   229

	template <typename T> FixedT operator +(const T &value) const

   230

	{

   231

		return Raw(Raw::m_data + (Tstorage)FixedHelperT<T>::Raw(value, Tdec_bits));

   232

	}

   233

   234

	/**

   235

	 * Subtraction for Fixed-point data types

   236

	 * @param value The non-floating point number to subtract

   237

	 */

   238

	template <typename T> FixedT operator -(const T &value) const

   239

	{

   240

		return Raw(Raw::m_data - (Tstorage)FixedHelperT<T>::Raw(value, Tdec_bits));

   241

	}

   242

   243

	/**

   244

	 * A simple multiplication for non-floating data types

   245

	 * @param value The non floating-point factor

   246

	 */

   247

	template <typename T> FixedT operator *(const T &value) const

   248

	{

   249

		return Raw( (((int64)Raw::m_data) * FixedHelperT<T>::Raw(value)) >> FixedHelperT<T>::dec_bits);

   250

	}

   251

   252

	/**

   253

	 * A simple division for non-floating data types

   254

	 * @param value The non floating-point divisor

   255

	 */

   256

	template <typename T> FixedT operator /(const T &value) const

   257

	{

   258

		return Raw( (((int64)Raw::m_data) << FixedHelperT<T>::dec_bits) / FixedHelperT<T>::Raw(value));

   259

	}

   260

   261

	/**

   262

	 * Addition-assignment for Fixed-point data types

   263

	 * @param value The non-floating point number to add

   264

	 */

   265

	template <typename T> FixedT& operator +=(const T &value)

   266

	{

   267

		Raw::m_data = Raw::m_data + (Tstorage)FixedHelperT<T>::Raw(value, Tdec_bits);

   268

		return *this;

   269

	}

   270

   271

	/**

   272

	 * Subtract-assignment for Fixed-point data types

   273

	 * @param value The non-floating point number to subtract

   274

	 */

   275

	template <typename T> FixedT& operator -=(const T &value)

   276

	{

   277

		Raw::m_data = Raw::m_data - (Tstorage)FixedHelperT<T>::Raw(value, Tdec_bits);

   278

		return *this;

   279

	}

   280

   281

	/**

   282

	 * Multiply-assignment for Fixed-point data types

   283

	 * @param value The non-floating point number to multiply

   284

	 */

   285

	template <typename T> FixedT& operator *=(const T &value)

   286

	{

   287

		Raw::m_data = ((int64)Raw::m_data * (Tstorage)FixedHelperT<T>::Raw(value)) >> FixedHelperT<T>::dec_bits;

   288

		return *this;

   289

	}

   290

   291

	/**

   292

	 * Divide-assignment for Fixed-point data types

   293

	 * @param value The non-floating point number to use as divisor

   294

	 */

   295

	template <typename T> FixedT& operator /=(const T &value)

   296

	{

   297

		Raw::m_data = (Raw::m_data << FixedHelperT<T>::dec_bits) / (Tstorage)FixedHelperT<T>::Raw(value);

   298

		return *this;

   299

	}

   300

   301

	/**

   302

	 * Unary minus operator

   303

	 */

   304

	FixedT operator -() const { return Raw(-Raw::m_data); }

   305

   306

	/**

   307

	 * Stream operator, used for floating point output

   308

	 * @param os    The stream we are going to write to

   309

	 * @param value The Fixed-point variable we want to write in the stream

   310

	 */

   311

	friend std::ostream& operator << (std::ostream &os, const FixedT &value) { os << (double)value.m_data / (1ULL << Raw::dec_bits); return os; }

   312

};

   313

   314

/** The value of \f$\pi\f$ with ample precision for our computations */

   315

static const FixedT<int64, 16> PI(3141592, 1000000);

   316

/** The number of elements used in Taylor approximations */

   317

static const int PRECISION = 5;

   318

   319

/**

   320

 * Computes a integral, positive power of a FixedT. Uses an optimized algorithm

   321

 * to keep computational requirement at bay, by decomposing the power into powers

   322

 * of two itself.

   323

 *

   324

 * @param arg   The number to compute the power for

   325

 * @param pow   the power

   326

 * @return      arg^pow

   327

 * @todo        Add a nice LaTeX forumla to the documentation

   328

 * @todo        Implement negative powers

   329

 */

   330

template <typename Tstorage, int Tdec_bits>

   331

FixedT<Tstorage, Tdec_bits> pow(const FixedT<Tstorage, Tdec_bits> &arg, int pow)

   332

{

   333

	FixedT<int64, Tdec_bits> temp = 1;

   334

	FixedT<int64, Tdec_bits> mul = arg;

   335

	if (pow == 0) return (FixedT<Tstorage, Tdec_bits>)1;

   336

	for (int i = 1; pow > 0; pow >>= 1, i <<= 1, mul *= mul) if (pow & 1) temp *= mul;

   337

   338

	return temp;

   339

}

   340

   341

/**

   342

 * Computes the n-th root of a number using an iterative scheme

   343

 * @param arg  the number to compute the root from

   344

 * @param root The "n" in the n-th root

   345

 * @return     the root

   346

 * @todo       Abort not after a fixed number of iterations but dynamically by step size

   347

 * @todo       Make some better first-guess.

   348

 */

   349

template <typename Tstorage, int Tdec_bits>

   350

FixedT<Tstorage, Tdec_bits> nroot(const FixedT<Tstorage, Tdec_bits> &arg, int root)

   351

{

   352

	FixedT<int64, 20> one = 1;

   353

	FixedT<int64, 20> guess(6, 5);

   354

	for (int i = 0; i < 150; i++) guess = (one / root) * ( guess * (root - 1) + arg / pow(guess, root - 1));

   355

   356

	return guess;

   357

}

   358

   359

/**

   360

 * Computes a single element of the Taylor series for the approximation of the cosine.

   361

 * A single element of the Taylor series computes as:

   362

 * \f[

   363

 *    \frac{(-1)^nx^{2n}}{(2n)!}

   364

 * \f]

   365

 * In order to prevent overflows, underflows and other nasty stuff, we do not compute

   366

 * the numerator and denominator separately in the formula but build a product of

   367

 * fractions:

   368

 * \f[

   369

 *    \frac{x}{2n} * \frac{x}{2n-1} * \frac{x}{2n-2} ...

   370

 * \f]

   371

 * @param arg The value of "x" in above formula

   372

 * @param pow The value of "2n" in above formula

   373

 * @return    The value of the element of the Taylor series

   374

 */

   375

template <typename Tstorage, int Tdec_bits>

   376

FixedT<Tstorage, Tdec_bits> trigbase(const FixedT<Tstorage, Tdec_bits> &arg, int pow)

   377

{

   378

	/* Shortcut for the very first element */

   379

	if (pow == 0) return (FixedT<Tstorage, Tdec_bits>)1;

   380

   381

	/* find out the (-1)^n part, element should be negative if n is non-even

   382

	 * but "pow" already holds 2n, account for that

   383

	 */

   384

	bool neg = ( (pow / 2) % 2 > 0);

   385

   386

	FixedT<Tstorage, Tdec_bits> element = arg / pow;

   387

	for (--pow; pow > 0; pow --) element *= arg / pow;

   388

   389

	return (neg) ? -element : element;

   390

}

   391

   392

/**

   393

 * Computes the cosine of an argument using a Taylor series.

   394

 * @param arg The number to compute the cosine from, assumes arg to be in radians.

   395

 * @return    The cosine

   396

 * @todo      Optimize the two while loops.

   397

 * @note      It is possible to factor out the series to optimize even further, but is it

   398

 *            worth the effort? It could reduce the opcount a little.

   399

 * @see       trigbase, PRECISION

   400

 */

   401

template <typename Tstorage, int Tdec_bits>

   402

FixedT<Tstorage, Tdec_bits> cos(const FixedT<Tstorage, Tdec_bits> &arg)

   403

{

   404

	FixedT<Tstorage, Tdec_bits> l_arg = arg;

   405

	FixedT<Tstorage, Tdec_bits> cosine = 0;

   406

   407

	/* We've got to adjust the argument around zero, otherwise we will need

   408

	 * (literally) hundreds of elements

   409

	 */

   410

	while (l_arg < -PI) l_arg += PI * 2;

   411

	while (l_arg >  PI) l_arg -= PI * 2;

   412

   413

	/* sum up the Taylor elements */

   414

	for(int i = 0; i < PRECISION; i++) cosine += trigbase(l_arg, 2 * i);

   415

   416

	return cosine;

   417

}

   418

   419

/**

   420

 * Computes the sine of an argument using a Taylor series, uses the cosine

   421

 * computation in fact, as \f$sin(x) = cos(x - \frac{\pi}{2}))\f$

   422

 * @param arg The number to compute the cosine from, assumes arg to be in radians.

   423

 * @return    The sine

   424

 */

   425

template <typename Tstorage, int Tdec_bits>

   426

FixedT<Tstorage, Tdec_bits> sin(const FixedT<Tstorage, Tdec_bits> &arg)

   427

{

   428

	return cos(arg - PI / 2);

   429

}

   430