#ifndef __JMATHTOOLKIT__
#define __JMATHTOOLKIT__

#include <limits>
#include <cmath>

#include "JLang/JException.hh"
#include "JMath/JZero.hh"


/**
 * \file
 * Auxiliary methods for geometrical methods.
 * In this, the action of each global method is transferred to 
 * the corresponding member method of the leading (first) argument, e.g:
 *
 * <tt>getXXX(first, second, ...)</tt>
 *
 * becomes
 *
 * <tt>first.getXXX(second, ...)</tt>
 *\author mdejong
 */

/**
 * Auxiliary classes and methods for mathematical operations.
 */
namespace JMATH {}
namespace JPP { using namespace JMATH; }

namespace JMATH {

  using JLANG::JValueOutOfRange;


  /**
   * Determine factorial.
   *
   * \param  n        number
   * \return          factorial (= n!)
   */
  inline long long int factorial(const long long int n)
  {
    if (n < 0) {
      THROW(JValueOutOfRange, "JMATH::factorial(): invalid argument " << n);
    }

    if (n != 1)
      return n * factorial(n - 1);
    else
      return 1;
  }


  /**
   * Determine combinatorics.
   *
   * \param  n        number
   * \param  m        number
   * \return          n!/(m!*(n-m)!)
   */
  inline long long int factorial(const long long int n, const long long int m)
  {
    if (n < 0 || m < 0 || n < m) {
      THROW(JValueOutOfRange, "JMATH::factorial(): invalid argument " << n << ' ' << m);
    }

    if      (n == m)
      return 1;
    else if (m == 0)
      return 1;
    else
      return factorial(n - 1, m - 1) + factorial(n - 1, m);
  }


  /**
   * Check equality.
   *
   * \param  first          first  object
   * \param  second         second object
   * \param  precision      precision
   * \return                true if two objects are equals; else false
   */
  template<class JFirst_t, class JSecond_t>
  inline bool equals(const JFirst_t&  first,
		     const JSecond_t& second,
		     const double     precision = std::numeric_limits<double>::min())
  {
    return first.equals(second, precision);
  }


  /**
   * Get square of distance between objects.
   *
   * \param  first          first  object
   * \param  second         second object
   * \return                square of distance
   */
  template<class JFirst_t, class JSecond_t>
  inline double getDistanceSquared(const JFirst_t&  first,
				   const JSecond_t& second)
  {
    return first.getDistanceSquared(second);    
  }


  /**
   * Get distance between objects.
   *
   * \param  first          first  object
   * \param  second         second object
   * \return                distance
   */
  template<class JFirst_t, class JSecond_t>
  inline double getDistance(const JFirst_t&  first,
			    const JSecond_t& second)
  {
    return first.getDistance(second);
  }


  /**
   * Get dot product of objects.
   *
   * \param  first          first  object
   * \param  second         second object
   * \return                dot product
   */
  template<class JFirst_t, class JSecond_t>
  inline double getDot(const JFirst_t&  first,
		       const JSecond_t& second)
  {
    return first.getDot(second);
  }


  /**
   * Get space angle between objects.
   *
   * \param  first          first  object
   * \param  second         second object
   * \return                angle [deg]
   */
  template<class JFirst_t, class JSecond_t>
  inline double getAngle(const JFirst_t&  first,
			 const JSecond_t& second)
  {
    const double dot = getDot(first,second);

    if      (dot >= +1.0)
      return   0.0;
    else if (dot <= -1.0)
      return 180.0;
    else
      return acos(dot) * 180.0 / acos(-1.0);
  }


  /**
   * Get perpendicular dot product of objects.
   *
   * \param  first          first  object
   * \param  second         second object
   * \return                perpendicular dot product
   */
  template<class JFirst_t, class JSecond_t>
  inline double getPerpDot(const JFirst_t&  first,
			   const JSecond_t& second)
  {
    return first.getPerpDot(second);
  }


  /**
   * Get cross product of objects.
   *
   * \param  first          first  object
   * \param  second         second object
   * \return                cross product
   */
  template<class T>
  inline T getCross(const T& first,
		    const T& second)
  {
    return T().getCross(first, second);
  }


  /**
   * Convolute data with given kernel.
   *
   * \param   input         input data
   * \param   kernel        convolution kernel
   * \return                convoluted data
   */
  template<class T>
  inline std::vector<T> convolve(const std::vector<T>& input, const std::vector<T>& kernel)
  {
    const size_t k = kernel.size();
    const size_t n = input .size() - k + 1;
    
    std::vector<T> out(n, getZero<T>());
     
    for (size_t i = 0; i != n; ++i) {
      for (size_t j = 0; j != k; ++j) {
	out[i] += input[i + j] * kernel[k - j - 1]; 
      }
    }

    return out;
  }
}

#endif