X-Git-Url: http://git.treefish.org/~alex/phys/su2clebsch.git/blobdiff_plain/e2c2ab8026276b7ebe343b1ebe86f3248547940a..refs/heads/master:/su2clebsch.cpp

diff --git a/su2clebsch.cpp b/su2clebsch.cpp
index 49a8047..831c9da 100644
--- a/su2clebsch.cpp
+++ b/su2clebsch.cpp
@@ -7,35 +7,91 @@
 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;
-    if ( ! (0 <= factor1_n && factor1_n < factor1_S.dimension()) ) {
-      //std::cerr << "clebsch warning: 0 <= factor1_state && factor1_state < factor1_dimension" << std::endl;
-      return 0.0;
-    }
-
     clebsch::weight factor2_S(2);
     factor2_S(1) = factor2_2j;
-    if ( ! (0 <= factor2_n && factor2_n < factor2_S.dimension()) ) {
-      //std::cerr << "clebsch warning: 0 <= factor2_state && factor2_state < factor2_dimension" << std::endl;
-      return 0.0;
-    }
-
     clebsch::weight irrep_S(2);
     irrep_S(1) = irrep_2j;
-    if ( ! (0 <= irrep_n && irrep_n < irrep_S.dimension()) ) {
-      //std::cerr << "clebsch warning: 0 <= irrep_state && irrep_state < irrep_dimension" << std::endl;
-      return 0.0;
-    }
-
+   
     clebsch::decomposition decomp(factor1_S, factor2_S);
 
     bool b = true;
@@ -46,8 +102,6 @@ namespace su2clebsch {
       }
     }
     if (b) {
-      cerr << "WARNING: 2j=" << irrep_2j 
-	   << " does not occour in decomposition!" << endl;
       return 0.0;
     }