share/contrib/sarag/intervalArithmetic.mac

   1 /* ------------------------------------------------------------------- */
   2 /* SARAG - Interval Arithmetic                                         */
   3 /* by Fabrizio Caruso                                                  */
   4
   5
   6 intervalProduct(lhs,rhs) :=
   7   [min(lhs[1]*rhs[1],lhs[1]*rhs[2],
   8        lhs[2]*rhs[1],lhs[2]*rhs[2]),
   9    max(lhs[1]*rhs[1],lhs[1]*rhs[2],
  10        lhs[2]*rhs[1],lhs[2]*rhs[2])];
  11
  12
  13 squaredIntervalAt(interval) :=
  14   squaredIntervalAtAux(interval[1],interval[2]);
  15
  16 squaredIntervalAtAux(leftEnd,rightEnd) :=
  17   if leftEnd*rightEnd<1 then
  18     [0,max(leftEnd^2,rightEnd^2)]
  19   else
  20     [min(leftEnd^2,rightEnd^2),max(leftEnd^2,rightEnd^2)];
  21
  22
  23 powerIntervalAt(interval,p) :=
  24   if oddp(p) then
  25      [interval[1]^p,interval[2]^p]
  26   else
  27      if interval[1]*interval[2]<1 then
  28          [0,max(interval[1],interval[2])^p]
  29      else
  30         if interval[1]>0 then
  31           [interval[1]^p,interval[2]^p]
  32         else
  33           [interval[2]^p,interval[1]^p];
  34
  35
  36 constantByInterval(c,interval) :=
  37   if c = 0 then
  38     [0]
  39   else
  40     if c>0 then
  41       [c*interval[1],c*interval[2]]
  42     else
  43       [c*interval[2],c*interval[1]];
  44
  45 constantPlusInterval(d,interval) :=
  46   [interval[1]+d,interval[2]+d];
  47
  48 evaluatePolAt(pol,var,interval) :=
  49   if degree(pol,var) < 1 then
  50     [pol]
  51   else
  52     evaluatePolAtAux(rest(reverse(poly2list(pol,var)),1),
  53                      constantByInterval(coeff(pol,var,degree(pol,var)),interval));
  54
  55 evaluatePolAtAux(revCoeffSeq,interval) :=
  56   if revCoeffSeq=[] then
  57     interval
  58   else
  59     evaluatePolAtAux(rest(revCoeffSeq,1),
  60                      constantPlusInterval(first(revCoeffSeq),
  61                                           interval)
  62                     );
  63
  64 0;