module jsymbols;
% Author: Matthew Rebbeck, ZIB.
% Advisor: Rene' Grognard.
% Date: March 1994
% Version 1 (experimental code for the symbolic 6j symbol)
% by Winfried Neun (ZIB) email: neun@zib-berlin.de
% ThreeJSymbol with reasoning on the input for optimal computing;
% Note: the code is 'pedestrian' but should be transparent and it
% seems to work !
% It reflects strongly the exploratory code I used to orientate myself
% It should be rewritten in a more 'professional' style;
load!-package 'matrix; % Needed for matrix input.
load!-package 'solve;
load!-package 'ineq;
symbolic procedure clean_up_sqrt(input);
%
% Takes input and factorises out all sqrt's.
%
begin
scalar num,denom,answer;
if not pairp input or car input neq 'quotient then answer := input
else <<
num := cadr input;
denom := caddr input;
num := collect_sqrt(num);
denom := collect_sqrt(denom);
answer := {'quotient,num,denom}; >>;
return answer;
end;
flag('(clean_up_sqrt),'opfn);
symbolic procedure collect_sqrt(input);
%
% Cleans up a series of multiplied sqrt's.
% eg: sqrt(x)*sqrt(y)->sqrt_of(x*y).
%
begin
scalar sqrt_args,extra_bit,minust,answer;
sqrt_args := {}; extra_bit := {};
if not pairp input then answer := input
else <<
if car input = 'minus then
<<minust := t;
if pairp cadr input then input := cadr input
else return input
>>;
if car input='times then input := cdr input else input := {input};
for each elt in input do
<< if eqcar(elt,'sqrt)
then sqrt_args := append(sqrt_args,{cadr elt})
else extra_bit := append(extra_bit,{elt}); >>;
if sqrt_args = {} then <<answer := extra_bit;>>
else <<
sqrt_args:=reval append({'sqrt_of},{append({'times},sqrt_args)});
answer := append({sqrt_args},extra_bit); >>;
answer := append({'times},answer);
if minust then answer := {'minus,answer}; >>;
return answer;
end;
symbolic operator listp;
algebraic operator sqrt_of, oddexpt;
algebraic << let sqrt_of(~x) => sqrt (x) when numberp x >>;
algebraic procedure ThreeJSymbol (u1,u2,u3);
if listp u1 and listp u2 and listp u3 then
begin
scalar j1,j2,j3,m1,m2,m3,tmp,lower,upper,better,best;
matrix range(6,2);
j1:=first u1; m1:=second u1;
j2:=first u2; m2:=second u2;
j3:=first u3; m3:=second u3;
% Vanishing conditions;
if numberp(tmp:=m1+m2+m3) and tmp neq 0 then return 0
else if numberp(tmp:=j1+j2+j3) and tmp < -1 then return 0
else if numberp(tmp:=j1+j2-j3) and tmp < 0 then return 0
else if numberp(tmp:=j1-j2+j3) and tmp < 0 then return 0
else if numberp(tmp:=j2+j3-j1) and tmp < 0 then return 0
else if numberp(tmp:=j1+m1) and tmp < 0 then return 0
else if numberp(tmp:=j1-m1) and tmp < 0 then return 0
else if numberp(tmp:=j2+m2) and tmp < 0 then return 0
else if numberp(tmp:=j2-m2) and tmp < 0 then return 0
else if numberp(tmp:=j3+m3) and tmp < 0 then return 0
else if numberp(tmp:=j3-m3) and tmp < 0 then return 0
else
% Restrictions on k
<<range:=mat((0,infinity),
(0,infinity),
(0,infinity),
(0,infinity),
(0,infinity),
(0,infinity));
% as is;
lower:=range(1,1);upper:=range(1,2);
if numberp(tmp:=j1+j2-j3) then upper:=tmp;
if numberp(tmp:=j1-m1) then
if upper = infinity then upper:=tmp
else upper:=min(upper,tmp);
if numberp(tmp:=j2+m2) then
if upper = infinity then upper:=tmp
else upper:=min(upper,tmp);
if numberp(tmp:=j3+j2+m1) then lower = max(lower,-tmp);
if numberp(tmp:=j2-j1-m2) then lower = max(lower,-tmp);
range(1,1):=lower;range(1,2):=upper;
if upper = infinity then <<better:=upper;best:=0>>
else <<better:=upper-lower+1;best:=1>>;
% {j1,m1} <-> {j2,m2};
lower:=range(2,1);upper:=range(2,2);
if numberp(tmp:=j2+j1-j3) then upper:=tmp;
if numberp(tmp:=j2-m2) then
if upper = infinity then
upper:=tmp else upper:=min(upper,tmp);
if numberp(tmp:=j1+m1) then
if upper = infinity then
upper:=tmp else upper:=min(upper,tmp);
if numberp(tmp:=j3+j1+m2) then lower = max(lower,-tmp);
if numberp(tmp:=j1-j2-m1) then lower = max(lower,-tmp);
range(2,1):=lower;range(2,2):=upper;
if upper neq infinity then
<< tmp:=upper-lower+1;
if (better = infinity) or (better > tmp) then
<< better:=tmp;best:=2 >> >>;
% {j1,m1} <-> {j3,m3};
lower:=range(3,1);upper:=range(3,2);
if numberp(tmp:=j3+j2-j1) then upper:=tmp;
if numberp(tmp:=j3-m3) then
if upper = infinity then
upper:=tmp else upper:=min(upper,tmp);
if numberp(tmp:=j2+m2) then
if upper = infinity then
upper:=tmp else upper:=min(upper,tmp);
if numberp(tmp:=j1+j2+m3) then lower = max(lower,-tmp);
if numberp(tmp:=j2-j3-m2) then lower = max(lower,-tmp);
range(3,1):=lower;range(3,2):=upper;
if upper neq infinity then
<< tmp:=upper-lower+1;
if (better = infinity) or (better > tmp) then
<< better:=tmp;best:=3 >> >>;
% {j2,m2} <-> {j3,m3};
lower:=range(4,1);upper:=range(4,2);
if numberp(tmp:=j1+j3-j2) then upper:=tmp;
if numberp(tmp:=j1-m1) then
if upper = infinity then
upper:=tmp else upper:=min(upper,tmp);
if numberp(tmp:=j3+m3) then
if upper = infinity then
upper:=tmp else upper:=min(upper,tmp);
if numberp(tmp:=j2+j3+m1) then lower = max(lower,-tmp);
if numberp(tmp:=j3-j1-m3) then lower = max(lower,-tmp);
range(4,1):=lower;range(4,2):=upper;
if upper neq infinity then
<< tmp:=upper-lower+1;
if (better = infinity) or (better > tmp) then
<< better:=tmp;best:=4 >> >>;
% {j1,m1} -> {j2,m2} -> {j3,m3} -> {j1,m1};
lower:=range(5,1);upper:=range(5,2);
if numberp(tmp:=j2+j3-j1) then upper:=tmp;
if numberp(tmp:=j2-m2) then
if upper = infinity then
upper:=tmp else upper:=min(upper,tmp);
if numberp(tmp:=j3+m3) then
if upper = infinity
then upper:=tmp else upper:=min(upper,tmp);
if numberp(tmp:=j1+j3+m2) then lower = max(lower,-tmp);
if numberp(tmp:=j3-j2-m3) then lower = max(lower,-tmp);
range(5,1):=lower;range(5,2):=upper;
if upper neq infinity then
<< tmp:=upper-lower+1;
if (better = infinity) or (better > tmp) then
<< better:=tmp;best:=5 >> >>;
% {j1,m1} -> {j3,m3} -> {j2,m2} -> {j1,m1};
lower:=range(6,1);upper:=range(6,2);
if numberp(tmp:=j3+j1-j2) then upper:=tmp;
if numberp(tmp:=j3-m3) then
if upper = infinity then
upper:=tmp else upper:=min(upper,tmp);
if numberp(tmp:=j1+m1) then
if upper = infinity then
upper:=tmp else upper:=min(upper,tmp);
if numberp(tmp:=j2+j1+m3) then lower = max(lower,-tmp);
if numberp(tmp:=j1-j3-m1) then lower = max(lower,-tmp);
range(6,1):=lower;range(6,2):=upper;
if upper neq infinity then
<< tmp:=upper-lower+1;
if (better = infinity) or (better > tmp) then
<< better:=tmp;best:=6 >> >>;
if best = 1 then return !*3j!-symbol!:sign!*(j1-j2-m3)
* clean_up_sqrt(
for k:=range(best,1):range(best,2) sum
if evenp(k)
then !*3j!-symbol!:one!-term!*(k,j1,j2,j3,m1,m2,m3)
else - !*3j!-symbol!:one!-term!*(k,j1,j2,j3,m1,m2,m3))
else if best = 2 then
return !*3j!-symbol!:sign!*((j1+j2+j3)+j2-j1-m3)
* clean_up_sqrt(
for k:=range(best,1):range(best,2) sum
if evenp(k)
then !*3j!-symbol!:one!-term!*(k,j2,j1,j3,m2,m1,m3)
else - !*3j!-symbol!:one!-term!*(k,j2,j1,j3,m2,m1,m3))
else if best = 3 then
return !*3j!-symbol!:sign!*((j1+j2+j3)+j3-j2-m1)
* clean_up_sqrt(
for k:=range(best,1):range(best,2) sum
if evenp(k)
then !*3j!-symbol!:one!-term!*(k,j3,j2,j1,m3,m2,m1)
else - !*3j!-symbol!:one!-term!*(k,j3,j2,j1,m3,m2,m1))
else if best = 4 then
return !*3j!-symbol!:sign!*((j1+j2+j3)+j1-j3-m2)
* clean_up_sqrt(
for k:=range(best,1):range(best,2) sum
if evenp(k)
then !*3j!-symbol!:one!-term!*(k,j1,j3,j2,m1,m3,m2)
else - !*3j!-symbol!:one!-term!*(k,j1,j3,j2,m1,m3,m2))
else if best = 5 then
return !*3j!-symbol!:sign!*(j2-j3-m1) * clean_up_sqrt(
for k:=range(best,1):range(best,2) sum
if evenp(k)
then !*3j!-symbol!:one!-term!*(k,j2,j3,j1,m2,m3,m1)
else - !*3j!-symbol!:one!-term!*(k,j2,j3,j1,m2,m3,m1))
else if best = 6 then
return !*3j!-symbol!:sign!*(j3-j1-m2) * clean_up_sqrt(
for k:=range(best,1):range(best,2) sum
if evenp(k)
then !*3j!-symbol!:one!-term!*(k,j3,j1,j2,m3,m1,m2)
else - !*3j!-symbol!:one!-term!*(k,j3,j1,j2,m3,m1,m2))
else return lisp lpri list("ThreeJSymbol({",
second u1,",",third u1,
"},{",
second u2,",",third u2,
"},{",
second u3,",",third u3,
"}) % symbol best left as is;")
>>
end
else
" the argument must be of the form: {j1,m1},{j2,m2},{j3,m3}"$
algebraic procedure !*3j!-symbol!:sign!* u;
if rounded then (-1)^u else (-1)^(remainder(num u,2*den u)/(den u))$
algebraic procedure !*3j!-symbol!:one!-term!*(k,j1,j2,j3,m1,m2,m3);
sqrt(simplify_factorial(
( factorial(j1+j2-j3)
*factorial(j3+j1-j2)
*factorial(j2+j3-j1)
*factorial(j1+m1)
*factorial(j1-m1)
*factorial(j2+m2)
*factorial(j2-m2)
*factorial(j3+m3)
*factorial(j3-m3)
)/(factorial(j1+j2+j3+1)
*(factorial(k)
*factorial(j1+j2-j3-k)
*factorial(j3-j2+m1+k)
*factorial(j3-j1-m2+k)
*factorial(j1-m1-k)
*factorial(j2+m2-k)
)^2
)
))$
algebraic <<
operator clebsch_gordan;
let clebsch_gordan({~j1,~m1},{~j2,~m2},{~j3,~m3}) =>
ThreeJSymbol ({~j1,~m1},{~j2,~m2},{~j3,~m3}) *
(2*j3+1)^(1/2) * (-1)^(-(j1-j2-m3));
% The 6 J symbol
% The naming of the functions follows Landolt-Boernstein
operator SixJSymbol;
let { SixJSymbol({~j1,~j2,~j3} , {~l1,~l2,~l3}) =>
(begin scalar nnn,mmm,!*factor,!*exp,signum;
% necessary conditions for existence
if (necess_6j (j1,j2,j3,l1,l2,l3) = nil) then return nil;
on factor;
mmm := Racah_W(j1,j2,j3,l1,l2,l3);
nnn :=
square_Racah_delta(j1,j2,j3) * square_Racah_delta(j1,l2,l3)*
square_Racah_delta(l1,j2,l3) * square_Racah_delta(l1,l2,j3)
* mmm^2;
nnn := simplify_factorial (nnn);
signum := sign mmm;
if numberp signum then return (signum * sqrt nnn)
else return (sign nnn * sqrt nnn); end)};
procedure square_Racah_delta(a,b,c);
simplify_factorial(factorial(a+b-c) *factorial(a-b+c)
* factorial(-a +b +c) / factorial (a + b + c +1));
procedure Racah_W(j1,j2,j3,l1,l2,l3);
begin scalar mini,maxi,interv;
mini := min(j1+j2+l1+l2,j2+j3+l2+l3,j3+j1+l3+l1);
maxi := max(0,j1+j2+j3,j1+l2+l3,l1+j2+l3,l1+l2+j3);
aaa:
if numberp (mini - maxi) then
return for k:=maxi:mini sum
((-1)^k*simplify_factorial (factorial (k+1) /
(factorial (k-j1-j2-j3)*
factorial (k-j1-l2-l3) *
factorial (k-l1-j2-l3) *
factorial (k-l1-l2-j3)*
factorial (j1+j2+l1+l2-k)*
factorial (j2+j3+l2+l3-k)*
factorial (j3+j1+l3+l1-k))))
else <<
interv :=findinterval (j1,j2,j3,l1,l2,l3);
if interv = {} then
return
sum((-1)^k* simplify_factorial
(factorial (k+1) /
(factorial (k-j1-j2-j3)*
factorial (k-j1-l2-l3) *
factorial (k-l1-j2-l3) *
factorial (k-l1-l2-j3)*
factorial (j1+j2+l1+l2-k)*
factorial (j2+j3+l2+l3-k)*
factorial (j3+j1+l3+l1-k))),k,maxi,mini)
else <<
maxi := first first interv;
mini := second first interv + maxi;
go to aaa >>;
>>;
end;
% conditions for non vanishing 6j symbol see Landolt/Boernstein
procedure necess_6j(j1,j2,j3,l1,l2,l3);
begin scalar nnn, !*rounded, dmode!*;
off rounded;
nnn := j1 + j2 + j3;
if (numberp nnn) then if not fixp nnn then return nil;
nnn := l1 + l2 + j3;
if (numberp nnn) then if not fixp nnn then return nil;
nnn := j1 + l2 +l3;
if (numberp nnn) then if not fixp nnn then return nil;
nnn := l1 + j2 +l3;
if (numberp nnn) then if not fixp nnn then return nil;
if not numberp j1 or not numberp j2 or not numberp j3
or not numberp l1 or not numberp l2 or not numberp l3
then return t; % don't know
% Triangle condtions
if j1 + j2 - j3 >=0 and j1 - j2 + j3 >=0 and
-j1 + j2 + j3 >=0 and l1 + l2 - j3 >=0 and
l1 - l2 + j3 >=0 and -l1 + l2 + j3 >=0 and
j1 + l2 - l3 >=0 and j1 - l2 + l3 >=0 and
-j1 + l2 + l3 >=0 and l1 + j2 - l3 >=0 and
l1 - j2 + l3 >=0 and -l1 + j2 + l3 >=0 then return t;
end;
procedure aconds!-6j(j1,j2,j3,l1,l2,l3);
{ % Triangle condtions
j1 + j2 - j3 >=0, j1 - j2 + j3 >=0,
-j1 + j2 + j3 >=0,
l1 + l2 - j3 >=0, l1 - l2 + j3 >=0,
-l1 + l2 + j3 >=0,
j1 + l2 - l3 >=0, j1 - l2 + l3 >=0,
-j1 + l2 + l3 >=0,
l1 + j2 - l3 >=0, l1 - j2 + l3 >=0,
-l1 + j2 + l3 >=0,
% condtions for the summation index
!=k +1 >= 0,
!=k-j1-j2-j3 >=0,
!=k-j1-l2-l3 >=0,
!=k-l1-j2-l3 >=0,
!=k-l1-l2-j3 >=0,
j1+j2+l1+l2-!=k >=0,
j2+j3+l2+l3-!=k >=0,
j3+j1+l3+l1-!=k >=0};
% same in Edmonds formulation
procedure conds!-6j(j1,j2,j3,l1,l2,l3);
{ % Triangle condtions
j1 + j2 - j3 >=0, j1 - j2 + j3 >=0,
-j1 + j2 + j3 >=0,
l1 + l2 - j3 >=0, l1 - l2 + j3 >=0,
-l1 + l2 + j3 >=0,
j1 + l2 - l3 >=0, j1 - l2 + l3 >=0,
-j1 + l2 + l3 >=0,
l1 + j2 - l3 >=0, l1 - j2 + l3 >=0,
-l1 + j2 + l3 >=0,
% condtions for the summation index
!=k >= 0,
j1 + j2 + l1 + l2 + 1 -!=k >=0,
j1 + j2 -j3 -!=k >=0,
l1 + l2 - j3 -!=k >=0,
j1 + l2 - l3 -!=k >=0,
l1 + j2 -l3 -!=k >=0,
-j1 -l1 + j3 + l3 + !=k >=0,
-j2 -l2 + j3 + l3 + !=k >=0};
% conditions for non vanishing 3j symbol see Landolt/Boernstein
procedure conds!-3j (j1,j2,j3,m1,m2,m3);
{ m1 + m2 + m3 =0,
j1 + j2 - j3 >=0,
j1 - j2 + j3 >=0,
-j1 + j2 + j3 >=0,
!=k >=0,
j1 + j2 -j3 -!=k >=0,
j1 - m1 -!=k >=0,
j3 + m2 -!=k >=0,
j3 -j2 + m1 + !=k >=0,
j3 - j1 -m2 + !=k >=0};
% Trying to find the approroate intervals for the summation in the
% 6j symbol computation using the ineq package by Herbert Melenk
procedure findinterval(j1,j2,j3,l1,l2,l3);
begin scalar interv,svars;
svars := lisp( 'list .
solvevars list(simp j1,simp j2,simp j3,
simp l1, simp l2,simp l3));
interv :=reverse
ineq_solve(aconds!-6j(j1,j2,j3,l1,l2,l3)
,!=k . svars,record=t);
return findinterval1(part(first interv),0);
end;
>>;
symbolic procedure findinterval1(maxis,minis);
(<<
if eqcar(maxis,'equal) and cadr maxis eq '!=k then
maxis := caddr maxis;
if not eqcar(maxis,'!*interval!*) then
list('list,list('list,maxis , 0))
else <<
minis := caddr maxis;
maxis := cadr maxis;
if eqcar(maxis,'max) then
maxis := cdr maxis else maxis := list maxis;
if eqcar(minis,'min) then
minis := cdr minis else minis := list minis;
foreach xx in maxis do
foreach yy in minis do
if numberp (difff :=reval list('difference,yy,xx)) then
fixedints := {'list,xx,difff} . fixedints;
'list . fixedints
>>
>>) where fixedints = nil , difff = nil;
flag('(findinterval1),'opfn);
symbolic procedure fancy!-clebsch_gordon(u);
begin scalar a,j1,m1,j2,m2,j,m;
u:=cdr u; j1:=cadr car u; m1:=caddr car u;
u:=cdr u; j2:=cadr car u; m2:=caddr car u;
u:=cdr u; j :=cadr car u; m :=caddr car u;
a:={j1,m1,j2,m2,'!|,j1,j2,j,m};
return fancy!-in!-brackets(
{'fancy!-inprint, mkquote 'times,0,mkquote a},
'!(,'!));
end;
put('clebsch_gordon,'fancy!-prifn, 'fancy!-clebsch_gordon);
symbolic procedure fancy!-ThreeJSymbol(u);
fancy!-matpri2({cdr cadr u,cdr caddr u},nil,nil);
put('ThreeJSymbol,'fancy!-prifn, 'fancy!-ThreeJSymbol);
symbolic procedure fancy!-SixJSymbol(u);
fancy!-matpri2({cdr cadr u,cdr caddr u},nil,'("{" . "}"));
put('SixJSymbol,'fancy!-prifn, 'fancy!-SixJSymbol);
symbolic procedure fancy!-NineJSymbol(u);
fancy!-matpri2({cdr cadr u,cdr caddr u, cdr cadddr u},nil,
'("{" . "}"));
put('NineJSymbol,'fancy!-prifn, 'fancy!-NineJSymbol);
endmodule;
end;