#MAPLE routines for computing Schlafli functions. #By Warren D. Smith Sept 2008. # #NS(n,sig,X) is the normalized Schlafli function. #Meaning X is an nXn matrix giving the n vertex coordinates (as its rows; these # must lie on the unit sphere or unit hyperboloid; the first coordinate is the special one) #and then NS(n,sig,X) yields the nonEuclidean (n-1)-volume of the simplex with those n vertices #NORMALIZED by dividing by SphereSurf(n). Sig is +1 for spherical and -1 for hyperbolic. #NonEucVol(n,sig,X) gives the un-normalized volume. # #Note these routines are designed to be used with exact arithmetic and may exhibit poor #numerical behavior (which can be overcome by using extended precision arithmetic). #Also, they might be slow, especially if run symbolically (faster if used on floats as input). # #I've tested it on about 100 examples using MAPLE9 #but my tests are not exhaustive so there might still be bugs. #The most likely kind of bug is an algebraic-branch error. #Please report bugs to warren.wds AT gmail.com with(LinearAlgebra); #example: spheresurf(2) = 2*Pi. SphereSurf := (n) -> if(n=0) then 0 else 2*Pi^(n/2) / GAMMA(n/2) fi; #example SigMat(3, -1) = diag(1, -1, -1). SigMat := proc(n,sig) local M; M := sign(sig)*Matrix(n,n,shape=identity); M[1,1] := 1; return(M); end proc; NF4 := proc(A) local L, U, G, j,k, Up, Um, q, sdG, ft, zp, zm, res, x1,xm,xp; L := (x) -> polylog(2, x); #Li2, use L for brevity U := (Z,z) -> ( #this (my defn) is a factor 2 off Murakami-Yano defn L(z) + L(Z[1,2]*Z[1,3]*Z[3,4]*Z[2,4]*z) + L(Z[1,2]*Z[2,3]*Z[3,4]*Z[1,4]*z) + L(Z[1,3]*Z[2,3]*Z[2,4]*Z[1,4]*z) - L(-Z[1,2]*Z[1,3]*Z[2,3]*z) - L(-Z[1,2]*Z[2,4]*Z[1,4]*z) - L(-Z[1,3]*Z[3,4]*Z[1,4]*z) - L(-Z[2,3]*Z[3,4]*Z[2,4]*z) ); Z := Matrix( 4,4, shape=symmetric ); G := Matrix( 4,4, shape=symmetric ); for j from 1 to 4 do G[j,j] := 1; Z[j,j] := 1; for k from 1 to j-1 do if(A[j,k] != A[k,j]) then print("nonsymm matrix in NF4"); return(FAIL); fi; G[j,k] := -cos(A[j,k]); #grammian matrix Z[j,k] := exp(I*A[j,k]); #exp(I*dihedral) od; od; sdG := sqrt(Determinant(G)); q := Z[1,2]*Z[3,4]+Z[1,3]*Z[2,4]+Z[2,3]*Z[1,4] + Z[1,2]*Z[1,3]*Z[1,4] + Z[1,2]*Z[2,3]*Z[2,4] + Z[1,3]*Z[2,3]*Z[3,4] + Z[1,4]*Z[2,4]*Z[3,4] + Z[1,2]*Z[1,3]*Z[1,4]*Z[2,3]*Z[2,4]*Z[3,4]; ft := sin(A[1,2])*sin(A[3,4])+sin(A[1,3])*sin(A[2,4])+sin(A[1,4])*sin(A[2,3]); zp := (-2/q) * ( ft + sdG ); zm := (-2/q) * ( ft - sdG ); Up := U(Z, zp); Um := U(Z, zm); x1 := z*diff(U(Z,z), z); xm := subs(z=zm, x1); xp := subs(z=zp, x1); res := (Up - xp*ln(zp) - Um + xm*ln(zm)) / (4 * SphereSurf(4)); return(res); end proc; #example NS(5, 1, Matrix(5,5,shape=identity) ) = 1/32. NS := proc( n, sig, X ) local L, U, DP, C, A, G, Z, j, k, Up, Um, q, sdG, ft, zp, zm, coldm, res,r2, x1,xm,xp; L := (x) -> polylog(2, x); #Li2, using L for brevity U := (Z,z) -> ( #this (my defn) is a factor 2 off Murakami-Yano defn L(z) + L(Z[1,2]*Z[1,3]*Z[3,4]*Z[2,4]*z) + L(Z[1,2]*Z[2,3]*Z[3,4]*Z[1,4]*z) + L(Z[1,3]*Z[2,3]*Z[2,4]*Z[1,4]*z) - L(-Z[1,2]*Z[1,3]*Z[2,3]*z) - L(-Z[1,2]*Z[2,4]*Z[1,4]*z) - L(-Z[1,3]*Z[3,4]*Z[1,4]*z) - L(-Z[2,3]*Z[3,4]*Z[2,4]*z) ); if(n=0) then 1 else if(n=1) then 1/SphereSurf(1) else if(2<=n and n<=5) then coldm := ColumnDimension(X); DP := evalm( MatrixMatrixMultiply(X, MatrixMatrixMultiply(SigMat(coldm, sig), Transpose(X) ))); C := evalm( DP^(-1) ); A := Matrix( n,n, shape=symmetric ); G := Matrix( n,n, shape=symmetric ); Z := Matrix( n,n, shape=symmetric ); for j from 1 to n do G[j,j] := 1; A[j,j] := 0; Z[j,j] := 1; if(abs(evalf(DP[j,j] - 1)) > 0.01) then print("DP=", DP, evalf(DP)); print("erroneous NS input, norm not 1"); return(FAIL); fi; for k from 1 to j-1 do G[j,k] := sign(sig)*C[j,k]/sqrt(C[k,k]*C[j,j]); #grammian matrix A[j,k] := Pi - arccos(G[j,k]); #internal dihedrals Z[j,k] := exp(I*A[j,k]); #exp(I*dihedral) if(sig<0 and evalf(DP[j,k]) < -0.01) then print("Erroneous NS input, Minkowski inner products not all same sign"); return(FAIL); fi; od; od; if(n=2) then if(sig>0) then arccos( DP[1,2] ) / SphereSurf(2); else arccosh( DP[1,2] ) / SphereSurf(2); fi; else if(n=3) then ( A[1,2] + A[2,3] + A[1,3] - Pi )*sign(sig) / SphereSurf(3); else if(n=4) then sdG := sqrt(Determinant(G)); q := Z[1,2]*Z[3,4]+Z[1,3]*Z[2,4]+Z[2,3]*Z[1,4] + Z[1,2]*Z[1,3]*Z[1,4] + Z[1,2]*Z[2,3]*Z[2,4] + Z[1,3]*Z[2,3]*Z[3,4] + Z[1,4]*Z[2,4]*Z[3,4] + Z[1,2]*Z[1,3]*Z[1,4]*Z[2,3]*Z[2,4]*Z[3,4]; ft := sin(A[1,2])*sin(A[3,4])+sin(A[1,3])*sin(A[2,4])+sin(A[1,4])*sin(A[2,3]); zp := (-2/q) * ( ft + sdG ); zm := (-2/q) * ( ft - sdG ); Up := U(Z, zp); Um := U(Z, zm); x1 := z*diff(U(Z,z), z); xm := subs(z=zm, x1); xp := subs(z=zp, x1); res := (Up - xp*ln(zp) - Um + xm*ln(zm)) * sqrt(sign(sig)) / (4 * SphereSurf(4)); return(res); #simplifier often makes it worse, so not simplifying here else res := (1/2) * ( #1 row+col missing NF4( Matrix(4,4, [ [A[2,2],A[2,3],A[2,4],A[2,5]], [A[3,2],A[3,3],A[3,4],A[3,5]], [A[4,2],A[4,3],A[4,4],A[4,5]], [A[5,2],A[5,3],A[5,4],A[5,5]] ])) + NF4( Matrix(4,4, [ [A[1,1],A[1,3],A[1,4],A[1,5]], [A[3,1],A[3,3],A[3,4],A[3,5]], [A[4,1],A[4,3],A[4,4],A[4,5]], [A[5,1],A[5,3],A[5,4],A[5,5]] ])) + NF4( Matrix(4,4, [ [A[1,1],A[1,2],A[1,4],A[1,5]], [A[2,1],A[2,2],A[2,4],A[2,5]], [A[4,1],A[4,2],A[4,4],A[4,5]], [A[5,1],A[5,2],A[5,4],A[5,5]] ])) + NF4( Matrix(4,4, [ [A[1,1],A[1,2],A[1,3],A[1,5]], [A[2,1],A[2,2],A[2,3],A[2,5]], [A[3,1],A[3,2],A[3,3],A[3,5]], [A[5,1],A[5,2],A[5,3],A[5,5]] ])) + NF4( Matrix(4,4, [ [A[1,1],A[1,2],A[1,3],A[1,4]], [A[2,1],A[2,2],A[2,3],A[2,4]], [A[3,1],A[3,2],A[3,3],A[3,4]], [A[4,1],A[4,2],A[4,3],A[4,4]] ] ))) - (1/4) * ( #sum of normalized dihedrals A[1,2] + A[1,3] + A[1,4] + A[1,5] + A[2,3] + A[2,4] + A[2,5] + A[3,4] + A[3,5] + A[4,5] ) / SphereSurf(2) + (1/2); #all rows missing return(res); fi; fi; fi; else print("Erroneous NS input, only support n in range 0..5"); return(FAIL); fi; fi; fi; end proc; NonEucVol := (n, sig, X) -> NS(n, sig, X) * SphereSurf(n); #example PureLinearWeightedNS( 6, Matrix(6,6,shape=identity), Vector( [6,5,4,3,2,1] ) ) = 21/32; PureLinearWeightedNS := proc (n, X, Lvec) local j,k,g,f,H,M,res; res := 0; H := Transpose(X)^(-1); M := Matrix(n-1, n); for j from 1 to n do f := Multiply( Normalize(Row(H,j)), Lvec ); for k from 1 to j-1 do for g from 1 to n do M[k,g] := X[k,g]; od; od; for k from j+1 to n do for g from 1 to n do M[k-1,g] := X[k,g]; od; od; res := res + f * NS(n-1, 1, M); od; return(res); end proc; #end---------------------------------------------