One Hat Cyber Team
Your IP :
216.73.216.115
Server IP :
194.44.31.54
Server :
Linux zen.imath.kiev.ua 4.18.0-553.77.1.el8_10.x86_64 #1 SMP Fri Oct 3 14:30:23 UTC 2025 x86_64
Server Software :
Apache/2.4.37 (Rocky Linux) OpenSSL/1.1.1k
PHP Version :
5.6.40
Buat File
|
Buat Folder
Eksekusi
Dir :
~
/
usr
/
share
/
Macaulay2
/
View File Name :
PHCpack.m2
phcPresent := run ("type phc >/dev/null 2>&1") === 0 phcVersion := if phcPresent then replace("PHCv([0-9.]+) .*\n","\\1",get "! phc --version") phcVersionNeeded := "2.4.77" phcPresentAndModern := phcPresent and match("^[0-9.]+$",phcVersion) and phcVersion >= phcVersionNeeded newPackage( "PHCpack", Version => "1.8", Date => "25 May 2016", Authors => { {Name => "Elizabeth Gross", Email => "egross7@uic.edu", HomePage => "http://www.math.uic.edu/~lizgross"}, {Name => "Sonja Petrovic", Email => "Sonja.Petrovic@iit.edu", HomePage => "http://mypages.iit.edu/~spetrov1/"}, {Name => "Jan Verschelde", Email => "jan@math.uic.edu", HomePage => "http://www.math.uic.edu/~jan"} }, Headline => "interface to PHCpack", Keywords => {"Numerical Algebraic Geometry", "Interfaces"}, Certification => { "journal name" => "The Journal of Software for Algebra and Geometry", "journal URI" => "https://msp.org/jsag/", "article title" => "Interfacing with PHCpack", "acceptance date" => "2013-07-11", "published article URI" => "https://msp.org/jsag/2013/5-1/p04.xhtml", "published article DOI" => "10.2140/jsag.2013.5.20", "published code URI" => "https://msp.org/jsag/2013/5-1/jsag-v5-n1-x04-code.zip", "release at publication" => "48d32ceffeeb807b23eee20a5c9c243ce085b3ef", "version at publication" => "1.6", "volume number" => "5", "volume URI" => "https://msp.org/jsag/2013/5-1/" }, DebuggingMode => false, AuxiliaryFiles => true, CacheExampleOutput => true, PackageExports => {"NAGtypes"}, OptionalComponentsPresent => phcPresentAndModern ) checkIsRunnable = () -> ( if not phcPresent then error "phc not present"; if not phcPresentAndModern then error ("phc present but not modern enough; need version ", phcVersionNeeded, " but found version ",phcVersion); ) --Copyright 2013 Elizabeth Gross, Sonja Petrovic, Jan Verschelde. -- You may redistribute this file under the terms of the GNU General -- Public License as published by the Free Software Foundation, -- either version 2 of the License, or any later version. export { "cascade", "computingPrecision", "constructEmbedding", "gamma", "factorWitnessSet", "interactive", "intermediateSolutions", "isCoordinateZero", "isWitnessSetMember", "loadSettingsPath", "mixedVolume", "nonZeroFilter", "numericalIrreducibleDecomposition", "numThreads", "randomSeed", "refineSolutions", "saveSettingsPath", "seeProgress", "solveRationalSystem", "solveSystem", "StableMixedVolume", "StartDimension", "StartSystem", "tDegree", "toLaurentPolynomial", "topWitnessSet", "trackPaths", "zeroFilter", "intersectSlice", "searchDelta", "searchNpoints", "searchTolerance", "realSlice1D", "realSlice2D", "versionNumber" } protect ErrorTolerance, protect Iterations, protect Bits, protect ResidualTolerance, protect Append --##########################################################################-- -- GLOBAL VARIABLES --##########################################################################-- PHCDBG = 0; -- debug level (10=keep temp files) -- We used to allow the user to set this in the "Configuration" of the package, but we need -- to know whether the program is present before "newPackage" runs, and thus there is no -- good way to get the option early enough. PHCexe = "phc "|(if member("--no-randomize",commandLine) then "-0123456 " else "") -- this is the executable string that make sures that calls to PHCpack run: -- NOTE: the absolute path should be put into the init-PHCpack.m2 file --##########################################################################-- -- INTERNAL METHODS --##########################################################################-- ---------------------------------- -- NumericalVariety ---------------------------------- --the following defines the type NumericalVariety --if the user is using an old version of NAGtypes. if not((class(NumericalVariety))===Type) then (--export {generalEquations, "IsIrreducible"}; --protect generalEquations; protect IsIrreducible; NumericalVariety = new Type of MutableHashTable; NumericalVariety.synonym = "numerical variety"; dim NumericalVariety := V -> max select(keys V, k->class k === ZZ); degree NumericalVariety := V -> ( d := dim V; sum(keys V, k->if k =!= d then 0 else sum(V#k,degree)) ); numericalVariety = method(TypicalValue=>NumericalVariety); numericalVariety List := Ws -> ( V := new NumericalVariety; scan(Ws, W->( d := dim W; if V#?d then V#d = V#d | {W} else V#d = {W}; )); check V; V ); check NumericalVariety := o-> V -> ( if any(keys V, k->(class k =!= ZZ or k<0)) then error "the keys of a NumericalVariety should be nonnegative integers"; scan(keys V, k->if class k === ZZ then scan(V#k, W->( if dim W != k then error "dimension of a witness set does not match the key in NumericalVariety"; ))); ); net NumericalVariety := V -> ( out := "A variety of dimension " | net dim V |" with components in"; scan(keys V, k->if class k === ZZ then ( row := "dim "|net k|": "; scan(V#k, W->row = row|" "|net W); out = out || row; )); out ) ) ---------------------------------- --- File read/write operations --- ---------------------------------- getFilename = () -> ( filename := temporaryFileName(); while fileExists(filename) or fileExists(filename|"PHCinput") or fileExists(filename|"PHCoutput") do filename = temporaryFileName(); filename ) parseSolutions = method(TypicalValue => Sequence, Options => {Bits => 53}) parseSolutions (String,Ring) := o -> (s,R) -> ( -- parses solutions in PHCpack format -- IN: s = string of solutions in PHCmaple format -- V = list of variable names -- OUT: List of solutions, each of type Point, -- carrying also other diagnostic information about each. oldprec := defaultPrecision; defaultPrecision = o.Bits; L := get s; L = replace("=", "=>", L); L = replace("I", "ii", L); L = replace("E\\+","e",L); L = replace("E", "e", L); L = replace("time", "\"time\"", L); L = replace("rco", "\"rco\"", L); L = replace("multiplicity", "\"mult\"", L); L = replace("\\bres\\b", "\"residual\"", L); L = replace("\\bresolution\\b", "\"residual\"", L); -- because M2 automatically thinks "res"=resolution sols := toList apply(value L, sol->new HashTable from toList sol); defaultPrecision = oldprec; apply(sols, sol->point( {apply(gens R, v->sol#v)} | outputToPoint sol )) ) ------------------------------------- -- The below method doesn't work unless you do simultaneous tracking, but a bug in the output of phc -p (option 2 for intermediate points) prevents that from working. ------------------------------------- -- parseIntermediateSolutionsOLD = method() -- parseIntermediateSolutionsOLD (String,ZZ,Ring) := (output,numsols,R) -> ( -- L := get output; -- rgx := "\\*{5} +path([[:digit:]]) +\\*+$"; -- start := (regex(rgx,L))#0#0; -- linesL := lines substring(start,L); -- solsize := 0; -- for i from 1 to #linesL-1 do ( -- if match("^\\*",linesL#i) then ( -- solsize = i; -- break; -- );); -- << "got solution size : " << solsize << endl; -- chunksize := numsols*solsize; -- chunks := {}; -- while match(rgx,linesL#0) do ( -- chunks = append(chunks,take(linesL,chunksize)); -- linesL = drop(linesL,chunksize); -- ); -- << #chunks << endl; -- << "first chunks "<< netList chunks#0 << endl << netList chunks#10 << endl; -- chunksNew := for chunk in chunks list ( -- for l in chunk list ( -- replace(rgx,"solution \\1 :",l); -- ); -- ); -- chunks = chunksNew; -- << #chunks << endl; -- << "first chunk again" << netList chunks#0 << endl; -- for chunk in chunks list ( -- --fname := getFilename(); -- fname := "this.that"; -- f := openOut fname; -- for sol in chunk do ( -- for l in sol do ( -- f << toString(l) << endl; -- ); -- ); -- close f; -- foutname := fname | ".sols"; -- if checkIsRunnable() then run(PHCexe|" -z "|fname|" "|foutname); -- parseSolutions(get foutname,R) -- ) -- ) parseIntermediateSolutions = method() parseIntermediateSolutions (String,Ring) := (output,R) -> ( linesL := lines get output; start := 0; (loc,len) := (regex(".*([[:digit:]])+$",linesL#0))#1; numvars := value(substring(loc,loc+len,linesL#0)); for i from 0 to #linesL-1 do ( if "OUTPUT INFORMATION DURING CONTINUATION :" == linesL#i then ( start = i; break; ); ); solsize := 0; for i from start+3 to #linesL-1 do ( if match("^== err ",linesL#i) then ( solsize = i-start-3 + 1; break; ); ); << "got solution size : " << solsize << endl; filename := getFilename(); << filename << endl; f := openOut filename; i := start+3; parsefiles := {}; numlines := 0; while(i<#linesL-1) do ( line := linesL#i; if match("^t :", line) then (f << "solution 1 : " << endl;); if match("== [[:digit:]] ", line) then ( i = i + solsize + 1; close f; foutname := filename | ".sols"; parsefiles = append(parsefiles, (filename,numlines // solsize)); filename = getFilename(); f = openOut filename; numlines=0; continue; ); f << line << endl; i = i + 1; numlines = numlines + 1; ); close f; results := {}; for pf in parsefiles do ( oldf := openIn pf#0; f = openOut(pf#0 | ".final"); f << "THE SOLUTIONS :" << endl; f << pf#1 << " " << numvars << endl; f << "===========================================================================" << endl; f << get oldf; close f; checkIsRunnable(); run(PHCexe|"-z "|pf#0|".final "|pf#0|".sols"); results = append(results,parseSolutions(pf#0|".sols",R)); ); results ) pointsToFile = method(TypicalValue => Nothing, Options => {Append => false}) pointsToFile (List,Ring,String) := o -> (S,R,name) -> ( -- writes list of points to file in PHCpack format -- IN: S, list of points, e.g.: obtained as points(WitnessSet); -- R, ring (with symbols for the variables); -- name, string with file name. file := if o.Append then openOutAppend name else openOut name; if o.Append then file << endl << "THE SOLUTIONS :" << endl; file << #S << " " << numgens R << endl << "===========================================================" << endl; scan(#S, i->( file << "solution " << i+1 << " :" << endl << "t : 0.00000000000000E+00 0.00000000000000E+00" << endl << "m : 1" << endl << "the solution for t :" << endl; scan(numgens R, v->( L := " "|toString R_v|" : "| format(0,-1,9,9,realPart toCC S#i#v)|" "| format(0,-1,9,9,imaginaryPart toCC S#i#v); file << L << endl; )); file << "== err : 0 = rco : 1 = res : 0 ==" << endl; )); close file; ) solutionsToFile = method(TypicalValue => Nothing, Options => {Append => false}) solutionsToFile (List,Ring,String) := o -> (S,R,name) -> ( -- writes solutions to file in PHCpack format -- IN: S, list of solutions; -- R, ring (with symbols for the variables); -- name, string with file name. file := if o.Append then openOutAppend name else openOut name; if o.Append then file << endl << "THE SOLUTIONS :" << endl; file << #S << " " << numgens R << endl << "===========================================================" << endl; scan(#S, i->( file << "solution " << i << " :" << endl << "t : 0.00000000000000E+00 0.00000000000000E+00" << endl << "m : 1" << endl << "the solution for t :" << endl; scan(numgens R, v->( L := " "|toString R_v|" : "| format(0,-1,9,9,realPart toCC S#i#Coordinates#v)|" "| format(0,-1,9,9,imaginaryPart toCC S#i#Coordinates#v); file << L << endl; )); file << "== err : 0 = rco : 1 = res : 0 ==" << endl; )); close file; ) startSystemFromFile = method(TypicalValue => List) startSystemFromFile (String) := (name) -> ( -- IN: file name starting with a random coefficient start system -- OUT: list of polynomials in a ring with coefficients in CC -- REQUIRED: the format of the system on file is a random coefficient -- system produced by phc -m, every term starts on a separate line -- with the "+" sign. The coefficient field must be CC. s := get name; s = replace("i","ii",s); s = replace("E","e",s); s = replace("e\\+00","",s); L := lines(s); n := value L_0; result := {}; i := 0; j := 1; local stop; local term; local p; while i < n do ( stop = false; p = 0; while not stop do ( if #L_j != 0 then ( -- we have to bite off the first "+" sign of the term if (L_j_(#L_j-1) != ";") then ( term = value substring(1,#L_j-1,L_j); p = p + term; ) else ( -- in this case (L_j_(#L_j-1) == ";") holds term = value substring(1,#L_j-2,L_j); p = p + term; stop = true; result = result | {p} ); ); j = j + 1; stop = stop or (j >= #L); ); i = i + 1; ); result ) systemFromFile = method(TypicalValue => List) systemFromFile (String) := (name) -> ( -- IN: file name starting with a polynomial system -- OUT: list of polynomials in a ring with coefficients in CC -- NOTE: the format of the system on file may be such that the -- first term starts with +1*x which cannot be digested by M2. -- In contrast to the startSystemFromFile, the first term could -- also be "-1*x" so we must be a bit more careful... -- Another problem are constants as "3.0e+00", at least on a Mac. s := get name; s = replace("i","ii",s); s = replace("E","e",s); s = replace("e\\+","e",s); -- M2 does not like 3.0e+00 as constant L := lines(s); dimL0 := separate(" ", L_0); -- deal with case of nonsquare systems n := value dimL0_0; -- first is always number of equations if not instance(n,ZZ) then n = value dimL0_1; -- deal with leading spaces result := {}; i := 0; j := 1; local stop; local term; local p; while i < n do ( stop = false; p = 0; while not stop do ( if #L_j != 0 then ( if (L_j_(#L_j-1) != ";") then ( -- we have to bite off the first "+" sign of the term -- but check if the first character is a plus! if (L_j_0 == "+") or (L_j_0 == "-") then term = value substring(1,#L_j, L_j) -- substring(1,#L_j-1,L_j) else term = value substring(0,#L_j, L_j); -- substring(0,#L_j-1,L_j); if (L_j_0 == "+") then p = p + term else p = p - term; ) else ( -- in this case (L_j_(#L_j-1) == ";") holds if (L_j_0 == "+") or (L_j_0 == "-") then term = value substring(1,#L_j-1,L_j) -- substring(1,#L_j-2,L_j) else term = value substring(0,#L_j-1,L_j); -- substring(0,#L_j-2,L_j) if (#L_j > 1) then -- take care of a lonely ";" ( if (L_j_0 == "+") then p = p + term else p = p - term; ); stop = true; result = result | {p} ); ); j = j + 1; stop = stop or (j >= #L); ); i = i + 1; ); result ) systemToFile = method(TypicalValue => Nothing) systemToFile (List,String) := (F,name) -> ( file := openOut name; if (#F == numgens ring F#0) then file << #F << endl else file << #F << " " << numgens ring F#0 << endl; scan(F, f->( L := toExternalString f; L = replace("ii", "I", L); L = replace("e", "E", L); L = replace("p53","",L); L = replace("p[0-9]+","",L); file << L << ";" << endl; )); close file; ) witnessSetFromFile = method(TypicalValue => WitnessSet) witnessSetFromFile (String) := (name) -> ( -- IN: file name which contains a witness set in PHCpack format -- OUT: a witness set e := systemFromFile(name); d := dimEmbedding(e); witnessPointsFile := temporaryFileName() | "PHCwitnessPoints"; if fileExists witnessPointsFile then removeFile witnessPointsFile; checkIsRunnable(); run(PHCexe|" -z " | name | " "|witnessPointsFile); eR := ring first e; g := parseSolutions(witnessPointsFile,eR); w := witnessSet(ideal(take(e,{0,#e-d-1})),ideal(take(e,{#e-d,#e-1})),g); return w; ) witnessSetToFile = method() witnessSetToFile (WitnessSet,String) := (witset,name) -> ( -- IN: name, a file name where to write a witness set to, -- witset, a witness set. R := ring ideal(witset); s := equations(witset)|slice(witset); p := points(witset); systemToFile(s,name); pointsToFile(p/coordinates,R,name,Append=>true); ) ---------------------------- -- embedding dimension ----- ---------------------------- dimEmbedding = method(TypicalValue => ZZ) dimEmbedding (List) := (system) -> ( -- IN: embedded system with slack variables. -- OUT: returns the number of slack variables = the dimension. eR := ring first system; v := gens eR; slack := v_(#v-1); -- slack is the last variable zz := toString(slack); -- zz is the name of the last variable if substring(0,2,zz) != "zz" then -- check if slack starts with zz ( return 0; ) else ( ds := substring(2,#zz-1,zz); dimension := if (value(ds)===null) then 0 else value(ds); ); return dimension; ) ----------------------------- --- conversion to Point --- ----------------------------- outputToPoint = method() outputToPoint HashTable := (H)->{ SolutionStatus => if H#"mult" == 1 then Regular else Singular, ConditionNumber => (H#"rco")^(-1), LastT => H#"time" } ---------------------------------- ----add slack variables----------- ---------------------------------- addSlackVars = method() addSlackVars WitnessSet := (W) -> ( -- creates a new system of polynomials, in variables: -- old set of variables, and zz1, ..., zzd, where -- d is the dimension of W. R := ring W; n := numgens R; d := dim W; -- this will be the number of slack variables to add W1 := generalEquations W; -- Add in new variables zz1, ..., zzd, -- this changes the equations, the slice, and the points slackvars := apply(d, i->getSymbol("zz"|toString (i+1))); newR := (coefficientRing R)[gens R, slackvars]; newvars := (vars newR)_{n..n+d-1}; -- new slice: newSlice := apply(d, i -> sub(W1.Slice#i,newR) + newR_(n + i)); -- add a linear matrix A := random(newR^(d),newR^(n-d)); AZ := newvars * A; newEqns := (sub(gens ideal W1, newR) + AZ) | newvars; -- new points zeros := toList apply(d, i -> 0_(coefficientRing R)); newPoints := apply(W1.Points, pt -> point({join(coordinates(pt),zeros)})); witnessSet(ideal newEqns, ideal newSlice, newPoints) ) ------------------------------------------------------- ------------ Witness Superset Filters --------------- ------------------------------------------------------- witnessSuperSetFilter = method() witnessSuperSetFilter (WitnessSet,List) := (witset,pts) -> ( -- This is an auxiliary procedure to witnessSuperSetsFilter, -- applying isWitnessSetMember to filter points in a witness superset. -- IN: witset, a witness set to represent a witness set; -- pts, a list of points to be tested for membership. -- OUT: a list of points in pts for which isWitnessSetMember is false. result := new MutableList from {}; for p in pts do ( if not isWitnessSetMember(witset,p) then result = append(result,p); ); return result; ) witnessSuperSetsFilter = method() witnessSuperSetsFilter (MutableList,List) := (witsets,pts) -> ( -- This is an auxiliary procedure to cascade, -- applying witnessSuperFilter to filter points in the -- witness supersets given in the list. -- IN: witsets, a list of tuples (dimension,witness set), -- pts, a list of points to be tested for membership. -- OUT: a list of points in pts for which isWitnessSetMember is false. result := new MutableList from {}; local c; for p in pts do ( found := false; if instance(p,AbstractPoint) then c = p else c = point{p}; for w in witsets when (not found) do ( ws := w#1; found = isWitnessSetMember(ws,c); ); if not found then result = append(result,p); ); return toList(result); ) --##########################################################################-- -- EXPORTED METHODS -- -- NOTE: -- trackPaths and refineSolutions are methods adapted -- from NumericalAlgebraicGEometry/PHCpack.interface.m2 --##########################################################################-- -------------- -- CASCADE -- -------------- cascade = method(TypicalValue => NumericalVariety, Options => {StartDimension => -1,Verbose => false}) cascade (List) := o -> (system) -> ( -- IN: system, a polynomial system; -- dimension, top dimension of the solution set. -- OUT: a hash table with keys the dimension of each component, -- values are witness sets for positive dimensions, -- or a list of isolated solutions for key equal to zero. if not(class coefficientRing ring ideal system===ComplexField) then error "coefficient ring is not complex"; R := ring ideal system; if # system > numgens R then error "the system is overdetermined"; if o.StartDimension==-1 then startdim := (numgens R)-1 else startdim = o.StartDimension; PHCinputFile := temporaryFileName() | "PHCinput"; PHCoutputFile := temporaryFileName() | "PHCoutput"; PHCbatchFile := temporaryFileName() | "PHCbatch"; PHCsolsFile := temporaryFileName() | "PHCsols"; PHCsessionFile := temporaryFileName() | "PHCsession"; for f in {PHCinputFile, PHCoutputFile, PHCbatchFile, PHCsolsFile, PHCsessionFile} do if fileExists f then removeFile f; toList (0..startdim) / (i->if fileExists (PHCoutputFile | "_sw" | i) then removeFile (PHCoutputFile | "_sw" | i) ); if o.Verbose then stdio << "writing output to file " << PHCoutputFile << endl; systemToFile(system,PHCinputFile); s := concatenate("0\n0\ny\n",PHCinputFile); -- option 0, system on file -- extra 0 for double precision s = concatenate(s,"\n",PHCoutputFile); -- add name of the output file s = concatenate(s,"\n"); s = concatenate(s,toString(startdim)); -- add top dimension s = concatenate(s,"\nn\n"); -- do not restrict slices bat := openOut PHCbatchFile; bat << s; close bat; if o.Verbose then ( stdio << "calling phc -c < " << PHCbatchFile; stdio << " > " << PHCsessionFile << endl ); checkIsRunnable(); run(PHCexe|" -c < " | PHCbatchFile | " > " | PHCsessionFile); if o.Verbose then ( stdio << "output of phc -c is in file " << PHCoutputFile << endl; stdio << "... constructing witness sets ... " << endl ); local slackvars; local RwithSlack; --get solutions result := new MutableList from {}; dims:=select(toList (0..startdim),j->(fileExists (PHCoutputFile | "_sw" | j) and match("THE SOLUTIONS",get (PHCoutputFile | "_sw" | j)))); topdimension := max dims; if o.Verbose then stdio << "the top dimension is " << topdimension << endl; i := topdimension; while i>=0 do ( if member(i,dims) then ( fil := (PHCoutputFile | "_sw" | i); if o.Verbose then stdio << "processing super witness set file " << fil << endl; if i > 0 then ( slackvars = apply(i, k->getSymbol("zz"|toString(k+1))); RwithSlack = (coefficientRing R)monoid(gens R | slackvars); use RwithSlack; supwit := witnessSetFromFile(fil); if i == topdimension then ( result = append(result,(i,supwit)); ) else ( supsols := points(supwit); genpts := witnessSuperSetsFilter(result,supsols); g := toList(apply(genpts,x->if instance(x,AbstractPoint) then point x else point{x})); ws := witnessSet(ideal(equations(supwit)),ideal(slice(supwit)),g); if #g!=0 then result = append(result,(i,ws)); ); ) else ( checkIsRunnable(); run(PHCexe | " -z " | fil | " " | PHCsolsFile); use R; supwit = witnessSetFromFile(fil); if o.Verbose then stdio << "the super witness set from file :\n" << supwit << endl; psols := parseSolutions(PHCsolsFile,R); isols := witnessSuperSetsFilter(result,psols); ws = witnessSet(ideal(equations(supwit)),ideal(slice(supwit)),isols); if #isols!=0 then result = append(result,(0,ws)); ); ); i = i-1 ); use R; numericalVariety(toList (apply(result,i-> last i))) ) ------------------------- -- CONSTRUCT EMBEDDING -- ------------------------- constructEmbedding = method(TypicalValue => List,Options => {Verbose => false}) constructEmbedding (List, ZZ) := o-> (system, dimension) -> ( -- IN: system, a polynomial system with complex coefficients; -- dimension, expected dimension of the solution set. -- OUT: system with as many random hyperplanes at the end -- as the value of dimension. if not(class coefficientRing ring first system===ComplexField) then error "coefficient ring of system is not complex"; PHCinputFile := temporaryFileName() | "PHCinput"; PHCoutputFile := temporaryFileName() | "PHCoutput"; PHCbatchFile := temporaryFileName() | "PHCbatch"; PHCsessionFile := temporaryFileName() | "PHCsession"; for f in {PHCinputFile, PHCoutputFile, PHCbatchFile, PHCsessionFile} do if fileExists f then removeFile f; systemToFile(system,PHCinputFile); s := concatenate("1\n0\ny\n",PHCinputFile); s = concatenate(s,"\n",PHCoutputFile); s = concatenate(s,"\n"); s = concatenate(s,toString(dimension)); s = concatenate(s,"\nn\n"); bat := openOut PHCbatchFile; bat << s; close bat; if o.Verbose then ( stdio << "calling phc -c < " << PHCbatchFile; stdio << " > " << PHCsessionFile << endl ); checkIsRunnable(); run(PHCexe|" -c < " | PHCbatchFile | " > " | PHCsessionFile); if o.Verbose then stdio << "output of phc -c is in file " << PHCoutputFile << endl; -- extending the ring with slack variables slackvars := apply(dimension, i->getSymbol("zz"|toString(i+1))); R := ring ideal system; -- surplus variables are used when the initial system is overconstrained nv := numgens R; nq := #system; nbss := nq - nv; local RwithSlack; if (nbss > 0) then ( surplusvars := apply(nbss, i->getSymbol("ss"|toString(i+1))); RwithSlack = (coefficientRing R)[gens R, surplusvars, slackvars]; ) else ( RwithSlack = (coefficientRing R)[gens R, slackvars]; ); use RwithSlack; return systemFromFile(PHCoutputFile); ) ------------------------ -- FACTOR WITNESS SET -- ------------------------ factorWitnessSet = method(TypicalValue=>List, Options => {Verbose => false}) factorWitnessSet (WitnessSet ) := o-> w -> ( -- IN: a witness set properly embedded with slack variables. -- OUT: a list of witness sets, every element is irreducible. PHCinputFile := temporaryFileName() | "PHCinput"; PHCoutputFile := temporaryFileName() | "PHCoutput"; PHCbatchFile := temporaryFileName() | "PHCbatch"; PHCsessionFile := temporaryFileName() | "PHCsession"; if o.Verbose then stdio << "preparing input file to " << PHCinputFile << endl; system := equations(w) | slice(w); systemToFile(system,PHCinputFile); R := ring first system; use R; L := toList(points(w)) / coordinates; pointsToFile(L,R,PHCinputFile,Append=>true); if o.Verbose then stdio << "preparing batch file to " << PHCbatchFile << endl; s := concatenate("2\n",PHCinputFile); -- option 2 of phc -f s = concatenate(s,"\n",PHCoutputFile); s = concatenate(s,"\n1\n"); -- use monodromy to factor s = concatenate(s,"0\n"); -- default settings of path trackers bat := openOut PHCbatchFile; bat << s; close bat; if o.Verbose then stdio << "... calling monodromy breakup ..." << endl; checkIsRunnable(); run(PHCexe|" -f < " | PHCbatchFile | " > " | PHCsessionFile); if o.Verbose then (stdio << "session information of phc -f is in " << PHCsessionFile << endl; stdio << "output of phc -f is in file " << PHCoutputFile << endl ); -- counting the number of factors count := 0; result := new MutableList from {}; name := PHCinputFile | "_f" | toString(count+1); while (fileExists name) do ( result = append(result,witnessSetFromFile(name)); count = count + 1; name = PHCinputFile | "_f" | toString(count+1); ); stdio << "found " << count << " irreducible factors " << endl; for ws in result do ws.cache.IsIrreducible=true; return numericalVariety(toList(result)); ) ------------------------ -- IS COORDINATE ZERO -- ------------------------ isCoordinateZero = method(TypicalValue => Boolean) isCoordinateZero (AbstractPoint,ZZ,RR) := (sol,k,tol) -> ( -- IN: sol, a solution of a polynomial system; -- k, index to a coordinate must be within range; -- tol, tolerance for the absolute value of the k-th coordinate. -- OUT: true, if the k-th coordinate has abs less than tol; -- false, otherwise. L := sol#Coordinates; return abs(L_k)<=tol; ) --------------------------- -- IS WITNESS SET MEMBER -- --------------------------- isWitnessSetMember = method(TypicalValue => Boolean, Options => {Verbose => false}) isWitnessSetMember (WitnessSet,AbstractPoint) := o-> (witset,testpoint) -> ( -- IN: witset, a witness set for a positive dimensional solution set, -- testpoint, does it belong to the solution set? -- OUT: true if testpoint is a member of the solution set, -- false otherwise. PHCwitnessFile := temporaryFileName() | "PHCwitset"; PHCtestpointFile := temporaryFileName() | "PHCtestpoint"; PHCbatchFile := temporaryFileName() | "PHCbatch"; PHCoutputFile := temporaryFileName() | "PHCoutput"; PHCsessionFile := temporaryFileName() | "PHCsession"; if o.Verbose then stdio << "writing witness set to file " << PHCwitnessFile << endl; witnessSetToFile(witset,PHCwitnessFile); d := dim(witset); -- pad test point with zero values for slack variables dzeros := toList(apply(0..d,i->0)); L := {coordinates(testpoint)|dzeros}; R := ring ideal(witset); if o.Verbose then stdio << "writing test point to file " << PHCtestpointFile << endl; pointsToFile(L,R,PHCtestpointFile); s := concatenate("1\n0\n",PHCwitnessFile); s = concatenate(s,"\n",PHCtestpointFile); s = concatenate(s,"\n",PHCoutputFile); s = concatenate(s,"\n0\n"); bat := openOut PHCbatchFile; bat << s; close bat; if o.Verbose then ( stdio << "calling phc -f < " << PHCbatchFile; stdio << " > " << PHCsessionFile << endl; ); checkIsRunnable(); run(PHCexe|" -f < " | PHCbatchFile | " > " | PHCsessionFile); if o.Verbose then stdio << "output of phc -f is in file " << PHCoutputFile << endl; -- if the point does not belong to the witness set, -- then the output file contains the word "not" r := get PHCoutputFile; result := not match("not",r); return result; ) ------------------ -- MIXED VOLUME -- ------------------ mixedVolume = method(Options => {StableMixedVolume => false, StartSystem => false, Verbose => false, numThreads => 0, interactive=>false}) mixedVolume List := Sequence => opt -> system -> ( -- IN: system = list of polynomials in the system -- OUT: mixed volume of the system. if optional inputs specified, then output is -- a sequence containing a subset of: mixed volume, stable mixed volume, -- start system, and solutions to start system. -- Calls an Ada translation of ACM TOMS Algorithm 846: -- "MixedVol: a software package for mixed-volume computation" -- by Tangan Gao, T. Y. Li, Mengnien Wu, ACM TOMS 31(4):555-560, 2005. -- With the introduction of double double and quad double arithmetic, -- the menu options after version 2.3.90 changed. -- Fixed in the distribution of 2.3.97 of PHCpack. R := ring ideal system; n := #system; if n < numgens R then error "the system is underdetermined"; if n > numgens R then error "the system is overdetermined"; if not(class coefficientRing R===ComplexField) then error "coefficient ring is not complex"; filename := getFilename(); if opt.Verbose then stdio << "using temporary files " << filename|"PHCinput" << " and " << filename|"PHCoutput" << endl; infile := filename|"PHCinput"; outfile := filename|"PHCoutput"; cmdfile := filename|"PHCcommands"; sesfile := filename|"PHCsession"; startfile := filename|"PHCstart"; -- writing data to the corresponding files file := openOut cmdfile; file << "4" << endl; -- call MixedVol in PHCpack if opt.StartSystem then (file << "1" << endl) -- random coefficient start system wanted else (file << "0" << endl); -- no random coefficient start system if opt.StableMixedVolume then (file << "y" << endl) -- stable mixed volume wanted else (file << "n" << endl); -- no stable mixed volume file << "n" << endl; -- no mixed-cell configuration on file if opt.StartSystem then ( -- file and options for start system file << startfile << endl; file << "0" << endl << "1" << endl; ); close file; systemToFile(system,infile); if opt.interactive then ( << endl << "If you need a start system, the filename MUST be " << endl << endl << startfile << endl << endl << endl; checkIsRunnable(); run(PHCexe|" -m "|infile|" "|outfile); ) else ( -- launching mixed volume calculator : execstr := PHCexe|" -m "|(if opt.numThreads > 1 then ("-t"|opt.numThreads|" ") else "")|infile|" "|outfile|" < "|cmdfile|" > "|sesfile; checkIsRunnable(); ret := run(execstr); if ret =!= 0 then error "error occurred while executing PHCpack command: phc -m"; ); F := get outfile; -- search lines of outfile for: " mixed volume : " -- once found, extract just the number and return its value: local mixvol; scanLines(line -> if substring(0,21,line) == "common mixed volume :" then ( mixvol = value replace("common mixed volume : ","",line); break ), outfile); local stabmv; if opt.interactive or opt.StableMixedVolume then ( scanLines(line -> if substring(0,21,line) == "stable mixed volume :" then ( stabmv = value replace("stable mixed volume : ","",line); break ), outfile); ); local result; if not fileExists(startfile) then ( if opt.StableMixedVolume or (opt.interactive and class(stabmv) =!= Nothing) then result = (mixvol, stabmv) else result = mixvol; ) else ( solsfile := startfile | ".sols"; p := startSystemFromFile(startfile); execstr = PHCexe|" -z "|startfile|" "|solsfile; checkIsRunnable(); ret = run(execstr); if ret =!= 0 then error "error occurred while executing PHCpack command: phc -z"; sols := parseSolutions(solsfile, ring ideal system); if class(stabmv)=!=Nothing then result = (mixvol,stabmv,p,sols) else result = (mixvol,p,sols); ); if opt.interactive then result = (mixvol,stabmv,p,sols); result ) --------------------- -- NON ZERO FILTER -- --------------------- nonZeroFilter = method(TypicalValue => List) nonZeroFilter (List,ZZ,RR) := (sols,k,tol) -> ( -- IN: sols, solutions of a polynomial system; -- k, index to a coordinate of a solution -- tol, tolerance for the absolute value of k-th coordinate -- OUT: list of solutions in sols where the k-th coordinate -- is less than the tolerance. return select(sols,t->(not isCoordinateZero(t,k,tol))); ) ----------------------------------------- -- NUMERICAL IRREDUCIBLE DECOMPOSITION -- ----------------------------------------- numericalIrreducibleDecomposition=method(TypicalValue=>NumericalVariety, Options=>{StartDimension=>-1, Verbose=>false}) numericalIrreducibleDecomposition (List) := o -> (L) -> ( --IN: an ideal, top dimension --OUT: a NumericalVariety setRandomSeed(random ZZ); if o.Verbose then stdio << "starting cascade of homotopies ..." << endl; W := cascade(L, StartDimension=>o.StartDimension, Verbose=>o.Verbose); witsets := apply(keys W, i->if i!=0 then (factorWitnessSet((W#i)_0, Verbose=>o.Verbose))#i else W#i); numericalVariety(flatten witsets) ) ------------------------ -- REFINING SOLUTIONS -- ------------------------ refineSolutions = method(TypicalValue=>List, Options => {Verbose => false}) refineSolutions (List,List,ZZ) := o-> (f,sols,dp) -> ( -- IN: f, a polynomial system with complex coefficients; -- sols, list of points, initial approximations for the solutions of f; -- dp, number of decimal places in the working precision. -- OUT: a list of refined solutions. if not(class coefficientRing ring first f===ComplexField) then error "coefficient ring of system is not complex"; PHCinputFile := temporaryFileName() | "PHCinput"; PHCoutputFile := temporaryFileName() | "PHCoutput"; PHCbatchFile := temporaryFileName() | "PHCbatch"; PHCsessionFile := temporaryFileName() | "PHCsession"; PHCsolutions := temporaryFileName() | "PHCsolutions"; if o.Verbose then stdio << "writing input system to " << PHCinputFile << endl; systemToFile(f,PHCinputFile); if o.Verbose then stdio << "appending solutions to " << PHCinputFile << endl; R := ring first f; solutionsToFile(sols, R, PHCinputFile, Append=>true); if o.Verbose then stdio << "preparing input data for phc -v in " << PHCbatchFile << endl; s := concatenate("3\n",PHCinputFile); s = concatenate(s,"\n"); s = concatenate(s,PHCoutputFile); s = concatenate(s,"\n3\n1.0E-"); s = concatenate(s,toString(dp-4)); -- tolerance for correction term s = concatenate(s,"\n4\n1.0E-"); s = concatenate(s,toString(dp+16)); -- tolerance for residual s = concatenate(s,"\n6\n"); nit := ceiling(dp/10.0); s = concatenate(s,toString(nit)); -- number of Newton iterations s = concatenate(s,"\n7\n"); s = concatenate(s,toString(dp)); -- decimal places in working precision s = concatenate(s,"\n0\n"); bat := openOut PHCbatchFile; bat << s; close bat; -- stdio << "running phc -v, writing output to " << PHCsessionFile << endl; checkIsRunnable(); run(PHCexe|" -v < " | PHCbatchFile | " > " | PHCsessionFile); if o.Verbose then ( stdio << "using temporary file " << PHCoutputFile; stdio << " for storing refined solutions " << endl; stdio << "solutions in Maple format in " << PHCsolutions << endl ); checkIsRunnable(); run(PHCexe|" -z " | PHCoutputFile | " " | PHCsolutions); b := ceiling(log_2(10^dp)); result := parseSolutions(PHCsolutions,R,Bits=>b); result ) ------------------ -- SOLVE SYSTEM -- ------------------ solveSystem = method(TypicalValue => List, Options => {Verbose => false, numThreads=>0, randomSeed => -1, computingPrecision => 1}) solveSystem List := List => o->system -> ( -- IN: system = list of polynomials with complex coefficients, -- i.e. the system to solved -- OUT: solutions to the system, a list of Points -- fixed removing of nonzero slack variables -- for overdetermined systems (JV 2015/05/27) if instance(ring ideal system, FractionField) then error "ring is a fraction field, use solveRationalSystem"; if not(class coefficientRing ring first system===ComplexField) then error "coefficient ring is not complex"; if not member(o.computingPrecision,{1,2,4}) then error "Precision must be set to 1, 2, or 4."; filename := getFilename(); if o.Verbose then stdio << "using temporary files " << filename|"PHCinput" << " and " << filename|"PHCoutput" << endl; infile := filename|"PHCinput"; outfile := filename|"PHCoutput"; solnsfile := filename|"PHCsolns"; local newR; R := ring ideal system; -- stdio << "*** variables in the ring : " << gens R << " ***" << endl; n := #system; if n < numgens R then error "the system is underdetermined, positive dimensional"; -- add slack variables if needed (i.e. if system is overdetermined) if n > numgens R then ( nSlacks := n - numgens R; if o.Verbose then stdio << "adding " << nSlacks << " slack variables to overdetermined system" << endl; slackVars := apply(nSlacks, i->getSymbol("OOOO"|toString i)); newR = (coefficientRing R)(gens R | slackVars); rM := random(CC^n,CC^nSlacks); system = apply(#system, i->sub(system#i,newR) +(rM^{i}*transpose submatrix'(vars newR,toList(0..numgens R - 1)))_(0,0)) ) else newR=R; -- needed for parsing the solutions -- writing data to the corresponding files: systemToFile(system,infile); -- launching blackbox solver: execstr := PHCexe|" -b" |(if o.computingPrecision == 2 then "2 " else if o.computingPrecision == 4 then "4 " else " ") |(if o.numThreads > 1 then ("-t"|o.numThreads|" ") else "") |(if o.randomSeed > -1 then ("-0"|o.randomSeed|" ") else "") |infile|" "|outfile; checkIsRunnable(); ret := run(execstr); if ret =!= 0 then error "error occurred while executing PHCpack command: phc -b"; if o.Verbose then stdio << "solutions are in the file " << solnsfile << endl; execstr = PHCexe|" -z "|infile|" " |solnsfile; checkIsRunnable(); ret = run(execstr); if ret =!= 0 then error "error occurred while executing PHCpack command: phc -z"; -- parse and output the solutions: local result; if n == numgens R then ( result = parseSolutions(solnsfile, R) ) else ( slackRing := (coefficientRing R)(gens R | slackVars); result = parseSolutions(solnsfile, slackRing); stdio << "*** after parseSolutions, ring has " << gens R << " ***" << endl; if o.Verbose then stdio << "computed " << #result << " solutions of system with slack variables" << endl; result = zeroFilter(result, n-1, 1.0e-8); if o.Verbose then stdio << "after filtering nonsolutions : " << #result << " solutions left" << endl; result = apply(result, sol -> point(take(coordinates sol, numgens R), sol.cache)); newR = (coefficientRing R)(gens R) -- put variables back in original ring ); result ) --------------------------- -- SOLVE RATIONAL SYSTEM -- --------------------------- solveRationalSystem = method(TypicalValue => List, Options => {Verbose => false}) solveRationalSystem List := o-> system -> ( -- IN: system = list of rational equations with complex coefficients, -- i.e. the system to solved -- OUT: solutions to the system, a list of Points origRing := ring ideal system; --convert to laurent polynomial if instance(ring ideal system, FractionField) then system = toLaurentPolynomial(system, getSymbol "LLLL"); filename := getFilename(); if o.Verbose then stdio << "using temporary files " << filename|"PHCinput" << " and " << filename|"PHCoutput" << endl; infile := filename|"PHCinput"; outfile := filename|"PHCoutput"; solnsfile := filename|"PHCsolns"; R := ring ideal system; vars R; n := #system; if n < numgens R then error "the system is underdetermined, positive dimensional"; -- add slack variables if needed (i.e. if system is overdetermined) if n > numgens R then ( nSlacks := n - numgens R; slackVars := apply(nSlacks, i->getSymbol("OOOO"|toString i)); newR := QQ(monoid(gens R | slackVars, Inverses=>true, MonomialOrder=>RevLex)); rM := random(QQ^n,QQ^nSlacks); system = apply(#system, i->sub(system#i,newR) +(rM^{i}*transpose submatrix'(vars newR,toList(0..numgens R - 1)))_(0,0)) ) else newR=R; -- writing data to the corresponding files: systemToFile(system,infile); -- launching blackbox solver: execstr := PHCexe|" -b " |infile|" "|outfile; checkIsRunnable(); ret := run(execstr); if ret =!= 0 then error "error occurred while executing PHCpack command: phc -b"; execstr = PHCexe|" -z "|infile|" " |solnsfile; checkIsRunnable(); ret = run(execstr); if ret =!= 0 then error "error occurred while executing PHCpack command: phc -z"; -- parse and output the solutions: result := parseSolutions(solnsfile, QQ[gens origRing]); result ) --------------------------- -- TO LAURENT POLYNOMIAL -- --------------------------- toLaurentPolynomial = method(TypicalValue => List) toLaurentPolynomial (List, Symbol) := (system, var) -> ( -- IN: system (or a polynomial) which is rational -- (i.e. lives in some field of fractions of a polynomial ring) -- and a symbol to be used for new indexed variables (e.g. x) -- OUT: same system converted to a Laurent polynomial system, -- where denominators are replaced with new indexed variables. R := ring ideal system; P := R.baseRings_(#R.baseRings-1); -- P is the polynomial ring whose field of fractions the system lives in counter := 0; scan(system, f-> ( if instance(class f, FractionField) then --if f is already polynomial, do not do anything! -- if liftable(f,P) then --* --if it can be lifted to P, then do so and update the system -- system = system-set{f} | {lift(f,P)} --* -- else --* ( -- add one new variable "var_counter", and define the -- appropriate Laurent polynomial ring: P = (coefficientRing P) monoid(flatten entries vars P | {var_counter}, Inverses=>true, MonomialOrder=>RevLex); -- add the new Laurent polynomial to replace the rational equation: newvar := P_(numgens P - 1); system = system - set{f} | {sub(numerator(f),P)*newvar^(-1)}; system = system | {newvar - sub(denominator(f),P)}; counter = counter+1; ) ) ); -- For consistency, we make sure that everyone lives in the Laurent -- polynomial ring. This is not necessary for PHCpack, -- but for any further M2 calculations for the system, it is. system=apply(system,f-> sub(f,P)); system ) --------------------- -- TOP WITNESS SET -- --------------------- topWitnessSet = method( Options => {Verbose => false}) topWitnessSet (List,ZZ) := o->(system,dimension) -> ( -- IN: system, a polynomial system; -- dimension, top dimension of the solution set. -- OUT: a witness set for the top dimensional component, -- a list of nonsolutions if o.Verbose then stdio << "... calling constructEmbedding ..." << endl; e := constructEmbedding(system,dimension); if o.Verbose then stdio << "... calling solveSystem ..." << endl; s := solveSystem(e); g := zeroFilter(s,#e-1,1.0e-10); ns := nonZeroFilter(s,#e-1,1.0e-10); if o.Verbose then stdio << "... constructing a witness set ... " << endl; w := witnessSet(ideal(take(e,{0,#e-dimension-1})), ideal(take(e,{#e-dimension,#e-1})),g); return (w,ns); ) ----------------- -- TRACK PATHS -- ----------------- trackPaths = method(TypicalValue => List, Options=>{gamma=>0, tDegree=>2, Verbose => false, numThreads=>0, seeProgress=>false, interactive => false, saveSettingsPath => "", loadSettingsPath => "", intermediateSolutions => false}) trackPaths (List,List,List) := List => o -> (T,S,Ssols) -> ( -- IN: T, target system to be solved; -- S, start system with solutions in Ssols; -- Ssols, solutions at the start of the paths. -- OPT INPUTS: gamma, constant for gamma trick -- tDegree, degree of continuation parameter -- OUT: Tsols, solutions at the end of the paths. if not(class coefficientRing ring first T===ComplexField) then error "coefficient ring of target system is not complex"; if not(class coefficientRing ring first S===ComplexField) then error "coefficient ring of start system is not complex"; if (o.loadSettingsPath != "") and o.interactive then error "You cannot both load settings and be in interactive mode. Please reset your options."; R := ring first T; n := #T; targetfile := temporaryFileName() | "PHCtarget"; systemToFile(T,targetfile); outfile := temporaryFileName() | "PHCoutput"; if not (o.numThreads > 1) then ( startfile := temporaryFileName() | "PHCstart"; systemToFile(S,startfile); Ssolsfile := temporaryFileName() | "PHCstartsols"; solutionsToFile(Ssols,R,Ssolsfile); ) else ( startandsolutionfile := temporaryFileName() | "PHCstartandsols"; systemToFile(S,startandsolutionfile); solutionsToFile(Ssols, R, startandsolutionfile, Append=>true); ); Tsolsfile := temporaryFileName() | "PHCtargetsols"; batchfile := temporaryFileName() | "PHCbat"; if o.Verbose then stdio << "using temporary files " << outfile << " and " << Tsolsfile << endl; if n < numgens R then error "the system is underdetermined"; if n > numgens R then error "the system is overdetermined"; bat := openOut batchfile; if (o.loadSettingsPath != "") then ( if o.numThreads > 1 then ( bat << targetfile << endl << outfile << endl << startandsolutionfile << endl; ) else ( bat << targetfile << endl << outfile << endl <<"n"<< endl << startfile << endl << Ssolsfile << endl; ); optionFileLines := lines get o.loadSettingsPath; for i from 0 to #optionFileLines - 1 do ( if o.intermediateSolutions == true and i == #optionFileLines-1 then ( bat << "2" << endl; ) else ( bat << optionFileLines#i << endl; ) ); close bat; << batchfile << endl; checkIsRunnable(); run(PHCexe|" -p "|(if o.numThreads > 1 then ("-t"|o.numThreads) else "")|"<"|batchfile|" >phc_session.log"); checkIsRunnable(); run(PHCexe|" -z "|outfile|" "|Tsolsfile); ) -- making batch file else if o.interactive then ( bat << targetfile << endl << outfile << endl <<"n"<< endl << startfile << endl << Ssolsfile << endl; close bat; << "running (cat "|batchfile|"; cat) | PHCexe -p "<< endl; -- (cat batch; cat) | phc -p run("(cat "|batchfile|"; cat) | "|PHCexe|" -p "); checkIsRunnable(); run(PHCexe|" -z "|outfile|" "|Tsolsfile); ) else ( if not (o.numThreads > 1) then ( bat << targetfile << endl << outfile << endl <<"n"<< endl << startfile << endl << Ssolsfile << endl; -- first menu with settings of the construction of the homotopy bat << "k" << endl << o.tDegree << endl; if o.gamma != 0 then ( bat << "a" << endl << realPart o.gamma << endl; bat << imaginaryPart o.gamma << endl; ); bat << "0" << endl; -- second menu bat << "0" << endl; -- exit for now -- third menu if o.intermediateSolutions then ( bat << "2" << endl; ) else ( bat << "0" << endl; -- exit for now ); -- fourth menu bat << "0" << endl; -- exit for now close bat; ); if o.numThreads > 1 then ( bat << targetfile << endl << outfile << endl << startandsolutionfile << endl; -- first menu with settings of the construction of the homotopy bat << "k" << endl << o.tDegree << endl; if o.gamma != 0 then ( bat << "a" << endl << realPart o.gamma << endl; bat << imaginaryPart o.gamma << endl; ); bat << "0" << endl; -- second menu bat << "0" << endl; -- exit for now -- third menu if o.seeProgress then (bat << "y" << endl) else (bat << "n" << endl); close bat; ); checkIsRunnable(); run(PHCexe|" -p "|(if o.numThreads > 1 then ("-t"|o.numThreads) else "")|"<"|batchfile|" >phc_session.log"); checkIsRunnable(); run(PHCexe|" -z "|outfile|" "|Tsolsfile); ); -- parse and output the solutions if o.intermediateSolutions then ( return parseIntermediateSolutions(outfile,R); ); result := parseSolutions(Tsolsfile, R); if n > numgens R then ( result = apply(result, s->( if any(drop(first s, numgens R), x->abs x > 0.01) then error "slack value is nonzero"; {take(first s, numgens R)}|drop(s,1) )); totalN := #result; scan(result, s->( if s#1#"mult">1 then error "multiple root encountered"; if s#1#"mult"<0 then error "negative multiplicity"; )); result = select(result, s-> max(s#0/abs)<10000 -- path failed and/or diverged ); if PHCDBG>0 and #result < totalN then -- error "discarded!" << "track[PHCpack]: discarded "<< totalN-#result << " out of " << totalN << " solutions" << endl; ); if o.saveSettingsPath != "" then saveTrackPathsOptions(outfile, o.saveSettingsPath); return result; ) stripSpaces = method(); stripSpaces (String) := S -> ( -- Helper function to make saveTrackPaths more readable return replace(" ", "", S); ); addToFile = method(); addToFile (File, ZZ, String, ZZ) := (f, OptionNumber, S,StartIndex) -> ( -- Helper function to make saveTrackPaths more readable f << OptionNumber << endl << stripSpaces(substring(S,StartIndex,11)) << endl; ); saveTrackPathsOptions = method() saveTrackPathsOptions (String, String) := (outputFileName, saveFileName) -> ( saveFile := openOut saveFileName; fileLines := lines get outputFileName; for i from 0 to #fileLines - 1 do ( -- Start by saving the homotopy parameters. if fileLines#i == "HOMOTOPY PARAMETERS :" then( saveFile << "d" << endl << stripSpaces(substring(fileLines#(i+1),5,4)) << endl; saveFile << "k" << endl << stripSpaces(substring(fileLines#(i+2),5,4)) << endl; saveFile << "a" << endl << substring(fileLines#(i+3),6,21) << endl << substring(fileLines#(i+3),29,21) << endl; saveFile << "t" << endl << substring(fileLines#(i+4),6,21) << endl << substring(fileLines#(i+4),29,21) << endl; if #select(".no projective.",fileLines#(i+5)) == 1 then ( saveFile << "p" << endl << "n" << endl; ) else ( saveFile << "p" << endl << "y" << endl; ); saveFile << 0 << endl; ); -- Save chunk of 34 parameters. if fileLines#i == "GLOBAL MONITOR : " then ( addToFile(saveFile,1,fileLines#(i+1),46); addToFile(saveFile,2,fileLines#(i+2),46); addToFile(saveFile,3,fileLines#(i+3),46); addToFile(saveFile,4,fileLines#(i+4),46); addToFile(saveFile,5,fileLines#(i+5),46); addToFile(saveFile,6,fileLines#(i+6),46); addToFile(saveFile,7,fileLines#(i+8),46); addToFile(saveFile,8,fileLines#(i+8),58); addToFile(saveFile,9,fileLines#(i+9),46); addToFile(saveFile,10,fileLines#(i+9),58); addToFile(saveFile,11,fileLines#(i+10),46); addToFile(saveFile,12,fileLines#(i+10),58); addToFile(saveFile,13,fileLines#(i+11),46); addToFile(saveFile,14,fileLines#(i+11),58); addToFile(saveFile,15,fileLines#(i+12),46); addToFile(saveFile,16,fileLines#(i+12),58); addToFile(saveFile,17,fileLines#(i+13),46); addToFile(saveFile,18,fileLines#(i+13),58); addToFile(saveFile,19,fileLines#(i+15),46); addToFile(saveFile,20,fileLines#(i+15),58); addToFile(saveFile,21,fileLines#(i+16),46); addToFile(saveFile,22,fileLines#(i+16),58); addToFile(saveFile,23,fileLines#(i+17),46); addToFile(saveFile,24,fileLines#(i+17),58); addToFile(saveFile,25,fileLines#(i+18),46); addToFile(saveFile,26,fileLines#(i+18),58); addToFile(saveFile,27,fileLines#(i+19),46); addToFile(saveFile,28,fileLines#(i+19),58); addToFile(saveFile,29,fileLines#(i+21),46); addToFile(saveFile,30,fileLines#(i+21),58); addToFile(saveFile,31,fileLines#(i+22),46); addToFile(saveFile,32,fileLines#(i+22),58); addToFile(saveFile,33,fileLines#(i+23),46); addToFile(saveFile,34,fileLines#(i+23),58); saveFile << 0 << endl; ); if fileLines#i == "OUTPUT INFORMATION DURING CONTINUATION :" then ( saveFile << stripSpaces(substring(fileLines#(i+1),0,4)) << endl; ); ); close saveFile; ) ----------------- -- ZERO FILTER -- ----------------- zeroFilter = method(TypicalValue => List) zeroFilter (List,ZZ,RR) := (sols,k,tol) -> ( -- IN: sols, solutions of a polynomial system; -- k, index to a coordinate of a solution -- tol, tolerance for the absolute value of k-th coordinate -- OUT: list of solutions in sols where the k-th coordinate -- is less than the tolerance. return select(sols,t->isCoordinateZero(t,k,tol)); ) -------------------- -- intersectSlice -- -------------------- intersectSlice = method(TypicalValue => List) intersectSlice (WitnessSet, List) := (w, slcRR) -> ( -- IN: w, a witness set; -- slcRR, a list of linear equations. -- OUT: solutions of the equations of the witness set w -- which satisfy the list of linear equations. startSys:=join(equations(w),slice(w)); targetSys := equations(w) | slcRR; trackPaths(targetSys,startSys,w.Points) ) ----------------------------------- -- line search with golden ratio -- ----------------------------------- lineSearch = (F,a,b,tol,npoints) -> ( delta := (b-a)/(2*npoints); xmin := discretization1D (F, a, b, npoints ); a = max( xmin - delta, a); b = min( xmin + delta, b); xmin = goldenSearch (F, a, b, tol); return xmin; ) goldenSearch = (F,a,b,tol) -> ( -- Applies the golden section search method -- to minimize a unimodal function F over [a,b]. -- IN: F, a function in one variable; -- a, the left bound of the search interval; -- b, the right bound of the search interval; -- tol, tolerance on the approximate minimum. -- OUT: if F is unimodal (it has a unique minimum), -- then the value x on return approximates -- the minimum with respect to the tolerance tol. -- EXAMPLE: myFunction = (x) -> (1-x)*sin((3*x)^3) -- goldenSearch( myFunction, .4, .7, 1.0e-4) gr := (sqrt(5) - 1)/2; c := b - gr * (b - a); d := a + gr * (b - a); while abs(c - d) > tol do ( Fc := F(c); Fd := F(d); if Fc < Fd then ( b = d;) else ( a = c;); c = b - gr * (b - a); d = a + gr * (b - a); ); return (b + a) / 2; ) ----------------------------------------------------------- -- convert coefficients of matrix to list of polynomials -- ----------------------------------------------------------- matrix2slice = (slcmat, w) -> ( -- Given a witness set w and a matrix of coefficients -- for a set of hyperplanes, uses the ring of w.Equations -- to return the list representation of the hyperplanes -- with coefficients in slcmat. -- IN: slcmat, a matrix of coefficients of hyperplanes; -- w, a witness set. -- OUT: a list of linear equations with coefficient -- from sclmat and variables from w.Equations. R2 := ring w.Equations; X := transpose (vars(R2) | 1); slc := flatten entries (promote(slcmat,R2) * X); return slc; ) ------------------------- -- the cost of a slice -- ------------------------- sliceCost = (slcmat, w) -> ( -- Returns the norm of the imaginary parts of the solutions -- that satisfy w.Equations on the slices with coefficient -- in the matrix slcmat. -- IN: sclmat, a matrix of coefficients of hyperplanes; -- w, a witness set. -- OUT: sum of the squares of the imaginary parts of the coordinates -- of the witness points on the hyperplanes defined by sclmat. slc := matrix2slice(slcmat,w); fsols := intersectSlice(w,slc); cost := sum ( coordinates(fsols_0) / imaginaryPart / (x->x^2) ); return cost; ) realPartMatrix = (m) -> matrix applyTable (entries m, x->1_CC*realPart x) rotationOfSlice = (t,startSlice) -> ( -- Returns a rotation of the coefficients of startSlice -- about the angle t. -- IN: t, an angle; -- startSlice, a matrix with coefficients of the slice. -- OUT: new slice, rotated from startSlice with angle t. c := numColumns(startSlice); M1 := id_(CC^c); M2 := mutableMatrix M1; M2_(0,0) = cos(t); M2_(0,1) = -sin(t); M2_(1,0) = sin(t); M2_(1,1) = cos(t); M3 := matrix M2; return startSlice*M3; ) rotationMatrix2D = (c,i,j,t)-> ( -- Returns a general rotation matrix in dimension c, -- involving variables i and j and angle t. -- IN: c, dimension of the matrix; -- i, first variable involved in the rotation matrix; -- j, second variable involved in the rotation matrix; -- t, angle in the rotation matrix. -- OUT: a rotation matrix of dimension c about angle t, -- which involves variables i and j. M1 := id_(CC^c); M2 := mutableMatrix M1; M2_(i,i) = cos(t); M2_(i,j) = -sin(t); M2_(j,i) = sin(t); M2_(j,j) = cos(t); M3 := matrix M2; return M3 ) changeOfSlice2D = (t1,t2,startSlice)-> ( -- Applies a rotation matrix using angles t1 and t2 on startSlice. -- IN: t1, first angle for the first hyperplane in startSlice; -- t2, second angle for the second hyperplane in startSlice; -- startSlice, coefficients of two hyperplanes. -- OUT: coefficients of a rotated startSlice. c := numColumns(startSlice); rotMatrix1 := rotationMatrix2D(c,0,1,t1); rotMatrix2 := rotationMatrix2D(c,1,2,t2); row1 := startSlice^{0}*rotMatrix1; row2 := startSlice^{1}*rotMatrix2; return row1||row2; ) -- discretization1D = (F,a,b,n) -> ( -- Evaluates the function F over n+1 equidistant points -- in the interval [a,b], including a and b, and returns -- the point where F takes the minimal value. -- IN: F, a function in one variable; -- a, left bound of the interval [a,b]; -- b, right bound of the interval [a,b]; -- n, number of points in the interval [a,b]. -- OUT: the point in the n equidistant points in [a,b] -- where F takes its minimal value. range := for i to n list a+(b-a)*i/n; functionValues := for x in range list F(x); minValue := min(functionValues); minPos := position(functionValues,a->(a==minValue)); return range_minPos; ) discretization2D = (F,a1,b1,a2,b2,n) -> ( -- Returns the slice with the largest number of real roots. range1 := for i to n-1 list a1+(b1-a1)*i/n; range2 := for i to n-1 list a2+(b2-a2)*i/n; functionValues := flatten for x in range1 list for y in range2 list F(x,y); minValue := min(functionValues); posInList := position(functionValues,a->(a==minValue)); minPos := (posInList//n,posInList%n); return (range1_(minPos#0),range2_(minPos#1)); ) alternatingMinimization = (F,a1,b1,a2,b2,tol) -> ( -- Applies the method of alternating minimization. -- EXAMPLE: myF = (x,y) -> x^2 + y^2 -- minF = alternatingMinimization(myF,-1,1,-1,1,1.0e-4) cOld := a1; dOld := a2; c:=a1+random(RR)*(b1-a1); d:=a2+random(RR)*(b2-a2); while abs(c - cOld) > tol and abs(d - dOld) > tol do ( Fc:=(y)->F(c,y); dOld=d; d=goldenSearch(Fc,a2,b2,tol); Fd:=(x)->F(x,d); cOld=c; c=goldenSearch(Fd,a1,b1,tol); ); return (c,d); ) ----------------- -- realSlice1D -- ----------------- realSlice1D = method(TypicalValue => List, Options => {searchNpoints => 20, searchDelta => 0.1, searchTolerance => 1.0e-4}) realSlice1D(WitnessSet) := o -> (w) -> ( -- Starting from the given one dimensional witness set, -- applies line search to find a real slice. -- IN: w, a witness set; -- searchNpoints (optional), number of points in the discretization; -- searchDelta (optional), 2*searchDelta is width of search interval; -- searchTolerance (optional), tolerance for the line search method. -- OUT: a slice where the number of real solutions was maximal. startSlice := realPartMatrix(w.Slice); costfun := (a) -> sliceCost(rotationOfSlice(a,startSlice), w); amin := lineSearch (costfun, 0, 2*pi, o.searchTolerance, o.searchNpoints); slcmin := rotationOfSlice(amin,startSlice); return matrix2slice(slcmin,w) ) ----------------- -- realSlice2D -- ----------------- realSlice2D = method(TypicalValue => List, Options => {searchNpoints => 5, searchDelta => 0.1, searchTolerance => 1.0e-4}) realSlice2D(WitnessSet) := o -> (w) -> ( startSlice := realPartMatrix(w.Slice); costfun := (a,b) -> sliceCost(changeOfSlice2D(a,b,startSlice),w); (min1,min2) := discretization2D(costfun,0,2*pi,0,2*pi,o.searchNpoints); a1 := min1 - o.searchDelta; b1 := min1 + o.searchDelta; a2 := min2 - o.searchDelta; b2 := min2 + o.searchDelta; tol := o.searchTolerance; (min1,min2) = alternatingMinimization(costfun,a1,b1,a2,b2,tol); slcmin := changeOfSlice2D(min1,min2,startSlice); return matrix2slice(slcmin,w) ) ------------------- -- versionNumber -- ------------------- versionNumber = method(TypicalValue => Nothing, Options => {Verbose => false}) versionNumber(Nothing) := o -> (Nothing) -> ( -- Calling versionNumber(null) returns a tuple of two strings, -- with the version number and release date. -- IN: if the option Verbose is true, versionNumber(Verbose=>true), -- then the output of phc --version is printed to screen. -- OUT: information about the current version of phc. filename := temporaryFileName() | "PHCversion"; run("phc --version > "|filename); data := get filename; if o.Verbose then stdio << data << endl; if #data < 31 then ( stdio << "Which version of phc is in your execution path?" << endl; stdio << data << endl; return ("", ""); ) else ( vnbr := substring(4,6,data); date := substring(#data-11,10,data); return (vnbr, date); ); ) --##########################################################################-- -- DOCUMENTATION --##########################################################################-- beginDocumentation() load "./PHCpack/PHCpackDoc.m2"; --##########################################################################-- -- TESTS --########################################################################## ----------------------------------- -- test 0: cascade ----------------------------------- TEST/// R=CC[x11,x22,x21,x12,x23,x13,x14,x24] L={x11*x22-x21*x12,x12*x23-x22*x13,x13*x24-x23*x14} -- assert( # cascade L == 1 )--there is one component of dim.5. C=cascade(L, Verbose=>true) assert(#C == 1) ///; ----------------------------------- --test 1: constructEmbedding ----------------------------------- TEST/// R=CC[x,y,z] L={x^2+y^2, x+y+z} Emb=constructEmbedding(L,1) assert( # Emb == 4) -- there are four equations assert ( member(zz1, Emb) == true ) --one of the eqns is zz1 ///; ----------------------------------- -- test 2: factorWitnessSet ----------------------------------- TEST/// R = CC[x,y]; system = {x*y}; (w, ns)=topWitnessSet(system,1); V=factorWitnessSet(w); assert ((# V#1) == 2) ///; ----------------------------------- -- test 3: isCoordinateZero ----------------------------------- TEST/// P=point({{0,1.0e-12}}) assert (isCoordinateZero(P,1,1.0e-10) == true) ///; ----------------------------------- -- test 4: isWitnessSetMember ----------------------------------- TEST/// R = CC [x,y] system = {2*y+2*x, 4*y + 4*x} (W, ns)=topWitnessSet(system, 1) assert isWitnessSetMember(W,point{{0,0}}) ///; ----------------------------------- -- test 5: mixedVolume ----------------------------------- TEST/// R=CC[x,y,z] S={y-x^2,z-x^3,x+y+z-1} m=mixedVolume(S) assert(m==3) R=CC[x,y,z] S1={y^3+z^2+3,z^2+x^4+x^4*z^2+4,x^4+y^3+x^4*y^3+5} M=mixedVolume(S1) assert(M==48)--testing output against by-hand calculation R=CC[x,y] S2={x^2+x*y+y^2+x+y+1,x^5+x^4*y+x^3*y^2+x^2*y^3+x*y^4+y^5} --another example M=mixedVolume(S2) assert(M==10)--testing output against by-hand calculation ///; ----------------------------------- -- test 6: nonZeroFilter ----------------------------------- TEST/// R = CC[x,y]; f = { x^3*y^5 + y^2 + x^2*y, x*y + x^2 - 1}; fSols = solveSystem(f); nonzeroSols = nonZeroFilter(fSols,0,1.0e-10); assert(# nonzeroSols ==10) -- this is just complement of zeroFilter which is tested in more detail. Here we are just counting the number of solutions to make sure the method ran. ///; ------------------------------------- -- test 7: numericalIrreducibleDecomposition ------------------------------------- TEST/// R = CC [x,y] system={x^2*(y^2-1)} V=numericalIrreducibleDecomposition(system) assert all(components V, c->dim c==1) --there are three components of dimension 1 assert (#components V == 3) ///; ----------------------------------- -- test 8: refineSolutions ----------------------------------- TEST/// R = CC[x,y]; S = {x^2 - 1/3, x*y - 1}; oldRoots = solveSystem(S); r0 = oldRoots#0#Coordinates#1 newRoots = refineSolutions(S,oldRoots,64) --recall that solutions are of type Point. --check if precision increased: assert(precision newRoots#0#Coordinates#1 > precision oldRoots#0#Coordinates#1) --check if input number of decimal places, 64, used correctly: assert(precision newRoots#0#Coordinates#1 == ceiling(log_2(10^64))) ///; ----------------------------------- -- test 9: solveRationalSystem ----------------------------------- TEST/// QQ[x,y,z]; sys = {y-x^2, z-x^3, (x+y+z-1)/x}; sols = solveRationalSystem(sys); assert(# sols == 3); --there are 3 solutions real = realPoints(sols); assert(# real ==1); --one solution is real ///; ----------------------------------- -- test 10: solveSystem ----------------------------------- TEST/// R=CC[x,y,z] S={x^2-y*z-3,y^2-x*z-4,z^2-x*y-5} L=solveSystem(S) n=# L assert(n==2) sol1={11/6.,-1/6.,-13/6.} sol2={-11/6.,1/6.,13/6.} assert((abs((sol1-L_0#Coordinates)_0)<.00000000001 and abs((sol1-L_0#Coordinates)_1)<.00000000001 and abs((sol1-L_0#Coordinates)_2)<.00000000001) or (abs((sol1-L_1#Coordinates)_0)<.00000000001 and abs((sol1-L_1#Coordinates)_1)<.00000000001 and abs((sol1-L_1#Coordinates)_2)<.00000000001)) assert((abs((sol2-L_0#Coordinates)_0)<.00000000001 and abs((sol2-L_0#Coordinates)_1)<.00000000001 and abs((sol2-L_0#Coordinates)_2)<.00000000001) or (abs((sol2-L_1#Coordinates)_0)<.00000000001 and abs((sol2-L_1#Coordinates)_1)<.00000000001 and abs((sol2-L_1#Coordinates)_2)<.00000000001)) ///; ----------------------------------- -- test 11: toLaurentPolynomial ----------------------------------- TEST/// QQ[x,y,z]; sys = {y-x^2, z-x^3, (x+y+z-1)/x}; convertedSys = toLaurentPolynomial(sys,w); R=ring ideal convertedSys assert(isPolynomialRing(R)) --make sure it is a poly ring O=options(R); assert(O#Inverses) -- want to make sure this is a Laurent ring ///; ----------------------------------- -- test 12: topWitnessSet ----------------------------------- TEST/// R = CC[x,y]; system = {x*y}; (w, ns)=topWitnessSet(system,1); assert (dim w == 1) assert (degree w == 2) ///; ----------------------------------- -- test 13: trackPaths ----------------------------------- TEST/// R = CC[x,y]; f = { x^3*y^5 + y^2 + x^2*y, x*y + x^2 - 1}; (m,q,qsols) = mixedVolume(f,StartSystem=>true); fsols = trackPaths(f,q,qsols) assert(# fsols == 8) ///; ----------------------------------- -- test 14: zeroFilter ----------------------------------- TEST/// R = CC[x,y]; f = { x^3*y^5 + y^2 + x^2*y, x*y + x^2 - 1}; fSols = solveSystem(f); zeroSols = zeroFilter(fSols,1,1.0e-10); assert( sort {zeroSols_0#Coordinates,zeroSols_1#Coordinates} == {{-1, 0}, {1, 0}} ) ///; --##########################################################################-- end -- terminate reading ... --##########################################################################-- restart installPackage("PHCpack",RemakeAllDocumentation=>true) installPackage ("PHCpack",RerunExamples=>true) check "PHCpack" viewHelp PHCpack --many errors seem to be due to the fact some things are not exported: --error: mutable unexported unset symbol(s) in package PHCpack: 'generalEquations', 'NumericalVariety', 'numericalVariety', 'IsIrreducible' /Users/sxp61/m2svn/trunk/M2/Macaulay2/packages/PHCpack.m2:329:12-329:28: here is the first use of 'generalEquations' /Users/sxp61/m2svn/trunk/M2/Macaulay2/packages/PHCpack.m2:401:34-401:50: here is the first use of 'NumericalVariety' /Users/sxp61/m2svn/trunk/M2/Macaulay2/packages/PHCpack.m2:488:3-488:19: here is the first use of 'numericalVariety' /Users/sxp61/m2svn/trunk/M2/Macaulay2/packages/PHCpack.m2:588:26-588:39: here is the first use of 'IsIrreducible' restart needsPackage "PHCpack"; var0 = {x_11,x_12,x_16,x_22,x_23,x_33,x_34,x_44,x_45,x_55,x_56,x_66,y_13,y_14,y_15, y_24,y_25,y_26,y_35,y_36,y_46}; QQ[var0]; rationalSystem = { (22/3)*x_11+(8/7)*x_12+2*x_16-1, x_23*y_13+(22/3)*x_12+(8/7)*x_22, x_33*y_13+x_34*y_14+(8/7)*x_23, x_34*y_13+x_44*y_14+x_45*y_15, x_45*y_14+x_55*y_15+2*x_56, x_56*y_15+(22/3)*x_16+2*x_66, (8/7)*x_12+(14/11)*x_22+(12/5)*x_23-1, x_34*y_24+(14/11)*x_23+(12/5)*x_33, x_44*y_24+x_45*y_25+(12/5)*x_34, x_45*y_24+x_55*y_25+x_56*y_26, (12/5)*x_23+(28/51)*x_33+(102/144)*x_34-1, x_56*y_25+x_66*y_26+(8/7)*x_16, x_45*y_35+(28/51)*x_34+(102/144)*x_44, x_55*y_35+x_56*y_36+(102/144)*x_45, (102/144)*x_34+(205/162)*x_44+(3/2)*x_45-1, x_16*y_13+x_56*y_35+x_66*y_36, x_56*y_46+(205/162)*x_45+(3/2)*x_55, x_16*y_14+x_66*y_46+(3/2)*x_56, (3/2)*x_45+(517/784)*x_55+(8/3)*x_56-1, x_16*y_15+(517/784)*x_56+(8/3)*x_66, 2*x_16+(8/3)*x_56+(29/196)*x_66-1}; system = (sub(ideal rationalSystem, CC[var0]))_* solutions = solveSystem system; # solutions solutions_0 newSystem = { (22531/300)*x_11+(821/70)*x_12+(4507/210)*x_16-1, x_23*y_13+(22531/300)*x_12+(821/70)*x_22, x_33*y_13+x_34*y_14+(821/70)*x_23, x_34*y_13+x_44*y_14+x_45*y_15, x_45*y_14+x_55*y_15+(4507/210)*x_56, x_56*y_15+(22531/300)*x_16+(4507/210)*x_66, (821/70)*x_12+(140953/11025)*x_22+(12325/504)*x_23-1, x_34*y_24+(140953/11025)*x_23+(12325/504)*x_33, x_44*y_24+x_45*y_25+(12325/504)*x_34, x_45*y_24+x_55*y_25+x_56*y_26, x_56*y_25+x_66*y_26+(821/70)*x_16, (12325/504)*x_23+(282013/5184)*x_33+(10231/1440)*x_34-1, x_45*y_35+(282013/5184)*x_34+(10231/1440)*x_44, x_55*y_35+x_56*y_36+(10231/1440)*x_45, x_16*y_13+x_56*y_35+x_66*y_36, (10231/1440)*x_34+(205697/16200)*x_44+(30529/2520)*x_45-1, x_56*y_46+(205697/16200)*x_45+(30529/2520)*x_55, x_16*y_14+x_66*y_46+(30529/2520)*x_56, (30529/2520)*x_45+(5175321/78400)*x_55+(897/35)*x_56-1, x_16*y_15+(5175321/78400)*x_56+(897/35)*x_66, (4507/210)*x_16+(897/35)*x_56+(293581/19600)*x_66-1}; newSolutions = solveSystem newSystem; # newSolutions smallSolution = zeroFilter(newSolutions, 11, 1.0e-18) smallerSolution = refineSolutions(newSystem, smallSolution, 64) mixedVolume system time degree ideal rationalSystem var1 = {x11, x22, x21, x12, x23, x13, x14, x24}; R = QQ[var1]; rationalSystem = { x11*x22-x21*x12, x12*x23-x22*x13, x13*x24-x23*x14}; system = (sub(ideal rationalSystem, CC[var1]))_* V = numericalIrreducibleDecomposition system WitSets = V#5; w = first WitSets w.cache.IsIrreducible R = ring rationalSystem_0 PD = primaryDecomposition ideal rationalSystem for I in PD list << "(dim=" << dim I << ", deg=" << degree I << ") " restart loadPackage "PHCpack" r = CC[x,y] solveSystem({2*x+y+5,5*y^2+3*x}) restart loadPackage "PHCpack" R = CC[x,y] f = {x^3*y^5 + y^2 + x^2*y, x*y + x^2 - 1}; I = ideal f; --(mv,smv,q,qsols) = mixedVolume(f,StableMixedVolume=>true,StartSystem=>true) (mv,q,qsols) = mixedVolume(f,StartSystem=>true) fsols = trackPaths(f,q,qsols, interactive=>true, saveSettingsPath=>"settings.phc") fsols = trackPaths(f,q,qsols, loadSettingsPath=>"settings.phc") --m = mixedVolume(f) --(mv,sv) = mixedVolume(f,StableMixedVolume => true) --mv = mixedVolume(f,interactive=>true) --(mv,smv,q,qsols) = mixedVolume(f,interactive=>true) --mixedVolume(f,interactive=>true) --fsols = trackPaths(f,q,qsols)