Artifact d097d5c6705a55f28e3872b46e6cac73afc281026b1aae72994853f31344f97f:
- Executable file
r37/packages/specfn/jsymbols.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 17282) [annotate] [blame] [check-ins using] [more...]
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 (sign mmm * sqrt nnn) else return 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;