/* This is a Sollya wrapper for checking the positivity of a polynomial */
/* by a sum-of-squares technique. */
/* Note: the GP tool (http://pari.math.u-bordeaux.fr/) is required. */
/* This is the script (written in GP) that performs the SOS decomposition */
gpSOSScript = "unisos.gp";
/* This is the file where unisos.gp expects to find the data of the problem */
gpInputFile = "sosInput.gp";
/* This is the name of the file where the GP scripts outputs its result. */
/* This result is simply a Sollya file defining a polynomial q. */
gpOutputFile = "sosOutput.sollya";
holSOSPath = "./holsos/";
holSOSBasename = "sosinstance";
/* sosInstances is a list L=[|l1,...,ln|], where the elements li are */
/* of the form li=[|p,I|] such that one wants to prove that p>0 over I. */
/* */
/* This procedure returns [|okay, L|], where okay is true if and only */
/* if all the instances have been proved positive, and L is a list of */
/* instances li=[|q, I|] where q is written as a sum of squares. */
procedure sos(sosInstances) {
var okay, res, i;
var sosInst, mySosInstances, returnInstances, myHolSOSPath;
var p, I, a, b, diffPQ;
var oldDisplay, oldAutosimplify, oldRationalmode, oldFullParentheses;
var sosInstanceName;
myHolSOSPath = holSOSPath;
if (myHolSOSPath[length(myHolSOSPath)-1] != "/")
then myHolSOSPath = myHolSOSPath @ "/";
bashexecute("if [ ! -d " @ holSOSPath @ " ]; then mkdir " @ holSOSPath @ ";fi");
bashexecute("rm -f " @ holSOSPath @ "*");
i = 1;
oldDisplay = display; display = powers!;
okay = true;
mySosInstances = sosInstances;
returnInstances = [||];
while (okay && (mySosInstances != [||])) do {
sosInst = head(mySosInstances);
mySosInstances = tail(mySosInstances);
p = simplifysafe(sosInst[0]);
I = sosInst[1];
a = inf(I);
b = sup(I);
write("p = ", p, ";\n") > gpInputFile;
write("a = ", a, ";\n") >> gpInputFile;
write("b = ", b, ";\n") >> gpInputFile;
bashexecute("gp -q " @ gpSOSScript @ " | head -n -1 > " @ gpOutputFile);
oldAutosimplify = autosimplify;
autosimplify = off!;
oldRationalmode = rationalmode;
rationalmode = off!;
q = error;
execute(gpOutputFile);
autosimplify = oldAutosimplify!;
rationalmode = oldRationalmode!;
if (q == q) then {
oldAutosimplify = autosimplify;
autosimplify = on!;
oldRationalmode = rationalmode;
rationalmode = on!;
diffPQ = simplifysafe(horner(simplifysafe(horner(q) - horner(p))));
rationalmode = oldRationalmode!;
autosimplify = oldAutosimplify!;
if (diffPQ != 0) then {
write("Warning: SOS returned a different polynomial for the following instance\n");
write("p = ", p, ";\n");
write("a = ", a, ";\n");
write("b = ", b, ";\n");
write("q = ", q, ";\n");
okay = false;
} else {
oldAutosimplify = autosimplify;
autosimplify = off!;
oldRationalmode = rationalmode;
rationalmode = off!;
returnInstances = returnInstances :. [| q, [a, b] |];
display = default!;
sosInstanceName = myHolSOSPath@holSOSBasename @ i;
display = powers!;
oldFullParentheses = fullparentheses;
fullparentheses = on!;
write("p = ",p,";\n") > sosInstanceName;
write("q = ",q,";\n") >> sosInstanceName;
write("a = ",a,";\n") >> sosInstanceName;
write("b = ",b,";\n") >> sosInstanceName;
fullparentheses = oldFullParentheses!;
i = i + 1;
autosimplify = oldAutosimplify!;
rationalmode = oldRationalmode!;
};
} else {
write("Warning: SOS did not work for the following instance\n");
write("p = ", p, ";\n");
write("a = ", a, ";\n");
write("b = ", b, ";\n");
okay = false;
};
};
display = oldDisplay!;
oldAutosimplify = autosimplify;
autosimplify = off!;
oldRationalmode = rationalmode;
rationalmode = off!;
res = [|okay, returnInstances |];
autosimplify = oldAutosimplify!;
rationalmode = oldRationalmode!;
bashexecute("rm -rf " @ gpInputFile @ " " @ gpOutputFile);
return res;
};