################################################################################################## # SPECTRA : A Maple package for solving semidefinite programs in exact arithmetic # ################################################################################################## #------------------------------------------------------------------------------------------------# SPECTRA := module() description "Semidefinite Programming solved Exactly with Computational Tools of Real Algebra": export SolveLMI: local mygbasis, mygbasiselim, AffineCut, CharPol, CharPolPar, CheckDimension, CheckHankelness, CheckPSD, CheckRegularity, ComputeDimensionIntrinsec, ComputeMinors, ComputeMinors2, ComputeSolutionsDimZero, ComputeTheRank, elim_factors, IncidenceVarietyHankel, IncidenceVarietySym, LagSyst, luckyprime, myindex, myisolate_gb, myisolate_rur, MyRand, myverb, myzerotest, NumericalCheckPSD, NumericalRankOfMatrix, PreprocessInput, PrintVerbose, PutBracketIfAlone, ReduceToList, ReduceToSequence, RestoreInputVariables, SampleLowRankHankel, SampleLowRankSym, SequenceOfSigns, SplitParametrization, SplitParametrization_old, subs_interval, Vectorize, zerochar: option package: uses LinearAlgebra, linalg, combinat, FGb, VectorCalculus, Statistics, fgbrs: ## parameters myzerotest := 10^(-8): luckyprime := 65521: zerochar := 65521: ## cambiato 0 -> 65521 -> 9001 myverb := 0: myindex := 5000000: ## standard 5000000 #myindex := 15000000: ## for Elliptope(5) rank 3,4 #### TO DO LIST # 1. Risolvere problema cambio interno variabili # 2. Preprocessing della matrice: se ha 0 sulla diagonale, o (X1,-X1) sulla diagonale ... # 3. rank and psd #Encapsulation of functions for computing GB mygbasis:=proc(F, char, vars1, vars2, opts:={"verb"=0}) local mysize,mygb; try mygb:=fgb_gbasis(F, char, vars1, vars2,opts); catch: mysize := 1000000; # mysize:=5000000; try mygb:=fgb_gbasis(F, char, vars1, vars2,{op(opts),"index"=mysize}); catch: mygb:=fgb_gbasis(F, char, vars1, vars2,{op(opts),"index"=2*mysize}); end: end: return mygb; end: mygbasiselim:=proc(F, char, vars1, vars2, opts:={"verb"=0}) local mysize,mygb; try mygb:=fgb_gbasis_elim(F, char, vars1, vars2,opts); catch: mysize := 1000000; # mysize:=5000000; try mygb:=fgb_gbasis_elim(F, char, vars1, vars2,{op(opts),"index"=mysize}); catch: mygb:=fgb_gbasis_elim(F, char, vars1, vars2,{op(opts),"index"=2*mysize}); end: end: return mygb; end: ################################################################################################## # Feasibility problem of Semidefinite Programing :: Linear Matrix Inequalities # ################################################################################################## # Main routine : SolveLMI # # Main subroutines : SampleLowRankSym, SampleLowRankHankel, LagSyst, ComputeSolutionsDimZero # ################################################################################################## #------------------------------------------------------------------------------------------------# SolveLMI := proc() ############################################################################### local M, m, r, ll, i, j, degg, checkdim, minmin, GB, PAR, ISOL, clm1, clm2, smpl1, relat, smpl2, toadd, algdeg, mytotaloutput, avoid_further_dimension_computations, fusion, mylist, fusiondeg, parsol, mysplit, myboolean, listparsol, glob_out_rank_r, tt, tba, constr, constr0, ddenom: global list_of_ranks, set_of_options, regularity, glob_out, invars, workingvars, hankelness, char_pol, optimizing_vars: regularity := 0: M := args[1]: if whattype(M) <> Matrix then error "invalid input: the first argument must be a symmetric linear matrix" fi: m := RowDimension(M): ## Think about preprocessing! #M,relat := PreprocessInput(M): #if relat <> 0 then # SolveLMI(M): #fi: char_pol := CharPol(M): #print(char_pol); ## substitution of variables (NOT ENABLED BECAUSE THERE IS SOME PROBLEM) invars := [op(indets(M))]: #for i from 1 to nops(invars) do # M := subs({invars[i] = cat('WW',i)}, M): #od: #workingvars := [op(indets(M))]: workingvars := invars: glob_out := []: list_of_ranks := [seq(i,i=0..m-1)]: set_of_options := {}: mytotaloutput := []: optimizing_vars := invars: ## 19 Apr 2017 ## control on the dimensions if m <> ColumnDimension(M) then error "invalid input: the matrix is not square" fi: ## organizing input arguments if nargs > 4 then ### 19 Apr 2017 error "too many input arguments": elif nargs = 2 then if whattype(args[2]) <> set then error "invalid input: the second argument must be a set": fi: set_of_options := args[2]: elif nargs = 3 then if whattype(args[2]) <> set then error "invalid input: the second argument must be a set": fi: if whattype(args[3]) <> list then error "invalid input: the third argument must be a list": fi: set_of_options := args[2]: list_of_ranks := [op({op(args[3])})]: elif nargs = 4 then ## 19 Apr 2017 if whattype(args[2]) <> set then error "invalid input: the second argument must be a set": fi: if whattype(args[3]) <> list then error "invalid input: the third argument must be a list": fi: if whattype(args[4]) <> list then error "invalid input: the fourth argument must be a list": fi: set_of_options := args[2]: list_of_ranks := [op({op(args[3])})]: optimizing_vars := [op({op(args[4])})]: fi: ## default list_of_ranks if op(list_of_ranks) = NULL then list_of_ranks := [seq(i,i=0..m-1)]: fi: ## eliminate non-admissible ranks list_of_ranks := [op({op(list_of_ranks)} intersect {seq(i,i=0..m-1)})]: if op(list_of_ranks) = NULL then error "invalid input: expected rank must be strictly less than matrix size": fi: ## eliminate non-admissible variables for linear forms optimizing_vars := [op({op(optimizing_vars)} intersect {op(invars)})]: ## eliminate non-admissible options set_of_options := set_of_options intersect {rnk,deg,all,reg,par,dbg2810,ver,nodim}: ## a debug function to keep secret if dbg2810 in set_of_options then debug(SampleLowRankSym,SampleLowRankHankel,LagSyst): fi: ## control on the symmetry for i from 2 to m do for j from 1 to i-1 do if M[i,j] <> M[j,i] then error "invalid input: the matrix is not symmetric" fi: od: od: ## control on the non triviality if nops(indets(M)) = 0 then error "invalid input: constant matrix": fi: ### checking the Hankelness hankelness := CheckHankelness(M): ## rank zero case (check if A(x) = 0 has a solution) if 0 in list_of_ranks then ll := solve(map(i->i=0,[Vectorize(M)])): if ll <> NULL then if {rnk,deg} subset set_of_options then glob_out := [op(glob_out), [op(ll), rnk = 0, deg = 1]]: elif rnk in set_of_options then glob_out := [op(glob_out), [op(ll), rnk = 0]]: elif deg in set_of_options then glob_out := [op(glob_out), [op(ll), deg = 1]]: else glob_out := [op(glob_out), [op(ll)]]: fi: if {all} union set_of_options <> set_of_options then #glob_out := RestoreInputVariables(glob_out, invars, workingvars): return([[ReduceToSequence(glob_out)]]): fi: fi: fi: avoid_further_dimension_computations := 0: if nodim in set_of_options then avoid_further_dimension_computations := 1: fi: ## higher rank defects for r in [op({op(list_of_ranks)} minus {0})] do ## compute the dimension of rank defect loci if dbg2810 in set_of_options then printf("Computing minors BEGIN\n"): fi: minmin := [ComputeMinors2(r+1,M)]: if dbg2810 in set_of_options then printf("Computing minors END\n"): fi: if avoid_further_dimension_computations = 0 and nops(indets(M)) > 1 then if dbg2810 in set_of_options then printf("Computing dimension low rank locus BEGIN\n"): fi: checkdim := CheckDimension(minmin): if dbg2810 in set_of_options then printf("Computing dimension low rank locus END --> "): lprint(checkdim): fi: if checkdim = false then avoid_further_dimension_computations := 1: fi: fi: ## Zero-dimensional case (I do not build incidence varieties) if nops(indets(M)) = 1 or ( avoid_further_dimension_computations = 0 and checkdim = 0 and op(indets(minmin)) = op(indets(M)) ) then if dbg2810 in set_of_options then printf("Zerodimensional case\n"): fi: GB := mygbasis(minmin,zerochar,[],[op(indets(minmin))]): constr0 := CharPol(M): ddenom := denom(constr0): constr := [seq(abs(ddenom[i])*constr0[i], i=1..nops(ddenom))]: PAR := []: if GB <> [] and GB <> [1] then if dbg2810 in set_of_options then printf("Non-trivial Groebner basis\n"): fi: PAR := [fgbrs:-rs_rur_gb(GB, [op(indets(minmin))], constraints=constr)]: fi: if PAR <> [] then myboolean := irreduc(PAR[1]): ## in case the polynomial q is reducible (zero-dim) if myboolean = false then mysplit := [SplitParametrization(PAR)]: mysplit := ReduceToList(mysplit): for i from 1 to nops(mysplit) do if whattype(mysplit[i][1]) <> integer then constr := mysplit[i][4]: parsol := ComputeSolutionsDimZero(M,mysplit[i][1..3],constr,[],[]): listparsol := ReduceToList(parsol): if ReduceToSequence(listparsol) <> NULL then if whattype(listparsol[1]) = list then tba := op(listparsol): else tba := listparsol: fi: if glob_out = [] then glob_out := [tba]: else tt := ReduceToList(glob_out): if tt <> [] then if whattype(tt[1]) = list then glob_out := [op(tt), tba]: else glob_out := [tt, tba]: fi: fi: fi: fi: ## this paragraph is to return [] when no solutions are computed if {all} union set_of_options <> set_of_options and glob_out <> [] then #glob_out := RestoreInputVariables(glob_out, invars, workingvars): if ReduceToSequence(glob_out) = NULL then return([]): fi: return([[ReduceToSequence(glob_out)]]): fi: fi: od: ## in case the polynomial q is IRreducible (zero-dim) else mysplit := PAR: mysplit := ReduceToList(mysplit): constr := mysplit[4]: parsol := ComputeSolutionsDimZero(M,mysplit[1..3],constr,[],[]): listparsol := ReduceToList(parsol): if ReduceToSequence(listparsol) <> NULL then if whattype(listparsol[1]) = list then tba := op(listparsol): else tba := listparsol: fi: if glob_out = [] then glob_out := [tba]; else tt := ReduceToList(glob_out): if tt <> [] then if whattype(tt[1]) = list then glob_out := [op(tt), tba]: else glob_out := [tt, tba]: fi: fi: fi: fi: ## this paragraph is to return [] when no solutions are computed if {all} union set_of_options <> set_of_options and glob_out <> [] then ## aggiunto questo :: rifletterci #glob_out := RestoreInputVariables(glob_out, invars, workingvars): if ReduceToSequence(glob_out) = NULL then return([]): fi: return([[ReduceToSequence(glob_out)]]): fi: fi: fi: fi: ## end of if avoid_further_dimension_computations = 0 ## Positive-dimensional case if avoid_further_dimension_computations = 1 then if dbg2810 in set_of_options then printf("Positive dimensional case\n"); fi: if hankelness = false then SampleLowRankSym(M,r): else SampleLowRankHankel(M,r): fi: ## this paragraph is to return [] when no solutions are computed if {all} union set_of_options <> set_of_options and glob_out <> [] then #glob_out := RestoreInputVariables(glob_out, invars, workingvars): if ReduceToSequence(glob_out) = NULL then return([]): fi: return([[ReduceToSequence(glob_out)]]): fi: fi: od: #glob_out := RestoreInputVariables(glob_out, invars, workingvars): glob_out := PutBracketIfAlone(glob_out): if ReduceToSequence(glob_out) = NULL then return([]): fi: return(glob_out): end proc: ######################################################################################## #------------------------------------------------------------------------------------------------# SampleLowRankSym := proc(oldM, oldr) ############################################################# local k, n, listaminori, varsX, L, locM, varsXY, HYP, SOLS, algdeg, PAR, clmi1, clmi2, s, EQS, var, varsY, myboolean, mysplit, i, ISOL, parsol, listparsol, tt, tba, listoflifts, liftVAR, nlift, j, l, constr; global list_of_ranks, set_of_options, regularity, glob_out, invars, workingvars, hankelness, char_pol, optimizing_vars: k := LinearAlgebra:-RowDimension(oldM): n := nops(indets(oldM)): listaminori := [op(combinat:-choose([seq(i, i=1..k)],k-oldr))]: for l from 1 to nops(listaminori) do ###varsX := [op(indets(oldM))]: ################# GOOOOOOOOD FIXED BUG! THIS COMMAND STAYS HERE! varsX := optimizing_vars: ## 25 aprile 2017 L := listaminori[l]: locM := oldM: EQS := IncidenceVarietySym(oldM, L, oldr): varsXY := [op(indets(EQS))]: varsY := [op({op(varsXY)} minus {op(varsX)})]: if ver in set_of_options then PrintVerbose(l,nops(listaminori),oldr): fi: ## checking regularity if reg in set_of_options then if CheckRegularity(EQS,varsXY) = false then regularity := regularity+1: WARNING("regularity assumptions not satisfied\n"): fi: fi: HYP := AffineCut(varsX): constr := CharPol(locM): PAR := LagSyst(EQS,HYP,varsX,varsY,n,constr): if PAR <> [] then myboolean := irreduc(PAR[1]): ## in case the polynomial q is reducible (positive-dim) if myboolean = false then mysplit := [SplitParametrization(PAR)]: mysplit := ReduceToList(mysplit): for i from 1 to nops(mysplit) do if whattype(mysplit[i][1]) <> integer then constr := mysplit[i][4]: parsol := ComputeSolutionsDimZero(locM,mysplit[i][1..3],constr,[],[]): listparsol := ReduceToList(parsol): if ReduceToSequence(listparsol) <> NULL then if whattype(listparsol[1]) = list then tba := op(listparsol): else tba := listparsol: fi: if glob_out = [] then glob_out := [tba]; else tt := ReduceToList(glob_out): if tt <> [] then if whattype(tt[1]) = list then glob_out := [op(tt), tba]: else glob_out := [tt, tba]: fi: fi: fi: fi: ## this paragraph is to return [] when no solutions are computed if {all} union set_of_options <> set_of_options and glob_out <> [] then #glob_out := RestoreInputVariables(glob_out, invars, workingvars): if ReduceToSequence(glob_out) = NULL then return([]): fi: return([[ReduceToSequence(glob_out)]]): fi: fi: od: ## in case the polynomial q is IRreducible (positive-dim) else mysplit := PAR: mysplit := ReduceToList(mysplit): constr := mysplit[4]: parsol := ComputeSolutionsDimZero(locM,mysplit[1..3],constr,[],[]): listparsol := ReduceToList(parsol): if ReduceToSequence(listparsol) <> NULL then if whattype(listparsol[1]) = list then tba := op(listparsol): else tba := listparsol: fi: if glob_out = [] then glob_out := [tba]; else tt := ReduceToList(glob_out): if tt <> [] then if whattype(tt[1]) = list then glob_out := [op(tt), tba]: else glob_out := [tt, tba]: fi: fi: fi: fi: ## this paragraph is to return [] when no solutions are computed if {all} union set_of_options <> set_of_options and glob_out <> [] then ## aggiunto questo :: rifletterci #glob_out := RestoreInputVariables(glob_out, invars, workingvars): if ReduceToSequence(glob_out) = NULL then return([]): fi: return([[ReduceToSequence(glob_out)]]): fi: fi: ## end of if(myboolean) = false fi: ## end of if PAR <> [] ## induction / recursion on variables listoflifts := []: liftVAR:=[]: while nops(varsX) > 0 do s := solve(HYP,varsX[1]): listoflifts := [op(listoflifts), s]: locM := subs({varsX[1] = s}, locM): EQS := subs(varsX[1]=s, EQS): ## here it has nops(varsX)-1 variables var := varsX[1]: liftVAR:=[op(liftVAR), var]: varsX := varsX[2..nops(varsX)]: constr := CharPol(locM): if nops(varsX) > 0 then HYP := AffineCut(varsX): PAR := LagSyst(EQS, HYP, varsX, varsY, nops(varsX), constr): if PAR <> [] then myboolean := irreduc(PAR[1]): ## in case the polynomial q is reducible (positive-dim-induction) if myboolean = false then mysplit := [SplitParametrization(PAR)]: mysplit := ReduceToList(mysplit): for i from 1 to nops(mysplit) do if whattype(mysplit[i][1]) <> integer then constr := mysplit[i][4]: parsol := ComputeSolutionsDimZero(oldM, mysplit[i][1..3],constr,liftVAR,listoflifts): listparsol := ReduceToList(parsol): if ReduceToSequence(listparsol) <> NULL then if whattype(listparsol[1]) = list then tba := op(listparsol): else tba := listparsol: fi: if glob_out = [] then glob_out := [tba]; else tt := ReduceToList(glob_out): if tt <> [] then if whattype(tt[1]) = list then glob_out := [op(tt), tba]: else glob_out := [tt, tba]: fi: fi: fi: fi: ## this paragraph is to return [] when no solutions are computed if {all} union set_of_options <> set_of_options and glob_out <> [] then ## aggiunto questo :: rifletterci #glob_out := RestoreInputVariables(glob_out, invars, workingvars): if ReduceToSequence(glob_out) = NULL then return([]): fi: return([[ReduceToSequence(glob_out)]]): fi: fi: ## ne rimane aperti 2 di if a questo punto (if nops(varsX) e if myboolean) od: ## buono fin qui ## in case the polynomial q is IRreducible (positive-dim-induction) else mysplit := PAR: mysplit := ReduceToList(mysplit): constr := mysplit[4]: parsol := ComputeSolutionsDimZero(oldM,mysplit[1..3],constr,liftVAR,listoflifts): ## errore qui da qualche parte listparsol := ReduceToList(parsol): if ReduceToSequence(listparsol) <> NULL then if whattype(listparsol[1]) = list then tba := op(listparsol): else tba := listparsol: fi: if glob_out = [] then glob_out := [tba]; else tt := ReduceToList(glob_out): if tt <> [] then if whattype(tt[1]) = list then glob_out := [op(tt), tba]: else glob_out := [tt, tba]: fi: fi: fi: fi: ## this paragraph is to return [] when no solutions are computed if {all} union set_of_options <> set_of_options and glob_out <> [] then ## aggiunto questo :: rifletterci #glob_out := RestoreInputVariables(glob_out, invars, workingvars): if ReduceToSequence(glob_out) = NULL then return([]): fi: return([[ReduceToSequence(glob_out)]]): fi: fi: ## end of if myboolean fi: ## end of if PAR <> [] fi: ## end of if nops(varsX) od: ## end of while(nops(varsX)) > 0 od: ## end of for L in listaminori return(glob_out): end proc: ######################################################################################## #------------------------------------------------------------------------------------------------# PreprocessInput := proc(MMM) ##################################################################### local MM, myvars, n, s, relat, count, solut, i, j, newM: global list_of_ranks, set_of_options, regularity, glob_out, invars, workingvars, hankelness, char_pol, optimizing_vars: MM := MMM: myvars := [op(indets(MM))]: n := nops(myvars): ## output specifications s := RowDimension(MM): ## 0 : if there is a zero on the diagonal, and a non zero integer in the same line ## (M,relat) : where M is the new matrix and relat are the relations if n = 0 then return(MM): fi: relat := []: count := 0: for i from 1 to s do if MM[i,i] = 0 then count := count + 1: for j from 1 to s do if whattype(MM[i,j]) = integer and MM[i,j] <> 0 then return(MM,0): elif whattype(MM[i,j]) <> integer then relat := [op(relat), MM[i,j] = 0]: fi: od: fi: od: if count <> 0 then WARNING("there are zeros on the diagonal, building a simpler LMI\n"): fi: solut := solve(relat,myvars): newM := subs(op(solut),MM): return(newM,relat): end proc: ######################################################################################## #------------------------------------------------------------------------------------------------# PrintVerbose := proc(ll, nn, rr) ################################################################# local p; global list_of_ranks, set_of_options, regularity, glob_out, invars, workingvars, hankelness, char_pol, optimizing_vars: p := trunc(((ll-1)/nn)*100): printf("Rank %d :: completed to %d%%\n", rr, p): #if ll = 1 then # printf("Rank %d :: completed %d%%", rr, p): #elif ll = nn then # printf(" 100%%\n"): #else # printf(" %d%%", p): #fi: return(0): end proc: ######################################################################################## #------------------------------------------------------------------------------------------------# SampleLowRankHankel := proc(oldM, oldr) ########################################################## local k, n, varsX, locM, EQS, varsXY, varsY, HYP, PAR, myboolean, mysplit, i, ISOL, parsol, listparsol, tba, tt, listoflifts, liftVAR, s, var, nlift, j, constr; global list_of_ranks, set_of_options, regularity, glob_out, invars, workingvars, hankelness, char_pol, optimizing_vars: k := LinearAlgebra:-RowDimension(oldM): n := nops(indets(oldM)): ###varsX := [op(indets(oldM))]: varsX := optimizing_vars: ## 25 aprile 2017 locM := oldM: EQS := IncidenceVarietyHankel(oldM,oldr): varsXY := [op(indets(EQS))]: varsY := [op({op(varsXY)} minus {op(varsX)})]: ## checking regularity if reg in set_of_options then if CheckRegularity(EQS,varsXY) = false then regularity := regularity+1: WARNING("regularity assumptions not satisfied\n"): fi: fi: HYP := AffineCut(varsX): constr := CharPol(locM): PAR := LagSyst(EQS,HYP,varsX,varsY,n, constr): if PAR <> [] then myboolean := irreduc(PAR[1]): ## in case the polynomial q is reducible (positive-dim) if myboolean = false then mysplit := [SplitParametrization(PAR)]: mysplit := ReduceToList(mysplit): for i from 1 to nops(mysplit) do if whattype(mysplit[i][1]) <> integer then constr := mysplit[i][4]: parsol := ComputeSolutionsDimZero(locM,mysplit[i][1..3],constr,[],[]): listparsol := ReduceToList(parsol): if ReduceToSequence(listparsol) <> NULL then if whattype(listparsol[1]) = list then tba := op(listparsol): else tba := listparsol: fi: if glob_out = [] then glob_out := [tba]; else tt := ReduceToList(glob_out): if tt <> [] then if whattype(tt[1]) = list then glob_out := [op(tt), tba]: else glob_out := [tt, tba]: fi: fi: fi: fi: ## this paragraph is to return [] when no solutions are computed if {all} union set_of_options <> set_of_options and glob_out <> [] then ## aggiunto questo :: rifletterci #glob_out := RestoreInputVariables(glob_out, invars, workingvars): if ReduceToSequence(glob_out) = NULL then return([]): fi: return([[ReduceToSequence(glob_out)]]): fi: fi: od: ## in case the polynomial q is IRreducible (positive-dim) else mysplit := PAR: mysplit := ReduceToList(mysplit): constr := mysplit[4]: parsol := ComputeSolutionsDimZero(locM,mysplit[1..3],constr,[],[]): listparsol := ReduceToList(parsol): if ReduceToSequence(listparsol) <> NULL then if whattype(listparsol[1]) = list then tba := op(listparsol): else tba := listparsol: fi: if glob_out = [] then glob_out := [tba]; else tt := ReduceToList(glob_out): if tt <> [] then if whattype(tt[1]) = list then glob_out := [op(tt), tba]: else glob_out := [tt, tba]: fi: fi: fi: fi: ## this paragraph is to return [] when no solutions are computed if {all} union set_of_options <> set_of_options and glob_out <> [] then ## aggiunto questo :: rifletterci #glob_out := RestoreInputVariables(glob_out, invars, workingvars): if ReduceToSequence(glob_out) = NULL then return([]): fi: return([[ReduceToSequence(glob_out)]]): fi: fi: ## end of if(myboolean) = false fi: ## end of if PAR <> [] ## induction / recursion on variables listoflifts := []: liftVAR:=[]: while nops(varsX) > 0 do s := solve(HYP,varsX[1]): listoflifts := [op(listoflifts), s]: locM := subs({varsX[1] = s}, locM): EQS := subs(varsX[1]=s, EQS): ## here it has nops(varsX)-1 variables var := varsX[1]: liftVAR:=[op(liftVAR), var]: varsX := varsX[2..nops(varsX)]: if nops(varsX) > 0 then HYP := AffineCut(varsX): constr := CharPol(locM): PAR := LagSyst(EQS, HYP, varsX, varsY, nops(varsX), constr): if PAR <> [] then myboolean := irreduc(PAR[1]): ## in case the polynomial q is reducible (positive-dim-induction) if myboolean = false then mysplit := [SplitParametrization(PAR)]: mysplit := ReduceToList(mysplit): for i from 1 to nops(mysplit) do if whattype(mysplit[i][1]) <> integer then constr := mysplit[i][4]: parsol := ComputeSolutionsDimZero(oldM, mysplit[i][1..3], constr,liftVAR,listoflifts): listparsol := ReduceToList(parsol): if ReduceToSequence(listparsol) <> NULL then if whattype(listparsol[1]) = list then tba := op(listparsol): else tba := listparsol: fi: if glob_out = [] then glob_out := [tba]; else tt := ReduceToList(glob_out): if tt <> [] then if whattype(tt[1]) = list then glob_out := [op(tt), tba]: else glob_out := [tt, tba]: fi: fi: fi: fi: ## this paragraph is to return [] when no solutions are computed if {all} union set_of_options <> set_of_options and glob_out <> [] then ## aggiunto questo :: rifletterci #glob_out := RestoreInputVariables(glob_out, invars, workingvars): if ReduceToSequence(glob_out) = NULL then return([]): fi: return([[ReduceToSequence(glob_out)]]): fi: fi: ## ne rimane aperti 2 di if a questo punto (if nops(varsX) e if myboolean) od: ## buono fin qui ## in case the polynomial q is IRreducible (positive-dim-induction) else mysplit := PAR: mysplit := ReduceToList(mysplit): constr := mysplit[4]: parsol := ComputeSolutionsDimZero(oldM,mysplit[1..3],constr,liftVAR,listoflifts): listparsol := ReduceToList(parsol): if ReduceToSequence(listparsol) <> NULL then if whattype(listparsol[1]) = list then tba := op(listparsol): else tba := listparsol: fi: if glob_out = [] then glob_out := [tba]; else tt := ReduceToList(glob_out): if tt <> [] then if whattype(tt[1]) = list then glob_out := [op(tt), tba]: else glob_out := [tt, tba]: fi: fi: fi: fi: ## this paragraph is to return [] when no solutions are computed if {all} union set_of_options <> set_of_options and glob_out <> [] then ## aggiunto questo :: rifletterci #glob_out := RestoreInputVariables(glob_out, invars, workingvars): if ReduceToSequence(glob_out) = NULL then return([]): fi: return([[ReduceToSequence(glob_out)]]): fi: fi: ## end of if myboolean fi: ## end of if PAR <> [] fi: ## end of if nops(varsX) od: ## end of while(nops(varsX)) > 0 return(glob_out): end proc: ######################################################################################## #------------------------------------------------------------------------------------------------# LagSyst := proc(oldeqs, oldhyp, oldvarsX, oldvarsY, oldn, oldconstr) ############################# local EQS, varsZ, HYP, varsXY, NEWEQS, GB, PAR, newconstr, i: global list_of_ranks, set_of_options, regularity, glob_out, invars, workingvars, hankelness, char_pol, optimizing_vars: EQS := oldeqs: HYP := oldhyp: varsZ := [seq(cat('LL',i),i=1..nops(EQS))]: varsXY := [op(oldvarsX), op(oldvarsY)]: NEWEQS := [op(EQS), seq(add(varsZ[i]*diff(EQS[i],varsXY[j]),i=1..nops(EQS))-diff(HYP,varsXY[j]), j=1..nops(varsXY))]: if ver in set_of_options then lprint(NEWEQS, zerochar, [op(varsZ), op(oldvarsY)], oldvarsX); fi: if ver in set_of_options then myverb := 1: printf("Beginning of Groebner basis computation -- Potentially costly ... \n"); fi: GB := mygbasiselim(NEWEQS, zerochar, [op(varsZ), op(oldvarsY)], oldvarsX, {"verb"=myverb}): if ver in set_of_options then printf("... End of Groebner basis computation\n"); fi: if GB = [1] or GB = [] then PAR := []: else newconstr := expand(oldconstr): for i from 1 to nops(newconstr) do newconstr[i] := abs(denom(newconstr[i]))*newconstr[i]: od: # printf("\n HERE THE ERROR \n"); if ver in set_of_options then printf("Beginning of Computation of Rational Parametrization -- Potentially costly ... \n"); fi: PAR := [fgbrs:-rs_rur_gb(GB, [op(indets(GB))], constraints=newconstr)]: if ver in set_of_options then printf("... End of Computation of Rational Parametrization\n"); fi: fi: return(PAR): end proc: ######################################################################################## #------------------------------------------------------------------------------------------------# IncidenceVarietySym := proc(oldM, oldL, oldr) ####################################################### local matY, locM, k, i, j, J, B, no_redund: global list_of_ranks, set_of_options, regularity, glob_out, invars, workingvars, hankelness, char_pol, optimizing_vars: locM := oldM: k := LinearAlgebra:-RowDimension(locM): matY := Matrix(k,k-oldr): for i from 1 to k do for j from 1 to k-oldr do matY[i,j] := cat('YY',i,j): od: od: matY[oldL,[1..k-oldr]] := IdentityMatrix(k-oldr): J := LinearAlgebra:-Multiply(locM,matY): B := J[oldL, [1..k-oldr]]: for i from 1 to k-oldr do for j from i+1 to k-oldr do B[i,j]:=0: od: od: J[oldL, [1..k-oldr]] := B: no_redund := []: for i from 1 to k do for j from 1 to k-oldr do if J[i,j] <> 0 then no_redund := [op(no_redund), J[i,j]]: fi: od: od: #if ver in set_of_options then # printf("Dimension Inc Var: %d\n", ComputeDimensionIntrinsec(no_redund)): #fi: return(no_redund); end proc: ######################################################################################## #------------------------------------------------------------------------------------------------# IncidenceVarietyHankel := proc(oldM, oldr) ####################################################### local locM, k, matY, s, i, no_redund, J: global list_of_ranks, set_of_options, regularity, glob_out, invars, workingvars, hankelness, char_pol, optimizing_vars: locM := oldM: k := LinearAlgebra:-RowDimension(locM): matY := Matrix(k,k-oldr): for s from 1 to oldr+1 do for i from 1 to k-oldr do matY[i+s-1,i] := cat('YY',s): od: od: J := LinearAlgebra:-Multiply(locM,matY): no_redund := [Vectorize(J)]: return(no_redund): end proc: ######################################################################################## #------------------------------------------------------------------------------------------------# CharPolPar := proc(MM, PAR) ###################################################################### local s, locM, mycharpol, st, ft, myvars, n, ddenom, numers, i: global list_of_ranks, set_of_options, regularity, glob_out, invars, workingvars, hankelness, char_pol, optimizing_vars: st := time(): mycharpol := 0: ddenom := PAR[2]: numers := PAR[3]: locM := MM: myvars := [op(indets(locM))]: n := nops(myvars): if nops(numers) <> n then WARNING("impossible to build characteristic polynomial\n"): else for i from 1 to n do locM := subs({myvars[i] = numers[i]/ddenom}, locM): od: ft := time()-st: printf("\n\n Computation of char pol in time : %.2f\n\n", ft): fi: return(0): end proc: ######################################################################################## #------------------------------------------------------------------------------------------------# CharPol := proc(MM) ############################################################################## local mycharpol, k, s, ll, i: global list_of_ranks, set_of_options, regularity, glob_out, invars, workingvars, hankelness, char_pol, optimizing_vars: k := RowDimension(MM): if k <> ColumnDimension(MM) then error "the characteristic polynomial cannot be computed": fi: mycharpol := simplify(det(s*IdentityMatrix(k)+MM)): ll := [subs({s=0},mycharpol)]: for i from 1 to k-1 do ll := [op(ll),coeff(mycharpol,s^i)]: od: ## it gives in output [p_k,p_{k-1},p_{k-2},...,p_1] coefficients of mycharpol ## except the first one which is one :: they are k where k is the size of the matrix ## charpol = s^k + p_1*s^{k-1} + ... + s*p_{k-1} + p_k return(expand(ll)): end proc: ######################################################################################## #------------------------------------------------------------------------------------------------# PutBracketIfAlone := proc(glo) ################################################################### local locglo, i: global list_of_ranks, set_of_options, regularity, glob_out, invars, workingvars, hankelness, char_pol, optimizing_vars: locglo := ReduceToList(glo): for i from 1 to nops(locglo) do if whattype(locglo[i]) <> list then return([locglo]): fi: od: return(locglo): end proc: ######################################################################################## #------------------------------------------------------------------------------------------------# ReduceToSequence := proc() ####################################################################### local mm: global list_of_ranks, set_of_options, regularity, glob_out, invars, workingvars, hankelness, char_pol, optimizing_vars: if nargs = 1 then mm:=args[1]: else return(args): fi: if whattype(mm) <> list then return(mm): fi: while whattype(op(mm)) = list do mm:=op(mm): od: return(op(mm)): end proc: ######################################################################################## #------------------------------------------------------------------------------------------------# ReduceToList := proc() ########################################################################### local mm: global list_of_ranks, set_of_options, regularity, glob_out, invars, workingvars, hankelness, char_pol, optimizing_vars: if nargs = 1 then mm:=args[1]: else return(args): fi: if whattype(mm) <> list then return(mm): fi: while whattype(op(mm)) = list do mm:=op(mm): od: return(mm): end proc: ######################################################################################## #------------------------------------------------------------------------------------------------# ComputeSolutionsDimZero := proc(oldM, oldpar, oldconstr, liftVAR, listoflifts) ################### local gg, i, j, k, varsX, locM, locQ, locdeg, locpar, v, s, a, b, myout, thisrank, myboole, oldQ, val_of_charpol, constr, nzero, npos, nneg, listsigns, l, valval, nlift, thissol, activevar: global list_of_ranks, set_of_options, regularity, glob_out, invars, workingvars, hankelness, char_pol, optimizing_vars: #oldQ := myisolate_rur(oldpar, [op(indets(oldM))]): varsX := [op(indets(oldM))]: locM := oldM: constr := oldconstr: activevar := [op(indets(oldM) minus {op(liftVAR)})]: locQ, val_of_charpol := myisolate_rur(oldpar, activevar, constr): if oldQ = [] then return([]): fi: locpar := oldpar: locdeg := degree(locpar[1]): myout := []: for i from 1 to nops(locQ) do locM := oldM: ### begin of psd test using coefficients of char_pol ########################################## nzero := 0: # npos := 0: # nneg := 0: # listsigns := []: # valval := val_of_charpol[i]: # # for j from 1 to nops(valval) do # if valval[j][1] = 0 and valval[j][2] = 0 then # nzero := nzero + 1: # listsigns := [op(listsigns), 0]: # elif sign(valval[j][1])*sign(valval[j][2]) = -1 then # nzero := nzero + 1: # listsigns := [op(listsigns), 0]: # elif sign(valval[j][1]) = 1 and sign(valval[j][2]) = 1 then # npos := npos + 1: # listsigns := [op(listsigns), 1]: # elif sign(valval[j][1]) = -1 and sign(valval[j][2]) = -1 then # nneg := nneg + 1: # listsigns := [op(listsigns), -1]: # fi: # od: # # if nzero + npos = RowDimension(oldM) then # myboole := true: # else # myboole := false: # fi: # # thisrank := RowDimension(oldM)-nzero: # ### end of psd test using coefficients of char_pol ############################################ ### begin of numerical psd test ############################################################### #for j from 1 to nops(locQ[i]) do # # locM := subs({varsX[j] = ([op(locQ[i][j])][2][1] + [op(locQ[i][j])][2][2])/2}, locM): # #od: # # #st := time(): # #myboole := NumericalCheckPSD(locM); # # #ft := time()-st: # # #printf("Time of num check psd :: %f\n", ft): # #### end of numerical psd test ################################################################ ## if locQ[i] is a feasible point if myboole = true then thissol := [locQ[i]]: ## recompute fiber variables for this solution nlift := nops(liftVAR): if nlift > 0 then for j from 1 to nlift do thissol := map(_point->[liftVAR[nlift-j+1] = subs_interval(_point, listoflifts[nlift-j+1]), op(_point)], thissol): od: fi: ## end of recompute fiber variables for this solution ## when rnk and deg are in set_of_options if {rnk,deg} subset set_of_options then if ver in set_of_options then printf("Found a solution of rank %d, and algebraic degree %d\n", thisrank, locdeg): fi: if par in set_of_options then locpar := elim_factors(locpar): myout := [op(myout), [ReduceToSequence(thissol), rnk = thisrank, deg = locdeg, par = [op(locpar)]]]: else myout := [op(myout), [ReduceToSequence(thissol), rnk = thisrank, deg = locdeg]]: fi: if {all} union set_of_options <> set_of_options then return(myout): fi: ## when rnk is in set_of_options elif rnk in set_of_options then if ver in set_of_options then printf("Found a solution of rank %d\n", thisrank): fi: if par in set_of_options then locpar := elim_factors(locpar): myout := [op(myout), [ReduceToSequence(thissol), rnk = thisrank, par = [op(locpar)]]]: else myout := [op(myout), [ReduceToSequence(thissol), rnk = thisrank]]: fi: if {all} union set_of_options <> set_of_options then return(myout): fi: ## when deg is in set_of_options elif deg in set_of_options then if ver in set_of_options then printf("Found a solution of algebraic degree %d\n", locdeg): fi: if par in set_of_options then locpar := elim_factors(locpar): myout := [op(myout), [ReduceToSequence(thissol), deg = locdeg, par = [op(locpar)]]]: else myout := [op(myout), [ReduceToSequence(thissol), deg = locdeg]]: fi: if {all} union set_of_options <> set_of_options then return(myout): fi: ## when par is in set_of_options elif par in set_of_options then locpar := elim_factors(locpar): if ver in set_of_options then printf("Found a solution parametrized by: "): lprint(locpar): fi: myout := [op(myout), [ReduceToSequence(thissol), par = [op(locpar)]]]: if {all} union set_of_options <> set_of_options then return(myout): fi: ## otherwise else myout := [op(myout), [ReduceToSequence(thissol)]]: if {all} union set_of_options <> set_of_options then return(myout): fi: fi: fi: od: return(myout): end proc: #------------------------------------------------------------------------------------------------# ComputeDimensionIntrinsec := proc(equations) ##################################################### local var, n, eqs, GB, d: global list_of_ranks, set_of_options, regularity, glob_out, invars, workingvars, hankelness, char_pol, optimizing_vars: var := [op(indets(equations))]: n := nops(var): eqs := equations: GB := mygbasis(eqs, luckyprime, var, [], {"verb"=myverb}): if GB = [1] then return(-1): fi: for d from 1 to n do eqs:=[op(eqs), AffineCut(var)]: GB:=mygbasis(eqs,luckyprime,var,[],{"verb"=myverb}): if GB=[1] then return(d-1): fi: od: end proc: ######################################################################################## #------------------------------------------------------------------------------------------------# MyRand := proc() ################################################################################# local X, t: global list_of_ranks, set_of_options, regularity, glob_out, invars, workingvars, hankelness, char_pol, optimizing_vars: X := RandomVariable(Uniform(1,1000)): ##t := ((-1)^(rand(-2..1)()))*floor(Sample(X,1)[1]): ## anche numeri negativi t:= abs(((-1)^(rand(-2..1)()))*floor(Sample(X,1)[1])): ## solo numeri positivi return(t): end proc: ######################################################################################## #------------------------------------------------------------------------------------------------# AffineCut := proc(varvar) ######################################################################## local t: global list_of_ranks, set_of_options, regularity, glob_out, invars, workingvars, hankelness, char_pol, optimizing_vars: t := MyRand()+add(MyRand()*varvar[i],i=1..nops(varvar)): return(t): end proc: ######################################################################################## #------------------------------------------------------------------------------------------------# CheckDimension := proc(equations) ################################################################ local GB, eqs, hilb, T: global list_of_ranks, set_of_options, regularity, glob_out, invars, workingvars, hankelness, char_pol, optimizing_vars: ## Possible Outputs ## -1 : empty eqs := equations: ## 0 : zero-dimensional GB := mygbasis(eqs, 0, [op(indets(eqs))], [], {}): ## false : positive dim if GB = [1] then return(-1): else hilb := FGb:-fgb_hilbert(GB,0, [], [op(indets(eqs))], 'T'): fi: if hilb[2] = 0 then return(0): else return(false): fi: end proc: ######################################################################################## #------------------------------------------------------------------------------------------------# ComputeMinors := proc(s,K) ####################################################################### local cc, rr, rrlist, cclist, outlist, i, j: global list_of_ranks, set_of_options, regularity, glob_out, invars, workingvars, hankelness, char_pol, optimizing_vars: rr := RowDimension(K): cc := ColumnDimension(K): rrlist := [op(combinat:-choose([seq(i, i=1..rr)],s))]: cclist := [op(combinat:-choose([seq(i, i=1..cc)],s))]: outlist := {}: for i from 1 to nops(rrlist) do for j from 1 to nops(cclist) do outlist := outlist union {det(K(rrlist[i],cclist[j]))}: od: od: return(op(outlist)): end proc: ######################################################################################## #------------------------------------------------------------------------------------------------# ComputeMinors2 := proc(s,K) ###################################################################### local nn, S, KK: global list_of_ranks, set_of_options, regularity, glob_out, invars, workingvars, hankelness, char_pol, optimizing_vars: KK := K: nn := RowDimension(K): S := combinat[choose](nn,s); return(seq(seq(Determinant(SubMatrix(KK, r, c)), r = S), c = S)); end proc: ######################################################################################## #------------------------------------------------------------------------------------------------# CheckPSD := proc(QQ,varsX) ####################################################################### local llinf, llsup, s, i, j, sgninf, sgnsup; global list_of_ranks, set_of_options, regularity, glob_out, invars, workingvars, hankelness, char_pol, optimizing_vars: llinf := char_pol: llsup := char_pol: for j from 1 to nops(QQ) do llinf := subs({varsX[j] = [op(QQ[j])][2][1]}, llinf): llsup := subs({varsX[j] = [op(QQ[j])][2][2]}, llsup): od: sgninf := SequenceOfSigns(llinf): sgnsup := SequenceOfSigns(llsup): if ver in set_of_options then printf("CHECKING PSD\n"): printf("values of minors on inf: "): lprint(evalf(llinf)): printf("signs of minors on inf: "): lprint(sgninf): printf("values of minors on sup: "): lprint(evalf(llsup)): printf("signs of minors on sup: "): lprint(sgnsup): fi: for i from 1 to nops(sgninf) do if sgninf[i] = sgnsup[i] and llinf[i] < 0 then if ver in set_of_options then printf("Solution discarded\n\n"): fi: return(false): fi: od: if ver in set_of_options then printf("Solution accepted\n\n"): fi: return(true): end proc: ######################################################################################## #------------------------------------------------------------------------------------------------# SequenceOfSigns := proc(inlist) ################################################################## local ll, i, mysigns; global list_of_ranks, set_of_options, regularity, glob_out, invars, workingvars, hankelness, char_pol, optimizing_vars: ll := inlist: mysigns := []: for i from 1 to nops(ll) do mysigns := [op(mysigns), sign(ll[i])]: od: return(mysigns) end proc: ######################################################################################## #------------------------------------------------------------------------------------------------# ComputeTheRank := proc(QQ,varsX) ################################################################# local llinf, llsup, i, j, sgninf, sgnsup, rr; global list_of_ranks, set_of_options, regularity, glob_out, invars, workingvars, hankelness, char_pol, optimizing_vars: llinf := char_pol: llsup := char_pol: for j from 1 to nops(QQ) do llinf := subs({varsX[j] = [op(QQ[j])][2][1]}, llinf): llsup := subs({varsX[j] = [op(QQ[j])][2][2]}, llsup): od: sgninf := SequenceOfSigns(llinf): sgnsup := SequenceOfSigns(llsup): rr := 0: for i from 1 to nops(sgninf) do if sgninf[i] = sgnsup[i] then rr := rr+1: fi: od: return(rr): end proc: ######################################################################################## #------------------------------------------------------------------------------------------------# NumericalCheckPSD := proc(MMM) ################################################################### local rr, s, i, eige, MM, v, k: global list_of_ranks, set_of_options, regularity, glob_out, invars, workingvars, hankelness, char_pol, optimizing_vars: MM := MMM: eige := [evalf(eigenvalues(MM))]: s := nops(eige): for i from 1 to s do if sign(eige[i]) = -1 and abs(eige[i]) > myzerotest then return(false): fi: od: return(true): end proc: ######################################################################################## #------------------------------------------------------------------------------------------------# NumericalRankOfMatrix := proc(MMM) ############################################################### local rr,s, c, i, eige, MM, v, k: global list_of_ranks, set_of_options, regularity, glob_out, invars, workingvars, hankelness, char_pol, optimizing_vars: MM := MMM: eige := [evalf(eigenvalues(MM))]: s := RowDimension(MM): c:=0: for i from 1 to s do if abs(eige[i]) < myzerotest then c:=c+1: fi: od: return(s-c): end proc: ######################################################################################## #------------------------------------------------------------------------------------------------# CheckHankelness := proc(H) ####################################################################### local m, i, j, c, s: global list_of_ranks, set_of_options, regularity, glob_out, invars, workingvars, hankelness, char_pol, optimizing_vars: m := LinearAlgebra:-RowDimension(H): if m = 1 then return(true): fi: for i from 2 to m do for j from 1 to m-1 do if H[i,j]<>H[i-1,j+1] then return(false): fi: od: od: return(true): end proc: ######################################################################################## #------------------------------------------------------------------------------------------------# Vectorize := proc(MMM) ########################################################################### local s, t, vec, i, j: global list_of_ranks, set_of_options, regularity, glob_out, invars, workingvars, hankelness, char_pol, optimizing_vars: s := RowDimension(MMM): t := ColumnDimension(MMM): vec := []: for i from 1 to t do for j from 1 to s do vec := [op(vec),expand(MMM[j,i])]: od: od: return(op(vec)): end proc: ######################################################################################## #------------------------------------------------------------------------------------------------# SplitParametrization_old := proc(PAR) ############################################################ local q, k, i, listofpars, listoffactors, ll, j, nll: global list_of_ranks, set_of_options, regularity, glob_out, invars, workingvars, hankelness, char_pol, optimizing_vars: q := factor(PAR[1]): if irreduc(q) = false then k := nops(q): if k = 1 then return(PAR): else ##### new stuff #ll := factors(q)[2]: #nll := nops(ll): #listofpars := []: #for j from 1 to nll do # listofpars := [op(listofpars), [ll[j][1], op(PAR[2..nops(PAR)])]]: #od: ##### end new stuff listoffactors := op(q): listofpars := []: for i in listoffactors do listofpars := [op(listofpars), [i, op(PAR[2..nops(PAR)])]]: od: return(listofpars): fi: else return(PAR): fi: end proc: ######################################################################################## #------------------------------------------------------------------------------------------------# SplitParametrization := proc(PAR) ################################################################ local q, k, i, listofpars, listoffactors, ll, j, nll: global list_of_ranks, set_of_options, regularity, glob_out, invars, workingvars, hankelness, char_pol, optimizing_vars: q := factor(PAR[1]): if irreduc(q) = false then k := nops(q): if k = 1 then return(PAR): else ##### new stuff ll := factors(q)[2]: nll := nops(ll): listofpars := []: for j from 1 to nll do listofpars := [op(listofpars), [ll[j][1], op(PAR[2..nops(PAR)])]]: od: ##### end new stuff #listoffactors := op(q): #listofpars := []: #for i in listoffactors do # listofpars := [op(listofpars), [i, op(PAR[2..nops(PAR)])]]: #od: return(listofpars): fi: else return(PAR): fi: end proc: ######################################################################################## #------------------------------------------------------------------------------------------------# myisolate_gb := proc(mygb, myvars) ############################################################### local t: global list_of_ranks, set_of_options, regularity, glob_out, invars, workingvars, hankelness, char_pol, optimizing_vars: t := map(_point->[seq(myvars[i]=_point[i], i=1..nops(myvars))], fgbrs:-rs_isolate_gb(mygb,myvars)): return(t): end proc: ######################################################################################## #------------------------------------------------------------------------------------------------# myisolate_rur := proc(myrur, myvars, constr) ##################################################### local locrur, isol1, isol2: global list_of_ranks, set_of_options, regularity, glob_out, invars, workingvars, hankelness, char_pol, optimizing_vars: locrur := myrur: if locrur = [] then return([]): fi: if whattype(locrur) <> list then locrur := [locrur]: fi: #if ver in set_of_options then # printf("This is the input of rs_isolate_rur : "): # lprint(locrur): #fi: locrur[1] := denom(locrur[1])*locrur[1]: isol1, isol2 := fgbrs:-rs_isolate_rur(op(locrur), constraints=constr): return(map(_point->[seq(myvars[i]=_point[i], i=1..nops(myvars))], isol1), isol2): end proc: ######################################################################################## #------------------------------------------------------------------------------------------------# RestoreInputVariables := proc(oldout, inputvars, outputvars) ##################################### local i, locout; global list_of_ranks, set_of_options, regularity, glob_out, invars, workingvars, hankelness, char_pol, optimizing_vars: locout := oldout: for i from 1 to nops(outputvars) do locout := subs({outputvars[i] = inputvars[i]}, locout): od: return(locout): end proc: ######################################################################################## #------------------------------------------------------------------------------------------------# subs_interval := proc(_point, lin) ############################################################### local myvars, _intervals, _l, i, a, b: global list_of_ranks, set_of_options, regularity, glob_out, invars, workingvars, hankelness, char_pol, optimizing_vars: # On suppose que tous les coefficients de lin sont positifs. (pas vrai) # printf("subsinterval with %d variables\n", nops(indets(lin))): # print(lin): myvars:=[op(indets(_point))]: _intervals:=map(_p-> rhs(_p), _point): for i from 1 to nops(_intervals) do _l:=_intervals[i]: if abs(_l[1])0 then #printf("positive\n"): #a:=subs(_l[1], lin): ## Mohab #b:=subs(_l[2], lin): ## Mohab a:=subs(myvars[i]=_l[1], a): ## Simone b:=subs(myvars[i]=_l[2], b): ## Simone #return [a, b]: ## Mohab else #printf("negative\n"): #b:=subs(_l[1], lin): ## Mohab #a:=subs(_l[2], lin): ## Mohab b:=subs(myvars[i]=_l[1], b): ## Simone a:=subs(myvars[i]=_l[2], a): ## Simone #return [a, b]: ## Mohab fi: od: return([a, b]): ## Simone end proc: ######################################################################################## #------------------------------------------------------------------------------------------------# elim_factors := proc(oldpar) ##################################################################### local locpar, d, den, di, i, ll, gg: global list_of_ranks, set_of_options, regularity, glob_out, invars, workingvars, hankelness, char_pol, optimizing_vars: locpar := oldpar: ## for the moment I eliminate integer factors den := locpar[3]: d := igcd(coeffs(expand(locpar[2]))): if d = 1 then return(oldpar): fi: ll := []: for i from 1 to nops(den) do di := igcd(coeffs(den[i])): if di = 1 then return(oldpar): else ll := [op(ll), di]: fi: od: gg := igcd(d,op(ll)): locpar[2] := locpar[2]/gg: locpar[3] := den/gg: return(locpar): end proc: ######################################################################################## end module: savelib(SPECTRA, "SPECTRA.mla"):