-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathexample3.8.m2
89 lines (81 loc) · 8.87 KB
/
example3.8.m2
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
-- ellipse X: x_1^2+4*x_2^2-1
-- objective function f_u = f_u1+f_u2,
-- where f_ui = -1/3a_i(u_i-x_i)^3-1/2b_i(u_i-x_i)^2-c_i(u_i-x_i)
restart
KK = QQ
R = KK[z,a_1,a_2,b_1,b_2,c_1,c_2,d_1,d_2,s,t_1,t_2,u_1,u_2,x_1,x_2]-- with symbolic coefficients
R = KK[z,d_1,d_2,s,t_1,t_2,u_1,u_2,x_1,x_2]-- with random coefficients
for i in 1..2 do a_i = random(KK);
for i in 1..2 do b_i = random(KK);
for i in 1..2 do c_i = random(KK);
E1 = d_1^2-b_1^2+4*a_1*c_1-4*a_1*s*(2*t_1^3+2*t_1*t_2^2+3*t_1^2+t_2^2)
E2 = d_2^2-b_2^2+4*a_2*c_2-4*a_2*s*t_2*(2*t_1^2+2*t_2^2+2*t_1-1)
G1 = x_1-t_1
G2 = x_2-t_2
G3 = 2*a_1*(u_1-t_1)+b_1-d_1
G4 = 2*a_2*(u_2-t_2)+b_2-d_2
H3 = a_1*(u_1-t_1)^2+b_1*(u_1-x_1)+c_1-s*(2*t_1^3+2*t_1*t_2^2+3*t_1^2+t_2^2)
H4 = a_2*(u_2-t_2)^2+b_2*(u_2-x_2)+c_2-s*t_2*(2*t_1^2+2*t_2^2+2*t_1-1)
Gz = z-1
F = (t_1^2+t_2^2+t_1)^2-(t_1^2+t_2^2)
-- all branches at once
I = ideal(E1,E2,G1,G2,H3,H4,Gz,F)
degree I, codim I --(80,8)
I' = eliminate({s,t_1,t_2},I)
degree I', codim I' --(44,5)
J = eliminate (z,I')
degree J, codim J --(44,4)
K = eliminate ({d_1,d_2},J)
degree K, codim K -- (11,2)
-- branch given by G3 and G4
I1 = ideal(E1,E2,G1,G2,G3,G4,Gz,F)
degree I1, codim I1 --(20,8)
time I1' = eliminate({s,t_1,t_2},I1)
degree I1', codim I1' --(11,5)
time J1 = eliminate (z,I1')
degree J1, codim J1 --(11,4)
time K1 = eliminate ({d_1,d_2},J1)
degree K1, codim K1 -- (11,2)
K==K1
-- check that V_P(X) is irreducible
restart
R = CC[u_1,u_2,x_1,x_2]
needsPackage "Bertini"
noiterations = 10
for i in 1..noiterations do (
for i in 1..2 do a_i = random(RR);
for i in 1..2 do b_i = random(RR);
for i in 1..2 do c_i = random(RR);
-- FK1 corresponds to the ideal K1 in the previous code
FK1 = {x_1^4+2*x_1^2*x_2^2+x_2^4+2*x_1^3+2*x_1*x_2^2-x_2^2,3*a_1*a_2^2*u_2^2*x_1^2-4*a_1^2*a_2*u_1^2*x_1*x_2+8*a_1^2*a_2*u_1*x_1^2*x_2-6*a_1*a_2^2*u_2*x_1^2*x_2-4*a_1^2*a_2*x_1^3*x_2-a_1*a_2^2*u_2^2*x_2^2+3*a_1*a_2^2*x_1^2*x_2^2+2*a_1*a_2^2*u_2*x_2^3-a_1*a_2^2*x_2^4+3*a_1*a_2*b_2*u_2*x_1^2+a_1^2*a_2*u_1^2*x_2-2*a_1^2*a_2*u_1*x_1*x_2-4*a_1*a_2*b_1*u_1*x_1*x_2+a_1^2*a_2*x_1^2*x_2+4*a_1*a_2*b_1*x_1^2*x_2-3*a_1*a_2*b_2*x_1^2*x_2-a_1*a_2*b_2*u_2*x_2^2+a_1*a_2*b_2*x_2^3+3*a_1*a_2*c_2*x_1^2+a_1*a_2*b_1*u_1*x_2-a_1*a_2*b_1*x_1*x_2-4*a_1*a_2*c_1*x_1*x_2-a_1*a_2*c_2*x_2^2+a_1*a_2*c_1*x_2,a_1^2*a_2*u_1^2*x_1^2-2*a_1^2*a_2*u_1*x_1^3+4*a_1*a_2^2*u_2^2*x_1*x_2-3*a_1^2*a_2*u_1^2*x_2^2+6*a_1^2*a_2*u_1*x_1*x_2^2-8*a_1*a_2^2*u_2*x_1*x_2^2-5*a_1^2*a_2*x_1^2*x_2^2+4*a_1*a_2^2*x_1*x_2^3-a_1^2*a_2*x_2^4+2*a_1^2*a_2*u_1^2*x_1-4*a_1^2*a_2*u_1*x_1^2+a_1*a_2*b_1*u_1*x_1^2-a_1*a_2*b_1*x_1^3+3*a_1*a_2^2*u_2^2*x_2+4*a_1*a_2*b_2*u_2*x_1*x_2-3*a_1*a_2*b_1*u_1*x_2^2-6*a_1*a_2^2*u_2*x_2^2-2*a_1^2*a_2*x_1*x_2^2+3*a_1*a_2*b_1*x_1*x_2^2-4*a_1*a_2*b_2*x_1*x_2^2+3*a_1*a_2^2*x_2^3+2*a_1*a_2*b_1*u_1*x_1-2*a_1*a_2*b_1*x_1^2+a_1*a_2*c_1*x_1^2+3*a_1*a_2*b_2*u_2*x_2+4*a_1*a_2*c_2*x_1*x_2+a_1^2*a_2*x_2^2-3*a_1*a_2*b_2*x_2^2-3*a_1*a_2*c_1*x_2^2+2*a_1*a_2*c_1*x_1+3*a_1*a_2*c_2*x_2,
24*a_1^2*a_2^2*u_1^2*u_2^2*x_1*x_2-36*a_1^3*a_2*u_1^4*x_2^2-12*a_1*a_2^3*u_2^4*x_2^2+80*a_1^3*a_2*u_1^3*x_1*x_2^2-48*a_1^2*a_2^2*u_1^2*u_2*x_1*x_2^2-96*a_1^3*a_2*u_1*x_1^3*x_2^2-16*a_1^2*a_2^2*u_1*u_2^2*x_2^3+48*a_1*a_2^3*u_2^3*x_2^3+24*a_1^2*a_2^2*u_1^2*x_1*x_2^3+232*a_1^2*a_2^2*u_2^2*x_1*x_2^3-168*a_1^3*a_2*u_1^2*x_2^4+32*a_1^2*a_2^2*u_1*u_2*x_2^4-72*a_1*a_2^3*u_2^2*x_2^4+336*a_1^3*a_2*u_1*x_1*x_2^4-464*a_1^2*a_2^2*u_2*x_1*x_2^4-272*a_1^3*a_2*x_1^2*x_2^4-16*a_1^2*a_2^2*u_1*x_2^5+48*a_1*a_2^3*u_2*x_2^5+232*a_1^2*a_2^2*x_1*x_2^5-52*a_1^3*a_2*x_2^6-12*a_1*a_2^3*x_2^6+3*a_1^3*a_2*u_1^4*x_1-27*a_1*a_2^3*u_2^4*x_1+42*a_1^2*a_2^2*u_1^2*u_2^2*x_2+24*a_1^2*a_2*b_2*u_1^2*u_2*x_1*x_2-36*a_1^2*a_2^2*u_1*u_2^2*x_1*x_2+24*a_1*a_2^2*b_1*u_1*u_2^2*x_1*x_2+108*a_1*a_2^3*u_2^3*x_1*x_2-20*a_1^3*a_2*u_1^3*x_2^2-72*a_1^2*a_2*b_1*u_1^3*x_2^2-84*a_1^2*a_2^2*u_1^2*u_2*x_2^2-24*a_1*a_2^2*b_2*u_2^3*x_2^2+214*a_1^3*a_2*u_1^2*x_1*x_2^2+120*a_1^2*a_2*b_1*u_1^2*x_1*x_2^2-24*a_1^2*a_2*b_2*u_1^2*x_1*x_2^2+72*a_1^2*a_2^2*u_1*u_2*x_1*x_2^2-48*a_1*a_2^2*b_1*u_1*u_2*x_1*x_2^2-162*a_1*a_2^3*u_2^2*x_1*x_2^2-344*a_1^3*a_2*u_1*x_1^2*x_2^2+52*a_1^3*a_2*x_1^3*x_2^2-48*a_1^2*a_2*b_1*x_1^3*x_2^2+42*a_1^2*a_2^2*u_1^2*x_2^3-16*a_1^2*a_2*b_2*u_1*u_2*x_2^3+190*a_1^2*a_2^2*u_2^2*x_2^3-8*a_1*a_2^2*b_1*u_2^2*x_2^3+72*a_1*a_2^2*b_2*u_2^2*x_2^3-36*a_1^2*a_2^2*u_1*x_1*x_2^3+24*a_1*a_2^2*b_1*u_1*x_1*x_2^3+108*a_1*a_2^3*u_2*x_1*x_2^3+232*a_1^2*a_2*b_2*u_2*x_1*x_2^3+12*a_1^3*a_2*u_1*x_2^4-168*a_1^2*a_2*b_1*u_1*x_2^4+16*a_1^2*a_2*b_2*u_1*x_2^4-380*a_1^2*a_2^2*u_2*x_2^4+16*a_1*a_2^2*b_1*u_2*x_2^4-72*a_1*a_2^2*b_2*u_2*x_2^4-113*a_1^3*a_2*x_1*x_2^4-27*a_1*a_2^3*x_1*x_2^4+168*a_1^2*a_2*b_1*x_1*x_2^4-232*a_1^2*a_2*b_2*x_1*x_2^4+190*a_1^2*a_2^2*x_2^5-
8*a_1*a_2^2*b_1*x_2^5+24*a_1*a_2^2*b_2*x_2^5+6*a_1^3*a_2*u_1^4+6*a_1^2*a_2*b_1*u_1^3*x_1-54*a_1*a_2^2*b_2*u_2^3*x_1+42*a_1^2*a_2*b_2*u_1^2*u_2*x_2+36*a_1^2*a_2^2*u_1*u_2^2*x_2+42*a_1*a_2^2*b_1*u_1*u_2^2*x_2+24*a_1^2*a_2*c_2*u_1^2*x_1*x_2-36*a_1^2*a_2*b_2*u_1*u_2*x_1*x_2+24*a_1*a_2*b_1*b_2*u_1*u_2*x_1*x_2+18*a_1^2*a_2^2*u_2^2*x_1*x_2-18*a_1*a_2^2*b_1*u_2^2*x_1*x_2+162*a_1*a_2^2*b_2*u_2^2*x_1*x_2+24*a_1*a_2^2*c_1*u_2^2*x_1*x_2-22*a_1^3*a_2*u_1^2*x_2^2-30*a_1^2*a_2*b_1*u_1^2*x_2^2-36*a_1*a_2*b_1^2*u_1^2*x_2^2-42*a_1^2*a_2*b_2*u_1^2*x_2^2-72*a_1^2*a_2*c_1*u_1^2*x_2^2-72*a_1^2*a_2^2*u_1*u_2*x_2^2-84*a_1*a_2^2*b_1*u_1*u_2*x_2^2-12*a_1*a_2*b_2^2*u_2^2*x_2^2-24*a_1*a_2^2*c_2*u_2^2*x_2^2+68*a_1^3*a_2*u_1*x_1*x_2^2+214*a_1^2*a_2*b_1*u_1*x_1*x_2^2+40*a_1*a_2*b_1^2*u_1*x_1*x_2^2+36*a_1^2*a_2*b_2*u_1*x_1*x_2^2-24*a_1*a_2*b_1*b_2*u_1*x_1*x_2^2+80*a_1^2*a_2*c_1*u_1*x_1*x_2^2-36*a_1^2*a_2^2*u_2*x_1*x_2^2+36*a_1*a_2^2*b_1*u_2*x_1*x_2^2-162*a_1*a_2^2*b_2*u_2*x_1*x_2^2-48*a_1*a_2^2*c_1*u_2*x_1*x_2^2-40*a_1^3*a_2*x_1^2*x_2^2-172*a_1^2*a_2*b_1*x_1^2*x_2^2-4*a_1*a_2*b_1^2*x_1^2*x_2^2+16*a_1^2*a_2*c_1*x_1^2*x_2^2+36*a_1^2*a_2^2*u_1*x_2^3+42*a_1*a_2^2*b_1*u_1*x_2^3-16*a_1^2*a_2*c_2*u_1*x_2^3+190*a_1^2*a_2*b_2*u_2*x_2^3-8*a_1*a_2*b_1*b_2*u_2*x_2^3+24*a_1*a_2*b_2^2*u_2*x_2^3+48*a_1*a_2^2*c_2*u_2*x_2^3+18*a_1^2*a_2^2*x_1*x_2^3-18*a_1*a_2^2*b_1*x_1*x_2^3+54*a_1*a_2^2*b_2*x_1*x_2^3+24*a_1*a_2^2*c_1*x_1*x_2^3+232*a_1^2*a_2*c_2*x_1*x_2^3+52*a_1^3*a_2*x_2^4+6*a_1^2*a_2*b_1*x_2^4-190*a_1^2*a_2*b_2*x_2^4+8*a_1*a_2*b_1*b_2*x_2^4-12*a_1*a_2*b_2^2*x_2^4-168*a_1^2*a_2*c_1*x_2^4-24*a_1*a_2^2*c_2*x_2^4+12*a_1^2*a_2*b_1*u_1^3+3*a_1*a_2*b_1^2*u_1^2*x_1+6*a_1^2*a_2*c_1*u_1^2*x_1-27*a_1*a_2*b_2^2*u_2^2*x_1-54*a_1*a_2^2*c_2*u_2^2*x_1-3*a_1*a_2*b_1^2*x_1^3+12*a_1^2*a_2*c_1*x_1^3+42*a_1^2*a_2*c_2*u_1^2*x_2+36*a_1^2*a_2*b_2*u_1*u_2*x_2+
42*a_1*a_2*b_1*b_2*u_1*u_2*x_2+18*a_1*a_2^2*b_1*u_2^2*x_2+42*a_1*a_2^2*c_1*u_2^2*x_2-36*a_1^2*a_2*c_2*u_1*x_1*x_2+24*a_1*a_2*b_1*c_2*u_1*x_1*x_2+18*a_1^2*a_2*b_2*u_2*x_1*x_2-18*a_1*a_2*b_1*b_2*u_2*x_1*x_2+54*a_1*a_2*b_2^2*u_2*x_1*x_2+24*a_1*a_2*b_2*c_1*u_2*x_1*x_2+108*a_1*a_2^2*c_2*u_2*x_1*x_2-12*a_1^3*a_2*u_1*x_2^2-22*a_1^2*a_2*b_1*u_1*x_2^2-10*a_1*a_2*b_1^2*u_1*x_2^2-36*a_1^2*a_2*b_2*u_1*x_2^2-42*a_1*a_2*b_1*b_2*u_1*x_2^2-20*a_1^2*a_2*c_1*u_1*x_2^2-72*a_1*a_2*b_1*c_1*u_1*x_2^2-36*a_1*a_2^2*b_1*u_2*x_2^2-84*a_1*a_2^2*c_1*u_2*x_2^2-24*a_1*a_2*b_2*c_2*u_2*x_2^2+9*a_1^3*a_2*x_1*x_2^2+34*a_1^2*a_2*b_1*x_1*x_2^2+10*a_1*a_2*b_1^2*x_1*x_2^2-18*a_1^2*a_2*b_2*x_1*x_2^2+18*a_1*a_2*b_1*b_2*x_1*x_2^2-27*a_1*a_2*b_2^2*x_1*x_2^2+174*a_1^2*a_2*c_1*x_1*x_2^2+40*a_1*a_2*b_1*c_1*x_1*x_2^2-24*a_1*a_2*b_2*c_1*x_1*x_2^2-54*a_1*a_2^2*c_2*x_1*x_2^2+18*a_1*a_2^2*b_1*x_2^3+42*a_1*a_2^2*c_1*x_2^3+190*a_1^2*a_2*c_2*x_2^3-8*a_1*a_2*b_1*c_2*x_2^3+24*a_1*a_2*b_2*c_2*x_2^3+6*a_1*a_2*b_1^2*u_1^2+12*a_1^2*a_2*c_1*u_1^2+6*a_1*a_2*b_1*c_1*u_1*x_1-54*a_1*a_2*b_2*c_2*u_2*x_1-6*a_1*a_2*b_1^2*x_1^2+24*a_1^2*a_2*c_1*x_1^2+36*a_1^2*a_2*c_2*u_1*x_2+42*a_1*a_2*b_1*c_2*u_1*x_2+18*a_1*a_2*b_1*b_2*u_2*x_2+42*a_1*a_2*b_2*c_1*u_2*x_2+18*a_1^2*a_2*c_2*x_1*x_2-18*a_1*a_2*b_1*c_2*x_1*x_2+54*a_1*a_2*b_2*c_2*x_1*x_2+24*a_1*a_2*c_1*c_2*x_1*x_2-6*a_1^2*a_2*b_1*x_2^2-18*a_1*a_2*b_1*b_2*x_2^2-22*a_1^2*a_2*c_1*x_2^2-10*a_1*a_2*b_1*c_1*x_2^2-42*a_1*a_2*b_2*c_1*x_2^2-36*a_1*a_2*c_1^2*x_2^2-12*a_1*a_2*c_2^2*x_2^2+12*a_1*a_2*b_1*c_1*u_1+3*a_1*a_2*c_1^2*x_1-27*a_1*a_2*c_2^2*x_1+18*a_1*a_2*b_1*c_2*x_2+42*a_1*a_2*c_1*c_2*x_2+6*a_1*a_2*c_1^2};
elapsedTime sol = bertiniPosDimSolve FK1; -- ~ 260 s per iteration
print(peek sol)
)
-- check that F_X is irreducible
restart
KK = QQ;
R = KK[x,y,u_1,u_2];
I = ideal((x^2+y^2+x)^2-(x^2+y^2)); -- cardioid
jacI = compress transpose jacobian I;
Ising = I + minors(codim(I), jacI);
for i in 1..2 do a_i = random(KK);
for i in 1..2 do b_i = random(KK);
for i in 1..2 do c_i = random(KK);
obj1 = -1/3*a_1*(u_1-x)^3-1/2*b_1*(u_1-x)^2-c_1*(u_1-x);
obj2 = -1/3*a_1*(u_1-y)^3-1/2*b_1*(u_1-y)^2-c_1*(u_1-y);
jacobj = diff(matrix{{x,y}},obj1+obj2);
Icrit = saturate(I+minors(codim(I)+1,jacI||jacobj),Ising);
degree Icrit, codim Icrit
-- this is a zero-dimensional variety, and the cardinality is the 3-degree of the variety.
S = CC[x,y,u_1,u_2];
gg = first entries gens Icrit;
F = apply(gg, i-> sub(i,S));
needsPackage "Bertini"
noiterations = 100
for i in 1..noiterations do (
elapsedTime sol = bertiniPosDimSolve F; -- ~ 1.5 s per iteration
print(peek sol)
)