-- PATCH for M2 version 1.7
needsPackage "NumericalAlgebraicGeometry"
solutionsWithMultiplicity List := o-> sols -> ( 
    sorted := sortSolutions(sols,o);
    i := 0; 
    while i<#sorted list (
	si := sorted#i;
	si.Multiplicity = 1;
	j := i + 1;
	while j < #sorted and areEqual(sorted#j,si,o) do (
	    si.Multiplicity = si.Multiplicity + 1;
	    j = j + 1;
	    );
	i = j;
	si
	) 
    )
end --------------------------------------------------------------------------

-- START pushing F11 here
restart
load "tensor-3x4x5-decompositions.m2"
---case of generic rank 6, 
R=CC[u_(1,1)..u_(6,2),v_(1,1)..v_(6,3),w_(1,1)..w_(6,5)] 
F = transpose sum(6,i->matrix{{u_(i+1,1),u_(i+1,2),1}}**matrix{{v_(i+1,1)..v_(i+1,3),1}}**matrix{{w_(i+1,1)..w_(i+1,5)}})

setRandomSeed 0
needsPackage "NumericalAlgebraicGeometry"

-- create a generic pair (s0,T0) = (decomposition,tensor)
s0 = random(CC^1,CC^60);
T0 = sub(F,s0)

-- we know one decomposition at the moment:
sols0 = {point s0};

{* outsource to an external solver:
setDefault(Software=>BERTINI)
setDefault(Software=>PHCPACK)
*}

--while #sols0 != 720 do (
while #sols0 != 721 do (
    T1 = random(CC^60,CC^1);
    T2 = random(CC^60,CC^1);
    elapsedTime sols1 = track(polySystem(F-T0),polySystem(F-T1),sols0);
    elapsedTime sols2 = track(polySystem(F-T1),polySystem(F-T2),sols1);
    elapsedTime sols0' = track(polySystem(F-T2),polySystem(F-T0),sols2);
    assert all(sols0', s->norm(T0 - sub(F,matrix s)) < 1e-6); -- check T0 = F(s)
    sols0 = solutionsWithMultiplicity(sols0 | sols0'); -- take the union
    << "found " << #sols0 << " decompositions so far" << endl;
    )

sols0
toList{x-2}

S=CC[x]
track(|x-2|,|x-1|,{2})