clebschgordan.cpp

   1 #include <iostream>
   2
   3 #include "clebsch.h"
   4
   5 using namespace std;
   6
   7 int main() {
   8     while (true) {
   9         int choice, N;
  10
  11         cout << "What would you like to do?" << endl;
  12         cout << "1) Translate an i-weight S to its index P(S)" << endl;
  13         cout << "2) Recover an i-weight S from its index P(S)" << endl;
  14         cout << "3) Translate a pattern M to its index Q(M)" << endl;
  15         cout << "4) Recover a pattern M from its index Q(M)" << endl;
  16         cout << "5) Calculate Clebsch-Gordan coefficients for S x S' -> S''" << endl;
  17         cout << "6) Calculate all Glebsch-Gordan coefficients for S x S'" << endl;
  18         cout << "0) Quit" << endl;
  19
  20         do {
  21             cin >> choice;
  22         } while (choice < 0 || choice > 6);
  23
  24         if (choice == 0) {
  25             break;
  26         }
  27
  28         cout << "N (e.g. 3): ";
  29         cin >> N;
  30
  31         switch (choice) {
  32             case 1: {
  33                 clebsch::weight S(N);
  34                 cout << "Irrep S: ";
  35                 for (int k = 1; k <= N; ++k) {
  36                     cin >> S(k);
  37                 }
  38                 cout << S.index() << endl;
  39                 break;
  40             }
  41             case 2: {
  42                 int P;
  43                 cout << "Index: ";
  44                 cin >> P;
  45                 clebsch::weight S(N, P);
  46                 cout << "I-weight:";
  47                 for (int k = 1; k <= N; ++k) {
  48                     cout << ' ' << S(k);
  49                 }
  50                 cout << endl;
  51                 break;
  52             }
  53             case 3: {
  54                 clebsch::pattern M(N);
  55                 for (int l = N; l >= 1; --l) {
  56                     cout << "Row l = " << l << ": ";
  57                     for (int k = 1; k <= l; ++k) {
  58                         cin >> M(k, l);
  59                     }
  60                 }
  61                 cout << "Index: " << M.index() + 1 << endl;
  62                 break;
  63             }
  64             case 4: {
  65                 clebsch::weight S(N);
  66                 cout << "Irrep S: ";
  67                 for (int i = 1; i <= N; ++i) {
  68                     cin >> S(i);
  69                 }
  70
  71                 int Q;
  72                 cout << "Index (1..dim(S)): ";
  73                 cin >> Q;
  74                 clebsch::pattern M(S, Q - 1);
  75                 for (int l = N; l >= 1; --l) {
  76                     for (int k = 1; k <= l; ++k) {
  77                         cout << M(k, l) << '\t';
  78                     }
  79                     cout << endl;
  80                 }
  81                 break;
  82             }
  83             case 5: {
  84                 clebsch::weight S(N);
  85                 cout << "Irrep S (e.g.";
  86                 for (int k = N - 1; k >= 0; --k) {
  87                     cout << ' ' << k;
  88                 }
  89                 cout << "): ";
  90                 for (int k = 1; k <= N; ++k) {
  91                     cin >> S(k);
  92                 }
  93
  94                 clebsch::weight Sprime(N);
  95                 cout << "Irrep S' (e.g.";
  96                 for (int k = N - 1; k >= 0; --k) {
  97                     cout << ' ' << k;
  98                 }
  99                 cout << "): ";
 100                 for (int k = 1; k <= N; ++k) {
 101                     cin >> Sprime(k);
 102                 }
 103
 104                 clebsch::decomposition decomp(S, Sprime);
 105                 cout << "Littlewood-Richardson decomposition S \\otimes S' = \\oplus S'':" << endl;
 106                 cout << "[irrep index] S'' (outer multiplicity) {dimension d_S}" << endl;
 107                 for (int i = 0; i < decomp.size(); ++i) {
 108                     cout << "[" << decomp(i).index() << "] ";
 109                     for (int k = 1; k <= N; ++k) {
 110                         cout << decomp(i)(k) << ' ';
 111                     }
 112                     cout << '(' << decomp.multiplicity(decomp(i)) << ") {"
 113                          << decomp(i).dimension() << "}" << endl;;
 114                 }
 115
 116                 clebsch::weight Sdoubleprime(N);
 117                 for (bool b = true; b; ) {
 118                     cout << "Irrep S'': ";
 119                     for (int k = 1; k <= N; ++k) {
 120                         cin >> Sdoubleprime(k);
 121                     }
 122                     for (int i = 0; i < decomp.size(); ++i) {
 123                         if (decomp(i) == Sdoubleprime) {
 124                             b = false;
 125                             break;
 126                         }
 127                     }
 128                     if (b) {
 129                         cout << "S'' does not occur in the decomposition" << endl;
 130                     }
 131                 }
 132
 133                 int alpha;
 134                 while (true) {
 135                     cout << "Outer multiplicity index: ";
 136                     cin >> alpha;
 137                     if (1 <= alpha && alpha <= decomp.multiplicity(Sdoubleprime)) {
 138                         break;
 139                     }
 140                     cout << "S'' does not have this multiplicity" << endl;
 141                 }
 142
 143                 string file_name;
 144                 cout << "Enter file name to write to file (leave blank for screen output): ";
 145                 cin.ignore(1234, '\n');
 146                 getline(cin, file_name);
 147
 148                 const clebsch::coefficients C(Sdoubleprime, S, Sprime);
 149                 int dimS = S.dimension(),
 150                     dimSprime = Sprime.dimension(),
 151                     dimSdoubleprime = Sdoubleprime.dimension();
 152
 153                 ofstream os(file_name.c_str());
 154                 (file_name.empty() ? cout : os).setf(ios::fixed);
 155                 (file_name.empty() ? cout : os).precision(15);
 156                 (file_name.empty() ? cout : os) << "List of nonzero CGCs for S x S' => S'', alpha" << endl;
 157                 (file_name.empty() ? cout : os) << "Q(M)\tQ(M')\tQ(M'')\tCGC" << endl;
 158                 for (int i = 0; i < dimSdoubleprime; ++i) {
 159                     for (int j = 0; j < dimS; ++j) {
 160                         for (int k = 0; k < dimSprime; ++k) {
 161                             double x = double(C(j, k, alpha - 1, i));
 162
 163                             if (fabs(x) > clebsch::EPS) {
 164                                 (file_name.empty() ? cout : os) << j + 1 << '\t'
 165                                     << k + 1 << '\t' << i + 1 << '\t' << x << endl;
 166                             }
 167                         }
 168                     }
 169                 }
 170
 171                 break;
 172             }
 173             case 6: {
 174                 clebsch::weight S(N);
 175                 cout << "Irrep S (e.g.";
 176                 for (int k = N - 1; k >= 0; --k) {
 177                     cout << ' ' << k;
 178                 }
 179                 cout << "): ";
 180                 for (int k = 1; k <= N; ++k) {
 181                     cin >> S(k);
 182                 }
 183
 184                 clebsch::weight Sprime(N);
 185                 cout << "Irrep S' (e.g.";
 186                 for (int k = N - 1; k >= 0; --k) {
 187                     cout << ' ' << k;
 188                 }
 189                 cout << "): ";
 190                 for (int k = 1; k <= N; ++k) {
 191                     cin >> Sprime(k);
 192                 }
 193
 194                 string file_name;
 195                 cout << "Enter file name to write to file (leave blank for screen output): ";
 196                 cin.ignore(1234, '\n');
 197                 getline(cin, file_name);
 198
 199                 ofstream os(file_name.c_str());
 200                 (file_name.empty() ? cout : os).setf(ios::fixed);
 201                 (file_name.empty() ? cout : os).precision(15);
 202
 203                 clebsch::decomposition decomp(S, Sprime);
 204                 (file_name.empty() ? cout : os) <<
 205                     "Littlewood-Richardson decomposition S \\otimes S' = \\oplus S'':" << endl;
 206                 (file_name.empty() ? cout : os) <<
 207                     "[irrep index] S'' (outer multiplicity) {dimension d_S}" << endl;
 208                 for (int i = 0; i < decomp.size(); ++i) {
 209                     (file_name.empty() ? cout : os) << "[" << decomp(i).index() << "] ";
 210                     for (int k = 1; k <= N; ++k) {
 211                         (file_name.empty() ? cout : os) << decomp(i)(k) << ' ';
 212                     }
 213                     (file_name.empty() ? cout : os) << '(' << decomp.multiplicity(decomp(i)) << ") {"
 214                          << decomp(i).dimension() << "}" << endl;;
 215                 }
 216
 217                 for (int i = 0; i < decomp.size(); ++i) {
 218                     const clebsch::coefficients C(decomp(i),S, Sprime);
 219                     int dimS = S.dimension(),
 220                         dimSprime = Sprime.dimension(),
 221                         dimSdoubleprime = decomp(i).dimension();
 222
 223                     for (int m = 0; m < C.multiplicity; ++m) {
 224                         (file_name.empty() ? cout : os) << "List of nonzero CGCs for S x S' => S'' = (";
 225                         for (int j = 1; j <= N; ++j) cout << decomp(i)(j) << (j < N ? ' ' : ')');
 226                         (file_name.empty() ? cout : os) << ", alpha = " << m + 1 << endl;
 227                         (file_name.empty() ? cout : os) << "Q(M)\tQ(M')\tQ(M'')\tCGC" << endl;
 228                         for (int i = 0; i < dimSdoubleprime; ++i) {
 229                             for (int j = 0; j < dimS; ++j) {
 230                                 for (int k = 0; k < dimSprime; ++k) {
 231                                     double x = double(C(j, k, m, i));
 232
 233                                     if (fabs(x) > clebsch::EPS) {
 234                                         (file_name.empty() ? cout : os) << j  + 1<< '\t'
 235                                             << k + 1 << '\t' << i + 1 << '\t' << x << endl;
 236                                     }
 237                                 }
 238                             }
 239                         }
 240
 241                         (file_name.empty() ? cout : os) << endl;
 242                     }
 243                 }
 244
 245                 break;
 246             }
 247         }
 248     }
 249
 250     return 0;
 251 }