Wronski Polynomial Systems with polymake and Singular (Mega2013 submission)

Here you find hints for the correct compilation and installation to use the new interface between polymake and Singular. First of all you need a new perpetual beta version (we used 20130307 in this case). And you need to install libsingular as a bundled extension.

To recreate the computations we did for our Mega2013 submission you can download two script files:

• random_solve.script
• mega2013.script

In the first script file two functions are defined: random_wronski and solve_random_wronski. Here is the help text of the two functions:

# Create a random (possible) Wronski System
# Usage: random_wronski(M,lambda,{seed=>r, triangulation=>T, param=>s, ring=>R});
# @param Matrix<Int> M points of the Triangulation (Homogenized)
# @param Vector<Int> lambda lifting function
# @option Int seed for randomizer
# @option topaz::GeometricSimplicialComplex triangulation the Triangulation
# @option Scalar<Rational> param the additional parameter for creating a wronski system
# @option Ring ring The Ring
sub random_wronski($$;$){
...
}


# Solves many random wronski systems and checks how many real solutions there are.
# the data will be stored in an output file
# All systems are seeded, so that one can check the individual systems
# Usage: solve_random_wronski(M,lambda,n,{seedoffset=>r, param=>s, upperparam=>s0, file=>"path", triangulation=>T});
# @param Matrix<Int> M points of the Triangulation (Homogenized)
# @param Vector<Int> lambda lifting function of the triangulation
# @param Int n Number of Systems to solve
# @option Int seedoffset the starting seed
# @option Scalar<Rational> param the additional parameter of the wronski system
# @option Scalar<Rational> upperparam the upper bound for the additional parameter
# @option File the output file of the analysis
# @option topaz::GeometricSimplicialComplex T Triangulation
sub solve_random_wronski($$$;$){
...
}

the second script file uses the first script to recreate the computations we did. Save them in the same directory and start polymake while you are in this directory. Because the second script uses the first one you can load both scripts at once using the command

polytope> script("mega2013");

Without talking about the mathematical background here is something one might do with the new interface. The goal is to check if the lower bound by Soprunova and Sottile for the number of real solutions of a wronski system is met. We start out with a foldable triangulation of a point configuration. We use the lattice points of a scaled simplex.

polytope > $M = new Matrix<Int>(simplex(2,3)->LATTICE_POINTS);
polymake: used package cddlib
  Implementation of the double description method of Motzkin et al.
  Copyright by Komei Fukuda.
  http://www.ifor.math.ethz.ch/~fukuda/cdd_home/cdd.html

polytope > $lambda = new Vector<Int>([12,3,0,0,8,1,0,9,5,15]);

polytope > $F = regular_subdivision($M, $lambda);

polytope > print $F;
{5 7 8}
{5 6 8}
{2 5 6}
{2 3 6}
{0 1 4}
{1 2 5}
{1 4 5}
{4 5 7}
{7 8 9}

polytope > $T = new topaz::GeometricSimplicialComplex(COORDINATES=>$M->minor(All,~[0]), FACETS=>$F);

polytope > print $T->FOLDABLE;
1
polytope > print $T->SIGNATURE;
3

The signature will be be our lower bound, if some conditions are satisfied. Next we need switch to the application 'ideal'

polytope > application 'ideal';

Application ideal currently uses following third-party software packages:
4ti2, Singular, cddlib, geomview, graphviz, jreality, latte, lrslib, nauty, normaliz2, permlib, povray, qhull, sketch, sympol, topcom, tosimplex
For more details:  show_credits;

Now we check one of the conditions. Therefore we use the new function 'wronski_center_ideal' which creates an ideal.

ideal > $I = wronski_center_ideal($M,$lambda);
ideal > print join "\n", @{$I->GENERATORS};
s^12 + x0^3*s^15 + x0*x1*s + x1^3
x1*s^3 + x0*x1^2 + x0^2*s^9
x0*s^8 + x1^2 + x0^2*x1*s^5

We want to look at the real solutions of this particular ideal. But it suffice to look at the 's' coordinate. We calculate a Groebner Basis, which is done with singular.

ideal > $e = $I->GROEBNER(ORDER=>'lex')->BASIS->[0];
polymake: used package Singular
  Singular is a computer algebra system for polynomial computations, 
  with special emphasis on commutative and non-commutative algebra, 
  algebraic geometry, and singularity theory. It is free and open-source 
  under the GNU General Public Licence.
  http://www.singular.uni-kl.de/

ideal > print $e;
s^33 -1*s^15 -1*s^30 + 3*s^23 -6*s^31 -3*s^32

A solution of this System must also induce a solution of the polynomial '$e'. Therefore we look at the solutions of '$e' since we are only interested in the 's' coordinate of the solutions. And because ideal and polynomial objects are linked together with a certain ring, the polynomial '$e' looks like an univariate polynomial but it is a polynomial in the ring of the ideal, which has 3 variables 'x0', 'x1' and 's'. So we need to do some transformation, to solve this polynomial with 'Singular'

ideal > $R = new Ring(["s"]);

ideal > $univariate_e = new Polynomial($e->monomials_as_matrix->minor(All,[2]),$e->coefficients_as_vector,$R);

ideal > print new Ideal(GENERATORS=>[$univariate_e], RING=>$R)->SOLVE;
(-0.8952189506 0)
(4.411470567 0)
...
(0 0)

Looking at the real solutions (the complex ones are left out) we see that there is no solution in the interval '(0,1]', so the condition we want to check is true. Next we will build the 'wronski_system' and solve it.

ideal > $coeffs = new Array<Array<Rational>>([[19,8,-19],[39,7,42]]);

ideal > $ws = wronski_system($M,$lambda,$coeffs,1,triangulation=>$T);

ideal > print $ws->SOLVE;
(1.471565025 0) (-1.474019364 0)
(-0.678417139 0) (-0.998334934 0)
(-1.001667843 0) (0.6795486323 0)
...

Again we left out the complex solutions. Here we've got exactly three real solutions, which is precisely the signature of the triangulation. Thus the bound of Soprunova and Sottile is met.

If you want to do your own tests, you just need random_solve.script. Here is a small example (in this example we solve 20000 systems):

polytope > $M = new Matrix<Int>([
polytope (2)> [1, 0, 0],
polytope (3)> [1, 1, 0],
polytope (4)> [1, 2, 0],
polytope (5)> [1, 3, 0],
polytope (6)> [1, 0, 1],
polytope (7)> [1, 1, 1],
polytope (8)> [1, 2, 1],
polytope (9)> [1, 0, 2],
polytope (10)> [1, 1, 2],
polytope (11)> [1, 0, 3]
polytope (12)> ]);

polytope > $lambda = new Vector<Int>([12, 8, 9, 15, 3, 1, 5, 0, 0, 0]);

polytope > $F = polytope::regular_subdivision($M, $lambda);
polymake: used package cddlib
  Implementation of the double description method of Motzkin et al.
  Copyright by Komei Fukuda.
  http://www.ifor.math.ethz.ch/~fukuda/cdd_home/cdd.html


polytope > $T = new topaz::GeometricSimplicialComplex(COORDINATES=>$M->minor(All,~[0]), FACETS=>$F);

polytope > solve_random_wronski($M,$lambda,20000,{file=>'mytest.csv', triangulation=>$T});

====================
==== Test Start ====
====================
polymake: used package Singular
  Singular is a computer algebra system for polynomial computations, 
  with special emphasis on commutative and non-commutative algebra, 
  algebraic geometry, and singularity theory. It is free and open-source 
  under the GNU General Public Licence.
  http://www.singular.uni-kl.de/

(7201) Seed: 7201 -- no wronski!  (expected: 9   got: 7)
(10000) Counter: 9999 -- No Wronski: 1
(15298) Seed: 15299 -- no wronski!  (expected: 9   got: 8)
(20000) Counter: 19998 -- No Wronski: 2
==== Test End ====

You get a notification every time the script found a system which is not a Wronski system (due to a bad random choice of the coefficients), and every 10000 successful iterations. After this test run the script created a .csv file in the directory of your choosing. In this case it is mytest.csv in the directory from where you started polymake. Here are few lines of this file:

"System";"Seed";"Parameter s=";"Number of real roots"
"1";"1";"0.69";"3"
"2";"2";"0.19";"9"
"3";"3";"0.68";"3"
"4";"4";"0.18";"9"
"5";"5";"0.67";"5"
"6";"6";"0.17";"7"
[...]

Now you can use your favorite spreadsheet software to analyze the data. Every column is separated with an semicolon (“;”). The first column is just a normal counter, which counts through all the systems. In the second column you see the seed which was used to create the system. The third column shows the parameter s which was used to solve the system. And the last column is the number of the real solutions of the system. If you choose to investigate a particular system more closely you look at the seed and the additional parameter. The following code shows the polynomials of the third system and its solutions:

polytope > $ws = random_wronski($M,$lambda,{seed=>3, param=>0.68, triangulation=>$T});

polytope > print join "\n\n", @{$ws->GENERATORS};
50712988728146283/144115188075855872 -125179915901626391/144115188075855872*x0 + 619533179139595/70368744177664*x1 + 31360863141670615/36028797018963968*x0^2 + 55124059439014875/2251799813685248*x0*x1 + 63783145887073983/576460752303423488*x0^3 -19*x1^2 -49764303530503557/18014398509481984*x0^2*x1 + 28*x0*x1^2 + 36*x1^3

5634776525349587/144115188075855872 -46118916384809723/144115188075855872*x0 -265514219631255/35184372088832*x1 -13440369917858835/18014398509481984*x0^2 + 6124895493223875/2251799813685248*x0*x1 + 7087016209674887/576460752303423488*x0^3 -7*x1^2 -18334217090185521/18014398509481984*x0^2*x1 -24*x0*x1^2 + 4*x1^3

polytope > print $ws->SOLVE;
(-24.61895356 0) (-0.6025803549 0)
(0.2492503671 0) (-0.01283577957 0)
(-0.5026977058 0) (1.007358785 0)
(0.05237435187 0.6635156581) (-0.183413567 0.4314160025)
(0.05237435187 -0.6635156581) (-0.183413567 -0.4314160025)
(5.406365466 -1.762730557) (-0.03037389905 0.4225105979)
(5.406365466 1.762730557) (-0.03037389905 -0.4225105979)
(-0.2950477646 -0.6816387146) (0.05130602948 -0.02018490937)
(-0.2950477646 0.6816387146) (0.05130602948 0.02018490937)

Every rows corresponds to one solution. Here you see the three real solutions on top and all the complex solutions below. For instance the first row corresponds to x0 = -24.61895356 and x1 = -0.6025803549. The last row corresponds to: x0 = -0.2950477646 + i*0.6816387146 and x1 = 0.05130602948 + i*0.02018490937.