Wigner 3j, 6j, 9j, Clebsch-Gordan, Racah and Gaunt coefficients

[jyr]

Hi,

I have coded some routines that calculate Wigner 3j, 6j, 9j, Clebsch-
Gordan, Racah and Gaunt coefficients (integrals over 3 spherical
harmonics) exactly. It is all in a single python file with doc tests.
If I am not mistaken then Sage currently does not have such
functionality and I think it would be very useful for atomic,
molecular, nuclear physicists as well as quantum chemists out there.

If there is interest I could post it here or get a trac account and
open a ticket for it for review.

regards,

Jens

[mabshoff]

On May 3, 3:07 am, jyr <jyr2...@googlemail.com> wrote:
> Hi,

Hi Jens,

Yes, do both.

> regards,
>
> Jens

Cheers,

Michael

[jyr]

Ok, here it is:

-----------------------------------------------------------------------------------------------
r"""
Calculate Wigner 3j, 6j, 9j, Clebsch-Gordan, Racah and Gaunt
coefficients

Collection of functions for calculating er 3j, 6j, 9j, Clebsch-Gordan,
Racah as well as Gaunt coefficients exactly, all evaluating to a
rational number times the square root of a rational number [Rasch03].

Please see the description of the individual functions for further
details and examples.


REFERENCES:

- [Rasch03] 'Efficient Storage Scheme for Precalculated Wigner 3j, 6j
and Gaunt Coefficients', J. Rasch and A. C. H. Yu, SIAM
J. Sci. Comput. Volume 25, Issue 4, pp. 1416-1428 (2003)


AUTHORS:

- Jens Rasch (2009-03-24): initial version for Sage

"""

#***********************************************************************
# Copyright (C) 2008 Jens Rasch <jyr...@gmail.com>
#
# Distributed under the terms of the GNU General Public License (GPL)
# http://www.gnu.org/licenses/
#***********************************************************************

from sage.all import *
# from sage.calculus.calculus import sqrt,factorial
# from sage.functions.constants import pi



# This list of precomputed factorials is needed to massively
# accelerate consequetive calculations of the various coefficients
_Factlist=[1]

def _calc_factlist(nn):
if nn>=len(_Factlist):
for ii in range(len(_Factlist),nn+1):
_Factlist.append(_Factlist[ii-1]*ii)
#_Factlist.append(factorial(ii))
#return _Factlist


def test_calc_factlist(nn):
r"""
Function calculates a list of precomputed factorials in order to
massively accelerate consequetive calculations of the various
coefficients.

INPUT:

- nn Highest factorial to be computed

OUTPUT:

integer list of factorials


EXAMPLES:

Calculate list of factorials:

sage: test_calc_factlist(10)
[1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800]
"""
_calc_factlist(nn)
return _Factlist



def Wigner3j(j_1,j_2,j_3,m_1,m_2,m_3,prec=None):
r"""
Calculate the Wigner 3j symbol

\left({j_1 \atop m_1} {j_2 \atop m_2} {j_3 \atop m_3} \right)


NOTES:

The Wigner 3j symbol obeys the following symmetry rules:

- invariant under any permutation of the columns (with the
exception of a sign change where $J:=j_1+j_2+j_3$):
\begin{eqnarray}
\left({j_1 \atop m_1} {j_2 \atop m_2} {j_3 \atop m_3}\right)
&=&
\left({j_3 \atop m_3} {j_1 \atop m_1} {j_2 \atop m_2}\right)
=\left({j_2 \atop m_2} {j_3 \atop m_3} {j_1 \atop m_1}\right)
\qquad \mbox{cyclic}
&=&
(-)^{J}\left({j_3\atop m_3} {j_2\atop m_2} {j_1\atop
m_1}\right)
=(-)^{J}\left({j_1\atop m_1} {j_3\atop m_3} {j_2\atop
m_2}\right)
=(-)^{J}\left({j_2\atop m_2} {j_1\atop m_1} {j_3\atop
m_3}\right)
\qquad\mbox{anti-cyclic}
\end{eqnarray}

- invariant under space inflection, i. e.
\begin{eqnarray}
\left({{j_1}\atop{m_1}} {{j_2}\atop{m_2}} {{j_3}\atop{m_3}}
\right)
=
(-)^{J}
\left({j_1 \atop -m_1} {j_2 \atop -m_2}{j_3 \atop -m_3}\right)
\end{eqnarray}

- symmetric with respect to the 72 additional symmetries based on
the work by [Regge58]

- zero for $j_1$, $j_2$, $j_3$ not fulfilling triangle relation

- zero for ${m_1}+{m_2}+{m_3}!= 0$

- zero for violating any one of the conditions
$j_1\ge|m_1|$, $j_2\ge|m_2|$, $j_3\ge|m_3|$


ALGORITHM:

This function uses the algorithm of [Edmonds74] to calculate the
value of the 3j symbol exactly. Note that the formula contains
alternating sums over large factorials and is therefore unsuitable
for finite precision arithmetic and only useful for a computer
algebra system [Rasch03].



INPUT:

j_1 - integer or half integer
j_2 - integer or half integer
j_3 - integer or half integer
m_1 - integer or half integer
m_2 - integer or half integer
m_3 - integer or half integer
prec - precission, default: None. Providing a precission can
drastically spead up the calculation.


OUTPUT:

The result will be a rational number times the square root of a
rational number, unless a precission is given.


EXAMPLES:

A couple of examples:

sage: Wigner3j(2,6,4,0,0,0)
sqrt(5)/sqrt(143)

sage: Wigner3j(2,6,4,0,0,1)
0

sage: Wigner3j(0.5,0.5,1,0.5,-0.5,0)
1/sqrt(6)

sage: Wigner3j(40,100,60,-10,60,-50)
95608*sqrt(21082735836735314343364163310)/(18702538494885*sqrt
(220491455010479533763))

sage: Wigner3j(2500,2500,5000,2488,2400,-4888 ,prec=64)
7.60424456883448589e-12


It is an error to have arguments that are not integer or half
integer values:

sage: Wigner3j(2.1,6,4,0,0,0)
Traceback (most recent call last):
...
ValueError: j values must be integer or half integer

sage: Wigner3j(2,6,4,1,0,-1.1)
Traceback (most recent call last):
...
ValueError: m values must be integer or half integer



REFERENCES:

- [Regge58] 'Symmetry Properties of Clebsch-Gordan Coefficients',
T. Regge, Nuovo Cimento, Volume 10, pp. 544 (1958)

- [Edmonds74] 'Angular Momentum in Quantum Mechanics',
A. R. Edmonds, Princeton University Press (1974)

- [Rasch03] 'Efficient Storage Scheme for Precalculated Wigner 3j,
6j and Gaunt Coefficients', J. Rasch and A. C. H. Yu, SIAM
J. Sci. Comput. Volume 25, Issue 4, pp. 1416-1428 (2003)


AUTHORS:

- Jens Rasch (2009-03-24): initial version

"""

if int(j_1*2)!=j_1*2 or int(j_2*2)!=j_2*2 or int(j_3*2)!=j_3*2:
raise ValueError("j values must be integer or half integer")
#return 0
if int(m_1*2)!=m_1*2 or int(m_2*2)!=m_2*2 or int(m_3*2)!=m_3*2:
raise ValueError("m values must be integer or half integer")
#return 0
if (m_1+m_2+m_3<>0):
return 0
prefid=Integer((-1)**(int(j_1-j_2-m_3)))
m_3=-m_3
a1=j_1+j_2-j_3
if (a1<0):
return 0
a2=j_1-j_2+j_3
if (a2<0):
return 0
a3=-j_1+j_2+j_3
if (a3<0):
return 0
if (abs(m_1)>j_1) or (abs(m_2)>j_2) or (abs(m_3)>j_3):
return 0

maxfact=max(j_1+j_2+j_3+1,j_1+abs(m_1),j_2+abs(m_2),j_3+abs(m_3))
_calc_factlist(maxfact)

#try:
argsqrt=Integer(_Factlist[int(j_1+j_2-j_3)]\
*_Factlist[int(j_1-j_2+j_3)]\
*_Factlist[int(-j_1+j_2+j_3)]\
*_Factlist[int(j_1-m_1)]*_Factlist[int(j_1+m_1)]\
*_Factlist[int(j_2-m_2)]*_Factlist[int(j_2+m_2)]\
*_Factlist[int(j_3-m_3)]*_Factlist[j_3+m_3])\
/_Factlist[int(j_1+j_2+j_3+1)]\
#except:
# return 0

ressqrt=sqrt(argsqrt,prec)
if type(ressqrt)==sage.rings.complex_number.ComplexNumber:
ressqrt=ressqrt.real()

imin=max(-j_3+j_1+m_2,-j_3+j_2-m_1,0)
imax=min(j_2+m_2,j_1-m_1,j_1+j_2-j_3)
sumres=0
for ii in range(imin,imax+1):
den=_Factlist[ii]*_Factlist[int(ii+j_3-j_1-m_2)]\
*_Factlist[int(j_2+m_2-ii)]*_Factlist[int(j_1-ii-m_1)]\
*_Factlist[int(ii+j_3-j_2+m_1)]\
*_Factlist[int(j_1+j_2-j_3-ii)]
sumres=sumres+Integer((-1)**ii)/den

res=ressqrt*sumres*prefid
return res


def ClebschGordan(j_1,j_2,j_3,m_1,m_2,m_3,prec=None):
r"""
Calculates the Clebsch-Gordan coefficient

$\left< j_1 m_1 \; j_2 m_2 | j_3 m_3 \right>$


NOTES:

The Clebsch-Gordan coefficient will be evaluated via its relation
to Wigner 3j symbols:

\begin{eqnarray}
\left< j_1 m_1 \; j_2 m_2 | j_3 m_3 \right>=
(-1)^(j_1-j_2+m_3) \sqrt(2j_3+1)
\left({j_1 \atop m_1} {j_2 \atop m_2} {j_3 \atop -m_3}\right)
\end{eqnarray}

See also the documentation on Wigner 3j symbols which exhibit much
higher symmetry releations that the Clebsch-Gordan coefficient.


INPUT:

j_1 - integer or half integer
j_2 - integer or half integer
j_3 - integer or half integer
m_1 - integer or half integer
m_2 - integer or half integer
m_3 - integer or half integer
prec - precission, default: None. Providing a precission can
drastically spead up the calculation.


OUTPUT:

The result will be a rational number times the square root of a
rational number, unless a precission is given.


EXAMPLES:

A couple of examples:

sage: ClebschGordan(3/2,1/2,2, 3/2,1/2,2)
1

sage: ClebschGordan(1.5,0.5,1, 1.5,-0.5,1)
sqrt(3)/2

sage: ClebschGordan(3/2,1/2,1, -1/2,1/2,0)
-sqrt(3)/sqrt(6)


REFERENCES:

- [Edmonds74] 'Angular Momentum in Quantum Mechanics',
A. R. Edmonds, Princeton University Press (1974)

- [Rasch03] 'Efficient Storage Scheme for Precalculated Wigner 3j,
6j and Gaunt Coefficients', J. Rasch and A. C. H. Yu, SIAM
J. Sci. Comput. Volume 25, Issue 4, pp. 1416-1428 (2003)


AUTHORS:

- Jens Rasch (2009-03-24): initial version

"""

res=(-1)**(int(j_1-j_2+m_3))*sqrt(2*j_3+1)\
*Wigner3j(j_1,j_2,j_3,m_1,m_2,-m_3,prec)
return res





def _bigDeltacoeff(aa,bb,cc,prec=None):
r"""
Calculates the Delta coefficient of the 3 angular momenta for
Racah symbols. Also checks that the differences are of integer
value.

INPUT:

- aa - first angular momentum, integer or half integer
- bb - second angular momentum, integer or half integer
- cc - third angular momentum, integer or half integer
- prec - precission of the sqrt() calculation

OUTPUT:

double - Value of the Delta coefficient

"""

if(int(aa+bb-cc)!=(aa+bb-cc)):
raise ValueError("j values must be integer or half integer and
fullfil the triangle relation")
#return 0
if(int(aa+cc-bb)!=(aa+cc-bb)):
raise ValueError("j values must be integer or half integer and
fullfil the triangle relation")
#return 0
if(int(bb+cc-aa)!=(bb+cc-aa)):
raise ValueError("j values must be integer or half integer and
fullfil the triangle relation")
#return 0
if(aa+bb-cc)<0:
return 0
if(aa+cc-bb)<0:
return 0
if(bb+cc-aa)<0:
return 0

maxfact=max(aa+bb-cc,aa+cc-bb,bb+cc-aa,aa+bb+cc+1)
_calc_factlist(maxfact)

argsqrt=Integer(_Factlist[int(aa+bb-cc)]*_Factlist[int(aa+cc-bb)]\
*_Factlist[int(bb+cc-aa)])\
/Integer(_Factlist[int(aa+bb+cc+1)])

ressqrt=sqrt(argsqrt,prec)
if type(ressqrt)==sage.rings.complex_number.ComplexNumber:
res=ressqrt.real()
else:
res=ressqrt
return res


def test_bigDeltacoeff(aa,bb,cc,prec=None):
r"""
Test for the Delta coefficient of the 3 angular momenta for
Racah symbols. Also checks that the differences are of integer
value.

INPUT:

- aa - first angular momentum, integer or half integer
- bb - second angular momentum, integer or half integer
- cc - third angular momentum, integer or half integer
- prec - precission of the sqrt() calculation

OUTPUT:

double - Value of the Delta coefficient

EXAMPLES:

Simple examples:

sage: test_bigDeltacoeff(1,1,1)
1/(2*sqrt(6))

"""

return _bigDeltacoeff(aa,bb,cc,prec)




def Racah(aa,bb,cc,dd,ee,ff,prec=None):
r"""
Calculate the Racah symbol

W(a,b,c,d;e,f)

NOTES:

The Racah symbol is related to the Wigner 6j symbol:
\begin{eqnarray}
\left({j_1 \atop j_4} {j_2 \atop j_5} {j_3 \atop j_6} \right)
=
(-1)^{j_1+j_2+j_4+j_5} W(j_1,j_2,j_5,j_4,j_3,j_6)
\end{eqnarray}

Please see the 6j symbol for its much richer symmetries and for
additional properties.


ALGORITHM:

This function uses the algorithm of [Edmonds74] to calculate the
value of the 6j symbol exactly. Note that the formula contains
alternating sums over large factorials and is therefore unsuitable
for finite precision arithmetic and only useful for a computer
algebra system [Rasch03].


INPUT:

a - integer or half integer
b - integer or half integer
c - integer or half integer
d - integer or half integer
e - integer or half integer
f - integer or half integer
prec - precission, default: None. Providing a precission can
drastically spead up the calculation.


OUTPUT:

The result will be a rational number times the square root of a
rational number, unless a precission is given.


EXAMPLES:

A couple of examples and testcases:

sage: Racah(3,3,3,3,3,3)
-1/14


REFERENCES:

- [Edmonds74] 'Angular Momentum in Quantum Mechanics',
A. R. Edmonds, Princeton University Press (1974)

- [Rasch03] 'Efficient Storage Scheme for Precalculated Wigner 3j,
6j and Gaunt Coefficients', J. Rasch and A. C. H. Yu, SIAM
J. Sci. Comput. Volume 25, Issue 4, pp. 1416-1428 (2003)


AUTHORS:

- Jens Rasch (2009-03-24): initial version

"""

prefac=_bigDeltacoeff(aa,bb,ee,prec)*_bigDeltacoeff(cc,dd,ee,prec)
\
*_bigDeltacoeff(aa,cc,ff,prec)*_bigDeltacoeff
(bb,dd,ff,prec)
if prefac==0:
return 0
imin=max(aa+bb+ee,cc+dd+ee,aa+cc+ff,bb+dd+ff)
imax=min(aa+bb+cc+dd,aa+dd+ee+ff,bb+cc+ee+ff)

maxfact=max(imax+1,aa+bb+cc+dd,aa+dd+ee+ff,bb+cc+ee+ff)
_calc_factlist(maxfact)

sumres=0
for kk in range(imin,imax+1):
den=_Factlist[int(kk-aa-bb-ee)]*_Factlist[int(kk-cc-dd-ee)]\
*_Factlist[int(kk-aa-cc-ff)]\
*_Factlist[int(kk-bb-dd-ff)]*_Factlist[int(aa+bb+cc+dd-kk)]
\
*_Factlist[int(aa+dd+ee+ff-kk)]\
*_Factlist[int(bb+cc+ee+ff-kk)]
sumres=sumres+Integer((-1)**kk*_Factlist[kk+1])/den

res=prefac*sumres*(-1)**(int(aa+bb+cc+dd))
return res





def Wigner6j(j_1,j_2,j_3,j_4,j_5,j_6,prec=None):
r"""
Calculate the Wigner 6j symbol

\left({j_1 \atop j_4} {j_2 \atop j_5} {j_3 \atop j_6} \right)


NOTES:

The Wigner 6j symbol is related to the Racah symbol but exhibits
more symmetries as detailed below.
\begin{eqnarray}
\left({j_1 \atop j_4} {j_2 \atop j_5} {j_3 \atop j_6} \right)
=
(-1)^{j_1+j_2+j_4+j_5} W(j_1,j_2,j_5,j_4,j_3,j_6)
\end{eqnarray}

The Wigner 6j symbol obeys the following symmetry rules:

- Wigner $6j$ symbols are left invariant under any permutation of
the columns:
\begin{eqnarray}
\left({j_1 \atop j_4} {j_2 \atop j_5} {j_3 \atop j_6} \right)
&=&
\left({j_3 \atop j_6} {j_1 \atop j_4} {j_2 \atop j_5} \right)
=\left({j_2 \atop j_5} {j_3 \atop j_6} {j_2 \atop j_4} \right)
\qquad \mbox{cyclic} \\
&=&
\left({j_3 \atop j_6} {j_2 \atop j_5} {j_1 \atop j_4} \right)
=\left({j_1 \atop j_4} {j_3 \atop j_6} {j_2 \atop j_5} \right)
=\left({j_2 \atop j_5} {j_1 \atop j_4} {j_3 \atop j_6} \right)
\qquad \mbox{anti-cyclic}
\end{eqnarray}

- They are invariant under the exchange of the upper and lower
arguments in each of any two columns, i. e.
\begin{eqnarray}
\left({j_1 \atop j_4} {j_2 \atop j_5} {j_3 \atop j_6} \right)
=
\left({j_1 \atop j_4} {j_5 \atop j_2} {j_6 \atop j_3} \right)
=
\left({j_4 \atop j_1} {j_2 \atop j_5} {j_6 \atop j_3} \right)
=
\left({j_4 \atop j_1} {j_5 \atop j_2} {j_3 \atop j_6} \right)
\end{eqnarray}

- additional 6 symmetries [Regge59] giving rise to 144 symmetries
in total

- only non-zero if any triple of $j$'s fulfil a triangle relation


ALGORITHM:

This function uses the algorithm of [Edmonds74] to calculate the
value of the 6j symbol exactly. Note that the formula contains
alternating sums over large factorials and is therefore unsuitable
for finite precision arithmetic and only useful for a computer
algebra system [Rasch03].


INPUT:

j_1 - integer or half integer
j_2 - integer or half integer
j_3 - integer or half integer
j_4 - integer or half integer
j_5 - integer or half integer
j_6 - integer or half integer
prec - precission, default: None. Providing a precission can
drastically spead up the calculation.


OUTPUT:

The result will be a rational number times the square root of a
rational number, unless a precission is given.


EXAMPLES:

A couple of examples and testcases:

sage: Wigner6j(3,3,3,3,3,3)
-1/14

sage: Wigner6j(5,5,5,5,5,5)
1/52

sage: Wigner6j(6,6,6,6,6,6)
309/10868

sage: Wigner6j(8,8,8,8,8,8)
-12219/965770

sage: Wigner6j(30,30,30,30,30,30)
36082186869033479581/87954851694828981714124

sage: Wigner6j(0.5,0.5,1,0.5,0.5,1)
1/6

sage: Wigner6j(200,200,200,200,200,200, prec=1000)*1.0
0.000155903212413242


It is an error to have arguments that are not integer or half
integer values or do not fulfill the triangle relation:

sage: Wigner6j(2.5,2.5,2.5,2.5,2.5,2.5)
Traceback (most recent call last):
...
ValueError: j values must be integer or half integer and
fullfil the triangle relation

sage: Wigner6j(0.5,0.5,1.1,0.5,0.5,1.1)
Traceback (most recent call last):
...
ValueError: j values must be integer or half integer and
fullfil the triangle relation


REFERENCES:

- [Regge59] 'Symmetry Properties of Racah Coefficients',
T. Regge, Nuovo Cimento, Volume 11, pp. 116 (1959)

- [Edmonds74] 'Angular Momentum in Quantum Mechanics',
A. R. Edmonds, Princeton University Press (1974)

- [Rasch03] 'Efficient Storage Scheme for Precalculated Wigner 3j,
6j and Gaunt Coefficients', J. Rasch and A. C. H. Yu, SIAM
J. Sci. Comput. Volume 25, Issue 4, pp. 1416-1428 (2003)

"""

res=(-1)**(int(j_1+j_2+j_4+j_5))*Racah
(j_1,j_2,j_5,j_4,j_3,j_6,prec)
return res




def Wigner9j(j_1,j_2,j_3,j_4,j_5,j_6,j_7,j_8,j_9,prec=None):
r"""
Calculate the Wigner 9j symbol

\left\{
\begin{matrix}
{j_1} & {j_2} & {j_3} \\
{j_4} & {j_5} & {j_6} \\
{j_7} & {j_8} & {j_9}
\end{matrix}
\right\}


ALGORITHM:

This function uses the algorithm of [Edmonds74] to calculate the
value of the 3j symbol exactly. Note that the formula contains
alternating sums over large factorials and is therefore unsuitable
for finite precision arithmetic and only useful for a computer
algebra system [Rasch03].


INPUT:

j_1 - integer or half integer
j_2 - integer or half integer
j_3 - integer or half integer
j_4 - integer or half integer
j_5 - integer or half integer
j_6 - integer or half integer
j_7 - integer or half integer
j_8 - integer or half integer
j_9 - integer or half integer
prec - precission, default: None. Providing a precission can
drastically spead up the calculation.


OUTPUT:

The result will be a rational number times the square root of a
rational number, unless a precission is given.


EXAMPLES:

A couple of examples and testcases, note that for speed reasons a
precission is given:

sage: Wigner9j(1,1,1, 1,1,1, 1,1,0 ,prec=64) # ==1/18
0.0555555555555555555

sage: Wigner9j(1,1,1, 1,1,1, 1,1,1)
0

sage: Wigner9j(1,1,1, 1,1,1, 1,1,2 ,prec=64) # ==1/18
0.0555555555555555556

sage: Wigner9j(1,2,1, 2,2,2, 1,2,1 ,prec=64) # ==-1/150
-0.00666666666666666667

sage: Wigner9j(3,3,2, 2,2,2, 3,3,2 ,prec=64) # ==157/14700
0.0106802721088435374

sage: Wigner9j(3,3,2, 3,3,2, 3,3,2 ,prec=64) # ==3221*sqrt(70)/
(246960*sqrt(105)) - 365/(3528*sqrt(70)*sqrt(105))
0.00944247746651111739

sage: Wigner9j(3,3,1, 3.5,3.5,2, 3.5,3.5,1 ,prec=64) #
==3221*sqrt(70)/(246960*sqrt(105)) - 365/(3528*sqrt(70)*sqrt(105))
0.0110216678544351364

sage: Wigner9j(100,80,50, 50,100,70, 60,50,100 ,prec=1000)*1.0
1.05597798065761e-7

sage: Wigner9j(30,30,10, 30.5,30.5,20, 30.5,30.5,10 ,prec=1000)
*1.0 # ==(80944680186359968990/95103769817469)*sqrt
(1/682288158959699477295)
0.0000325841699408828

sage: Wigner9j(64,62.5,114.5, 61.5,61,112.5, 113.5,110.5,60,
prec=1000)*1.0
-3.41407910055520e-39

sage: Wigner9j(15,15,15, 15,3,15, 15,18,10, prec=1000)*1.0
-0.0000778324615309539

sage: Wigner9j(1.5,1,1.5, 1,1,1, 1.5,1,1.5)
0


It is an error to have arguments that are not integer or half
integer values or do not fulfill the triangle relation:

sage: Wigner9j(0.5,0.5,0.5, 0.5,0.5,0.5, 0.5,0.5,0.5,prec=64)
Traceback (most recent call last):
...
ValueError: j values must be integer or half integer and
fullfil the triangle relation

sage: Wigner9j(1,1,1, 0.5,1,1.5, 0.5,1,2.5,prec=64)
Traceback (most recent call last):
...
ValueError: j values must be integer or half integer and
fullfil the triangle relation


REFERENCES:

- [Edmonds74] 'Angular Momentum in Quantum Mechanics',
A. R. Edmonds, Princeton University Press (1974)

- [Rasch03] 'Efficient Storage Scheme for Precalculated Wigner 3j,
6j and Gaunt Coefficients', J. Rasch and A. C. H. Yu, SIAM
J. Sci. Comput. Volume 25, Issue 4, pp. 1416-1428 (2003)

"""

imin=0
imax=min(j_1+j_9,j_2+j_6,j_4+j_8)

sumres=0
for kk in range(imin,imax+1):
sumres=sumres+(2*kk+1)*Racah(j_1,j_2,j_9,j_6,j_3,kk,prec)\
*Racah(j_4,j_6,j_8,j_2,j_5,kk,prec)\
*Racah(j_1,j_4,j_9,j_8,j_7,kk,prec)
return sumres




def Gaunt(l_1,l_2,l_3,m_1,m_2,m_3,prec=None):
r"""
Calculate the Gaunt coefficient which is defined as the integral
over three spherical harmonics:

\begin{eqnarray}
Y{j_1 \atop m_1} {j_2 \atop m_2} {j_3 \atop m_3}
&=&
\int Y_{l_1,m_1}(\Omega)
Y_{l_2,m_2}(\Omega) Y_{l_3,m_3}(\Omega)\; d\Omega \\
&=&
\sqrt{\frac{(2l_1+1)(2l_2+1)(2l_3+1)}{4\pi}}
\left({l_1 \atop 0 } {l_2 \atop 0 } {l_3 \atop 0} \right)
\left({l_1 \atop m_1} {l_2 \atop m_2} {l_3 \atop m_3} \right)
\end{eqnarray}


NOTES:

The Gaunt coefficient obeys the following symmetry rules:

- invariant under any permutation of the columns
\begin{eqnarray}
Y{j_1 \atop m_1} {j_2 \atop m_2} {j_3 \atop m_3}
&=&
Y{j_3 \atop m_3} {j_1 \atop m_1} {j_2 \atop m_2}
=Y{j_2 \atop m_2} {j_3 \atop m_3} {j_1 \atop m_1}
\qquad \mbox{cyclic} \\
&=&
Y{j_3 \atop m_3} {j_2 \atop m_2} {j_1 \atop m_1}
=Y{j_1 \atop m_1} {j_3 \atop m_3} {j_2 \atop m_2}
=Y{j_2 \atop m_2} {j_1 \atop m_1} {j_3 \atop m_3}
\qquad\mbox{anti-cyclic}
\end{eqnarray}

- invariant under space inflection, i.e.
\begin{eqnarray}
Y{j_1 \atop m_1} {j_2 \atop m_2} {j_3 \atop m_3}
=
Y{j_1 \atop -m_1} {j_2 \atop -m_2} {j_3 \atop -m_3}
\end{eqnarray}

- symmetric with respect to the 72 Regge symmetries as inherited
for the $3j$ symbols [Regge58]

- zero for $l_1$, $l_2$, $l_3$ not fulfilling triangle relation

- zero for violating any one of the conditions: $l_1\ge|m_1|$,
$l_2\ge|m_2|$, $l_3\ge|m_3|

- non-zero only for an even sum of the $l_i$, i. e.
$J=l_1+l_2+l_3=2n$ for $n$ in $\mathbf{N}$


ALGORITHM:

This function uses the algorithm of [Liberatodebrito82] to
calculate the value of the Gaunt coefficient exactly. Note that
the formula contains alternating sums over large factorials and is
therefore unsuitable for finite precision arithmetic and only
useful for a computer algebra system [Rasch03].


INPUT:

j_1 - integer
j_2 - integer
j_3 - integer
m_1 - integer
m_2 - integer
m_3 - integer
prec - precission, default: None. Providing a precission can
drastically spead up the calculation.


OUTPUT:

The result will be a rational number times the square root of a
rational number, unless a precission is given.


EXAMPLES:

sage: Gaunt(1,0,1,1,0,-1)
-1/(2*sqrt(pi))

sage: Gaunt(1,0,1,1,0,0)
0

sage: Gaunt(29,29,34,10,-5,-5)
1821867940156*sqrt(22134)/(215552371055153321*sqrt(pi))

sage: Gaunt(20,20,40,1,-1,0)
28384503878959800/(74029560764440771*sqrt(pi))

sage: Gaunt(12,15,5,2,3,-5)
91*sqrt(36890)/(124062*sqrt(pi))

sage: Gaunt(10,10,12,9,3,-12)
-98*sqrt(6279)/(62031*sqrt(pi))

sage: Gaunt(1000,1000,1200,9,3,-12).n(64)
0.00689500421922113448


It is an error to use non-integer values for $l$ and $m$:

sage: Gaunt(1.2,0,1.2,0,0,0)
Traceback (most recent call last):
...
ValueError: l values must be integer

sage: Gaunt(1,0,1,1.1,0,-1.1)
Traceback (most recent call last):
...
ValueError: m values must be integer



REFERENCES:

- [Regge58] 'Symmetry Properties of Clebsch-Gordan Coefficients',
T. Regge, Nuovo Cimento, Volume 10, pp. 544 (1958)

- [Liberatodebrito82] 'FORTRAN program for the integral of three
spherical harmonics', A. Liberato de Brito,
Comput. Phys. Commun., Volume 25, pp. 81-85 (1982)

- [Rasch03] 'Efficient Storage Scheme for Precalculated Wigner 3j,
6j and Gaunt Coefficients', J. Rasch and A. C. H. Yu, SIAM
J. Sci. Comput. Volume 25, Issue 4, pp. 1416-1428 (2003)


AUTHORS:

- Jens Rasch (2009-03-24): initial version for Sage

"""

if int(l_1)!=l_1 or int(l_2)!=l_2 or int(l_3)!=l_3:
raise ValueError("l values must be integer")
#return 0
if int(m_1)!=m_1 or int(m_2)!=m_2 or int(m_3)!=m_3:
raise ValueError("m values must be integer")

bigL=(l_1+l_2+l_3)/2
a1=l_1+l_2-l_3
if (a1<0):
return 0
a2=l_1-l_2+l_3
if (a2<0):
return 0
a3=-l_1+l_2+l_3
if (a3<0):
return 0
if Mod(2*bigL,2)<>0:
# if int(int(2*bigL)/2)<>int(bigL):
return 0
if (m_1+m_2+m_3)<>0:
return 0
if (abs(m_1)>l_1) or (abs(m_2)>l_2) or (abs(m_3)>l_3):
return 0

imin=max(-l_3+l_1+m_2,-l_3+l_2-m_1,0)
imax=min(l_2+m_2,l_1-m_1,l_1+l_2-l_3)

maxfact=max(l_1+l_2+l_3+1,imax+1)
_calc_factlist(maxfact)

argsqrt=(2*l_1+1)*(2*l_2+1)*(2*l_3+1)*_Factlist[l_1-m_1]\
*_Factlist[l_1+m_1]*_Factlist[l_2-m_2]*_Factlist
[l_2+m_2]\
*_Factlist[l_3-m_3]*_Factlist[l_3+m_3]/(4*pi)
ressqrt=sqrt(argsqrt,prec)
if type(ressqrt)==sage.rings.complex_number.ComplexNumber:
ressqrt=ressqrt.real()
prefac=Integer(_Factlist[bigL]*_Factlist[l_2-l_1+l_3]\
*_Factlist[l_1-l_2+l_3]\
*_Factlist[l_1+l_2-l_3])/_Factlist[2*bigL+1]\
/(_Factlist[bigL-l_1]*_Factlist[bigL-l_2]*_Factlist[bigL-l_3])

sumres=0
for ii in range(imin,imax+1):
den=_Factlist[ii]*_Factlist[ii+l_3-l_1-m_2]\
*_Factlist[l_2+m_2-ii]*_Factlist[l_1-ii-m_1]\
*_Factlist[ii+l_3-l_2+m_1]*_Factlist[l_1+l_2-l_3-ii]
sumres=sumres+Integer((-1)**ii)/den

res=ressqrt*prefac*sumres*(-1)**(bigL+l_3+m_1-m_2)
return res

[Ondrej Certik]

On Sun, May 3, 2009 at 3:59 AM, jyr <jyr...@googlemail.com> wrote:
>
> Ok, here it is:
>
> -----------------------------------------------------------------------------------------------
> r"""
> Calculate Wigner 3j, 6j, 9j, Clebsch-Gordan, Racah and Gaunt
> coefficients
>
> Collection of functions for calculating er 3j, 6j, 9j, Clebsch-Gordan,
> Racah as well as Gaunt coefficients exactly, all evaluating to a
> rational number times the square root of a rational number [Rasch03].
>
> Please see the description of the individual functions for further
> details and examples.
>
>
> REFERENCES:
>
> - [Rasch03] 'Efficient Storage Scheme for Precalculated Wigner 3j, 6j
>  and Gaunt Coefficients', J. Rasch and A. C. H. Yu, SIAM
>  J. Sci. Comput. Volume 25, Issue 4, pp. 1416-1428 (2003)

Nice article. Is the C program, that you use in there available somewhere?

I especially like the magic square, eq. 2-10. Indeed, the sum in each
row or column is j1+j2+j3, it never occured to me the 3j symbols have
the same symmetries at the magic square.

Ondrej

[jyr]

Hi Ondrej,

No, I have not come round to publishing it yet in Comput. Phys.
Commun. or some similar journal.

The thought has occured to me that one could use the index functions
for the 3j, 6j, and Gaunt coefficients for a much simpler storage
scheme in python by using the index as a key for a dictionary of
stored symbols. I could then extend the above published routines with
an optional 'remember' keyword that would store them from one call to
the next.
However, I am not sure how much memory dictionaries consume in python
and whether this is efficient. At the end of the day for any realistic
calculation you are usually talking about millions of them.

regards,

Jens

[jyr]

Hi Michael,

Could you open a trac account for me please, so that I can add the
code to a ticket.

I tried to email you earlier, but the email probably died in some spam
filter.

regards,

Jens

[Dan Christensen]

jyr <jyr...@googlemail.com> writes:

> The thought has occured to me that one could use the index functions
> for the 3j, 6j, and Gaunt coefficients for a much simpler storage
> scheme in python by using the index as a key for a dictionary of
> stored symbols. I could then extend the above published routines with
> an optional 'remember' keyword that would store them from one call to
> the next.
> However, I am not sure how much memory dictionaries consume in python
> and whether this is efficient. At the end of the day for any realistic
> calculation you are usually talking about millions of them.

In my experience doing extremely large-scale computations using 10^10 or
more 6j symbols, it is faster to evaluate a 6j symbol than to look it up
in a hash table. On the other hand, if you repeatedly need sequences of
6j symbols with one argument varying, it does pay off to store arrays of
6j symbols in a hash table with the key being the 5 fixed arguments.
That is all in C using floating point, so it may be worthwhile using
a hash table if you want exact rational results.

By the way, I have a C library, written with Igor Khavkine, that
computes 3j, 6j, 10j and other networks using your choice of exact
rational or floating point arithmetic. It also does the q-deformed
calculations, either for fixed complex q (returning a complex) or for
symbolic q (returning an exact rational function of q). It is
optimized for speed and even has a python interface.

The bad news is that it is complicated to build, not well-documented,
and not really in a state for others to use.

So it would be great if Jens' library could be included in sage!

Dan

PS: I'm currently working on using sage to evaluate 3j and 6j symbols
for SU(n) for n > 2 (the above being the n=2 case). I'll report to the
list if I come up with something appropriate for inclusion in sage.

[Dan Christensen]

jyr <jyr...@googlemail.com> writes:

> def test_calc_factlist(nn):
> r"""
> Function calculates a list of precomputed factorials in order to
> massively accelerate consequetive calculations of the various

^^^^^^^^^^^^

Typo. And maybe say "future" instead of "consecutive"?

> EXAMPLES:
>
> Calculate list of factorials:
>
> sage: test_calc_factlist(10)
> [1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800]
> """
> _calc_factlist(nn)
> return _Factlist

Maybe "return _Factlist[:nn+1]" to ensure that the test passes even
if other tests run which create a larger _Factlist?

> This function uses the algorithm of [Edmonds74] to calculate the
> value of the 3j symbol exactly. Note that the formula contains
> alternating sums over large factorials and is therefore unsuitable
> for finite precision arithmetic and only useful for a computer
> algebra system [Rasch03].

Do you have a specific example illustrating what you mean by the last
sentence? I've done lots of calculations of theta symbols which contain
a similar sum. For these, if I use long double's in C, I retain good
precision even for fairly large arguments.

And I've found that floating point is also fine for 6j symbols.

> higher symmetry releations that the Clebsch-Gordan coefficient.

Typo ^^^^^^^^^^

> prec - precission, default: None. Providing a precission can

Only one s in "precision" (twice here and also in other places).

Probably worth running the file through a spell-checker.

Dan

[Dan Christensen]

Jens,

I now see that you've written an article on this topic and tested
exactly the case I've also tested: floating point 6j symbols in a
compiled language. Since your conclusions are different from mine, I'm
curious whether your storage system is faster or your 6j routine is
slower. Is your code available somewhere?

Dan

[Ondrej Certik]

Jens already replied above that it is not, until he publishes an
article. But I am certainly interested too in the C code, once the
article is done.

Ondrej

[jyr]

Hi Dan,

Thanks for reviewing the code. I have included your suggestions and I
have done a spell check on the file now and will submit the updated
version to trac once I get an account. (Michael?)


As a side remark concerning the prec argument, it is passed straight
into the sqrt() functions used in the algorithms and not used in any
other way. I am not much of a sage expert to know whether this is
sensible or whether there are better options available by now.


Concerning other algorithms for evaluating these symbols using finite
precision arithmetic, it has been my experience that most algorithms
loose precision pretty quickly. In addition to the algorithms used in
the paper I have also evaluated a number of other routines given to
me privately by other physicists and the conclusion is: they all
suffer from the cancellation effects of alternating sums of large
numbers even for moderatly low j values from around 15. The only
exception is the Schulten and Gordon paper where they have gone to
extraordinary length to evaluate which recursion relations are stable
for which parameter values and even then did the results have to be
renormalised.

I have therefore compared speed only to those routines in the paper
and found that the speed improvement is up to 40 times for Gaunt
coefficients, but only about 3-5 for 6j symbols. The latter is likely
due to the fact that due to the number of triangle relations the
recursion relations are much shorter and therefore faster to evaluate.
However, I have not redone these tests recently and on 64bit machines,
so I don't know what modern machines might give.

More recently people have suggested evaluating the symbols using prime
factor decomposition, essentially doing exact rational number
arithmetic. The problem here is that the factorial represenation is
always limited and for the very large factorials approximations are
used which totally kills the result. A nice example here is the
calculator at Cambridge University:
http://www-stone.ch.cam.ac.uk/wigner.shtml
Try a $3j$ with (18,18,18,0,0,0)

I furthermore get the impression that these algorithms are not very
fast, but I have not done a detailed speed evalutation. For finite
precision symbol calculations the only code that I have found that is
reliable is the Schulten and Gordon one.

For specialised applications it might indeed be advisable to save only
subsets of the symbols or use precise formulas for special parameter
values. But these are optimisations and do not provide generic random
access values which I have tried to address in my paper.

I will try to get the storage codes published soon.


regards

Jens

[jyr]

Just to say that I have opened a ticket now with Dan's changes:

http://trac.sagemath.org/sage_trac/ticket/5996

Jens