#include <iostream>
#include <complex>
#include <boost/foreach.hpp>
using namespace std;
// use rational numbers
#include <boost/rational.hpp>
// use math constants
#include <boost/math/constants/constants.hpp>
// use math factorial
#include <boost/math/special_functions/factorials.hpp>
//using namespace boost::math
// // use units
// #include <boost/typeof/std/complex.hpp>
// #include <boost/units/cmath.hpp>
// #include <boost/units/systems/si.hpp>
// #include <boost/units/systems/si/prefixes.hpp>
// #include <boost/units/systems/si/codata/universal_constants.hpp>
// 
// using namespace boost::units;
// using namespace boost::units::si;
// using namespace boost::units::si::constants;
// //use units end

// use unlimited precision
//#include<gmpxx.h> // the mpz_class gives arbitrary precision. I need working factorials of 50, or even better 100!
//typedef mpz_class       REAL;


typedef double REAL;

vector<boost::rational<signed long>> hermite_polynomial_coefficients(int order, boost::rational<signed long> lambda_per_alpha)
{
    vector<boost::rational<signed long>> c_even; // can't use static, because lambda_per_alpha changes between calls
    vector<boost::rational<signed long>> c_odd;
    if(c_even.size() == 0)
    {
        c_even.push_back(1);
        c_even.push_back(0);
    }
    if(c_odd.size() == 0)
    {
        c_odd.push_back(0);
        c_odd.push_back(1);
    }
    vector<boost::rational<signed long>>& a( (order % 2 == 0) ?c_even:c_odd );
    
    if(order+1 < a.size())
        return a;

    for(int i = a.size() ; i <= order ; ++i)
    {
        boost::rational<signed long> next = a[i-2]*(2*(i-2)+1-lambda_per_alpha)/(i*(i-1));
        a.push_back(next);
    }
    
//    BOOST_FOREACH(boost::rational<signed long> v,a)
//        std::cout << v << " ";
//    std::cout << "\n";
    
    return a;
}

vector<boost::rational<signed long>> hermite_polynomial_scaled(int order, boost::rational<signed long> lambda_per_alpha)
{
    vector<boost::rational<signed long>> a(hermite_polynomial_coefficients(order,lambda_per_alpha));
    vector<boost::rational<signed long>> b;
    boost::rational<signed long> factor(std::pow(2,order)/a[order]);
    for(int i = 0 ; i<=order ; ++i)
        b.push_back(a[i]*factor);
    return b;
}

// I want some results, this is not compiling, so I rewrite this without units.
//
// /* this function returns 1/sqrt(meter) */ quantum_wavefunction_1D quantum_oscillator_wavefunction(
//       int n                             // n - order of wavefunction
//     , const quantity<length> x          // position
//     , const quantity<mass> mass         // mass - mass of particle
//     , const quantity<frequency>  omega  // oscillator frequency
// )
// {
//     typedef derived_dimension<length_base_dimension,-2>::type  reciprocal_area;
// //    typedef derived_dimension<length_base_dimension,root<2> >::type  quantum_wavefunction_1D;
// 
//     quantity<energy>            E      ( hbar*omega*(n+0.5)       );        // energy of oscillator
//     quantity<reciprocal_area>   lambda ( 2.0*mass*E / (hbar*hbar) );        // 
//     quantity<reciprocal_area>   alpha  ( mass*omega / hbar        );
//     auto                        alpha_root(root<2>(alpha));
//     REAL                      lambda_per_alpha ( 2*n+1 );
// 
//     // N_n should be 1/sqrt(meter) because it has a dimension of quantum_wavefunction_1D
//     auto                        N_n(root<2>( root<2>(alpha/constants::pi<REAL>() ) /(std::pow(2,n)*factorial(n) ) ));
// 
//     complex<REAL> i(0,1);
//     complex<REAL> result(0);
//     vector<boost::rational<signed long>> a(hermite_polynomial_scaled(n,lambda_per_alpha))
//     for(int v=0 ; v<=n ; v++)
//         result+=N_n*a[v]*pow<v>(alpha_root*x)*exp(-0.5*alpha*pow<2>(x));
//     retirn result;
// }
// 

REAL quantum_oscillator_wavefunction( // assume hbar=1, mass=1, frequency=1
      int n                             // n - order of wavefunction
    , REAL x          // position
    , REAL mass=1.0   // mass - mass of particle
    , REAL omega=1.0  // oscillator frequency
)
{
    REAL hbar=1.0;
    REAL                      E      ( hbar*omega*(n+0.5)       );        // energy of oscillator
    REAL                      lambda ( 2.0*mass*E / (hbar*hbar) );        // 
    REAL                      alpha  ( mass*omega / hbar        );
    REAL                      alpha_root(pow(alpha,0.5));
    //REAL                      lambda_per_alpha ( 2*n+1 );

    // N_n should be 1/sqrt(meter) because it has a dimension of quantum_wavefunction_1D
    REAL                      N_n(  pow( pow(alpha/boost::math::constants::pi<REAL>() , 0.5) /(std::pow(2,n)*boost::math::factorial<REAL>(n) ) , 0.5)  );

    complex<REAL> i(0,1);
    complex<REAL> result(0);
    vector<boost::rational<signed long>> a(hermite_polynomial_scaled(n,  2*n+1 ));
    for(int v=0 ; v<=n ; v++)
    {
        REAL a_v(boost::rational_cast<REAL>(a[v])); // Hermite Polynomial coefficient
        result+=N_n*a_v*pow(alpha_root*x , v) *exp(-0.5*alpha*pow(x,2));
    }
    return abs(result);
}

int main()
{
    // this code is checking if Hermite Polynomials are OK
    // it's correct so comment it out.
    //for(int n=0 ; n<12 ; n++)
    //{
    //    int pot(n);
    //    std::cout << n << "  :  ";
    //    vector<boost::rational<signed long>> a(hermite_polynomial_scaled(n,2*n+1));
    //    BOOST_REVERSE_FOREACH(boost::rational<signed long> v,a)
    //    {
    //        if(v!=0)
    //            std::cout << boost::rational_cast<signed long long>(v) << "*x^"<< pot << " ";
    //        --pot;
    //    }
    //    std::cout << "\n";
    //}

//    quantity<length> s(-1*meter);
//    quantity<length> ds(0.01*milli*meter);
//    for(quantity<length> x=-s ; x<=s ; x+=ds)
//    {
//        cout << x << " " << quantum_oscillator_wavefunction(...) << "\n";
//    }

// without units, for now.
    REAL s=15.0;
    REAL ds=0.01;
    for(REAL x=-s ; x<=s ; x+=ds)
    {
        cout << x << " " << quantum_oscillator_wavefunction(25,x) << "\n";
    }


    return true;
}