Quantenmechanik/ Drehimpuls-Kopplung/ CG-Skript

Skript zu Clebsch-Gordan-Koeffizienten

Das Skript listet Clebsch-Gordan-Koeffizienten zu Paaren von Drehimpulsen L und S. Sie sollen zum Gesamt-Drehimpuls mit dem Operator J=L+S kombiniert werden. Nur die nichtverschwindenden Werte für die möglichen Eigenvektoren $|J,J_{3}\rangle$ des Gesamt-Drehimpulses werden aufgezählt.

Alle Rechnerei ist Ganzzahl-Arithmetik, mit der "Rationale Wurzelzahlen" aufgebaut werden. Das sind hier Quintupel von Ganzzahlen [r,p,q,u,v] mit q>0, v>0 und der Bedeutung: ${\frac {(r+ip)}{q}}\cdot {\sqrt {\frac {u}{v}}}$.
Der Imaginärteil p wird in diesem Skript nicht benutzt, er ist immer Null. Die Funktion reformat() verwandelt ganze und einfach-rationale, d.h. Paare [p,q], in solche Quintupel. Rechenschritte erfolgen mit den Funktionen
multip,addit,sqroot,recip,scale,lincomb,quotient,iszero.
Die Addition geht nur, wenn beide Argument rationale Vielfache derselben Wurzel sind, sonst bricht das Programm ab.

Zu Eingabe-Paaren (L,S) berechnet das Skript alle Clebsch-Gordan-Koeffizienten. L und S sind Vielfache von (1/2). Also (2L)=1,2,3,... und (2S)=1,2,3,... Eine Clebsch-Gordan-Matrix enthält (2L+1)*(2S+1) zum Quadrat rationale Wurzelzahlen.
De Matrix ist hier eine Tabelle mit 4 Indizes $cgm[i][j][k][l]$.
Quantenzahl J ist die absteigende Folge (L+S),(L+S-1)...|L-S|.
J,L und S sind entweder halb- oder ganzzahlig.
i,j,k,l fangen bei Null an, meinen aber folgende physikalischen Werte:

Index i adressiert den Wert J, vom Maximum (L+S) abwärts.

Index j adressiert den Wert J_3, in absteigender Reihenfolge J,J-1,...,-J.

Index k adressiert den Wert L_3, Reihenfolge abfallend L,L-1,...,-L

Index l adressiert S_3, Reihenfolge S,S-1,...,-S.

In physikalischer Notation ist $cgm[i][j][k][l]=\langle L,S;L_{3},S_{3}|J,J_{3}\rangle$.
Die Elemente mit festem $(J,J_{3})$ und $(L_{3},S_{3})$ variabel nennen wir etwas willkürlich eine Zeile der Matrix, bei umgekehrter Variation eine Spalte. Für L=3,S=2 hat eine Zeile die Länge (2*3+1)(2*2+1)=35 und eine Spalte (2*1+1)+(2*2+1)+(2*3+1)+(2*4+1)+(2*5+1)=3+5+7+9+11=35, ganz wie geplant.

Der Algorithmus folgt sklavisch dem Standard-Rezept. Angefangen wird mit dem maximalen J,J_3 und einem Koeffizienten 1. Mit den Absteige-Operatoren werden die Zeilenvektoren für die Folge der J_3 errechnet. Beim nächst niedrigeren J wird angefangen mit dem orthogonalen Komplement zu den bis dahin bereits vorhandenen Vektoren zur Quantenzahl J_3=J. Der Orthogonalisierungs-Algorithmus bekommt n-tupel von Wurzelzahlen als Eingabe zu sehen, jeweils (n-1) davon. Ist $\{x_{i}\}$ eine Menge von schon orthogonalen Vektoren und y ein Kandidat für ihr Komplement, dann werden einfach alle seine Projektionen abgezogen: $y'=y-\sum _{i}{\frac {\langle y|x_{i}\rangle }{\langle x_{i}|x_{i}\rangle }}\cdot x_{i}$. Für $y$ werden brutal die Standard-Einheitsvektoren durchprobiert bis ein orthonormierbares $y'$ da ist. Es gibt mindestens eine der (n) Achsen, die von den angelieferten (n-1) Vektoren linear unabhängig ist, daher wird ein Komplement-Vektor gefunden.
Konvention: Der Eintrag zum maximalen Wert von J_3 und L_3 bei jeder neuen Treppenstufe von J wird positiv gemacht.

Halbzahlige Variablen werden vermieden zugunsten ihres Doppels. Wenn ein Variablenname wie j2 oder m2 mit der Ziffer 2 endet, enthält er zweimal die physikalische Größe j oder m.

Das Skript wurde durchprobiert mit den Quantenzahlen L und S = (1/2),1,(3/2),2,(5/2),3. Keine Garantie darüber hinaus. Mit L=S=3 hat die CG-Matrix 2401=(7*7)*(7*7) Elemente, davon 216 nicht Null!

Der Code hat (langsame) Testsequenzen dazu, dass alle Zeilen und Spalten der Clebsch-Gordan-Matrix orthonormal sind, dass die Symmetrie beim Umklappen von $\{l_{3},s_{3},j_{3}\}$ stimmt und dass die Kopfzeile jedes j-Blocks auch mit den Racah-Formeln genau gleich herauskommt.

Definierende Gleichung der CG-Matrix:

$|j,j_{3}>=\sum _{l_{3}=-l}^{l}\sum _{s_{3}=-s}^{s}[|l,l_{3}\rangle \otimes |s,s_{3}\rangle ]\cdot \langle l,s;l_{3},s_{3}|j,j_{3}\rangle$

Symmetrien der Koeffizienten:

• $\langle l,s;l_{3},s_{3}|j,j_{3}\rangle =(-1)^{l+s-j}\langle s,l;s_{3},l_{3}|j,j_{3}\rangle$
• $\langle l,s;-l_{3},-s_{3}|j,-j_{3}\rangle =(-1)^{l+s-j}\langle l,s;l_{3},s_{3}|j,j_{3}\rangle$

Auf- und Absteigeregeln für festes J:

$A(j,j_{3})\langle l,s;l_{3},s_{3}|j,j_{3}\rangle =A(l,l_{3})\langle l,s;l_{3}+1,s_{3}|j,j_{3}+1\rangle +A(s,s_{3})\langle l,s;l_{3},s_{3}+1|j,j_{3}+1\rangle$

$B(j,j_{3})\langle l,s;l_{3},s_{3}|j,j_{3}\rangle =B(l,l_{3})\langle l,s;l_{3}-1,s_{3}|j,j_{3}-1\rangle +B(s,s_{3})\langle l,s;l_{3},s_{3}-1|j,j_{3}-1\rangle$

$A(j,m):={\sqrt {j(j+1)-m(m+1)}};\;B(j,m):={\sqrt {j(j+1)-m(m-1)}}$

Beziehung zwischen J-Nachbarn bei festem $j_{3}=m$ :

• $A\cdot \langle l,s;l_{3},s_{3}|j,m\rangle =B\cdot \langle l,s;l_{3},s_{3}|(j+1),m\rangle +C\cdot \langle l,s;l_{3},s_{3}|(j-1),m\rangle$

$A:=(l_{3}-s_{3})+(l_{3}+s_{3}){\frac {[s(s+1)-l(l+1)]}{[j(j+1)]}};\;B:={\sqrt {f(j+1)}};\;C:={\sqrt {f(j)}}$

$f(j):=(j^{2}-m^{2}){\frac {[(l+s+1)^{2}-j^{2}][j^{2}-(l-s)^{2}]}{4j^{2}(2j-1)(2j+1)}}$

Racah-Formel für den Spezialfall $j_{3}=j;\;l_{3}+s_{3}=j$

• $\langle l,s;l_{3},s_{3}|j,j\rangle =(-1)^{l-l_{3}}{\sqrt {P}}{\sqrt {Q}}$

$P:={\frac {(2j+1)!(l+s-j)!}{(l+s+j+1)!(j+l-s)!(j+s-l)!}}$

$Q:={\frac {(l+l_{3})!(s+s_{3})!}{(l-l_{3})!(s-s_{3})!}}$

Es existiert eine alleswissende, abschreckend komplizierte Racah-Formel für die Gesamtheit der CG-Koeffizienten.

Wigner-3j-Symbole sind eine alternative Form der CG-Koeffizienten.

${\begin{pmatrix}s&l&j\\s_{3}&l_{3}&-j_{3}\end{pmatrix}}=\langle s,l;s_{3},l_{3}|j,j_{3}\rangle {\frac {(-1)^{s-l+j_{3}}}{\sqrt {2j+1}}}$.

Kostprobe der Ausgabe

Clebsch-Gordan L=(1/2) S=(1/2)
${\begin{array}{l}\langle {\tfrac {1}{2}},{\tfrac {1}{2}};{\tfrac {1}{2}},{\tfrac {1}{2}}|1,1\rangle =1\\\langle {\tfrac {1}{2}},{\tfrac {1}{2}};{\tfrac {1}{2}},-{\tfrac {1}{2}}|1,0\rangle ={\sqrt {\frac {1}{2}}}\\\langle {\tfrac {1}{2}},{\tfrac {1}{2}};-{\tfrac {1}{2}},{\tfrac {1}{2}}|1,0\rangle ={\sqrt {\frac {1}{2}}}\\\langle {\tfrac {1}{2}},{\tfrac {1}{2}};-{\tfrac {1}{2}},-{\tfrac {1}{2}}|1,-1\rangle =1\\\langle {\tfrac {1}{2}},{\tfrac {1}{2}};{\tfrac {1}{2}},-{\tfrac {1}{2}}|0,0\rangle ={\sqrt {\frac {1}{2}}}\\\langle {\tfrac {1}{2}},{\tfrac {1}{2}};-{\tfrac {1}{2}},{\tfrac {1}{2}}|0,0\rangle =-1{\sqrt {\frac {1}{2}}}\\\end{array}}$

Clebsch-Gordan L=(1/2) S=1
${\begin{array}{l}\langle {\tfrac {1}{2}},1;{\tfrac {1}{2}},1|{\tfrac {3}{2}},{\tfrac {3}{2}}\rangle =1\\\langle {\tfrac {1}{2}},1;{\tfrac {1}{2}},0|{\tfrac {3}{2}},{\tfrac {1}{2}}\rangle ={\sqrt {\frac {2}{3}}}\\\langle {\tfrac {1}{2}},1;-{\tfrac {1}{2}},1|{\tfrac {3}{2}},{\tfrac {1}{2}}\rangle ={\sqrt {\frac {1}{3}}}\\\langle {\tfrac {1}{2}},1;{\tfrac {1}{2}},-1|{\tfrac {3}{2}},-{\tfrac {1}{2}}\rangle ={\sqrt {\frac {1}{3}}}\\\langle {\tfrac {1}{2}},1;-{\tfrac {1}{2}},0|{\tfrac {3}{2}},-{\tfrac {1}{2}}\rangle ={\frac {1}{3}}{\sqrt {6}}\\\langle {\tfrac {1}{2}},1;-{\tfrac {1}{2}},-1|{\tfrac {3}{2}},-{\tfrac {3}{2}}\rangle =1\\\langle {\tfrac {1}{2}},1;{\tfrac {1}{2}},0|{\tfrac {1}{2}},{\tfrac {1}{2}}\rangle ={\sqrt {\frac {1}{3}}}\\\langle {\tfrac {1}{2}},1;-{\tfrac {1}{2}},1|{\tfrac {1}{2}},{\tfrac {1}{2}}\rangle =-1{\sqrt {\frac {2}{3}}}\\\langle {\tfrac {1}{2}},1;{\tfrac {1}{2}},-1|{\tfrac {1}{2}},-{\tfrac {1}{2}}\rangle ={\sqrt {\frac {2}{3}}}\\\langle {\tfrac {1}{2}},1;-{\tfrac {1}{2}},0|{\tfrac {1}{2}},-{\tfrac {1}{2}}\rangle ={\frac {-1}{3}}{\sqrt {3}}\\\end{array}}$

Clebsch-Gordan L=(1/2) S=(3/2)
${\begin{array}{l}\langle {\tfrac {1}{2}},{\tfrac {3}{2}};{\tfrac {1}{2}},{\tfrac {3}{2}}|2,2\rangle =1\\\langle {\tfrac {1}{2}},{\tfrac {3}{2}};{\tfrac {1}{2}},{\tfrac {1}{2}}|2,1\rangle ={\frac {1}{2}}{\sqrt {3}}\\\langle {\tfrac {1}{2}},{\tfrac {3}{2}};-{\tfrac {1}{2}},{\tfrac {3}{2}}|2,1\rangle ={\frac {1}{2}}\\\langle {\tfrac {1}{2}},{\tfrac {3}{2}};{\tfrac {1}{2}},-{\tfrac {1}{2}}|2,0\rangle ={\sqrt {\frac {1}{2}}}\\\langle {\tfrac {1}{2}},{\tfrac {3}{2}};-{\tfrac {1}{2}},{\tfrac {1}{2}}|2,0\rangle ={\frac {1}{2}}{\sqrt {2}}\\\langle {\tfrac {1}{2}},{\tfrac {3}{2}};{\tfrac {1}{2}},-{\tfrac {3}{2}}|2,-1\rangle ={\frac {1}{2}}\\\langle {\tfrac {1}{2}},{\tfrac {3}{2}};-{\tfrac {1}{2}},-{\tfrac {1}{2}}|2,-1\rangle ={\frac {1}{2}}{\sqrt {3}}\\\langle {\tfrac {1}{2}},{\tfrac {3}{2}};-{\tfrac {1}{2}},-{\tfrac {3}{2}}|2,-2\rangle =1\\\langle {\tfrac {1}{2}},{\tfrac {3}{2}};{\tfrac {1}{2}},{\tfrac {1}{2}}|1,1\rangle ={\frac {1}{2}}\\\langle {\tfrac {1}{2}},{\tfrac {3}{2}};-{\tfrac {1}{2}},{\tfrac {3}{2}}|1,1\rangle ={\frac {-1}{2}}{\sqrt {3}}\\\langle {\tfrac {1}{2}},{\tfrac {3}{2}};{\tfrac {1}{2}},-{\tfrac {1}{2}}|1,0\rangle ={\sqrt {\frac {1}{2}}}\\\langle {\tfrac {1}{2}},{\tfrac {3}{2}};-{\tfrac {1}{2}},{\tfrac {1}{2}}|1,0\rangle ={\frac {-1}{2}}{\sqrt {2}}\\\langle {\tfrac {1}{2}},{\tfrac {3}{2}};{\tfrac {1}{2}},-{\tfrac {3}{2}}|1,-1\rangle ={\frac {1}{2}}{\sqrt {3}}\\\langle {\tfrac {1}{2}},{\tfrac {3}{2}};-{\tfrac {1}{2}},-{\tfrac {1}{2}}|1,-1\rangle ={\frac {-1}{2}}\\\end{array}}$

Clebsch-Gordan L=(1/2) S=2
${\begin{array}{l}\langle {\tfrac {1}{2}},2;{\tfrac {1}{2}},2|{\tfrac {5}{2}},{\tfrac {5}{2}}\rangle =1\\\langle {\tfrac {1}{2}},2;{\tfrac {1}{2}},1|{\tfrac {5}{2}},{\tfrac {3}{2}}\rangle =2{\sqrt {\frac {1}{5}}}\\\langle {\tfrac {1}{2}},2;-{\tfrac {1}{2}},2|{\tfrac {5}{2}},{\tfrac {3}{2}}\rangle ={\sqrt {\frac {1}{5}}}\\\langle {\tfrac {1}{2}},2;{\tfrac {1}{2}},0|{\tfrac {5}{2}},{\tfrac {1}{2}}\rangle ={\sqrt {\frac {3}{5}}}\\\langle {\tfrac {1}{2}},2;-{\tfrac {1}{2}},1|{\tfrac {5}{2}},{\tfrac {1}{2}}\rangle ={\frac {1}{5}}{\sqrt {10}}\\\langle {\tfrac {1}{2}},2;{\tfrac {1}{2}},-1|{\tfrac {5}{2}},-{\tfrac {1}{2}}\rangle ={\sqrt {\frac {2}{5}}}\\\langle {\tfrac {1}{2}},2;-{\tfrac {1}{2}},0|{\tfrac {5}{2}},-{\tfrac {1}{2}}\rangle ={\sqrt {\frac {3}{5}}}\\\langle {\tfrac {1}{2}},2;{\tfrac {1}{2}},-2|{\tfrac {5}{2}},-{\tfrac {3}{2}}\rangle ={\sqrt {\frac {1}{5}}}\\\langle {\tfrac {1}{2}},2;-{\tfrac {1}{2}},-1|{\tfrac {5}{2}},-{\tfrac {3}{2}}\rangle ={\frac {2}{5}}{\sqrt {5}}\\\langle {\tfrac {1}{2}},2;-{\tfrac {1}{2}},-2|{\tfrac {5}{2}},-{\tfrac {5}{2}}\rangle =1\\\langle {\tfrac {1}{2}},2;{\tfrac {1}{2}},1|{\tfrac {3}{2}},{\tfrac {3}{2}}\rangle ={\sqrt {\frac {1}{5}}}\\\langle {\tfrac {1}{2}},2;-{\tfrac {1}{2}},2|{\tfrac {3}{2}},{\tfrac {3}{2}}\rangle =-2{\sqrt {\frac {1}{5}}}\\\langle {\tfrac {1}{2}},2;{\tfrac {1}{2}},0|{\tfrac {3}{2}},{\tfrac {1}{2}}\rangle ={\sqrt {\frac {2}{5}}}\\\langle {\tfrac {1}{2}},2;-{\tfrac {1}{2}},1|{\tfrac {3}{2}},{\tfrac {1}{2}}\rangle ={\frac {-1}{5}}{\sqrt {15}}\\\langle {\tfrac {1}{2}},2;{\tfrac {1}{2}},-1|{\tfrac {3}{2}},-{\tfrac {1}{2}}\rangle ={\sqrt {\frac {3}{5}}}\\\langle {\tfrac {1}{2}},2;-{\tfrac {1}{2}},0|{\tfrac {3}{2}},-{\tfrac {1}{2}}\rangle =-1{\sqrt {\frac {2}{5}}}\\\langle {\tfrac {1}{2}},2;{\tfrac {1}{2}},-2|{\tfrac {3}{2}},-{\tfrac {3}{2}}\rangle =2{\sqrt {\frac {1}{5}}}\\\langle {\tfrac {1}{2}},2;-{\tfrac {1}{2}},-1|{\tfrac {3}{2}},-{\tfrac {3}{2}}\rangle ={\frac {-1}{5}}{\sqrt {5}}\\\end{array}}$

Clebsch-Gordan L=1 S=(1/2)
${\begin{array}{l}\langle 1,{\tfrac {1}{2}};1,{\tfrac {1}{2}}|{\tfrac {3}{2}},{\tfrac {3}{2}}\rangle =1\\\langle 1,{\tfrac {1}{2}};1,-{\tfrac {1}{2}}|{\tfrac {3}{2}},{\tfrac {1}{2}}\rangle ={\sqrt {\frac {1}{3}}}\\\langle 1,{\tfrac {1}{2}};0,{\tfrac {1}{2}}|{\tfrac {3}{2}},{\tfrac {1}{2}}\rangle ={\sqrt {\frac {2}{3}}}\\\langle 1,{\tfrac {1}{2}};0,-{\tfrac {1}{2}}|{\tfrac {3}{2}},-{\tfrac {1}{2}}\rangle ={\frac {1}{3}}{\sqrt {6}}\\\langle 1,{\tfrac {1}{2}};-1,{\tfrac {1}{2}}|{\tfrac {3}{2}},-{\tfrac {1}{2}}\rangle ={\sqrt {\frac {1}{3}}}\\\langle 1,{\tfrac {1}{2}};-1,-{\tfrac {1}{2}}|{\tfrac {3}{2}},-{\tfrac {3}{2}}\rangle =1\\\langle 1,{\tfrac {1}{2}};1,-{\tfrac {1}{2}}|{\tfrac {1}{2}},{\tfrac {1}{2}}\rangle ={\sqrt {\frac {2}{3}}}\\\langle 1,{\tfrac {1}{2}};0,{\tfrac {1}{2}}|{\tfrac {1}{2}},{\tfrac {1}{2}}\rangle =-1{\sqrt {\frac {1}{3}}}\\\langle 1,{\tfrac {1}{2}};0,-{\tfrac {1}{2}}|{\tfrac {1}{2}},-{\tfrac {1}{2}}\rangle ={\frac {1}{3}}{\sqrt {3}}\\\langle 1,{\tfrac {1}{2}};-1,{\tfrac {1}{2}}|{\tfrac {1}{2}},-{\tfrac {1}{2}}\rangle =-1{\sqrt {\frac {2}{3}}}\\\end{array}}$

Clebsch-Gordan L=1 S=1
${\begin{array}{l}\langle 1,1;1,1|2,2\rangle =1\\\langle 1,1;1,0|2,1\rangle ={\frac {1}{2}}{\sqrt {2}}\\\langle 1,1;0,1|2,1\rangle ={\frac {1}{2}}{\sqrt {2}}\\\langle 1,1;1,-1|2,0\rangle ={\frac {1}{2}}{\sqrt {\frac {2}{3}}}\\\langle 1,1;0,0|2,0\rangle ={\frac {1}{3}}{\sqrt {6}}\\\langle 1,1;-1,1|2,0\rangle ={\frac {1}{2}}{\sqrt {\frac {2}{3}}}\\\langle 1,1;0,-1|2,-1\rangle ={\frac {1}{2}}{\sqrt {2}}\\\langle 1,1;-1,0|2,-1\rangle ={\frac {1}{2}}{\sqrt {2}}\\\langle 1,1;-1,-1|2,-2\rangle =1\\\langle 1,1;1,0|1,1\rangle ={\frac {1}{2}}{\sqrt {2}}\\\langle 1,1;0,1|1,1\rangle ={\frac {-1}{2}}{\sqrt {2}}\\\langle 1,1;1,-1|1,0\rangle ={\frac {1}{2}}{\sqrt {2}}\\\langle 1,1;-1,1|1,0\rangle ={\frac {-1}{2}}{\sqrt {2}}\\\langle 1,1;0,-1|1,-1\rangle ={\frac {1}{2}}{\sqrt {2}}\\\langle 1,1;-1,0|1,-1\rangle ={\frac {-1}{2}}{\sqrt {2}}\\\langle 1,1;1,-1|0,0\rangle ={\frac {1}{3}}{\sqrt {3}}\\\langle 1,1;0,0|0,0\rangle ={\frac {-1}{3}}{\sqrt {3}}\\\langle 1,1;-1,1|0,0\rangle ={\frac {1}{3}}{\sqrt {3}}\\\end{array}}$

Python-Listing

def colsalign(s) : # align string array s columns, strings start with ' '
  n=len(s); m=32; cols=[0]*m; kmax=0
  for i in range(n) :
    t=s[i]; h=0; j=t.find(' ',h+1); k=0
    while j>=0 : 
      if ((j-h)>cols[k]) : cols[k]=j-h-1
      h=j; j=t.find(' ',h+1); k+=1
    j=len(t) 
    if ((j-h)>cols[k]) : cols[k]=j-h-1
    if (k+1)>kmax : kmax=k+1
  # dbg print('Cols: '+str(cols[:kmax]))
  for i in range(n) :
    t=s[i]; h=0; j=t.find(' ',h+1); lt=len(t); k=0; v=''
    while j<=lt :
      u=t[h:j]+(' '*(cols[k]-(j-h-1))); v+=u  
      h=j; j= (lt+10) if (h==lt) else t.find(' ',h+1)
      k+=1; j= lt if (j<0) else j
    s[i]= v
  return s

def hasfact(x,z) : return (x>0) and ((x % z)==0)
def havefact(x,y,z) : # x,y have common factor z?
  fx=hasfact(x,z); fy=hasfact(y,z) 
  return (fx and fy) or (fx and (y==0)) or (fy and (x==0))

def normalize(x) : # rational-square root number x
  pn=[2,3,5,7,11,13,17,19,23]; ln= len(pn) # small prime nbrs
  ps=[0,1,4,9,16,25,36,49,64,81,100,121,144]; ls=len(ps) # small squares
  pr,pi,q,sp,sq = tuple(x)
  pr,sgr= (-pr,-1) if (pr<0) else (pr,1) 
  pi,sgi= (-pi,-1) if (pi<0) else (pi,1)
  for i in range(2,len(ps)) :
    z=ps[i]
    while hasfact(sp,z) : sp=div(sp,z); pr *= i; pi *= i
    while hasfact(sq,z) : sq=div(sq,z); q *=i
  for i in range(ln) :
    z=pn[i]
    while hasfact(sp,z) and hasfact(sq,z) : sp=div(sp,z); sq=div(sq,z)
    while havefact(pr,pi,z) and hasfact(q,z) :
      pr=div(pr,z); pi=div(pi,z); q=div(q,z)
  if (pr==0)and(pi==0) : q=1;sp=1;sq=1
  x[0],x[1],x[2],x[3],x[4] = (pr*sgr,pi*sgi,q,sp,sq)

def real(x) : # rational square root number, get real part
  pr,pi,q,sp,sq = tuple(x); return sqrt(sp/(sq+0.0))*pr/(q+0.0)

def encodenb(x) :
  pr,pi,q,sp,sq = tuple(x); s=''; w=''
  if (sp!=1)or(sq!=1) :
    w='\\sqrt{'; 
    w += str(sp) if(sq==1) else ('\\frac{'+str(sp)+'}{'+str(sq)+'}')
    w +='}'
  if pr != 0 :
    u=str(pr); s += u if (q==1) else ('\\frac{'+u+'}{'+str(q)+'}')
    s = w if((s=='1')and(w!='')) else (s if (s=='1') else (s+w))
  if pi != 0 : 
    if (pr!=0) and (pi>0) : s +='+i'
    elif (pi<0) : s += '-i'; pi=-pi
    else : s += 'i'
    u=str(pi); u='' if ((pi==1)and(q==1)) else u
    s += u if (q==1) else ('\\frac{'+u+'}{'+str(q)+'}')
    s = s+w
  if (pr==0)and(pi==0) : s='0' 
  return s

def finetune(x) : 
  pn=[2,3,5,7,11,13,17,19,23]; ln= len(pn) # small primes
  xx=reformat(x); pr,pi,q,sp,sq = tuple(xx)
  nmax=max(sp,sq);imax=0 # minimize values under root
  for i in range(2,5) : 
    tp=sp*i*i; tq=sq
    for k in range(ln) :
      z=pn[k]
      while hasfact(tp,z) and hasfact(tq,z) : tp=div(tp,z); tq=div(tq,z)
    if max(tp,tq)<=nmax: nmax=max(tp,tq);imax=i; pmax=tp;qmax=tq
  if imax>0 : q *=imax; sp=pmax; sq=qmax
  return [pr,pi,q,sp,sq]

def displaynb(x) : # x is 5-tuple. to readable string, root in square brackets
  # bug: ambiguous output rt(2/3)=2*rt(1/6)= (1/3)*sqrt(6) 
  # minimize max(p,q) in sqrt(p/q)? reduce {4,9,1/4,1/9)*p/q= p'/q'?
  # prefer sqrt(p) over sqrt(1/q)
  xx=reformat(x); pr,pi,q,sp,sq = tuple(xx); s=''; w=''
  if (sp!=1)or(sq!=1) :
    w='['; 
    w += str(sp) if(sq==1) else (str(sp)+'/'+str(sq))
    w +=']'
  if pr != 0 :
    u=str(pr); s += u if (q==1) else (u+'/'+str(q))
    s = w if((s=='1')and(w!='')) else (s if (s=='1') else (s+w))
  if pi != 0 : 
    if (pr!=0) and (pi>0) : s +='+i'
    elif (pi<0) : s += '-i'; pi=-pi
    else : s += 'i'
    u=str(pi); u='' if ((pi==1)and(q==1)) else u
    s += u if (q==1) else (u+'/'+str(q))
    s = s+w
  if (pr==0)and(pi==0) : s='0' 
  return s

def reformat(x) : # int, rational, rcomplex to root numbers
  if type(x)==type(1) : z=[x,0,1,1,1]; n=-1
  elif type(x)==type([]) : n=len(x)
  if n==1 : z=[x[0],0,1,1,1]
  if n==2 : z=[x[0],0,x[1],1,1]
  if n==3 : z=[x[0],x[1],x[2],1,1]
  if n==5 : z=x
  if (n>5)or(n==0)or(n==4): print('reformat '+str(x)+' not supported.'); exit()
  return z

def multip(x,y) : # format (u,v,w,p,q) = (u+iv)/w sqrt(p/q)
  xx=reformat(x); yy=reformat(y)
  xpr,xpi,xq,xsp,xsq= tuple(xx); ypr,ypi,yq,ysp,ysq= tuple(yy)
  pr= xpr*ypr-xpi*ypi; pi= xpr*ypi+ypr*xpi; q=xq*yq; sp=xsp*ysp;sq=xsq*ysq
  z= [pr,pi,q,sp,sq]; normalize(z); return z

def addit(x,y) :
  xx=reformat(x); yy=reformat(y)
  return lincomb(xx,yy,1,1)

def sqroot(x) : # input should be (p,0,q,1,1) 
  xx=reformat(x); pr,pi,q,sp,sq = tuple(xx)
  if (pi!=0)or(sp!=1)or(sq!=1) : print('sqroot undefined.'); exit()
  sp=pr;sq=q;pr=1;q=1; z=[pr,pi,q,sp,sq]; normalize(z); return z

def recip(x) :
  xx=reformat(x); pr,pi,q,sp,sq = tuple(xx); qq=pr*pr+pi*pi; 
  y= [pr*q,-pi*q,qq,sq,sp]; normalize(y); return y

def imagi(x) : pr,pi,q,sp,sq = tuple(x); return [-pi,pr,q,sp,sq]

def scale(x,f) :
  xx=reformat(x); pr,pi,q,sp,sq = tuple(xx); return [pr*f,pi*f,q,sp,sq]

def lincomb(x,y,f,g) : # f,g integer, x,y fractional root numbers
  # a sqrt(r/s) + b sqrt(p/q) ?= sqrt(rp/qs) ?
  # can only be combined if r/s and p/q differ by square factor T/U ?
  # first write rq/qs, ps/qs and test if rq,ps have square factor
  ps=[0,1,4,9,16,25,36,49,64,81,100,121,144]; ls=len(ps) # small squares
  xpr,xpi,xq,xsp,xsq= tuple(x); xpr *= f; xpi *=f
  ypr,ypi,yq,ysp,ysq= tuple(y); ypr *= g; ypi *=g
  if (xpr==0)and(xpi==0) :   z=[ypr,ypi,yq,ysp,ysq]
  elif (ypr==0)and(ypi==0) : z= [xpr,xpi,xq,xsp,xsq] 
  else :
    sq=xsq*ysq; xp=xsp*ysq; yp=ysp*xsq # test common sqroot-denom
    for i in range(2,len(ps)) :
      z=ps[i]
      while hasfact(xp,z) : xp=div(xp,z); xpr *= i; xpi *= i
      while hasfact(yp,z) : yp=div(yp,z); ypr *= i; ypi *= i 
    # here xp/sq new sqrt arg? if yp=n xp, rewrite x ?
    if xp != yp : print('lincomb:sqrt clash'+str([xp,yp,x,y])); exit()
    sp=xp; q=xq*yq; pr=xpr*yq+ypr*xq; pi=xpi*yq+ypi*xq
    z= [pr,pi,q,sp,sq]
  normalize(z); return z

def matshow(s,x,new=0) :
  # new= 1 simplified, columns. new=0 raw.
  # new= 2 use encodenb, latex-like, check square hermitian 
  n=len(x); m=len(x[0]); print(s); t=''; ok=True; tt=['']*n
  for i in range(n) :
    for k in range(m) :
      if (new==2) and (k>0) : t +=' &'
      if (new==1) : t +=' '+displaynb(x[i][k])
      if (new==2) : t +=' '+encodenb(x[i][k])
      if new==0 : t +='  '+str(x[i][k]) # beware about v, we need a copy!
      if new==2 :
        u=x[i][k]; v= [] + x[k][i]; v[1]= -v[1]; ok= ok and isequal(u,v)
    if (new==2) and (i<(n-1)) : t +=' \\\\';
    tt[i]=t; t=''
  if new==1 : colsalign(tt)
  for i in range(n) : print(tt[i])
  if new==2 : print('Test Hermitian: '+str(ok))
  return tt

def vecshow(s,x) :
  n=len(x); print(s); t=''
  for k in range(n) : t +=' '+displaynb(x[k])
  print(t)

def div(p,q) :
  if p>=0 : return int(p/q)
  else : return -int(-p/q)

def getprimes(n) : # want n prime numbers
  p=[2]; i=1; plast=2
  while i<=n :
    q=plast+1; found=False
    while not found:
      j=0
      while ((q%p[j])!=0) and ((p[j]*p[j])<=q) :  j +=1
      if (q%p[j] != 0) : p += [q]; plast=q; found=True; i+=1
      else : q+=1
  return p

def greatestdiv(x) : # greatest divisor, skip 0s in list
  n=len(x); y=[0]*n; w=[0]*n; h=0; ymin=-1; gcd=1
  pn=[2,3,5,7,11,13,17,19,23,29,31,37,41,43,47]; ln= len(pn) # small primes
  for i in range(n) :
    z=abs(x[i])
    if z>0 : ymin= z if(ymin<0) else min(ymin,z); w[h]=z; h+=1
  w=w[:h]; ip=0
  while (ip<ln)and(pn[ip]<=ymin) :
    z=pn[ip]; ok=True
    while ok :
      for i in range(h) : ok= ok and ((w[i]%z)==0)
      if ok : 
        for i in range(h) : w[i]= div(w[i],z)
        ymin= div(ymin,z); gcd *= z
    ip +=1 
  return gcd    

def smallestmult(x) : # smallest multiple of nonzero numbers
  n=len(x); y=abs(x[0])
  for i in range(1,n) : 
    z=abs(x[i]); y = y if (z==0) else div(y*z,greatestdiv([y,z]))
  return y

def divide(p,q) : # improve integer fractions
  def hasfact(x,z) : return (x>0) and ((x % z)==0)
  pn=[2,3,5,7,11,13,17,19,23,29,31,37]; ln= len(pn) # small prime nbrs
  p,q,s = (abs(p),abs(q),-1) if ((p*q)<0) else (abs(p),abs(q),1) 
  if p==0 : q=1
  for i in range(ln) :
    z=pn[i]
    while hasfact(p,z) and hasfact(q,z) : p=div(p,z); q=div(q,z)
  return [s*p,q]

def jj3(j,j3, l,l3, s,s3) : # get coeffs of vector (j,j3). numbers halfints
  j2=div(j[0]*2,j[1]); j32=div(j3[0]*2,j3[1]) 
  l2=div(l[0]*2,l[1]); l32=div(l3[0]*2,l3[1]) 
  s2=div(s[0]*2,s[1]); s32=div(s3[0]*2,s3[1]) 
  ok=True; jmax=abs(l2)+abs(s2); jmin=abs(l2-s2); z=[0,1]; sign=0
  if (s2 !=1) : print('spin s '+str(s)+' not supported'); ok=False
  if (j2>jmax)or(j2<jmin) : print('j='+str(j)+' outside |l-s|,(l+s)'); ok=False
  if (j32>j2)or(j32<-j2)or(l32>l2)or(l32<-l2)or(s32>s2)or(s32<-s2) :
    print('one of '+str([j3,l3,s3])+' out of range'); ok=False
  if (j2>l2) : # contribute l3,s3= (j3-0.5,+) and (j3+0.5,-)
    if (l32==(j32-1))and(s32==1) : z=[l2+j32+1,2*l2+2]; sign=1
    if (l32==(j32+1))and(s32==-1) : z=[l2-j32+1,2*l2+2]; sign=1
  else  : # contribute l3,s3= (j3+0.5,-) and (j3-0.5,+)
    if (l32==(j32+1))and(s32==-1) : z=[l2+j32+1,2*l2+2]; sign=1
    if (l32==(j32-1))and(s32==1) : z=[l2-j32+1,2*l2+2]; sign=-1
  if (not ok)or(z[0]==0) : z= [0,0,1,1,1]
  else : z=[sign,0,1,z[0],z[1]]
  return z

def testcleb(j,l,s=[1,2]) : # Clebsch-Gordan case spin 1/2 only
  # make a matrix of coeffs for (j,j3) in (l,l3)(x)(s,s3)
  # return matrix mx, rows in order j3max...j3min
  # mx rows ordered: (lmax,smax)...(lmax..smin) (lmax-1,smax)...
  noisy=2 # example l=[1,1]; j=[3,2]; s fixed for spin 1/2
  nl= div(2*l[0],l[1])+1
  ns= div(2*s[0],s[1])+1
  nj= div(2*j[0],j[1])+1; mx=[[]]*nj 
  for ij in range(nj) :
    j3=j if (ij==0) else [j3[0]-j3[1],j3[1]]; hj=0; mx[ij]=[[]]*(nl*ns) 
    if ij==0 : s3val=[[]]*(nl*ns); l3val=[[]]*(nl*ns)
    for il in range(nl) :
      l3=l if (il==0) else [l3[0]-l3[1],l3[1]] 
      for ks in range(ns) :
        s3=s if (ks==0) else [s3[0]-s3[1],s3[1]] 
        if ij==0 : s3val[hj]=s3; l3val[hj]=l3
        z= jj3(j,j3, l,l3,s,s3); normalize(z); mx[ij][hj]=z; hj+=1
        if noisy==1 : print(str([ij,il,ks,z]))
  if noisy==2 : matshow('',mx,1)
  return mx,s3val,l3val

def rootclass(x) : # x is list of real 5-tuple numbers, are they additive?
  noisy=False; # maybe root-incompatible subsets cancel?
  if noisy : print('  rootclass'+str(x))
  n=len(x); q=[0]*n; p=[0]*n; r=[0]*n; clas=[0]*n; pclas=[0]*n; ncl=0
  for i in range(n) : p[i]=x[i][3];q[i]=x[i][4];r[i]=x[i][0]
  qq= smallestmult(q); ls=12 
  for i in range(n) : 
    p[i]= div(qq*p[i],q[i]); j=0; found=False 
    for k in range(2,ls) : # take squares out of root
      z= k*k
      while hasfact(p[i],z) : p[i]=div(p[i],z); r[i] *= k
    while (not found) and (j<ncl) : found=(p[j]==pclas[j]); j+=1
    if not found : pclas[j]=p[i]; clas[i]=j; ncl +=1
    else : clas[i]= j-1 # equivalence class of same p factor
  if noisy : print('p-classes: '+str(pclas[:ncl])); 
  sp=[0]*ncl; sq=[1]*ncl
  for k in range(ncl) :
    for i in range(n) :
      if clas[i]==k : sp[k]= r[i]*sq[k]+ sp[k]*x[i][2]; sq[k]=sq[k]*x[i][2]
  if noisy : print('p-sums: '+str(sp)) # only one of them non-zero! 
  for i in range(n) : x[i]=[r[i],0,x[i][2],p[i],qq]
  return clas,pclas[:ncl]

def scalprod(x,y) : # square-root compatible subset
  n=len(x); m=len(y); w=[0]*n; z=0; bigbug=False
  if bigbug : print('x= '+str(x)+'\ny= '+str(y))
  if n != m : print('scalprod length mismatch.')
  for i in range(n) : w[i]=multip(x[i],y[i]); z=addit(z,w[i]) 
  if bigbug : print('w= '+str(w))
  return z

def quotient(x,y) : z=recip(y); return multip(z,x)
def iszero(x) :
  if type(x)==type(1) : return (x==0)
  else : return (x[0]==0)and(x[1]==0)

def orthogonize(y,x) : # make vector y orthonormal to vectorlist x[i]
  # as startvector y try maximal n linear indep ones, y=delta(i,k)
  n=len(y);m=len(x); ok=True;sok=True; z=[0]*n; norm=[[]]*m;pxy=[[]]*m
  for i in range(m) : ok= ok and (len(x[i])==n)
  if not ok : print('Vector Length Error')
  for i in range(m) :
    for j in range(i,m) :
      sp=scalprod(x[i],x[j]) # root-non-compatible here?
      if i==j : norm[i]=sp; pxy[i]=scalprod(y,x[i])
      if i != j : sok= sok and iszero(sp)
  if not sok : print('Orthogonality error')
  for k in range(n) : z[k]= reformat(y[k])
  for i in range(m) : 
    u= quotient(pxy[i],norm[i])
    for k in range(n) : v=multip(u,x[i][k]); z[k]=lincomb(z[k],v,1,-1)
  nz= scalprod(z,z); ok=True
  if iszero(nz) : print('Orthogonize result zero.'); ok=False
  else :
    for i in range(m) : sp=scalprod(x[i],z); ok= ok and iszero(sp)
    if not ok : print('Result z not orthogonal.')
    nz= sqroot(nz); qz= recip(nz)
    for k in range(n) : z[k] = multip(z[k],qz)
    sp=scalprod(z,z); sp=lincomb(sp,reformat(1),1,-1); nok= iszero(sp)
    if not nok : print('Result z not normalized.') 
  return z,ok

def signflip(z) : # make first entry positive
  n=len(z); flip= z[0][0]<0; nf= n if flip else 0
  for i in range(nf) : z[i][0]= -z[i][0]
  if flip : print('Signflip')

def orthoplement(x,noisy) : # orthogonal complement
  # if all of x root-compatible, forget the roots
  # solve linear equ on rationals, then normalize
  n=len(x); m=len(x[0]); ok=True
  if noisy : print(' orthoplem '+str(n)+' vectors of length '+str(m))
  if (n+1)!=m : print('Bad number of of vectors!') 
  if m==2 : return [scale(x[0][1],1),scale(x[0][0],-1)] # other sign?
  for i in range(n) : clas,pclas=rootclass(x[i]); ok= ok and (len(pclas)==1)
  if not ok : print('Roots not compatible.')
  if noisy : print('Orthoplement '+str(n))
  n2=n; n2=0 #dbg, n2=n is fatally bad
  for k in range(n2) :
    print(str(x[k]))
    for j in range(len(x[k])) :
      if len(x[k][j])>4 : x[k][j][3]=1;x[k][j][4]=1 # roots away is fatal
  y=[0]*m; found=False; i=0
  while (not found) and (i<m) : # try all unit vector startpoints
    y[i]=1; z,found=orthogonize(y,x); y[i]=0; i+=1
  if not found : print('Cannot orthoplement m='+str(m)+', exit.'); exit() 
  else : signflip(z)
  return z
  # faster? guess int-only startvector y; Fi=prod(k,k!=i,(xk,xk))(xi,y);
  # G=prod(k,(xk,xk)); y'= Gy- sum(i,Fi xi) . if zero, modify y. int-solution?

def initclebsch(ja2,jb2,noisy) : # ja,jb half-integral, ja2,jb2 their double
  # matrices quad-index: [n][d][ja2+1][jb2+1], d=(da+db-x)...|da-db|(2)
  max2=ja2+jb2; min2=abs(ja2-jb2); n=div(max2-min2,2)+1; da=ja2+1;db=jb2+1
  cg=[[]]*n; sd=0; d=max2+1; s='CG-Matrices Cols='+str(da*db)+' Rows=' 
  for i in range(n) :
    cg[i]=[[]]*d ; s+=' '+str(d); sd+=d
    for j in range(d) :
      cg[i][j]=[[]]*da
      for k in range(da) : cg[i][j][k]=[0]*db # defines cg[i,j,k,h]
    d -= 2    
  j2=ja2+jb2; i=0; cg[i][0][0][0]=1
  if noisy : print(s+' Sd='+str(sd))
  return cg # quadruple-index array: j,m, ma,mb

def gathervectors(cgm,jmax2,m2, ja2,jb2, noisy) : # cgm vectors filled at m2?
  # extract relevant row vectors from da*db-size submatrices
  # put together the row vectors for j2.3==m2 
  # to start (j2=m2), find an orthogonal vector to those found here
  n=len(cgm); mydim=m2+1; k=0; x=[[]]*n # all cgm with higher dims contribute
  terms=[0]*n ; nj=[0]*n; j2=m2; h=0; jm2=jmax2
  for i in range(n) :
    cg=cgm[i]; cdim=len(cg); okdim= cdim>mydim # if ok, add a vector x[k]
    if (noisy and okdim) : print('gather from submatrix dim='+str(cdim))
    # extract the da*db matrix entries to m2
    # cg submatrix i, row that counts belongs to m2: 
    ma2=jmax2;j=0; y=[0]*(jmax2+1) # items in cg start with max=ma2, min=mb2
    while ma2>ja2 : ma2 -=1
    k= div(jm2-m2,2) # relevant row in matrix cg
    while okdim and (ma2>=(-jmax2)) : # scan diagonal ma2+mb2=m2
      mb2=m2-ma2; ok= (abs(ma2)<=ja2)and(abs(mb2)<=jb2); z=0
      if ok : ia=div(ja2-ma2,2);ib=div(jb2-mb2,2); z=cg[k][ia][ib]; terms[i]+=1
      if noisy : print('     ma2,mb2,ok='+str([ma2,mb2,ok])) 
      if ok : y[j]=z; j+=1; dn= 1 if (not iszero(z)) else 0; nj[i]+=dn 
      ma2 -=2
    if (terms[i]>0) and (nj[i]>0) : x[h]=y[:j]; h+=1
    jm2 -= 2
  if noisy : print('Gather Terms: '+str(terms)+' NonZero: '+str(nj))
  return x[:h]

def downsweep(cgm, z, j2, ja2,jb2,noisy) :  
  # one diagonal z is ready: for m2=j2, all ma,mb with ma+mb=m2 are defined.
  # m index ordering is top down. input vector z = diagonal m2=j2.
  nz=len(z); da=ja2+1; db=jb2+1; ix= div(ja2+jb2-j2,2); cg=cgm[ix]; d=j2+1; k=0
  if noisy : print('downsweep ix,nz,j2,ja2,jb2='+str([ix,nz,j2,ja2,jb2]))
  if noisy : print(' input vector='+str(z));
  jmax2=ja2+jb2
  ma2=jmax2 # items in z ordered to start with maximal ma2, minimal mb2?
  while ma2>ja2 : ma2 -=1  # j2=2 ja2=jb2=1, must start ma2=1 not 2!
  while ma2>=(-jmax2) : # insert start-diagonal from vector z to cg[0] row
    mb2=j2-ma2; ok= (abs(ma2)<=ja2)and(abs(mb2)<=jb2); i=0;ia=-1;ib=-1
    if ok :
      ia=div(ja2-ma2,2); ib=div(jb2-mb2,2); cg[i][ia][ib]=z[k]; k+=1
      if noisy: print(' ins ma2,mb2, i,ia,ib,z='+str([ma2,mb2,i,ia,ib,z[k-1]]))
    ma2 -=2
  m2=j2; i=0; nput=0; imax=j2 # reach source row cg[2j]
  while m2>(-jmax2) : # fetch data from index m2
    ma2=jmax2
    if noisy : print(' Pass i,j2,m2,ja2='+str([i,j2,m2,ja2]))
    while ma2>ja2 : ma2 -=1
    while ma2>=(-jmax2) :
      mb2=m2-ma2; ok= (abs(ma2)<=ja2)and(abs(mb2)<=jb2)and(i<imax)
      if ok :  # add stuff to cg[ia+1][ib] and cg[ia][ib+1] if in range
        ia=div(ja2-ma2,2); ib=div(jb2-mb2,2); x=cg[i][ia][ib]
        q=j2*(j2+2)-m2*(m2-2) # divisor for norm
      if ok and ((ia+1)<da) :
        if noisy : print(' put1 i,ia,ib'+str([i,ia,ib, i+1,ia+1,ib,x]))
        y=ja2*(ja2+2)-ma2*(ma2-2); y=[y,q]; y= sqroot(y); nput+=1
        y=multip(y,x); z=cg[i+1][ia+1][ib]; cg[i+1][ia+1][ib]=addit(z,y)
      if ok and ((ib+1)<db) :
        if noisy : print(' put2 i,ia,ib'+str([i,ia,ib, i+1,ia,ib+1,x]))
        y=jb2*(jb2+2)-mb2*(mb2-2); y=[y,q]; y= sqroot(y); nput+=1
        y=multip(y,x); z=cg[i+1][ia][ib+1]; cg[i+1][ia][ib+1]=addit(z,y)
      ma2 -= 2
    m2 -=2; i+=1
  if noisy : print('Downsweep done, nput='+str(nput))

def cgshow(j2,cg) : 
  nj=len(cg); na=len(cg[0]); nb=len(cg[0][0]); ja2=na-1; jb2=nb-1; m2=j2
  print('CG-Matrix Nj,j2,ja2,jb2='+str([nj,j2,ja2,jb2])); s=''
  for i in range(nj) :
    cgi=cg[i]; na=len(cgi); nb=len(cgi[0]) 
    print('MatrixRow '+str(i)+' J2,M2='+str(j2)+','+str(m2))
    for j in range(na) :
      for k in range(nb) : z=finetune(cgi[j][k]); s+=' '+displaynb(z)
    print(s); s=''; m2 -=2
  
def clebschgordan(ja2,jb2,noisy=True) : # ja2,jb2 are double of j-quanta.
  # ja=(ja2/2) may be called 'L' and jb=(jb2/2) called 'S'
  cgm= initclebsch(ja2,jb2,noisy)
  j2=ja2+jb2;k2=j2; jmax2=j2; jmin2=abs(ja2-jb2);  i=0; z=[1]
  while j2>= jmin2 :
    if i>0 : x=gathervectors(cgm,jmax2,j2,ja2,jb2,noisy)
    if i>0 : z=orthoplement(x,noisy)
    downsweep(cgm, z, j2,ja2,jb2, noisy); j2-=2; i+=1
  if noisy :
    for i in range(len(cgm)) : cgshow(k2,cgm[i]); k2-=2
  return cgm

def wikidump(cg) : # list nonzero entries in latex style
  def wnum(n2) :
    s='' if (n2>=0) else '-'; na=abs(n2); whole= (na%2)==0; n=div(na,2)
    if whole : s+=str(n)
    else : s+='\\tfrac{'+str(na)+'}2'
    return s
  def label(l2,l32,s2,s32,j2,j32,x) :
    s='\\langle'+wnum(l2)+','+wnum(s2)+';'+wnum(l32)+','+wnum(s32)+'|'
    s+=wnum(j2)+','+wnum(j32)+'\\rangle='+encodenb(reformat(x))+'\\\\\n'
    return s 
  nj=len(cg); na=len(cg[0][0]); nb=len(cg[0][0][0]); s='\\begin{array}{l}\n'
  for ij in range(nj) : 
    mj=len(cg[ij]) 
    for im in range(mj) :
      for ia in range(na) :
        for ib in range(nb) :
          x=cg[ij][im][ia][ib]; ok= not iszero(x)
          j2=mj-1; l2=na-1; s2=nb-1; j32=j2-2*im; l32=l2-2*ia; s32=s2-2*ib
          if ok : s+= label(l2,l32,s2,s32,j2,j32,x)
  s+='\\end{array}\n'; return s

def bugcheck(cgm,mxa,mxb) :
  def subcheck(mx,my) :
    na=len(mx); nb=len(my[0]); nc=len(my[0][0]); bugs=0
    for i in range(na) :
      h=0
      for j in range(nb) :
        for k in range(nc) :
          x=reformat(mx[i][h]); y=reformat(my[i][j][k]); z=lincomb(x,y,1,-1)
          h+=1; bugs+= (0 if iszero(z) else 1) 
    return bugs
  n=len(cgm); ok=(n==2)
  na=len(cgm[0]); nb=len(cgm[1]); ma=cgm[0]; mb=cgm[1]
  ok= ok and (na==len(mxa)) and (nb==len(mxb))
  if not ok : print('Dimensions wrong, n,na,nb='+str([n,na,nb]))
  wa=len(ma[0])*len(ma[0][0]); va=len(mxa[0]); ok= ok and (va==wa)
  wb=len(mb[0])*len(mb[0][0]); vb=len(mxb[0]); ok= ok and (vb==wb)
  if not ok : print('Also wrong, wa,wb,va,vb='+str([wa,wb,va,vb]))
  bugs= subcheck(mxa,ma)+subcheck(mxb,mb)
  print('Nb of mismatch: '+str(bugs))

def cgtest(n2=2) : # compare cgm to mxa,mxb, for spin-(1/2) 2nd factors. ok.
  # n2=2L added to 2S=1, n2=1,2,3...
  cgm= clebschgordan(n2,1)
  mxa,s,t =testcleb([n2+1,2],[n2,2])
  mxb,s,t =testcleb([n2-1,2],[n2,2])
  bugcheck(cgm,mxa,mxb) # special check for spin-(1/2), ok no mismatch.
  #clebschgordan(1,2); print('')

def cg_getrow(cg,i) : # extract full row of CG matrix
  na=len(cg[0][0]); nb=len(cg[0][0][0]); n=na*nb; z=[[]]*n 
  j=0; h=i; ia=0;ib=0
  while h>=len(cg[j]) : h-=len(cg[j]);j+=1 # row is cg[j][h]
  for k in range(n) : 
    z[k]=cg[j][h][ia][ib]; ib+=1; (ia,ib)= (ia+1,0) if (ib>=nb) else (ia,ib) 
  return z

def cg_getcol(cg,i) :
  na=len(cg[0][0]); nb=len(cg[0][0][0]); n=na*nb; z=[[]]*n 
  ia= div(i,nb); ib=i-ia*nb; j=0;h=0
  for k in range(n) : 
    z[k]=cg[j][h][ia][ib]; h+=1; (j,h)= (j+1,0) if (h>=len(cg[j])) else (j,h)
  return z

def orthotest(cg,mode,noisy) : # is cg orthogon.matrix? mode=0 rows =1 columns.
  na=len(cg[0][0]); nb=len(cg[0][0][0]); m=na*nb; x=[[]]*m 
  hint=['Rows','Cols']; nerr=0;pairs=0 
  for i in range(m) : x[i]= cg_getrow(cg,i) if(mode==0) else cg_getcol(cg,i)
  for i in range(m) :
    for j in range(i,m) :
      sp=scalprod(x[i],x[j]); pairs+=1
      if i==j : sp=lincomb(sp,reformat(1),1,-1)
      sok= iszero(sp); nerr= (nerr+1) if (not sok) else nerr
  if noisy : print('Orthonorm '+hint[mode]+' errs: '+str(nerr)+'/'+str(pairs))
  return nerr

def mirrortest(cg,noisy) :
  nonz=0 # test symmetry : <l,s;-l3,-s3|j,-j3>=(-1)^{l+s-j}<l,s;l3,s3|j,j3>
  nj=len(cg); na=len(cg[0][0]); nb=len(cg[0][0][0]); n=na*nb; nerr=0;pairs=0
  for ij in range(nj) : 
    mj=len(cg[ij]) 
    for im in range(mj) :
      for ia in range(na) :
        for ib in range(nb) :
          x=cg[ij][im][ia][ib]; y=cg[ij][mj-1-im][na-1-ia][nb-1-ib]
          j2=mj-1; a2=na-1; b2=nb-1; z=a2+b2-j2; nonz+=(0 if iszero(x) else 1)
          if (z%2)==1 : print('Bug '+str([j2,a2,b2])); exit()
          z=div(abs(z),2); sign=1 if ((z%2)==0) else (-1); pairs+=1 
          w=lincomb(reformat(x),reformat(y),1,-sign); ok= iszero(w)
          if not ok : nerr +=1; print(' Err'+str([ij,im,ia,ib,x,y]))
  if noisy :
    print('Mirror errs/pairs/nonz '+str(nerr)+'/'+str(pairs)+'/'+str(nonz))
  return nerr
  
def racah(j2,l2,s2,l32) : # racah formula, alternative j sequence start
  def facto(n) : # n!
    y=1; i=n
    while i>1 : y*=i; i-=1
    return y
  def fac2(n) : # (n/2)!
    if (n<=2) : return 1
    if (n%2)!=0 : print('fac2 wrong argument '+str(n)+'. exit.'); exit()
    i=div(n,2); return facto(i)
  s32=j2-l32; pp=facto(j2+1)*fac2(l2+s2-j2)
  pq= fac2(l2+s2+j2+2)*fac2(j2+l2-s2)*fac2(j2+s2-l2)
  qp= fac2(l2+l32)*fac2(s2+s32); qq=fac2(l2-l32)*fac2(s2-s32)
  x=div(l2-l32,2); sign= 1 if (x%2)==0 else -1
  z=sqroot([pp*qp,pq*qq]); z[0]=sign*z[0]; return z 

def racahtest(cg,noisy) : # test all cg[j][0] lines against the formula
  n=len(cg); j2=len(cg[0])-1; count=0; errs=0
  for i in range(n) :
    row=cg[i][0]; nl=len(row); ns=len(row[0]); l2=nl-1; l32=l2
    for k in range(nl) :
      s2= ns-1; s32=s2
      for h in range(ns) : 
        ok= (l32+s32)==j2; z=row[k][h]
        if ok : y=racah(j2,l2,s2,l32); count+=1; w=lincomb(y,reformat(z),1,-1)
        if ok : errs = errs if iszero(w) else (errs+1)
        s32 -=2
      l32 -=2
    j2 -=2
  if noisy : print('Racahtest errors='+str(errs)+'/'+str(count))
  return errs

def cgexample(l2,s2,noisy) :
  cgm= clebschgordan(l2,s2,noisy); errs=0; s=''
  errs += orthotest(cgm,0,noisy)
  errs += orthotest(cgm,1,noisy)
  errs += mirrortest(cgm,noisy)
  errs += racahtest(cgm,noisy)
  if errs==0 : s=wikidump(cgm); print('<math>'+s+'</math>')
  return errs,s

def fulltest() :
  tit=['','(1/2)','1','(3/2)','2']
  for i in range(1,5) :
    for k in range(1,5) :
      title="\n\'\'\'Clebsch-Gordan L="+tit[i]+" S="+tit[k]+"\'\'\'<br>" 
      print(title); cgexample(i,k,False)
  #errs=cgexample(6,6,True) 

fulltest()