share/contrib/sarag/certificateOfPositivity.mac

   1
   2 /* -------------------------------------- */
   3
   4 isListCertLess(lhs,rhs) :=
   5   isListLess(first(lhs),first(rhs));
   6
   7
   8 /* Order Isolating with certificates List */
   9 isListCertSort(isListCert) :=
  10   sort(isListCert,isListCertLess);
  11
  12
  13 subInterval(checkSub,bigInt) :=
  14   is(first(checkSub)>=first(bigInt) and
  15      second(checkSub)<=second(bigInt));
  16
  17 subIntervalInSet(intSet,bigInt) :=
  18   block([found,i],
  19   found : false,
  20   if intSet=[] then
  21     return(false),
  22   i : 1,
  23
  24   while not(found) and i<= length(intSet) do
  25     (
  26     if subInterval(part(intSet,i),bigInt) then
  27       found : true
  28     else
  29       i : i + 1
  30     ),
  31   return(found)
  32
  33   );
  34
  35
  36
  37 /* Yun's square free decomposition algorithm */
  38 squareFreeFactor(pol,var) :=
  39   block([g,res,c],
  40   res : [],
  41
  42   g : gcd(pol,diff(pol,var),var),
  43   q : g,
  44
  45   c : ratexpand(pol/q),
  46   d : ratexpand(diff(pol,var)/g-diff(c,var)),
  47   while not(degree(c,var)=0) do
  48     (
  49     q : gcd(c,d),
  50     res : endcons(q,res),
  51     c : ratexpand(c/q),
  52     d : ratexpand(d/q-diff(c,var))
  53     ),
  54   return(res)
  55   );
  56
  57
  58
  59 bernsteinSplitWithZ(bcoeff,l,r,m) :=
  60   block([
  61         p : length(bcoeff)-1,
  62         bern_first,
  63         bern_second,
  64         reduced_l,
  65         reduced_r,
  66         b : make_array('any),
  67         i,j,nnorm,dnorm,ratnorm],
  68
  69
  70   ratnorm : (r-m)/(r-l),
  71   nnorm : num(ratnorm),
  72   dnorm : denom(ratnorm),
  73
  74   b : make_array('any,p+1),
  75   for i : 0 thru p do
  76     b[i] : bcoeff[i+1],
  77
  78   bern_first : [dnorm^p * b[0]],
  79   bern_second : [dnorm^p * b[p]],
  80
  81   for i : 1 thru p do
  82     (
  83
  84     for j: 0 thru p-i do
  85        b[j] : nnorm*b[j] + (dnorm-nnorm)*b[j+1],
  86
  87     bern_first : endcons(dnorm^(p-i) * b[0],bern_first),
  88
  89     bern_second : cons(dnorm^(p-i)*b[p-i],bern_second)
  90
  91     ),
  92
  93   return([bern_first,bern_second])
  94   );
  95
  96
  97
  98
  99 bernsteinMerge(lhs,rhs):=
 100   block([extB,normfact,lhsNorm],
 101
 102   extB:bernsteinSplitWithZ(lhs[3],lhs[1][1],lhs[1][2],rhs[1][2])[1],
 103
 104   degP: length(lhs[3])-1,
 105   leftEnd: lhs[1][1],
 106   rightEnd: lhs[1][2],
 107   middlePoint :rhs[1][2],
 108
 109   ratnorm : (rightEnd-middlePoint)/(rightEnd-leftEnd),
 110   dnorm : denom(ratnorm),
 111
 112   if signChanges(extB)=0 then
 113     return([[lhs[1][1],rhs[1][2]],lhs[2],extB/dnorm^degP])
 114   else
 115     return(false)
 116   );
 117
 118
 119
 120 /* It compresses a certificate by checking the positivity */
 121 /* on adjacient intervals */
 122 /* resList contains items of the form */
 123 /* [[leftEnd,rightEnd],normFact,bernCoeff]] */
 124
 125 compressCertificate(certList) :=
 126   block([isList,resList,extBCoeff,firstItem,secondItem,degP,lhsNorm,
 127          dnorm,ratnorm,leftEnd,rightEnd,middlePoint,merged],
 128     degP: length(certList[1][3]),
 129     resList : [],
 130     isList : certList,
 131     while (length(isList)>1) do
 132       (
 133       firstItem : first(isList),
 134       secondItem : second(isList),
 135
 136       merged:bernsteinMerge(firstItem,secondItem),
 137
 138       if not(merged=false) then
 139         (
 140         extBCoeff: merged[3],
 141         lhsNorm: lst_gcd(extBCoeff),
 142
 143         isList : rest(isList,2),
 144
 145         isList : cons(merged,isList)
 146         )
 147       else
 148         (
 149
 150         resList : append(resList,[first(isList)]),
 151         isList: rest(isList)
 152         )
 153
 154       ),
 155     if length(isList)=1 then
 156       resList : endcons(first(isList),resList),
 157   return(resList)
 158   );
 159
 160
 161
 162 /* Certification of positivity for a square free polynomial and */
 163 /* for a polynomial with no roots in the considered interval    */
 164 deCasteljauNoCheckCertificateBetween(pol,var,
 165                        search_interval) :=
 166   block([sfList,sfCert,i,j,intUnion,noCertificate,
 167          searchLeftEnd:search_interval[1],
 168          searchRightEnd:search_interval[2],
 169          sg,
 170          degP:degree(pol,var),
 171          res,bernIntSet,item,bernInt,noSubInt,
 172          bcoeff,normFact,bsplit,leftEnd,rightEnd,middlePoint,bCoeff,
 173 lhsSplit, lhsNorm, rhsSplit, rhsNorm, ratnorm, dnorm, certificate,sgnCh],
 174
 175   res : [],
 176
 177   bCoeff : bernsteinCoeffList(pol,var,searchLeftEnd,searchRightEnd),
 178   sg: sgn(bCoeff[1]),
 179   normFact : apply(sarag_lcm,map(denom,bCoeff)),
 180   bCoeff : normFact*bCoeff,
 181   bernIntSet : [[search_interval,normFact,bCoeff]],
 182
 183   certificate:true,
 184   while not(bernIntSet=[]) and certificate do
 185     (
 186     item : first(bernIntSet),
 187     bernIntSet : rest(bernIntSet),
 188
 189     bernInt : first(item),
 190     normFact : second(item),
 191     bCoeff : third(item),
 192
 193     leftEnd : first(bernInt),
 194     rightEnd : second(bernInt),
 195
 196     if first(bCoeff)=0 then
 197      (
 198      certificate:false,
 199      res:[0,[leftEnd]],
 200      return(res)
 201      )
 202     else
 203      (
 204      if last(bCoeff)=0 then
 205       (
 206       certificate:false,
 207       res:[0,[rightEnd]],
 208       return(res)
 209       )
 210      else
 211       (
 212       if first(bCoeff)*last(bCoeff)<0 then
 213        (
 214        certificate:false,
 215        res:[0,[leftEnd,rightEnd]],
 216        return(res)
 217        )
 218       )
 219      ),
 220
 221      sgnCh:signChanges(bCoeff),
 222      if sgnCh=0 and certificate then
 223         (
 224         res : endcons([[leftEnd,rightEnd],normFact,bCoeff],res)
 225         )
 226      else
 227         (
 228         middlePoint : (leftEnd+rightEnd)/2,
 229         bsplit : bernsteinSplitWithZ(bCoeff,leftEnd,rightEnd,middlePoint),
 230
 231         lhsSplit : first(bsplit),
 232         lhsNorm : lst_gcd(lhsSplit),
 233
 234         rhsSplit : second(bsplit),
 235         rhsNorm : lst_gcd(rhsSplit),
 236
 237         ratnorm : (rightEnd-middlePoint)/(rightEnd-leftEnd),
 238         dnorm : denom(ratnorm),
 239
 240         bernIntSet : endcons([[leftEnd,middlePoint],
 241                              normFact/lhsNorm*dnorm^degP,
 242                              lhsSplit/lhsNorm],bernIntSet),
 243         bernIntSet : endcons([[middlePoint,rightEnd],
 244                              normFact/rhsNorm*dnorm^degP,
 245                              rhsSplit/rhsNorm],bernIntSet)
 246
 247         ) /* else */
 248    ), /* while */
 249   if certificate then
 250     (
 251     if COMPRESS_CERTIFICATE then
 252       return([sg,compressCertificate(isListCertSort(res))])
 253     else
 254       return([sg,isListCertSort(res)])
 255     )
 256   else
 257     return(res)
 258 );
 259
 260
 261 /* Certification of positivity for a generic polynomial */
 262 deCasteljauCertificateBetween(pol,var,search_interval) :=
 263   block([sqFreeCert,sfCert,i,j,intUnion,noCertificate,
 264          res,sqFree,item,bernInt,noSubInt,
 265          bcoeff,normFact,bsplit,leftEnd,rightEnd,middlePoint,exPol, searchLeftEnd, searchRightEnd],
 266
 267   exPol:expandIf(pol),
 268   searchLeftEnd : first(search_interval),
 269   searchRightEnd : second(search_interval),
 270
 271   sqFree: SQUARE_FREE_ALGORITHM(exPol,var),
 272   sqFreeCert: deCasteljauNoCheckCertificateBetween(sqFree,var,
 273                                                     search_interval),
 274
 275   if first(sqFreeCert)=0 then
 276     (
 277     certificate: false,
 278     return([0,sqFreeCert[2],sqFree])
 279     )
 280   else
 281     return(
 282            deCasteljauNoCheckCertificateBetween(exPol,var,
 283            search_interval)
 284            )
 285  );
 286
 287
 288
 289
 290 localBernsteinProof(poly,var,local_cert,printFlag):=
 291   block([cert_sgn,sgnCh,bernCoeffList,i, verifySum],
 292
 293   sgnCh:signChanges(local_cert[3]),
 294   bernCoeffList: (1/local_cert[2])*(local_cert[3]),
 295
 296   if sgnCh=0 then
 297    cert_sgn:sgn(local_cert[3][1])
 298   else
 299     cert_sgn:0,
 300   if cert_sgn=0 then
 301     (
 302     if printFlag then
 303       (
 304       sprint("It has a zero in ",local_cert[1])
 305       ),
 306     print("because in this interval it has "),
 307     print("the following Bernstein coefficients ",
 308             bernCoeffList)
 309     )
 310   else
 311     (
 312     if cert_sgn= 1 then
 313       (
 314       if printFlag then
 315         (
 316         sprint("It is positive in ",local_cert[1])
 317         ),
 318       print("because in this interval it has "),
 319
 320       print("the following Bernstein coefficients ",
 321             bernCoeffList)
 322       )
 323     else
 324           (
 325       if cert_sgn=-1 then
 326         (
 327         if printFlag then
 328           (
 329           sprint("It is negative in ",local_cert[1])
 330           ),
 331         print("because in this interval it has "),
 332         print("the following Bernstein coefficients ",
 333               bernCoeffList)
 334         )
 335           ),
 336           print("as we see below: "),
 337           for i: 0 thru length(bernCoeffList) - 2 do
 338             (
 339                         print(bernCoeffList[i+1], bernstein(length(bernCoeffList)-1,i, local_cert[1][1], local_cert[1][2],var), "+")
 340             ),
 341           print(bernCoeffList[length(bernCoeffList)], bernstein(length(bernCoeffList) - 1, length(bernCoeffList) - 1, local_cert[1][1], local_cert[1][2],var), "="),
 342           for i: 0 thru length(bernCoeffList) - 2 do
 343             (
 344                         print(bernCoeffList[i+1] * ratsimp(bernstein(length(bernCoeffList)-1,i, local_cert[1][1], local_cert[1][2],var)), "+")
 345             ),
 346           print(bernCoeffList[length(bernCoeffList)] * ratsimp(bernstein(length(bernCoeffList) - 1, length(bernCoeffList) - 1, local_cert[1][1], local_cert[1][2],var)), "="),
 347       verifySum : 0,
 348           for i: 0 thru length(bernCoeffList) - 1 do
 349             (
 350                     verifySum : verifySum + bernCoeffList[i+1] * ratsimp(bernstein(length(bernCoeffList)-1,i, local_cert[1][1], local_cert[1][2],var))
 351                 ),
 352           print(ratsimp(verifySum))
 353     ),
 354   return(true)
 355   );
 356
 357
 358 bernsteinProof(poly,var,search_interval,cert):=
 359   block([i,exPol],
 360   exPol:expandIf(pol),
 361
 362   if cert[1]=0 then
 363      (
 364      sprint("The polynomial ", exPol, " has a zero in ",
 365             search_interval, "."),
 366      if length(cert[2])= 1 then /* root found */
 367        print("In fact it is zero for ", var, " = ", cert[2][1])
 368
 369      else
 370        (
 371        if expand(cert[3])=exPol then
 372          ( /* square-free, interval found */
 373
 374          print("In fact it has value ",
 375                 subst(cert[2][1], var, exPol), " for ", var, " = ", cert[2][1]),
 376          print("and it has value ",
 377                 subst(cert[2][2], var, exPol), " for ", var, " = ", cert[2][2])
 378          )
 379        else
 380          ( /* non-square-free, interval found */
 381
 382          print("In fact its square-free part ", cert[3]),
 383          print("has value ", subst(cert[2][1],var,cert[3]),
 384                " for ", var, " = ", cert[2][1]),
 385          print("and  has value ",
 386                subst(cert[2][2],var,cert[3])," for ", var, " = ",  cert[2][2])
 387          )
 388        )
 389      )
 390   else
 391     if cert[1]=1 then
 392       (
 393       sprint("The polynomial ", exPol, " is positive in ",
 394               search_interval),
 395       if length(cert[2])>1 then
 396         (
 397         sprint("."),
 398         print("In fact we have that:  "),
 399         for i: 1 thru length(cert[2]) do
 400           (
 401           sprint("(",i,")"),
 402           localBernsteinProof(poly,var,cert[2][i],true)
 403           )
 404         )
 405       else
 406         localBernsteinProof(poly,var,cert[2][1],false)
 407       )
 408     else
 409       (
 410       sprint("The polynomial ", exPol, " is negative in ",
 411               search_interval),
 412       if length(cert[2])>1 then
 413         (
 414         sprint("."),
 415         print("In fact we have: "),
 416         for i: 1 thru length(cert[2]) do
 417           (
 418           sprint("(",i,")"),
 419           localBernsteinProof(poly,var,cert[2][i],true)
 420           )
 421         )
 422       else
 423         localBernsteinProof(poly,var,cert[2][1],false)
 424       ),
 425   print(""),
 426
 427   if cert[1]=0 then
 428     return(1)
 429   else
 430     return(length(cert[2]))
 431  );
 432
 433
 434 certificateProofBetween(pol,var, search_interval):=
 435   bernsteinProof(pol,var,search_interval, certificateBetween(pol,var,search_interval));
 436
 437
 438
 439 certificateProof([args]):=
 440   if length(args)=2 then
 441     certificateProofBetween(args[1],args[2],
 442                             [DEFAULT_LEFT_END, DEFAULT_RIGHT_END])
 443   else
 444     certificateProofBetween(args[1],args[2], args[3]);