Compute implicit cartesian equation for curves given by a parametric equation including trigonometric functions (sage script and form script for output optimization)

implicit_curve.frm

#-
Symbol x,y,u;
Format float;
*Example curve, just paste the output of the sage script in here
Local F=x^24 - 13*x^22*y^2 - 9*x^20*y^4 + 285*x^18*y^6 + 970*x^16*y^8 + 1246*x^14*y^10 + 462*x^12*y^12 - 454*x^10*y^14 - 475*x^8*y^16 - 65*x^6*y^18 + 75*x^4*y^20 + 25*x^2*y^22 - 42*x^22 - 637*x^20*y^2 - 2660*x^18*y^4 - 6125*x^16*y^6 - 11340*x^14*y^8 - 18746*x^12*y^10 - 23296*x^10*y^12 - 18690*x^8*y^14 - 8890*x^6*y^16 - 2345*x^4*y^18 - 364*x^2*y^20 - 49*y^22 + 511*x^20 + 4585*x^18*y^2 + 22470*x^16*y^4 + 64260*x^14*y^6 + 111930*x^12*y^8 + 126126*x^10*y^10 + 98280*x^8*y^12 + 55860*x^6*y^14 + 22575*x^4*y^16 + 5425*x^2*y^18 + 490*y^20 - 1484*x^18 - 14231*x^16*y^2 - 53424*x^14*y^4 - 119756*x^12*y^6 - 184184*x^10*y^8 - 194194*x^8*y^10 - 132496*x^6*y^12 - 54684*x^4*y^14 - 12796*x^2*y^16 - 1519*y^18 + 1519*x^16 + 11277*x^14*y^2 + 43407*x^12*y^4 + 89089*x^10*y^6 + 105105*x^8*y^8 + 79079*x^6*y^10 + 40677*x^4*y^12 + 12747*x^2*y^14 + 1484*y^16 - 490*x^14 - 3955*x^12*y^2 - 9240*x^10*y^4 - 15785*x^8*y^6 - 19250*x^6*y^8 - 11781*x^4*y^10 - 3052*x^2*y^12 - 511*y^14 + 49*x^12 + 119*x^10*y^2 + 1260*x^8*y^4 + 910*x^6*y^6 + 105*x^4*y^8 + 427*x^2*y^10 + 42*y^12 - 25*x^8*y^2 + 100*x^6*y^4 - 110*x^4*y^6 + 20*x^2*y^8 - y^10;
Local Fx=24*x^23 - 286*x^21*y^2 - 180*x^19*y^4 + 5130*x^17*y^6 + 15520*x^15*y^8 + 17444*x^13*y^10 + 5544*x^11*y^12 - 4540*x^9*y^14 - 3800*x^7*y^16 - 390*x^5*y^18 + 300*x^3*y^20 + 50*x*y^22 - 924*x^21 - 12740*x^19*y^2 - 47880*x^17*y^4 - 98000*x^15*y^6 - 158760*x^13*y^8 - 224952*x^11*y^10 - 232960*x^9*y^12 - 149520*x^7*y^14 - 53340*x^5*y^16 - 9380*x^3*y^18 - 728*x*y^20 + 10220*x^19 + 82530*x^17*y^2 + 359520*x^15*y^4 + 899640*x^13*y^6 + 1343160*x^11*y^8 + 1261260*x^9*y^10 + 786240*x^7*y^12 + 335160*x^5*y^14 + 90300*x^3*y^16 + 10850*x*y^18 - 26712*x^17 - 227696*x^15*y^2 - 747936*x^13*y^4 - 1437072*x^11*y^6 - 1841840*x^9*y^8 - 1553552*x^7*y^10 - 794976*x^5*y^12 - 218736*x^3*y^14 - 25592*x*y^16 + 24304*x^15 + 157878*x^13*y^2 + 520884*x^11*y^4 + 890890*x^9*y^6 + 840840*x^7*y^8 + 474474*x^5*y^10 + 162708*x^3*y^12 + 25494*x*y^14 - 6860*x^13 - 47460*x^11*y^2 - 92400*x^9*y^4 - 126280*x^7*y^6 - 115500*x^5*y^8 - 47124*x^3*y^10 - 6104*x*y^12 + 588*x^11 + 1190*x^9*y^2 + 10080*x^7*y^4 + 5460*x^5*y^6 + 420*x^3*y^8 + 854*x*y^10 - 200*x^7*y^2 + 600*x^5*y^4 - 440*x^3*y^6 + 40*x*y^8;
Local Fy=-26*x^22*y - 36*x^20*y^3 + 1710*x^18*y^5 + 7760*x^16*y^7 + 12460*x^14*y^9 + 5544*x^12*y^11 - 6356*x^10*y^13 - 7600*x^8*y^15 - 1170*x^6*y^17 + 1500*x^4*y^19 + 550*x^2*y^21 - 1274*x^20*y - 10640*x^18*y^3 - 36750*x^16*y^5 - 90720*x^14*y^7 - 187460*x^12*y^9 - 279552*x^10*y^11 - 261660*x^8*y^13 - 142240*x^6*y^15 - 42210*x^4*y^17 - 7280*x^2*y^19 - 1078*y^21 + 9170*x^18*y + 89880*x^16*y^3 + 385560*x^14*y^5 + 895440*x^12*y^7 + 1261260*x^10*y^9 + 1179360*x^8*y^11 + 782040*x^6*y^13 + 361200*x^4*y^15 + 97650*x^2*y^17 + 9800*y^19 - 28462*x^16*y - 213696*x^14*y^3 - 718536*x^12*y^5 - 1473472*x^10*y^7 - 1941940*x^8*y^9 - 1589952*x^6*y^11 - 765576*x^4*y^13 - 204736*x^2*y^15 - 27342*y^17 + 22554*x^14*y + 173628*x^12*y^3 + 534534*x^10*y^5 + 840840*x^8*y^7 + 790790*x^6*y^9 + 488124*x^4*y^11 + 178458*x^2*y^13 + 23744*y^15 - 7910*x^12*y - 36960*x^10*y^3 - 94710*x^8*y^5 - 154000*x^6*y^7 - 117810*x^4*y^9 - 36624*x^2*y^11 - 7154*y^13 + 238*x^10*y + 5040*x^8*y^3 + 5460*x^6*y^5 + 840*x^4*y^7 + 4270*x^2*y^9 + 504*y^11 - 50*x^8*y + 400*x^6*y^3 - 660*x^4*y^5 + 160*x^2*y^7 - 10*y^9;
*Hack for Simultaneous optimization as described in the Paper (https://arxiv.org/pdf/1310.7007.pdf), F Fx and Fy can be easily extracted
Local G = u*F+u^2*Fx+u^3*Fy;
Bracket u;

*Simple Variant of Common Subexpression Elimination, High Op Count, Low Variable Usage
*Format O4,saIter=1000,stats=on,Method=CSE;

*Combined Variant, CSE First, then Greedy, Low Op Count, High Variable Usage (very near to Greedy)
*Format O4,saIter=1000,stats=on,Method=CSEGreedy;

*Greedy Variant, Lowest Op Count, Highest Variable Usage
Format O4,saIter=1000,stats=on,Method=Greedy;

*Optimize separately
*Print F, Fx, Fy;

*Optimize simultaneous
Print G;

.end

trigonometric_parametric_curve_implicitization.sage

#x(t) and y(t) give the parametric formula,
#d is the assumed degree of the implicit equation
#if degree is too low, the matrix will have no solution
#if degree is too high, it will take a lot of time and ram
#closed curves have even degree
#x_symmetric and y_symmetric should be true if the curve is known to be _exactly_ symmetric along these axes (performance optimization)

#uncomment out the example curve you want to try or write a new one

#circle
#x(t)=cos(t)
#y(t)=sin(t)
#x_symmetric=True
#y_symmetric=True
#d=2

#lemniscate of bernoulli
#x(t)=cos(t)/(sin(t)^2+1)
#y(t)=sin(t)*cos(t)/(sin(t)^2+1)
#x_symmetric=True
#y_symmetric=True
#d=4

#deltoid
#x(t)=2*cos(t)+cos(2*t)
#y(t)=2*sin(t)-sin(2*t)
#x_symmetric=False
#y_symmetric=True
#d=4

#cardioid
#x(t)=(1-cos(t))*cos(t)
#y(t)=(1-cos(t))*sin(t)
#x_symmetric=False
#y_symmetric=True
#d=4

#kampyle of eudoxus
#x(t)=sec(t)
#y(t)=tan(t)*sec(t)
#x_symmetric=True
#y_symmetric=True
#d=4

#bicorn
#x(t)=sin(t)
#y(t)=cos(t)^2*(2+cos(t))/(3+sin(t)^2)
#x_symmetric=True
#y_symmetric=False
#d=4

#astroid
#x(t)=cos(t)^3
#y(t)=sin(t)^3
#x_symmetric=True
#y_symmetric=True
#d=6

#quadrifolium
x(t)=cos(2*t)*cos(t)
y(t)=cos(2*t)*sin(t)
x_symmetric=True
y_symmetric=True
d=6

#x(t)=tan(3*t)*cos(5*t)
#y(t)=tan(3*t)*sin(5*t)
#x_symmetric=True
#y_symmetric=True
#d=16

#x(t)=17/20*cos(t)*((3/5)+2/5*cos(t*7))
#y(t)=17/20*cos(t+1)*(3/5+2/5*cos(t*7+1))
#x_symmetric=False
#y_symmetric=False
#d=16

#x(t)=tan(5*t)*cos(7*t)
#y(t)=tan(5*t)*sin(7*t)
#x_symmetric=True
#y_symmetric=True
#d=24

#x(t)=(1+3/4*cos(t*10))*cos(t)
#y(t)=(1-3/4*cos(t*10))*sin(t)
#x_symmetric=True
#y_symmetric=True
#d=22

#x(t)=(3/2+cos(10*t))*cos(3*t)
#y(t)=(3/2-cos(10*t))*sin(3*t)
#x_symmetric=True
#y_symmetric=True
#d=26

#pre_x(t)=17/20*cos(t)*((3/5)+2/5*cos(t*7))
#pre_y(t)=17/20*cos(t+1/3*pi)*(3/5+2/5*cos(t*7+1/3*pi))
#x(t)=1/2*sqrt(2)*(pre_x(t)+pre_y(t))
#y(t)=1/2*sqrt(2)*(pre_x(t)-pre_y(t))
#x_symmetric=False
#y_symmetric=True
#d=16

var("u v")

x1=x(t).trig_expand(full=True).trig_simplify().subs(cos(t)==u, sin(t)==v)
y1=y(t).trig_expand(full=True).trig_simplify().subs(cos(t)==u, sin(t)==v)

max_denom=100

#approximate potentially occuring constant sin, cos and sqrt terms with small rationals
#w0=SR.wild(0)
#x1=x1.substitution_delayed(sin(w0),lambda d:sin(d[w0]).n().nearby_rational(max_denominator=max_denom))
#x1=x1.substitution_delayed(cos(w0),lambda d:cos(d[w0]).n().nearby_rational(max_denominator=max_denom))
#x1=x1.substitution_delayed(sqrt(w0),lambda d:sqrt(d[w0]).n().nearby_rational(max_denominator=max_denom))

#y1=y1.substitution_delayed(sin(w0),lambda d:sin(d[w0]).n().nearby_rational(max_denominator=max_denom))
#y1=y1.substitution_delayed(cos(w0),lambda d:cos(d[w0]).n().nearby_rational(max_denominator=max_denom))
#y1=y1.substitution_delayed(sqrt(w0),lambda d:sqrt(d[w0]).n().nearby_rational(max_denominator=max_denom))

#a=None
#K=QQ

#approximate potentially occuring constant sin and cos terms with argument small rational to pi (leads to algebraic numbers)
#for some curves this variant is better, for others the first variant is better, for many it doesn't matter
w0=SR.wild(0)
x1=x1.substitution_delayed(sin(w0),lambda d:sin((d[w0]/pi).n().nearby_rational(max_denominator=max_denom)*pi))
x1=x1.substitution_delayed(cos(w0),lambda d:cos((d[w0]/pi).n().nearby_rational(max_denominator=max_denom)*pi))

y1=y1.substitution_delayed(sin(w0),lambda d:sin((d[w0]/pi).n().nearby_rational(max_denominator=max_denom)*pi))
y1=y1.substitution_delayed(cos(w0),lambda d:cos((d[w0]/pi).n().nearby_rational(max_denominator=max_denom)*pi))

K=AA

R.<u,v>=K[]

x1=R(x1)
y1=R(y1)

K=number_field_elements_from_algebraics(list(x1.dict().values())+list(y1.dict().values()),embedded=True)[0]
a=K.gen()

R.<u,v>=K[]
S=FractionField(R)

x1=S(x1)
y1=S(y1)

if x1.denominator() == 1 and y1.denominator() == 1:
	x1=R(x1)
	y1=R(y1)

	rational=False
else:
	l=R(lcm(x1.denominator(),y1.denominator()))

	x1=R(x1.numerator() * (l / x1.denominator()) )
	y1=R(y1.numerator() * (l / y1.denominator()) )

	rational=True

print("x(t)="+str(x(t)))
print("y(t)="+str(y(t)))
print("x_symmetric="+str(x_symmetric))
print("y_symmetric="+str(y_symmetric))
print("d="+str(d))

polys=[]

max_d=0

trig_identity=u^2+v^2-1

for i in range(d+1):
	for j in range(i+1):
		if (not x_symmetric or (j%2)==0) and (not y_symmetric or ((i-j)%2)==0):
			if rational:
				polys+=[(x1^j*y1^(i-j)*l^(d-i))%trig_identity]
			else:
				polys+=[(x1^j*y1^(i-j))%trig_identity]
			max_d=max(max_d, polys[-1].degree())
		else:
			polys+=[R(0)]

m=binomial(max_d+2,2)
n=binomial(d+2,2)

M=MatrixSpace(K, m, n, sparse=True)

print("max_d="+str(max_d))

A=copy(M.zero_matrix())
for (i,poly) in enumerate(polys):
	for key in poly.dict():
		r = key[0]
		s = key[1]
		j = binomial(r+s+1,2)+r
		A[j,i] = poly.dict()[key]

B=A.delete_columns([n-1])
V=A.column(n-1)

print("\nReady to solve linear equation system\n")

sol = B \ V

sol2 = list(sol) + [-1]

sol_poly=0

R.<u,v>=K["x, y"]

index=0
for i in range(d+1):
	for j in range(i+1):
		sol_poly+=sol2[index]*u^j*v^(i-j)
		index+=1

facs=factor(sol_poly)

sol_poly2=0

for fac in facs:
	if rational:
		if bool((fac[0](x=x1/l,y=y1/l)).numerator()%trig_identity==0):
			sol_poly2=fac[0]
	else:
		if bool(fac[0](x=x1,y=y1)%trig_identity==0):
			sol_poly2=fac[0]

if not a==None and not a.is_rational():
	print("a.n(100)="+str(a.n(100)))
	print("a.n(100).nearby_rational(max_denominator=10000)="+str(a.n(100).nearby_rational(max_denominator=10000)))
	print("a.minpoly()="+str(a.minpoly())+"\n")

print("Implicit cartesian equation:")
print(sol_poly2)
print("\nPartial derivative with respect to x:")
print(diff(sol_poly2,u))
print("Partial derivative with respect to y:")
print(diff(sol_poly2,v))