Given that you know that one $m_i$ value (let's say $m_3$) will have magnitude at most 2, which means that $m_2$ must be within 2 of $-m_1$. Also, the triangle inequality helps cut down on the storage, so we can cut it down to scaling as $j_{max}^4$.
Because of symmetries, we are free to stipulate that $j_2 \leq j_1$. And let's suppose that $j_1$ can go up to $j_{max} = 80$. We can also cut off about half of the values by using symmetry to stipulate that $m_1 \geq 0$. The naive way to do all this (shown below) includes some invalid indices. But that should be okay, since you'll just never access them.
from __future__ import division, print_function
import numpy as np
from sympy import init_printing, symbols, summation, horner
from sympy.physics.wigner import wigner_3j
init_printing()
try:
from numba import njit
from numba.utils import IS_PY3
except ImportError:
print("Couldn't find numba; that `length` function below will be slow...")
import sys
def _identity_decorator_outer(*args, **kwargs):
def _identity_decorator_inner(fn):
return fn
return _identity_decorator_inner
njit = _identity_decorator_outer
IS_PY3 = (sys.version_info[:2] >= (3, 0))
if IS_PY3:
xrange = range
else:
xrange = xrange
Construct an array of Wigner 3j values as follows:
jmax=10
@njit
def length(jmax):
i=0
for j1 in xrange(jmax+1):
for m1 in xrange(j1+1):
for j2 in xrange(j1+1):
for m2 in xrange(-m1-2,-m1+2+1):
for j3 in xrange((j1-j2),(j1+j2)+1): # Assume j1>=j2
i += 1
return i
print(length(jmax))
21780
%%time
w3j = np.empty((length(jmax),), dtype=float)
i=0
for j1 in range(jmax+1):
for m1 in range(j1+1):
for j2 in range(j1+1):
for m2 in range(-m1-2,-m1+2+1):
for j3 in range((j1-j2),(j1+j2)+1): # Assume j1>=j2
for m3 in [-m1-m2]:
tmp = wigner_3j(j1,j2,j3,m1,m2,m3)
try:
w3j[i] = float(tmp.evalf(n=18))
except AttributeError:
w3j[i] = float(tmp)
i += 1
if i!= w3j.size:
raise ValueError("Something went wrong in assigning the values!")
CPU times: user 14.9 s, sys: 194 ms, total: 15.1 s Wall time: 15 s
print('I expect the computation with jmax=80 will take at least {0:.2g} hours.'.format(15*length(80)/(3600*length(10))))
I expect the computation with jmax=80 will take at least 10 hours.
You can save that array to a text file with something like
np.savetxt('Wigner3j_jmax_80_m3_2.txt', w3j)
Then, load that array from the text file into whatever program you're using, and express the 3j values as a function with something like the following code:
w3j = np.loadtxt('Wigner3j_jmax_80_m3_2.txt')
def wigner3j(j1,j2,j3,m1,m2,m3):
if abs(m1)>j1 or abs(m2)>j2 or abs(m3)>j3 or abs(j1-j2)>j3 or j1+j2<j3 or m1+m2+m3!=0:
return 0.0
if abs(m3)>2: # Re-order so that abs(m3) is <=2, if possible
if abs(m2)<=2:
return wigner3j(j1,j3,j2,m1,m3,m2) * (-1)**(j1+j2+j3)
elif abs(m1)<=2:
return wigner3j(j3,j2,j1,m3,m2,m1) * (-1)**(j1+j2+j3)
else:
raise ValueError("wigner3j was not designed to handle indices {0}".format((j1,j2,j3,m1,m2,m3)))
if j2>j1: # Re-order so that j1>=j2
return wigner3j(j2,j1,j3,m2,m1,m3) * (-1)**(j1+j2+j3)
if m1<0: # Use symmetry to eliminate ~1/2 the storage
return wigner3j(j1,j2,j3,-m1,-m2,-m3) * (-1)**(j1+j2+j3)
# The // division by 2 is standard c-style integer division
return w3j[(j1*(j1*(j1*(5*j1 + 10) + 20*m1 + 5) + 40*m1 - 4) + j2*(20*j2 + 8*m1 + 8*m2 + 20) + 4*j3 + 24*m1 + 4*m2 + 8) // 4]
Here, you can do a simple loop to get the indices. It will include some invalid sets of indices, but that should be fine; you'll just never use them.
indices = [[j1,j2,j3,m1,m2,m3]
for j1 in range(jmax+1) for m1 in range(j1+1)
for j2 in range(j1+1) for m2 in range(-m1-2,-m1+2+1)
for j3 in range(j1-j2,j1+j2+1) for m3 in [-m1-m2]]
print("Total number of elements:",len(indices))
Total number of elements: 21780
The array is created by a nested loop, where the innermost loop goes over values of $j_3$, and the outermost loop goes over values of $j_1$. [Actually, the innermost loop as written is a trivial loop over one value of $m_3$, so we ignore it.] So to index it, you have to start with the number of iterations required to get the innermost loop from the initial value of $j_3$ to the desired value of $j_3$. For the second-to-innermost loop, you have to increase the index by the number of iterations of both loops required to increase the value of $m_2$ from its initial value to the desired value. And so on.
This is achieved by the following sum:
j1,j2,j3,m1,m2,m3,j_max = symbols('j1,j2,j3,m1,m2,m3,j_{max}', integer=True)
index_value = (
summation(1, (j3,j1-j2,j3-1))
+summation(summation(1, (j3,j1-j2,j1+j2)), (m2,-m1-2,m2-1))
+summation(summation(summation(1, (j3,j1-j2,j1+j2)), (m2,-m1-2,-m1+2)), (j2,0,j2-1))
+summation(summation(summation(summation(1, (j3,j1-j2,j1+j2)), (m2,-m1-2,-m1+2)), (j2,0,j1)), (m1,0,m1-1))
+summation(summation(summation(summation(summation(1, (j3,j1-j2,j1+j2)), (m2,-m1-2,-m1+2)), (j2,0,j1)), (m1,0,j1)), (j1,0,j1-1))
).expand().simplify()
index_value
We can check that this actually correctly indexes the indices
array like this:
# Insert the result into the `index` function and the `wigner3j` function above
print(horner(4*index_value))
j1*(j1*(j1*(5*j1 + 10) + 20*m1 + 5) + 40*m1 - 4) + j2*(20*j2 + 8*m1 + 8*m2 + 20) + 4*j3 + 24*m1 + 4*m2 + 8
def index(j1,j2,j3,m1,m2,m3):
return (j1*(j1*(j1*(5*j1 + 10) + 20*m1 + 5) + 40*m1 - 4) + j2*(20*j2 + 8*m1 + 8*m2 + 20) + 4*j3 + 24*m1 + 4*m2 + 8) // 4
def isvalid(j1,j2,j3,m1,m2,m3):
if abs(m1)>j1 or abs(m2)>j2 or abs(m3)>j3 or abs(j1-j2)>j3 or j3>j1+j2:
return False
return True
for i in range(len(indices)):
assert i == index(*indices[i]), '{0}, {1}, {2}'.format(i,indices[i],index(*indices[i]))
print("That seems to have worked")
That seems to have worked
As an additional check, we can calculate the total size of the array, as the index with $j_1 = j_{max}+1$ and everything else as the lowest possible indices for that value of $j_1$.
n_indices = index_value
n_indices = n_indices.subs({j2:0,j3:j1,m1:0,m2:-2,m3:2}).subs(j1,j_max+1).expand().simplify()
display(n_indices)
assert (n_indices.subs(j_max,jmax)) == len(indices), ((n_indices.subs(j_max,jmax)),len(indices))
print("The size {0} seems to be correct.".format(n_indices.subs(j_max,jmax)))
The size 21780 seems to be correct.
Just for fun, let's actually test that the result is the same as sympy's for every value you'll need up to $j_{1,2}=10$, and for every such value with $j_{1,2}=j_{max}$. (That limited range should be convincing enough, and the rest just take too long.)
for j1 in range(10) + [jmax]:
for m1 in range(-j1,j1+1):
for j2 in range(10) + [jmax]:
for m2 in range(-j2,j2+1):
for j3 in range(10) + [jmax]:
for m3 in range(-j3,j3+1):
if abs(m3)>2 and abs(m2)>2 and abs(m1)>2:
continue
a = wigner_3j(j1,j2,j3,m1,m2,m3)
try:
a = float(a.evalf(n=18))
except AttributeError:
a = float(a)
b = wigner3j(j1,j2,j3,m1,m2,m3)
assert a==b, 'sympy...wigner_3j(*{0})={1}; wigner3j(*{0})={2}'.format([j1,j2,j3,m1,m2,m3], a, b)
print("That seems to have worked")
That seems to have worked