I have written a function to compute the Macaulay resultant of 3
degree 2 homogenous polynomials via the determinant of a matrix
depending on the coefficients. It seems to take a very long time to
compute the determinant, I didn't actually get it to finish (which is
not true in Pari/gp). I was able to write another version that
converts that matrix to a gp object and then takes the determinant in
a reasonable amount of time. However, doing it this way seems like I
am missing something in how objects work in sage. Perhaps it is since
I am solving a system that results in free variables, and then using
those free variables as the symoblic entries in the matrix having
never defined them anywhere. Is converting to gp and then taking the
determinant the "right" way to speed up this calculation?
Here is the code with both versions of the resultant.
R.<X,Y,Z> = PolynomialRing(SR, 3)
def f(A,B):
return [A[0]*B[0]^2 + A[1]*B[1]^2 + A[2]*B[2]^2 +
A[3]*B[0]*B[1] + A[4]*B[0]*B[2] + A[5]*B[1]*B[2],A[6]*B[0]^2 +
A[7]*B[1]^2 + A[8]*B[2]^2 + A[9]*B[0]*B[1] + A[10]*B[0]*B[2] +
A[11]*B[1]*B[2],A[12]*B[0]^2 + A[13]*B[1]^2 + A[14]*B[2]^2 +
A[15]*B[0]*B[1] + A[16]*B[0]*B[2] + A[17]*B[1]*B[2]]
x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6,z1,z2,z3,z4,z5,z6=var('x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6,z1,z2,z3,z4,z5,z6')
Avar=[x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6,z1,z2,z3,z4,z5,z6]
def res222(A):
L=f(A,[X,Y,Z])
J=matrix([[derivative(L[0],X),derivative(L[0],Y),derivative(L[0],Z)],
[derivative(L[1],X),derivative(L[1],Y),derivative(L[1],Z)],
[derivative(L[2],X),derivative(L[2],Y),derivative(L[2],Z)]]).determinant()
Jx = derivative(J,X)
Jy = derivative(J,Y)
Jz = derivative(J,Z)
R=-1/512 * (matrix([[A[0],A[1],A[2],A[3],A[4],A[5]],
[A[6],A[7],A[8],A[9],A[10],A[11]],
[A[12],A[13],A[14],A[15],A[16],A[17]],[Jx.coefficient({X:2,Y:0,Z:
0}),Jx.coefficient({X:0,Y:2,Z:0}),Jx.coefficient({X:0,Y:0,Z:
2}),Jx.coefficient({X:1,Y:1,Z:0}),Jx.coefficient({X:1,Y:0,Z:
1}),Jx.coefficient({X:0,Y:1,Z:1})],[Jy.coefficient({X:2,Y:0,Z:
0}),Jy.coefficient({X:0,Y:2,Z:0}),Jy.coefficient({X:0,Y:0,Z:
2}),Jy.coefficient({X:1,Y:1,Z:0}),Jy.coefficient({X:1,Y:0,Z:
1}),Jy.coefficient({X:0,Y:1,Z:1})],[Jz.coefficient({X:2,Y:0,Z:
0}),Jz.coefficient({X:0,Y:2,Z:0}),Jz.coefficient({X:0,Y:0,Z:
2}),Jz.coefficient({X:1,Y:1,Z:0}),Jz.coefficient({X:1,Y:0,Z:
1}),Jz.coefficient({X:0,Y:1,Z:1})]]).determinant())
return R
def res222gp(A):
L=f(A,[X,Y,Z])
J=matrix([[derivative(L[0],X),derivative(L[0],Y),derivative(L[0],Z)],
[derivative(L[1],X),derivative(L[1],Y),derivative(L[1],Z)],
[derivative(L[2],X),derivative(L[2],Y),derivative(L[2],Z)]]).determinant()
Jx = derivative(J,X)
Jy = derivative(J,Y)
Jz = derivative(J,Z)
R=gp(matrix([[A[0],A[1],A[2],A[3],A[4],A[5]],
[A[6],A[7],A[8],A[9],A[10],A[11]],
[A[12],A[13],A[14],A[15],A[16],A[17]],[Jx.coefficient({X:2,Y:0,Z:
0}),Jx.coefficient({X:0,Y:2,Z:0}),Jx.coefficient({X:0,Y:0,Z:
2}),Jx.coefficient({X:1,Y:1,Z:0}),Jx.coefficient({X:1,Y:0,Z:
1}),Jx.coefficient({X:0,Y:1,Z:1})],[Jy.coefficient({X:2,Y:0,Z:
0}),Jy.coefficient({X:0,Y:2,Z:0}),Jy.coefficient({X:0,Y:0,Z:
2}),Jy.coefficient({X:1,Y:1,Z:0}),Jy.coefficient({X:1,Y:0,Z:
1}),Jy.coefficient({X:0,Y:1,Z:1})],[Jz.coefficient({X:2,Y:0,Z:
0}),Jz.coefficient({X:0,Y:2,Z:0}),Jz.coefficient({X:0,Y:0,Z:
2}),Jz.coefficient({X:1,Y:1,Z:0}),Jz.coefficient({X:1,Y:0,Z:
1}),Jz.coefficient({X:0,Y:1,Z:1})]]))
return -1/512 *R.matdet()
###########################################
M=matrix([[0,-1,0],[1,-1,0],[-1,-1,1]])
T=f(Avar,M*vector([X,Y,Z]))
ZZ=M.inverse()*vector(T)
s=solve([ZZ[0].coefficients()[0]==x1,ZZ[0].coefficients()
[2]==x2,ZZ[0].coefficients()[5]==x3,ZZ[0].coefficients()
[1]==x4,ZZ[0].coefficients()[3]==x5,ZZ[0].coefficients()
[4]==x6,ZZ[1].coefficients()[0]==y1,ZZ[1].coefficients()
[2]==y2,ZZ[1].coefficients()[5]==y3,ZZ[1].coefficients()
[1]==y4,ZZ[1].coefficients()[3]==y5,ZZ[1].coefficients()
[4]==y6,ZZ[2].coefficients()[0]==z1,ZZ[2].coefficients()
[2]==z2,ZZ[2].coefficients()[5]==z3,ZZ[2].coefficients()
[1]==z4,ZZ[2].coefficients()[3]==z5,ZZ[2].coefficients()
[4]==z6],x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6,z1,z2,z3,z4,z5,z6)
A= [s[0][i].rhs() for i in [0,..,17]]
if res222gp(A) ==0:
print("not a morphism")
else:
print(f(A,[X,Y,Z]))
--
To post to this group, send email to [email protected]
To unsubscribe from this group, send email to
[email protected]
For more options, visit this group at
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org
