with(combinat):with(ListTools): with(ArrayTools): with(LinearAlgebra):with(convex): # For a given monomial, the number of mons with degree difference being at most 1 is: # sum_{k=1}^floor(n/2) (2*k)! choose(n,2k) + 2 * sum_{k=1}^ceil(n/2) (2*k-1)! choose(n,2k-1) # something close to n/2 * factorial(n) monshom := proc(n,d) local l1; l1 := [seq(1,i=1..n)]; return map(_c -> _c - l1, composition(n+d, n)): end; # For both revmons functions, the bottlneck is the composition function from package combinat revmons := proc(n,d) if d = 0 then return [[seq(0,i=1..n)]]: else return [op(monshom(n,d)), op(revmons(n,d-1))]: fi: end; revmons2 := proc(n,d) local mons0; mons0 := [[seq(0,i=1..n)]]: return [seq(op(monshom(n,d-i+1)),i=1..d),mons0] end; monspt := proc(n,d) return Reverse(revmons(n,d)): end; mons := proc(n,d,X,ms::list := monspt(n,d)) return map(a->mul(X[i]^a[i],i=1..n),ms); # return [seq(mul(X^~a), a in ms)]; end; benchmons := proc(ms) local m, i; #ms := mons(n,d); for i from 1 to nops(ms) do m := ms[i]; Search(m,ms): od: return 0: end; support := proc (f,X,monslist) local cf, idxf,monsf; cf := coeffs(f,X,'monsf'); monsf := [monsf]; idxf := seq(Search(m,monslist),m in monsf); return [cf],[idxf]; end; dense_perturbation := proc(ms) return add(m^2,m in ms); end; sos2sdp := proc(r,X,relaxorder, mspt,msk,rc,ridx,precSVD::integer:=10,precSDP::integer:=200,epsStar::integer:=3,precRound::integer:=30,dataprec::integer:=10,gmp::boolean:=false,algo::integer:=1,g::list:=[],gc::list:=[],gidx::list:=[]) local start_idx, sub_idx,nvars,zero,d,k,n,ng,r0,msdp,nsdp,nblock,fd,rList,i,rfloat,j,lowerbnd,eigs,eivs,eigs2,eivs2,Y, Yrat, absorb_tbl,m, absorb_PP,nsdpg,nsdponeg,gic,giidx,ni,ig,igc,gifloat,mgi,eigsi,eivsi; k := relaxorder; d = 2*k: n := nops(X); ng := nops(g): zero := [seq(0,i=1..n)]: msdp := nops(mspt); nsdp := nops(msk); nblock := 1+ng: fd := fopen("in.dat-s",WRITE,TEXT); rList := Array([seq(0,i=1..msdp)]): for i from 1 to nops(rc) do rList[ridx[i]]:=rc[i]: od: if mspt[1] = zero then nvars := msdp - 1: start_idx := 2: r0 := rList[1]: else sub_idx := Search(2*mspt[1],mspt)-1: nvars := msdp - sub_idx: start_idx := 1: r0 := 0: fi: # test whether 0 belongs the multi-index list writeline(fd,convert(nvars,string)); # Number of moment variables - 1 (since y0 = 1) writeline(fd,convert(nblock,string)); # Number of blocks: 1 for the moment matrix fprintf(fd,"%d ",nsdp): # Size of the first block: the moment matrix has size nsdp nsdpg := locmatsizes(g,k,n): #lprint(nsdpg); for ni in nsdpg do fprintf(fd,"%d ",ni): od: # Size of the localizing matrices: for a polynomial gi of degree di, the size is binomial(n+k-floor(di/2),n) fprintf(fd,"\n"); rfloat := convert(rList,float,dataprec): #lprint(rfloat); if mspt[1] = zero then for i from start_idx to msdp do fprintf(fd,"%s ",convert(rfloat[i],string));od; else for i from sub_idx + 1 to msdp do fprintf(fd,"%s ",convert(rfloat[i],string));od; fi: fprintf(fd,"\n"); ### The case when y0 = 1 because 0 belongs to the multi-index list # starts with the first row of the moment matrix if mspt[1] = zero then writeline(fd,"0 1 1 1 -1"): # Moment variable y0 = 1 for j from 2 to nsdp do fprintf(fd, "%d 1 %d %d 1\n", j-1, 1, j); # Moment variable y_{j-1} at first row and column j od; # then the other rows of the moment matrix for i from 2 to nsdp do for j from i to nsdp do k := Search(mspt[i]+mspt[j],mspt): fprintf(fd, "%d 1 %d %d 1\n", k-1, i, j); # Moment variable y_{k-1} at row i and column j od: od: # then the localizing matrices for ig from 1 to ng do gic := gc[ig]: giidx := gidx[ig]: for igc from 1 to nops(gic) do mgi := mspt[giidx[igc]]: #lprint(mgi); gifloat := convert(gic[igc],float,dataprec): if mgi = zero then fprintf(fd, "0 %d 1 1 %s\n", ig+1,convert(-gifloat,string)); # Moment variable y_{j-1} at first row and column j for j from 2 to nsdpg[ig] do fprintf(fd, "%d %d 1 %d %s\n", j-1, ig+1, j, convert(gifloat,string)); # Moment variable y_{j-1} at first row and column j od; for i from 2 to nsdpg[ig] do for j from i to nsdpg[ig] do k := Search(mspt[i]+mspt[j]+mgi,mspt): fprintf(fd, "%d %d %d %d %s\n",k-1, ig+1, i, j, convert(gifloat,string)); # Moment variable y_{k-1} at row i and column j od: od: else for i from 1 to nsdpg[ig] do for j from i to nsdpg[ig] do k := Search(mspt[i]+mspt[j]+mgi,mspt): fprintf(fd, "%d %d %d %d %s\n",k-1, ig+1, i, j, convert(gifloat,string)); # Moment variable y_{k-1} at row i and column j od: od: fi: od: od: ### The case when y0 = 1 because 0 does not belong to the multi-index list # there is a bug in sdpa-gmp when not subtracting the "useless" moment variables with k - sub_idx else for i from 1 to nsdp do for j from i to nsdp do k := Search(mspt[i]+mspt[j],mspt): fprintf(fd, "%d 1 %d %d 1\n", k-sub_idx, i, j); # Moment variable y_{k-sub_idx} at row i and column j od; od; fi: #### Loop again to find the absorbers (this loop can be saved thanks to the previous one) absorb_tbl := table([seq(0,i=1..msdp)]): for i from 1 to nsdp do for j from i to nsdp do m := mspt[i]+mspt[j]: if not_even(m) then k:= Search(m,mspt): absorb_tbl[k]:=i: fi: od: od: #### In this table, for each gamma (appearing at index k in mspt), one indicates all monomials alpha (and implicitely beta) such that gamma = alpha + beta absorb_PP := table([seq([],i=1..msdp)]): for i from 1 to nsdp do k:= Search(2*mspt[i],mspt): absorb_PP[k]:=[op(absorb_PP[k]),i,i]: for j from i+1 to nsdp do m := mspt[i]+mspt[j]: k:= Search(m,mspt): absorb_PP[k]:=[op(absorb_PP[k]),i,j,j,i]: od: od: #lprint(msdp); lprint(mspt); lprint(absorb_PP); #### old version, works only when all monomials belong to the NP # for j from 2 to nsdp do # fprintf(fd, "%d 1 %d %d 1\n", j-1, 1, j); # Moment variable y_{j-1} at first row and column j # od; # for i from 2 to nsdp do # for j from i to nsdp do # k := Search(mspt[i]+mspt[j],mspt): # fprintf(fd, "%d 1 %d %d 1\n", k-1, i, j); # Moment variable y_{k-1} at row i and column j # od; # od; fclose(fd); system("sed -i 's/ \\./ 0\\./g' in.dat-s"); system("sed -i 's/^\\./0\\./g' in.dat-s"); system("sed -i 's/^-\\./-0\\./g' in.dat-s"); system("sed -i 's/ -\\./ -0\\./g' in.dat-s"); system("rm -f out.dat-s"); system("rm -f out.mm"); if not gmp then system("sdpa -ds in.dat-s -o out.dat-s -p param.sdpa > /dev/null"); else write_param(precSDP,epsStar,epsStar); system("sdpa_gmp -ds in.dat-s -o out.dat-s -p my_param_gmp.sdpa > /dev/null"); fi: system("echo $(grep objValPrimal out.dat-s) ';' 'yMat:='$(sed -n '/yMat/,/main/{//!p}' out.dat-s) ';' >> out.mm"); #system("echo 'yMat:=' $(sed -n '/yMat/,/main/{//!p}' out.dat-s) ';' >> out.mm"); system("sed -i 's/ =/ :=/g' out.mm"): system("sed -i 's/{/[/g' out.mm"): system("sed -i 's/}/]/g' out.mm"); system("sed -i 's/] \\[/], \\[/g' out.mm"); read "out.mm": if algo = 1 then lowerbnd := r0 + convert(objValPrimal,rational,exact) : else lowerbnd := 0: fi: printf("lower bound: "); lprint(evalf(lowerbnd)); eigs:=Array([]); eivs := Array([]); nsdponeg:=[nsdp,op(nsdpg)]: for i from 1 to ng+1 do Y := Matrix(yMat[i]): if algo = 1 then Yrat:=Matrix(Y): else # case when algo = 2 if precRound >= 30 then Y := convert(Y,rational,exact): else Y := convert(Y,rational,precRound) fi: Y := (Y + Transpose(Y))/2; if ng = 0 then Yrat := absorber_PP(Y,rList,absorb_PP,nsdp,mspt): else Yrat:=Y: fi: printf("err="); lprint(expand(Transpose(Vector(msk)).Yrat.Vector(msk)+lowerbnd-r)) fi: (eigsi,eivsi) := eigseivs(Yrat,X,msk[1..nsdponeg[i]],precSVD,precRound,algo); eigs:=Concatenate(1,eigs,eigsi): eivs:=Concatenate(1,eivs,eivsi): od: return convert(eigs,list), convert(eivs,list),lowerbnd,absorb_tbl,nsdponeg: end; absorber_PP := proc(Yrat,rList,absorb_PP,nsdp,mspt) local Yproj,i,j,i1,k,m,list_m,nopsm,maxY,minY,yij; # for i from 1 to nsdp do # for j from i to nsdp do # yij := Yrat[i,j]: if abs(yij) < 1/10^3 then Yrat[i,j] := 0: Yrat[j,i]:=0: fi: # od: #od: #Yproj := Matrix(Yrat,shape='symmetric'); Yproj := Matrix(nsdp,nsdp,shape='symmetric'); maxY := max(max(abs(Yproj))): minY := min(min(abs(Yproj))): #lprint(evalf(minY)); lprint(evalf(maxY)); #lprint(Y); #lprint(rList); for i from 1 to nsdp do for j from i to nsdp do m := mspt[i] + mspt[j]; k:= Search(m,mspt): list_m := [op(absorb_PP[k])]: nopsm := nops(list_m)/2: Yproj[i,j] := Yrat[i,j] - 1/nopsm*(add(Yrat[list_m[2*im-1],list_m[2*im]],im=1..nopsm) - rList[k]): od: od: return Yproj: end: write_param := proc(precSDP,epsStar,epsDash) local fd; fd := fopen("my_param_gmp.sdpa",WRITE,TEXT); fprintf(fd,"300 unsigned int maxIteration;\n"); fprintf(fd,"1.0E-%d \t double 0.0 < epsilonStar;\n",epsStar); fprintf(fd,"1.0E5 double 0.0 < lambdaStar;\n"); fprintf(fd,"2.0 double 1.0 < omegaStar;\n"); fprintf(fd,"-1.0E5 double lowerBound;\n"); fprintf(fd,"1.0E5 double upperBound;\n"); fprintf(fd,"0.1 double 0.0 <= betaStar < 1.0;\n"); fprintf(fd,"0.3 double 0.0 <= betaBar < 1.0, betaStar <= betaBar;\n"); fprintf(fd,"0.9 double 0.0 < gammaStar < 1.0;\n"); fprintf(fd,"1.0E-%d \t double 0.0 < epsilonDash;\n",epsDash); fprintf(fd,"%d precision\n",precSDP); fclose(fd); end; checkrational := proc(U) local v: for v in U do: if not type(convert(v,rational),realcons) then lprint(v): error "Non Rational Cholesky factor, retry with gmp = true":fi: od: return: end; eigseivs := proc(Yrat,X,ms,precSVD,precRound,algo) local v, e, msvec, eigs, eivs, U,S,V, Ysvd,SVD,ti,tcmp; SVD := false: msvec := Vector(ms); # if SVD then # Digits:=precSVD; # lprint("starting SVD"); # ti := time(): # (U,S,V) := MTM[svd](Yrat); # tcmp := time() - ti: lprint (tcmp); # lprint("ending SVD"); #else ti := time(): lprint("starting Cholesky"); if algo = 1 then Digits := precSVD:fi: U := LUDecomposition(Yrat,method='Cholesky'); checkrational(U): tcmp := time() - ti: #lprint(U); lprint (tcmp); lprint("ending Cholesky"); S := IdentityMatrix(nops(ms)); # fi: eigs := Diagonal(S); eivs := Transpose(Transpose(msvec).U); Digits:=10; Ysvd := U.S.V^%T; #lprint(Ysvd); #lprint(max(max(Y),-min(Y))); #lprint(max(max(Ysvd),-min(Ysvd))); #lprint(max(Y - Ysvd)); if precRound >= 30 then lprint("exact"); eigs := convert(eigs,rational,exact); eivs := map (_e -> convert(_e,rational,exact), eivs): else eigs := convert(eigs,rational,precRound); eivs := map (_e -> convert(_e,rational,precRound), eivs); fi: return (eigs, eivs): end; neighbours1 := proc(n) local i,nones,nmones,nzeros,nlist,n1; nlist := []; for i from 1 to floor(n/2) do nones:=i; nmones:=i; nzeros:=n-2*i; nlist:=[op(nlist),op(permute([1$nones, 0$nzeros, -1$nmones]))]; od: for i from 1 to ceil(n/2) do nones:=i; nmones:=i-1; nzeros:=n-2*i+1; n1 := permute([1$nones, 0$nzeros, -1\$nmones]); nlist := [op(nlist),op(n1),op(-n1)]: od: return nlist: end; count_n1 := proc(n) add(multinomial(n,i,i,n-2*i),i=1..floor(n/2)) + 2*add(multinomial(n,i,i-1,n-2*i+1),i=1..ceil(n/2)); end; even_n1 := proc(n,mspt) local n1,a,list_n1,cardn1,n1a,b,ab; list_n1 := []; n1 := neighbours1(n); cardn1 := nops(n1); for a in mspt do: n1a := []; for b from 1 to cardn1 do ab := a+n1[b]: if not has(ab,-1) then n1a := [op(n1a),ab]:fi: od: list_n1:=[op(list_n1),n1a]: od: return list_n1: end; # The following degree_diff1 procedure returns wrong results: for instance Y^2 does not belong to neighbours1(X^2) but degree_diff1(X^2,Y^2) = true # degree_diff1 := proc(a,b) # local c := b - a; # if abs(add(i, i in c))<= 1 then return true else return false:fi: #end; not_even := proc(m) return has(1,map(_c -> irem(_c,2), m)): end; ### This procedure sometimes fails to find absorbing polynomials since the resulting decomposition yields monomials whose squares do not necessarily belong to the NP decomp_mon := proc(a,n) local c,b,degb,cm,cnt,i; c := map(_c -> iquo(_c,2), a); b := a - 2*c; degb := ceil(add(i,i in b)/2); cm := [seq(0,i=1..n)]: cnt := 0; for i from 1 to n do: if b[i] = 1 then cm[i] := b[i]: cnt := cnt + 1: fi: if cnt >= degb then break: fi: od: return c,cm,b - cm; end; truncate_withNP := proc(f,X,mspt,msptk) local NP,msptkNP,msptNP,m,i; NP := newtonpolytope(f, X): msptkNP := []; for m in msptk do if contains(NP,2*m) then msptkNP := [op(msptkNP),m] : fi: od: msptNP := [op(msptkNP)]: for i from nops(msptk) + 1 to nops(mspt) do m := mspt[i]: if contains(NP,m) then msptNP := [op(msptNP),m] : fi: od: return msptNP,msptkNP: end: testNP := proc(f,X) local d,k,n,card_nk,mspt,msptk; d := degree(f): k := d/2: n := nops(X): card_nk := binomial(n+k,k): mspt := monspt(n,d): msptk := mspt[1..card_nk]: return truncate_withNP(f,X,mspt,msptk): end: absorber := proc(u,X,e,even_mons,ms,mspt,absorb_tbl) local i,j,k,ucoeffs,uidx,uc,err_list,err,m,bad_m,n,c,cm,cp,cfs,sqs,n1,m1,m2,k1,k2; #printf("\nu = "); lprint(evalf(u)); ucoeffs,uidx := support(u,X,ms); n := nops(X): #printf("\neven_mons = "); lprint(even_mons); err_list := Array([seq(e,i=1..nops(even_mons))]); cfs := []; sqs := []; for i from 1 to nops(ucoeffs) do uc := ucoeffs[i]; bad_m := mspt[uidx[i]]: if not_even(bad_m) then k := Search(bad_m,mspt):m1:=mspt[absorb_tbl[k]]:m2:=bad_m - m1; #printf("\n"); lprint(bad_m); lprint(m1); lprint(m2); cfs := [op(cfs),abs(uc)/2]; sqs := [op(sqs),mul(X[i]^m1[i],i=1..n)+sign(uc)*mul(X[i]^m2[i],i=1..n)]: k1:=Search(2*m1,even_mons): k2:=Search(2*m2,even_mons): err_list[k1] := err_list[k1]-1/2*abs(uc): err_list[k2] := err_list[k2]-1/2*abs(uc): else k := Search(bad_m,even_mons): err_list[k]:=err_list[k]+uc: fi: od: err_list := convert(err_list,list); # printf("\nerr list = "); lprint(evalf(err_list)); cfs := [op(cfs),op(err_list)]; sqs := [op(sqs),seq(mul(X[i]^(m[i]/2),i=1..n), m in even_mons)]; return cfs,sqs: end; old_absorber := proc(u,X,e,even_mons,ms,mspt,even_mons_n1,absorb_tbl) local i,j,k,ucoeffs,uidx,uc,err_list,err,m,bad_m,n,c,cm,cp,cfs,sqs,n1,m1,m2; printf("\nu = "); lprint(evalf(u)); ucoeffs,uidx := support(u,X,ms); n := nops(X): #bm := [seq(0,i=1..n)]: bp := bm: #for i from 1 to ceil(n/2) do bm[i]:=1:od:for i from ceil(n/2)+1 to n do bp[i]:=1:od: #card_absorbed := []; # counts how many monomials to absorbe for each even monomial with support in even_mons # printf("\neven_mons = "); lprint(even_mons); err_list := []; for j from 1 to nops(even_mons) do m := even_mons[j]; #lprint(m); n1 := even_mons_n1[j]; #lprint(n1); err := e; for i from 1 to nops(ucoeffs) do # printf("0 "); lprint(bad_m); lprint(err); # printf("1 "): lprint(bad_m); lprint(err); uc := ucoeffs[i]; bad_m := mspt[uidx[i]]: if m = bad_m then err := err + uc: else if has(bad_m,n1) then err := err - 1/2*abs(uc): fi:fi: od: err_list := [op(err_list),err]; od: printf("\nerr list = "); lprint(evalf(err_list)); cfs := [op(err_list)]; sqs := [seq(mul(X[i]^(m[i]/2),i=1..n), m in even_mons)]; for i from 1 to nops(ucoeffs) do uc := ucoeffs[i]; bad_m := mspt[uidx[i]]: if not_even(bad_m) then #c,cm,cp := decomp_mon(bad_m,n); #printf("\n"); lprint(bad_m); #lprint(c+cm); lprint(c+cp); #m1 := c+cm: m2 := c + cp: k := Search(bad_m,mspt):m1:=mspt[absorb_tbl[k]]:m2:=bad_m - m1; printf("\n"); lprint(bad_m); lprint(m1); lprint(m2); cfs := [op(cfs),abs(uc)/2]; sqs := [op(sqs),mul(X[i]^m1[i],i=1..n)+sign(uc)*mul(X[i]^m2[i],i=1..n)]: fi: od: return cfs,sqs: end; relaxordermin := proc(g) return ceil(max(seq(degree(gi),gi in g))/2): end: locmatsizes := proc(g,k,n) return [seq(binomial(n + k - ceil(degree(gi)/2),n),gi in g)]: end; perturbate := proc(p,prec) local X,d,k,n,mspt,card_nk,msptk,ms,msk; X:=[op(indets(p))]: d := degree(p): k := ceil(d/2); n := nops(X): mspt := monspt(n,d): card_nk := binomial(n+k,n); msptk := mspt[1..card_nk]; mspt,msptk := truncate_withNP(p,X,mspt,msptk): ms := mons(n,d,X,mspt): msk := mons(n,k,X,msptk): return p + 2^(-prec)*dense_perturbation(msk): end: multivsos1:=proc(f,prec::integer := 10, precSVD::integer := 10, precSDP::integer := 200, epsStar::integer := 30, precRound::integer:=30, dataprec::integer:=10, gmp::boolean:=false, algo::integer:=1, glist::list:=[],relaxorder::integer:=0) local p,d,mspt,ms,rc,ridx,ng,gc,gidx,gic,giidx,S,s,c, q,n,k,t,e,r,l,a,s1,s2,u,v,i,j,sqs,cfs,sos,rfloat, eigs, eigsg, eigsgi, soslist,soslistg, soslistgi, sumsos,cnd,maxq,obj_plus_r0, card_nk, even_mons,err_list,err,msptk,msk,absorb_tbl,rmin,nsdponeg,idx,oneg,idxi,g,X,ti,tf,gi,cg; ti := time[real](): c := max(map(_c -> abs(_c),coeffs(expand(f)))); p := 1/c*f; X := [op(indets(p))]: g := []; ng := nops(glist); for i from 1 to ng do gi := expand(glist[i]); cg := max(map(_c -> abs(_c),coeffs(expand(gi)))): g := [op(g),gi/cg] od: #g := expand(glist); if lcoeff(p) < 0 and ng = 0 then lprint(p): error "There is no decomposition into sum of squares for this polynomial"; fi; # Let d = 2 k be the degree; fail if this is odd (can't be pos def) # If it's a constant (i.e. the polynomial is a constant multiple of # a perfect square) then return almost immediately. if ng = 0 then d := degree(p): else rmin := relaxordermin([p,op(g)]): d := 2*max(relaxorder,rmin): fi: k := ceil(d/2); if (2 * k <> d and ng = 0) then lprint(p): error "There is no decomposition into sum of squares for this polynomial"; fi; if (d = 0 and ng = 0) then lprint(p, " * (", s, ")^2"); fi; n := nops(X): mspt := monspt(n,d): card_nk := binomial(n+k,n); msptk := mspt[1..card_nk]; if ng = 0 then mspt,msptk := truncate_withNP(p,X,mspt,msptk):fi: # if no constraints then compute NP ms := mons(n,d,X,mspt): msk := mons(n,k,X,msptk): # let r = p - e * t be a safe perturbation if algo = 1 then t := dense_perturbation(msk); # printf("\nperturber = "); lprint(t): even_mons := 2*msptk; # even_mons_n1 := even_n1(n,even_mons); e := 1/2^prec: else e := 0: t := 0: fi: r := expand(p - e*t); rc,ridx := support(r,X,ms): gc := []: gidx := []: for i from 1 to ng do gic,giidx := support(g[i],X,ms): gc:=[op(gc),gic]; gidx:=[op(gidx),giidx]; od: (eigs,soslist,obj_plus_r0,absorb_tbl,nsdponeg) := sos2sdp(r,X,k,mspt,msk,rc,ridx,precSVD,precSDP,epsStar,precRound,dataprec,gmp,algo,g,gc,gidx); idx := 0: sumsos := obj_plus_r0: oneg := [1,op(g)]; for i from 1 to ng+1 do idxi := idx + nsdponeg[i]: sumsos := sumsos + oneg[i]*sum(eigs[j]*soslist[j]^2,j=idx+1..idxi); idx := idxi: od: # sumsos := obj_plus_r0 + sum(eigs[j]*soslist[j]^2,j=1..nops(soslist)); # for i from 1 to ng do # eigsgi := eigsg[i]: soslistgi := soslistg[i]: # sumsos := sumsos + sum(eigsgi[j]*soslistgi[j]^2,j=1..nops(soslistgi))*g[i]; # od: #lprint(evalf(expand(sumsos + obj_plus_r0))); if algo = 1 then ## if ng > 0 then ## t := dense_perturbation(msk); ## printf("\nperturber = "); lprint(t): ## even_mons := 2*msptk;e := 1/2^prec: r := expand(p - e*t); ## fi: u := expand(r - sumsos+obj_plus_r0); lprint(soslist); lprint(u); cfs,sqs := absorber(u,X,e,even_mons,ms,mspt,absorb_tbl); nsdponeg[1]:=nsdponeg[1]+nops(sqs): else u:=0; cfs := []: sqs := []: fi: err := expand(u+e*t - add(cfs[i]*sqs[i]^2,i=1..nops(sqs))); if algo = 1 then printf("\nerr = "); lprint(evalf(err)); fi: cfs := [op(cfs), op(eigs)]: sqs := [op(sqs),op(soslist)]: sos := []; for i from 1 to nops(sqs) do sos := [op(sos), c*cfs[i],sqs[i]]: od; if soscheck1(f,sos,nsdponeg,g) = 0 then printf("\n Certified SOS Decomposition\n"): fi: tf := time[real]()-ti; printf("bitsize= %d\n",BitSizePolSeq(sos,X)); printf("time= %esecs\n",tf); return sos,nsdponeg: end; soscheck1:=proc(f, sos, nsdponeg, g::list:=[]) local s,i,j,idx,idxi,oneg; oneg := [1,op(g)]: s := 0; idx:=0; for i from 1 to nops(g)+1 do idxi:=idx+nsdponeg[i]: for j from idx+1 to idxi do if sos[2*j-1] < 0 then lprint(evalf(sos[2*j-1])): error "Negative number => invalid sum of squares decomposition"; else s := s + oneg[i]*sos[2*j-1]*sos[2*j]^2 fi: od: idx:=idxi: od: if not expand(f - s) = 0 then lprint(evalf(expand(f - s))); error "Inexact sum of squares decomposition"; else return 0: fi: end; BitSizeSos := proc (sos,X) return add(BitSizePol(p[1], X) + BitSizePolQuadr(p[2], X), p in sos); end; BitSizePolQuadr := proc(q,X) return BitRat(q[1]) + BitSizePol(q[2], X) + BitRat(q[3]); end; BitSizePolSeq3 := proc(listpol,X) return add(BitSizePol(p[1], X) + BitSizePolSeq(p[2], X), p in listpol); end; BitSizePolSeq := proc(listpol,X) return add(BitSizePol(p, X), p in listpol); end; BitSizePol:=proc(p, X) local res; res := [coeffs(expand(p),X)]; return BitSizeSeq(res); end; BitSizeSeq:=proc(l) return add(BitRat(c), c in l); end; BitRat := proc(r) local n, d, res,rs; if type(r,rational) then rs :=r : else rs:=r^2: fi: if rs = 0 then return 1; fi; (n, d) := (abs(numer(rs)), abs(denom(rs))); if d = 1 then res := ilog2(n) + 1 else res := ilog2(n) + ilog2(d) + 2 fi; return res; end; benchRAGLib := proc(f,g::list:=[]) local sys,i,ti,tf,sol; ti := time[real](): sys := [expand(f) < 0]; for i from 1 to nops(g) do sys := [op(sys), expand(-g[i]) < 0] od: sol := HasRealSolutions(sys); tf := time[real]()-ti; printf("time= %esecs\n",tf); return sol; end; benchSamplePoints := proc(f,g::list:=[]) local sys,vars,R,P,i,ti,tf; ti := time[real](): sys := [f < 0]; for i from 1 to nops(g) do sys := [op(sys), g[i] >= 0] od: vars:=[op(indets(f))]: R:=PolynomialRing(vars): P:=SamplePoints(sys,R); tf := time[real]()-ti; printf("time= %esecs\n",tf); return P; end;