ABSOLUTEINFAUX = 0; SHOWPOSITIVITY = 1; FINDPOLYWITHGIVENERROR = 2; SUPNORMABSOLUTE = 3; SUPNORMRELATIVENOPOLE = 4; /* Vectors of verbosity. If verb[SOMEPROCEDURE] > 0, it displays information */ verb = [|0, 0, 0, 0, 0|]; /* Level of recursion in the procedure that finds F. */ recLevelInAbsoluteInf = 0; /* List of removable discontinuities of f (provided by the user) */ removableDiscontinuitiesOfF = [||]; /****************************************************************************/ /* Procedure absoluteInf: ====================== Return a value F such that forall x in I, |f(x)|>=F The procedure may return NaN, in case of failure */ procedure absoluteInfAux(f, I, recLevel) { var J, res, F1, F2; if (verb[ABSOLUTEINFAUX]>0) then { print("absoluteInfAux called:"); if (verb[ABSOLUTEINFAUX]>1) then { print("f =", f); print("I =", I); }; print("recLevel =", recLevel); }; J = evaluate(f, I); if (inf(J)*sup(J) > 0) then res = inf(abs(J)) else { if (recLevel == 0) then res = @NaN@ else { F1 = absoluteInfAux(f, [inf(I), mid(I)], recLevel-1); F2 = absoluteInfAux(f, [mid(I), sup(I)], recLevel-1); if ( (!(F1 == F1)) || (!(F2 == F2)) ) then res = @NaN@ else { if (F1 < F2) then res = F1 else res = F2; }; }; }; if ((res==res) && (verb[ABSOLUTEINFAUX]>0)) then { if (verb[ABSOLUTEINFAUX]==0) then print("absoluteInfAux found a suitable value") else print("absoluteInfAux found a suitable value (F =", round(res, 15, RN), ")"); }; return res; }; procedure absoluteInf(f, I) { return absoluteInfAux(f, I, recLevelInAbsoluteInf); }; /****************************************************************************/ /* Procedure computeLowerBound: ============================ Given a polynomial p, a function f, an interval I, a real number gamma and a mode mode in {absolute, relative} compute a real number ell such that ||eps||_infty^I >= ell where eps is the function defined by * eps = p - f if mode is absolute * eps = p/f - 1 if mode is relative Additionally, make sure that if things are pretty that ell <= ||eps|| <= ell*(1+gamma) */ procedure computeLowerBound(p,f,I,gamma,mode) { var res; var g, gPrime; var zeros; var oldPrec; var oldPoints; var z, y; // Construct approximation error function, differentiate it // (resp. take the numerator of its derivative), find the zeros // of the derivative and do an interval evaluation of these // zeros plus the interval endpoints. The lower bounds // of the image bounding intervals are near the real // supremum norm but certainly below it. if (mode == absolute) then { g = p - f; gPrime = diff(p) - diff(f); } else { g = p/f-1; gPrime = diff(p) * f - p * diff(f); }; oldPrec = prec; oldPoints = points; prec = 10 + (ceil(-log2(abs(gamma))/2))!; points = 3*degree(p)!; zeros = dirtyfindzeros(gPrime,I); zeros = [|inf(I),sup(I)|] @ zeros; prec = 10 + ceil(-log2(abs(gamma)))!; res = 0; for z in zeros do { y = abs(evaluate(g,[z])); if (inf(y) > res) then res = inf(y); }; prec = oldPrec!; points = oldPoints!; return res; }; /****************************************************************************/ /* Procedure showPositivity ======================== Given a polynomial p, an interval I and a boolean silent, compute a boolean res such that res is true implies forall x in I p(x) > 0. Additionally, make sure that if things are pretty, forall x in I p(x) > 0 implies res is true */ procedure showPositivity(p,I) { var res; var zeros; if (verb[SHOWPOSITIVITY]>2) then print("Showing positivity of:",p); /* A polynomial is positive if it is positive somewhere and does not cross x-axis. */ res = (numberroots(p,I) == 0) && (inf(evaluate(p,[mid(I)])) > 0); if (!res) then { if (verb[SHOWPOSITIVITY]==1) then print("Warning: could not show positivity") else { if (verb[SHOWPOSITIVITY]>1) then print("Warning: could not show positivity of", p,"over",I); }; }; return res; }; /****************************************************************************/ /* Procedure upperBoundPolyCoeffRests ================================== Given a list l of intervals l = [|c_0,...,c_n|] and an interval I, compute an interval bound such that for all x in I, forall a_i in c_i, sum a_i * x^i in bound */ procedure upperBoundPolyCoeffRests(l,I) { var bound; var myList; /* Details: just a simple Horner in Interval Arithmetic */ if (l == [||]) then { bound = [0]; } else { myList = revert(l); bound = head(myList); myList = tail(myList); while (myList != [||]) do { bound = bound * I + head(myList); myList = tail(myList); }; }; return bound; }; /* procedure findPolyWithGivenError ================================ Computes a polynomial T such that ||T-f|| < errMax over I. T is computed using Taylor Models. discontinuities is a list that contains at most one element. if discontinuities is empty, T is computed with mode=absolute around mid(I) if singulariy==[|x0|], then T is computed with mode=relative around the point x0. The good degree for T is first searched by bisection, neglecting roundoff errors. Once the good degree has been found, all roundoff errors are taken into account for computing effectively a polynomial T with floating-point coefficients with a safe bound. */ procedure findPolyWithGivenError(f,I,errMax, discontinuities) { var mode, x0, n, done; var TM, TMmax, DeltaLagrange, DeltaCoeffs; var T, TRounded, TDiff, DeltaRound, Delta; var nmin, nmax, p, oldRationalMode; if (discontinuities==[||]) then { mode=absolute; x0 = mid(I); } else { /* discontinuities is supposed to contain at most one point */ x0 = head(discontinuities); mode = relative; }; if (verb[FINDPOLYWITHGIVENERROR] > 1) then { if (verb[FINDPOLYWITHGIVENERROR] > 2) then { print("findPolyWithGivenError: f =", f, "over", I); }; print("findPolyWithGivenError: errMax = ",round(errMax,15,RN)); print("findPolyWithGivenError: mode =",mode,"and x0 = ",x0); }; n = 1; done = false; while (!done) do { TM = taylorform(f,n,x0,I,mode); if (mode == absolute) then { DeltaLagrange = TM[2]; } else { DeltaLagrange = evaluate(x^n,I-x0) * TM[2]; }; if (verb[FINDPOLYWITHGIVENERROR] > 1) then { print("findPolyWithGivenError: trying degree ",n, "(current bound:", round(sup(abs(DeltaLagrange)),15,RN),")"); }; if (sup(abs(DeltaLagrange))>=errMax) then n = n*2 else done=true; }; nmin = n/2; nmax = n; /* TMmax and DeltaLagrangeMax are used to remember the TM and DeltaLagrange corresponding to nmax */ TMmax = TM; DeltaLagrangeMax = DeltaLagrange; while (nmax-nmin > 1) do { n = floor((nmin+nmax)/2); TM = taylorform(f,n,x0,I,mode); if (mode == absolute) then { DeltaLagrange = TM[2]; } else { DeltaLagrange = evaluate(x^(n),I-x0) * TM[2]; }; if (verb[FINDPOLYWITHGIVENERROR] > 1) then { print("findPolyWithGivenError: trying degree ",n, "(current bound:", round(sup(abs(DeltaLagrange)),15,RN)); }; if (sup(abs(DeltaLagrange))>=errMax) then nmin = n else { nmax = n; TMmax = TM; DeltaLagrangeMax = DeltaLagrange; }; }; n = nmax; if (verb[FINDPOLYWITHGIVENERROR] > 0) then { print("\nfindPolyWithGivenError: deg(T) =", n); }; TM = TMmax; DeltaLagrange = DeltaLagrangeMax; p = TM[0]; /* Bound the error due to the coefficient rests. This means evaluating the family of the polynomials having the rests as coefficients over the interval. */ DeltaCoeffs = upperBoundPolyCoeffRests(TM[1],I-x0); /* Now shift the Taylorform so that T is a polynomial approximating f */ oldRationalMode = rationalmode; rationalmode = on!; T = horner(simplifysafe(p(x - x0))); rationalmode = oldRationalMode!; /* Now round the coefficients of T into small floating-point numbers. This must be done after the translation. */ TRounded = roundcoefficients(T,[| prec... |]); TDiff = horner(T - TRounded); /* Bound the extra rounding error due to the rounding of the coefficient. This means bounding (evaluating in Interval Arithmetic) of the round-off polynomial. */ DeltaRound = evaluate(TDiff,I); /* Add up all 3 sources of error */ Delta = DeltaCoeffs + DeltaLagrange + DeltaRound; /* This test is very unlikely to fail. If we want to do really clean things, we should however think about what to do in this case (increase n? increase prec?) */ if (sup(abs(Delta)) > errMax) then print("ERROR in findPolyWithGivenError: roundoff errors made it fail"); return TRounded; }; /* The proof of this procedure can be found in the article. */ procedure supnormAbsolute(p, f, I, quality, discontinuities) { var ell, errMax, T, bound, s1, s2; var timeL, timeT, timePositivity; var oldRoundingwarnings, checked, sosInstances, res; oldRoundingwarnings = roundingwarnings; roundingwarnings = off!; if (verb[SUPNORMABSOLUTE]>0) then print("supnormAbsolute called with discontinuities =", discontinuities); timeL = time(ell = computeLowerBound(p, f, I, quality/32, absolute)); /* Note that errMax is surely smaller than ell*quality*(15/32) */ errMax = inf([ell*quality*(15/32)]); timeT = time(T = findPolyWithGivenError(f, I, errMax, discontinuities)); /* Same trick: bound <= ell*(1+quality/2) We actually prove ||p-T|| <= bound which proves that ||p-T|| <= ell*(1+quality/2) */ bound = inf([ell*(1+quality/2)]); s1 = bound-(p-T); s2 = bound-(T-p); timePositivity = time(checked = (showPositivity(s1,I) && showPositivity(s2,I))); if(!checked) then { print("ERROR: the computed lower bound for ||p-f|| does not achieve the required quality"); res = @NaN@; } else res = [ell, ell*(1+31*quality/32)]; sosInstances = [| [| s1, I |], [| s2, I |] |]; if (verb[SUPNORMABSOLUTE]>0) then { print(""); print("Result of timings: "); print("==================="); print("Dirty lower bound : ", floor(timeL*1000), "ms"); print("Polynomial T : ", floor(timeT*1000), "ms"); print("Showing positivity: ", floor(timePositivity*1000), "ms"); print(""); }; roundingwarnings = oldRoundingwarnings!; return [| res, sosInstances |]; }; /* The proof of this procedure can be found in the article. */ procedure supnormRelativeNoPole(p,f,I,quality, discontinuities) { var F, ell, u, errMax, T, signT, bound, s1, s2; var timeL, timeT, timeF, timePositivity; var oldRoundingwarnings, checked, sosInstances, res; oldRoundingwarnings = roundingwarnings; roundingwarnings = off!; if (verb[SUPNORMRELATIVENOPOLE]>0) then print("supnormRelativeNoPole called with discontinuities =", discontinuities); timeF = time(F = absoluteInf(f, I)); if (!(F==F)) then { print("ERROR: failed to prove that",f,"does not vanish on",I); res = @NaN@; sosInstances = [||]; } else { timeL = time(ell = computeLowerBound(p, f, I, quality/32, relative)); /* Note that errMax is surely smaller than ell*quality*(15/32) */ u = sup(ell*(1+31*quality/32)); errMax = inf( [ell * quality * (15/32) * 1/(1+u) * F/(1+15*quality/32)] ); timeT = time(T = findPolyWithGivenError(f, I, errMax, discontinuities)); if (T(mid(I)) > 0) then signT = 1 else signT = -1; /* Same trick: bound <= ell*(1+quality/2) We actually prove ||p/T-1|| <= bound which proves that ||p/T-1|| <= ell*(1+quality/2) */ bound = inf([ell*(1+quality/2)]); s1 = bound*T*signT-(p-T); s2 = bound*T*signT-(T-p); timePositivity = time(checked = (showPositivity(s1,I) && showPositivity(s2,I))); if(!checked) then { print("ERROR: the computed lower bound for ||p/f-1|| does not achieve the required quality"); res = @NaN@; } else res = [ell, u]; sosInstances = [| [| s1, I |], [| s2, I |] |]; if (verb[SUPNORMRELATIVENOPOLE]>0) then { print(""); print("Result of timings: "); print("==================="); print("Bound F : ", floor(timeF*1000), "ms"); print("Dirty lower bound : ", floor(timeL*1000), "ms"); print("Polynomial T : ", floor(timeT*1000), "ms"); print("Showing positivity: ", floor(timePositivity*1000), "ms"); }; }; roundingwarnings = oldRoundingwarnings!; return [| res, sosInstances |]; }; /* Procedure determineOrderOfZero ============================== Given a function f, a real number x0 and an integer n, compute an integer k such that k is an upper bound for the order of the zero of f in x0, provided that the zero is of order smaller than n. If x0 is a zero of order >= n, the procedure fails with @NaN. */ procedure determineOrderOfZero(f,x0,n) { var k, TM, T, coeffBounds, test; TM = taylorform(f,n,x0,relative); T = TM[0]; coeffBounds = TM[1]; k = 0; test = true; while ( (k <= n-1) && test ) do { if ( (coeff(T,k)==0) && coeffBounds[k] == [0] ) then k = k + 1 else test = false; }; if (test) then k = @NaN@; return k; }; /* Procedure findpoles =================== Given a polynomial p, a function f and an interval I, computes a list res of real numbers x_i such that if things are pretty, forall x_i in res, either f or p vanishes in or near x_i. */ procedure findpoles(p,f,I) { var res, oldPoints; oldPoints = points; points = 20!; /* We do not expect more than a few discontinuities */ /* They should be all found by sampling with 20 points */ res = dirtyfindzeros(f,I); points = oldPoints!; return res; }; /* The wrapper =========== In the Sollya language, we cannot detect the removable discontinuities *inside* the expression of f. In order to allow testing with such functions, we use the following workaround: the user may supply a list of removable discontinuities of f in the list removableDiscontinuitiesOfF. TODO: this code currently does not manage the case when there are several removable discontinuities. */ procedure supnorm(p,f,I,quality,mode) { var res, poles, x0, k, pTilde, fTilde, singularities, a; singularities = removableDiscontinuitiesOfF; res = [||]; if (mode==absolute) then res = supnormAbsolute(p,f,I,quality, singularities) else { poles = findpoles(p,f,I); if (poles==[||]) then res = supnormRelativeNoPole(p,f,I,quality, singularities) else { if ((length(poles)==1) && (singularities == [||])) then { x0 = head(poles); k = determineOrderOfZero(f, x0, degree(p)); if (!(k == k)) then { print("ERROR: unable to determine the order of ",x0,"as a zero of",f); res = error; } else { /* The next terrible line is a Sollya trick to divide the polynomial p by (x - x0)^k Sollya actually can divide by x^k but not by (x - x0)^k */ pTilde = horner((horner(p(x + x0))/(x^k))(x - x0)); /* Check if the division worked */ if (degree(pTilde) == degree(p) - k) then { fTilde = f/((x - x0)^k); res = supnormRelativeNoPole(pTilde,fTilde,I,quality, [|x0|]); }; }; } else print("ERROR: the function has several poles. We do not currently manage this case"); }; }; return res; };