restart --Water. {*Macaulay2: software system devoted to supporting research in algebraic geometry and commutative algebra*} --The Authors: Dan Grayson and Mike Stillman --Goal: Describe algebraic varieties. R=QQ[x,y,z]--This is a ring I=ideal(x^2+y^2-1,x-z)--This is an ideal generated by two polynomials. degree I dim I codim I --Let's eliminate! --Twisted cubic 1 R=QQ[x0,x1,x2,x3]**QQ[s,t] I1=ideal({x0,x1,x2,x3}-{s^3,s^2*t,s*t^2,t^3}) I2=eliminate({s,t},I1) --saturate --Twisted cubic 2 R=QQ[x0,x1,x2,x3] M=matrix{{x0,x1,x2},{x1,x2,x3}} I=minors(2,M) --hilbertPolynomial I --resI=resolution I class I theGenerators=I_* J=ideal drop(theGenerators,1) numgens I numgens J decJ=decompose J netList decJ decJ/degree --degree of each element in the list. help apply I I2 ----- {* Cut for time aHyperplane=ideal random({1},R) theIntersection=J+aHyperplane E=eliminate({x2,x3},theIntersection) factor E_0 gens ring E *} {* Cut for time. --Let's do a function to do this all in one step. eliminateB=(decomposeList,aHyperplane)->( varsEliminated:=drop(gens ring first decomposeList,2); print varsEliminated; apply(decomposeList,i->( eliminate(varsEliminated,i+aHyperplane)))) eliminateB(decJ,aHyperplane) *} {* --We can also work over finite fields and multigraded rings. kk=ZZ/30103 isPrime 30103 R=kk[x0,x1]**kk[y0,y1] f= x0^2*y1-x1^2*y0 degree f multidegree ideal f *} viewHelp "Macaulay2" viewHelp "Package" --Mention Taylor B. and Chris E. --Click on the list of packages link. --file:///Applications/Macaulay2-1.10/share/doc/Macaulay2/Macaulay2Doc/html/_packages_spprovided_spwith_sp__Macaulay2.html loadPackage"Resultants"; F = genericPolynomials(QQ,{4,4}) first F time resF = Resultant F degree resF ---This didn't give an answer 'now' --F = genericPolynomials(QQ,{4,4,4}) --time resF = Resultant F --degree resF --------------------------------------------------- --Numerical algebraic geometry --------------------------------------------------- kk=CC R=kk[x1,x2] --Anton Leykin loadPackage("NumericalAlgebraicGeometry",Reload=>true) help solveSystem theSols=solveSystem({x1^2+x2^2-1,x1-x2}) (sqrt 2)/2 onePoint= first theSols class onePoint keys onePoint onePoint#SolutionStatus netList pairs onePoint --How well does this implementation work? testTiming=(numVars,theDegree)->( R=kk[apply(numVars,i->"x"|i+1)]; -- F:=apply(numVars,i->random(theDegree,R)+random CC); F=apply(numVars,i->sum apply(theDegree+1,d->random(d,R))); startTime:=currentTime(); win=solveSystem(F); endTime:=currentTime(); totalTime:=(endTime-startTime); timePerPath:=totalTime/#win*1.0; ("Total Time: ",totalTime, "Time per path", timePerPath, "Bezout Bound",theDegree^numVars, "Found solutions",#win)) testTiming(2,5) testResiduals=(k,aSol)->( thePoly:=F_k; theSubs:=apply(gens ring F_k,aSol#Coordinates,(i,j)->i=>j); sub(F_k,theSubs) ) testResiduals(1,first win) theResids=apply(win,i->abs testResiduals(0,i)); max theResids min theResids --Resultant degrees. restart loadPackage"Bertini"--I am directly involved. loadPackage"Resultants" F = genericPolynomials(QQ,{4,4}) R= ring first F theCoefficients=gens coefficientRing R S=CC[s,t] coefSubs=apply(theCoefficients,i->i=>random(1,S)) makeB'InputFile(storeBM2Files, HomVariableGroup=>{gens S,gens R}, B'Polynomials=>F, B'Functions=>coefSubs ) readFile(storeBM2Files,"input",10000) runBertini(storeBM2Files) readFile(storeBM2Files) readFile(storeBM2Files,"nonsingular_solutions",1000) theWitnessPointsCoordinates=importSolutionsFile(storeBM2Files); first theWitnessPointsCoordinates readFile(storeBM2Files,"main_data",1000) witnessPoints=importMainDataFile(storeBM2Files); aWitnessPoint=first witnessPoints keys aWitnessPoint netList pairs aWitnessPoint --Bertini.m2 kk=CC loadPackage("Bertini",Reload=>true) testBertiniTiming=(numVars,theDegree)->( R=kk[apply(numVars,i->"x"|i+1)]; -- F:=apply(numVars,i->random(theDegree,R)+random CC); F=apply(numVars,i->sum apply(theDegree+1,d->random(d,R))); makeB'InputFile(storeBM2Files, AffVariableGroup=>{gens R}, B'Configs=>{UseRegeneration=>0},--Change to 1 to use regenereation. B'Polynomials=>F ); startTime:=currentTime(); runBertini(storeBM2Files); endTime:=currentTime(); win=importMainDataFile(storeBM2Files); totalTime:=(endTime-startTime); timePerPath:=totalTime/#win*1.0; ("Total Time: ",totalTime, "Time per path", timePerPath, "Bezout Bound",theDegree^numVars, "Found solutions",#win)) print testBertiniTiming(4,2) 2^10 startTime:=currentTime(); win=solveSystem(F); endTime:=currentTime(); totalTime:=(endTime-startTime); timePerPath:=totalTime/#win*1.0; ("Total Time: ",totalTime, "Time per path", timePerPath, "Bezout Bound",theDegree^numVars, "Found solutions",#win)) testResiduals=(k,aSol)->( thePoly:=F_k; theSubs:=apply(gens ring F_k,aSol#Coordinates,(i,j)->i=>j); sub(F_k,theSubs) ) {* --Sparse polynomial example. sparseThePolynomials=(aPoly,monomialExponents)->( theVars:=gens ring aPoly; monomialSupport:=apply(monomialExponents, i->product apply(i,theVars,(e,x)->x^e )); print monomialSupport; (monPoly,coefPoly):=(coefficients aPoly)/entries/flatten//toSequence; sum apply(monPoly,coefPoly, (i,j)->if member(i,monomialSupport) then i*j else 0 )) monomialExponents={{0,6,0},{1,4,1},{2,2,2},{3,0,3}} F = genericPolynomials(QQ,{6,6,6}) sparseThePolynomials(first F,monomialExponents) F=apply(F,f->sparseThePolynomials(f,monomialExponents)) *}