#include <iostream>

#include "clebsch.h"

#include "su2clebsch.hpp"

using namespace std;

namespace su2clebsch {
  int ipow(int base, int exp)
  {
    int result = 1;
    while (exp)
      {
        if (exp & 1)
	  result *= base;
        exp >>= 1;
        base *= base;
      }

    return result;
  }

  prestore::~prestore ()
  {
    delete[] cgcdata;
  }

  prestore::prestore(unsigned int _max2j) 
  {
    max2j = _max2j;
    cgcdata = new double[ipow(max2j+1,3)*ipow(2*max2j+1,3)];
    
    for (int factor1_2j=0; factor1_2j<=max2j; factor1_2j++)
      for (int factor1_mp=0; factor1_mp<=2*max2j; factor1_mp++)
	for (int factor2_2j=0; factor2_2j<=max2j; factor2_2j++)
	  for (int factor2_mp=0; factor2_mp<=2*max2j; factor2_mp++)
	    for (int irrep_2j=0; irrep_2j<=max2j; irrep_2j++)
	      for (int irrep_mp=0; irrep_mp<=2*max2j; irrep_mp++)
		{
		  cgcdata[ factor1_2j*ipow(max2j+1,2)*ipow(2*max2j+1,3) + 
			   factor1_mp*ipow(max2j+1,2)*ipow(2*max2j+1,2) +
			   factor2_2j*(max2j+1)*ipow(2*max2j+1,2) + 
			   factor2_mp*(max2j+1)*(2*max2j+1) +
			   irrep_2j*(2*max2j+1) + 
			   irrep_mp ] 
		    =
		    cgc( factor1_2j, factor1_mp-max2j, factor2_2j, factor2_mp-max2j,
			 irrep_2j, irrep_mp-max2j );
		}
  }

  bool stateValid(int state_2j, int state_2m) 
  {
    if ( abs(state_2m) > state_2j || (state_2j + state_2m) % 2 != 0 )
      return false;
    else
      return true;
  }

  double cgc (int factor1_2j, int factor1_2m, int factor2_2j, int factor2_2m,
	      int irrep_2j, int irrep_2m, const prestore& Prestore) 
  {   
    if ( ! ( stateValid(factor1_2j, factor1_2m) && stateValid(factor2_2j, factor2_2m) && stateValid(irrep_2j, irrep_2m) ) )
      return 0.0;
  
    if ( factor1_2j > Prestore.max2j || factor2_2j > Prestore.max2j || irrep_2j > Prestore.max2j )
      return cgc(factor1_2j, factor1_2m, factor2_2j, factor2_2m, irrep_2j, irrep_2m);
    
    return Prestore.cgcdata[ factor1_2j*ipow(Prestore.max2j+1,2)*ipow(2*Prestore.max2j+1,3) + 
			     (factor1_2m+Prestore.max2j) * ipow(Prestore.max2j+1,2)*ipow(2*Prestore.max2j+1,2) +
			     factor2_2j*(Prestore.max2j+1)*ipow(2*Prestore.max2j+1,2) + 
			     (factor2_2m+Prestore.max2j) * (Prestore.max2j+1)*(2*Prestore.max2j+1) +
			     irrep_2j*(2*Prestore.max2j+1) + 
			     (irrep_2m+Prestore.max2j) ];
  } 

  double cgc (int factor1_2j, int factor1_2m, int factor2_2j, int factor2_2m,
	      int irrep_2j, int irrep_2m) 
  {
    if ( ! ( stateValid(factor1_2j, factor1_2m) && stateValid(factor2_2j, factor2_2m) && stateValid(irrep_2j, irrep_2m) ) )
      return 0.0;

    const int factor1_n = (factor1_2j + factor1_2m) / 2;
    const int factor2_n = (factor2_2j + factor2_2m) / 2;
    const int irrep_n =   (irrep_2j + irrep_2m) / 2;

    clebsch::weight factor1_S(2);
    factor1_S(1) = factor1_2j;
    clebsch::weight factor2_S(2);
    factor2_S(1) = factor2_2j;
    clebsch::weight irrep_S(2);
    irrep_S(1) = irrep_2j;
   
    clebsch::decomposition decomp(factor1_S, factor2_S);

    bool b = true;
    for (int i = 0; i < decomp.size(); ++i) {
      if (decomp(i) == irrep_S) {
	b = false;
	break;
      }
    }
    if (b) {
      return 0.0;
    }

    const clebsch::coefficients C(irrep_S, factor1_S, factor2_S);

    return C( factor1_n, factor2_n, 0, irrep_n );
  }

};