/**
 * @file   helper.hpp
 * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
 * @date   Tue Apr  5 21:03:47 2011
 * 
 * @brief  Provides various helper functions.
 * 
 * 
 */

#ifndef HELPER_HPP_
#define HELPER_HPP_

#include <cassert>
#include <cstdio>
#include <limits>

namespace VMG
{

namespace Helper
{

  /**
   * Checks a number for validity, i.e. it is neither nan nor inf.
   * 
   */
  inline bool CheckNumber(const vmg_float& number)
  {
    bool valid = true;

    if (std::numeric_limits<vmg_float>::has_quiet_NaN) {
      valid &= number != std::numeric_limits<vmg_float>::quiet_NaN();
      assert(number != std::numeric_limits<vmg_float>::quiet_NaN());
    }

    if (std::numeric_limits<vmg_float>::has_signaling_NaN) {
      valid &= number != std::numeric_limits<vmg_float>::signaling_NaN();
      assert(number != std::numeric_limits<vmg_float>::signaling_NaN());
    }

    if (std::numeric_limits<vmg_float>::has_infinity) {
      valid &= number != std::numeric_limits<vmg_float>::infinity();
      valid &= number != -1 * std::numeric_limits<vmg_float>::infinity();
      assert(number != std::numeric_limits<vmg_float>::infinity());
      assert(number != -1 * std::numeric_limits<vmg_float>::infinity());
    }

    return valid;
  }

  inline int intpow(unsigned int base, unsigned int power)
  {
    assert (power >= 0);
    unsigned int result = 1;
    while (power != 0) {
      if (power & 1)
	result *= base;
      base *= base;
      power >>= 1;
    }
    return result;
  }

  /**
   * Tests two arbitrary objects for equality and prints
   * a warning if they differ.
   */
  template <class T>
  bool IsEq(const T& val1, const T& val2, const char name[])
  {
    bool rval = (val1 == val2);

    if (!rval)
      printf("WARNING: Values are not equal (%s)\n", name);

    assert(rval);

    return rval;
  }

  inline void MM(int nn, vmg_float u, vmg_float *tbl)
  {
    int i, n, m;

    assert(u >= 0.0 && u <= 1.0);

    /* n = 1: */
    /* i = 0: */
    tbl[0] = 1.0;

    for (n = 2; n <= nn; n++) {
      m = n - 1;

      /* i = n - 1: */
      tbl[n-1] =  (1 - u) * tbl[n-2] / m;

      for (i = n - 2; i >= 1; i--)
	tbl[i] = ((u + i) * tbl[i] + (n - (u + i)) * tbl[i-1]) / m;

      /* i = 0: */
      tbl[0] *= u / m;
    }
  }

  inline void MdM(int nn, vmg_float u,
		  vmg_float * __restrict__ tbl, vmg_float * __restrict__ dtbl)
  {
    int i, n, m;

    assert(nn >= 2);
    assert(u >= 0.0 && u <= 1.0);

    /* n = 1: */
    /* i = 0: */
    tbl[0] = 1.0;

    for (n = 2; n <= nn; n++) {
      m = n - 1;

      /* differentiation */
      if (n == nn) {
	dtbl[n-1] = - tbl[n-2];

	for (i = n - 2; i >= 1; i--)
	  dtbl[i] = tbl[i] - tbl[i-1];

	dtbl[0] = tbl[0];
      }

      /* i = n - 1: */
      tbl[n-1] =  (1 - u) * tbl[n-2] / m;

      for (i = n - 2; i >= 1; i--)
	tbl[i] = ((u + i) * tbl[i] + (n - (u + i)) * tbl[i-1]) / m;

      /* i = 0: */
      tbl[0] *= u / m;

    }
  }

}

}

#endif /* HELPER_HPP_ */