// Copyright (c) 2003
// Kevin Atkinson
//
// Permission to use, copy, modify, distribute and sell this software
// and its documentation for any purpose is hereby granted without
// fee, provided that the above copyright notice appear in all copies
// and that both that copyright notice and this permission notice
// appear in supporting documentation.  Kevin Atkinson makes no
// representations about the suitability of this software for any
// purpose.  It is provided "as is" without express or implied
// warranty.

#include <cmath>

template <typename T>
static inline T pow2(T v, int e)
{
  return e < 0 ? v >> -e : v << e;
}

template <int EXP, typename T = int>
class FixedPoint
{
public:
  T val;
  static const int exponent = EXP;
  typedef T BaseType;
  FixedPoint() : val(0) {}

  FixedPoint(T v)
    : val(pow2(v,-EXP)) {}
  FixedPoint(double v0) {
    int e;
    double v = std::frexp(v0,&e);
    val = static_cast<T>(std::ldexp(v,e-EXP));
  }
  FixedPoint(T v, int e)
    : val(pow2(v, e - EXP)) {}
  FixedPoint(const FixedPoint & v)
    : val(v.val) {}
  template <int E2,typename T2> FixedPoint(const FixedPoint<E2,T2> & v)
    : val(pow2(v.val, E2 - EXP)) {}

  operator T () const {return pow2(val,EXP);}
  operator double () const {return std::ldexp(static_cast<double>(val), EXP);}

  FixedPoint operator- () const {return FixedPoint(-val, EXP);}
  static FixedPoint frac(T x, T y)
    {return FixedPoint(pow2(x, -EXP) / y, EXP);}

  // Use these static members to force the return type to a particular
  //   fixed point type
  template<int E1, typename T1, int E2, typename T2>
  static FixedPoint add(FixedPoint<E1,T1> x, FixedPoint<E2,T2> y)
    {return FixedPoint(pow2(x.val, E1-EXP) + pow2(y.val, E2-EXP), EXP);}
  template<int E1, typename T1, int E2, typename T2>
  static FixedPoint sub(FixedPoint<E1,T1> x, FixedPoint<E2,T2> y)
    {return FixedPoint(pow2(x.val, E1-EXP) - pow2(y.val, E2-EXP), EXP);}
  template<int E1, typename T1, int E2, typename T2>
  static FixedPoint mul(FixedPoint<E1,T1> x, FixedPoint<E2,T2> y)
    {return FixedPoint(x.val * y.val, E1 + E2);}
  template<int E1, typename T1, int E2, typename T2>
  static FixedPoint div(FixedPoint<E1,T1> x, FixedPoint<E2,T2> y)
    {return FixedPoint(x.val / y.val, E1 - E2);}
};

//
// Comparision
//

template <int E, typename T>
static inline bool operator== (FixedPoint<E,T> x, FixedPoint<E,T> y)
{
  return x.val == y.val;
}

template <int E, typename T>
static inline bool operator!= (FixedPoint<E,T> x, FixedPoint<E,T> y)
{
  return x.val != y.val;
}

template <int E, typename T>
static inline bool operator< (FixedPoint<E,T> x, FixedPoint<E,T> y)
{
  return x.val < y.val;
}

template <int E, typename T>
static inline bool operator> (FixedPoint<E,T> x, FixedPoint<E,T> y)
{
  return x.val > y.val;
}

template <int E, typename T>
static inline bool operator<= (FixedPoint<E,T> x, FixedPoint<E,T> y)
{
  return x.val <= y.val;
}

template <int E, typename T>
static inline bool operator>= (FixedPoint<E,T> x, FixedPoint<E,T> y)
{
  return x.val >= y.val;
}

//
// Arithmetic
//

template <int E, typename T>
static inline FixedPoint<E,T> operator+ (FixedPoint<E,T> x, FixedPoint<E,T> y)
{
  return FixedPoint<E,T>(x.val + y.val, E);
}

template <int E, typename T>
static inline FixedPoint<E,T> operator- (FixedPoint<E,T> x, FixedPoint<E,T> y)
{
  return FixedPoint<E,T>(x.val - y.val, E);
}

template <int E, typename T>
static inline FixedPoint<E,T> operator* (FixedPoint<E,T> x, FixedPoint<E,T> y)
{
  return FixedPoint<E,T>(x.val * y.val, 2*E);
}

template <int E, typename T>
static inline FixedPoint<E,T> operator/ (FixedPoint<E,T> x, FixedPoint<E,T> y)
{
  return FixedPoint<E,T>(x.val / y.val, -E);
}

template <int E, typename T>
static inline FixedPoint<E,T> operator<< (FixedPoint<E,T> x, int e)
{
  return FixedPoint<E,T>(x.val << e, E);
}

template <int E, typename T>
static inline FixedPoint<E,T> operator>> (FixedPoint<E,T> x, int e)
{
  return FixedPoint<E,T>(x.val >> e, E);
}

//
// Diffrent Type
//

//
// For addition and subtraction the the best exponent for the return
//   type is non-obvious.  So I just used the first one.  Use the
//   static members add and sub for control over the return type
//

template<int E1, int E2, typename T>
static inline FixedPoint<E1,T> operator+ (FixedPoint<E1,T> x, FixedPoint<E2,T> y)
{
  return FixedPoint<E1,T>(x.val + pow2(y.val, E2 - E1), E1);
}

template<int E1, int E2, typename T>
static inline FixedPoint<E1,T> operator- (FixedPoint<E1,T> x, FixedPoint<E2,T> y)
{
  return FixedPoint<E1,T>(x.val - pow2(y.val, E2 - E1), E1);
}

template<int E1, int E2, typename T>
static inline FixedPoint<E1+E2,T> operator* (FixedPoint<E1,T> x, FixedPoint<E2,T> y)
{
  return FixedPoint<E1+E2,T>(x.val * y.val, E1 + E2);
}

template<int E1, int E2, typename T>
static inline FixedPoint<E1-E2,T> operator/ (FixedPoint<E1,T> x, FixedPoint<E2,T> y)
{
  return FixedPoint<E1-E2,T>(x.val / y.val, E1-E2);
}

// Implementing arithmatic when the underlying type is not the same is
// rather tricky.  What should the return type be.  It should be
// typeof(T1 <operation> T2) but typeof is GNU extension.
//
// If you want to avoid the loss of information than it is very difficult
// the answer no only depends on T1 and T2 but also E1 and E2