share/contrib/sarag/rootIsolation.mac

   1 /* ------------------------------------------------------------------- */
   2 /* SARAG - Root Isolation                                              */
   3 /* by Fabrizio Caruso modified by Alexandre Le Meur and Marie-Françoise Roy */
   4
   5
   6 /* i-th Bernstein polynomial of degree p for c and d */
   7 bernstein(p,i,c,d,var) :=
   8   binomial(p,i) * (var-c)^i * (d-var)^(p-i) / (d-c)^p;
   9
  10
  11 /* -------------------------------------- */
  12 /* Bounds on real roots                   */
  13
  14 /* Cauchy Root Upper Bound */
  15 /* Lemma 10.2 */
  16 /* (Defined for non-zero polynomials) */
  17 cauchyRootUpperBound(pol,var) :=
  18   block([p,q,i,den],
  19     p : expandIf(p),
  20     p : degree(pol,var),
  21     q : 0,
  22     while coeff(pol,var,q)=0 do
  23        q : q + 1,
  24     den : abs(coeff(pol,var,p)),
  25
  26     return(sum(abs(coeff(pol,var,i)),i,q,p)/den)
  27     ); /* end proc */
  28
  29
  30 /* Cauchy Root Lower Bound for non-zero roots */
  31 cauchyRootLowerBound(pol,var) :=
  32    block([p,q,i,den],
  33     p : expandIf(p),
  34     p : degree(pol,var),
  35     q : 0,
  36     while coeff(pol,var,q)=0 do
  37        q : q + 1,
  38     den : abs(coeff(pol,var,q)),
  39     return((sum(abs(coeff(pol,var,i)),i,q,p)/den)^(-1))
  40     ); /* end proc */
  41
  42
  43
  44 /* Prime Cauchy Root Upper Bound */
  45 /* Lemma 10.5 */
  46 primeCauchyRootUpperBound(pol,var) :=
  47   block([p,q,i,den],
  48     p : expandIf(p),
  49     p : degree(pol,var),
  50     q : 0,
  51     while coeff(pol,var,q)=0 do
  52        q : q + 1,
  53     den : coeff(pol,var,p)^2,
  54     return((p+1)* sum(coeff(pol,var,i)^2,i,q,p)/den)
  55     ); /* end proc */
  56
  57
  58 /* Prime Cauchy Root Lower Bound for non-zero roots */
  59 primeCauchyRootLowerBound(pol,var) :=
  60   block([p,q,i,den],
  61     p : expandIf(p),
  62     p : degree(pol,var),
  63     q : 0,
  64     while coeff(pol,var,q)=0 do
  65        q : q + 1,
  66     den : coeff(pol,var,q)^2,
  67     return(((p+1)* sum(coeff(pol,var,i)^2,i,q,p)/den)^(-1))
  68     ); /* end proc */
  69
  70
  71 /* ----------------------------------------- */
  72 /* Auxiliary routines for sign determination */
  73
  74 /* It determines the sign */
  75 /* to the right of the left end of open interval  */
  76 signToTheLeft(bern) :=
  77   block([i],
  78     i : 1,
  79     while(sgn(bern[i])=0) do
  80       i : i+1,
  81     return(sgn(bern[i]))
  82     ); /* end proc */
  83
  84
  85 /* It determines the sign */
  86 /* to the left of the right end of open interval  */
  87 signToTheLeft(bern) :=
  88   block([i],
  89     i : length(bern),
  90     while(sgn(bern[i])=0) do
  91       i : i-1,
  92     return(sgn(bern[i]))
  93     ); /* end proc */
  94
  95
  96
  97
  98 /* ------------------------------------- */
  99 /* Splitting of Bezier curves */
 100
 101 /* Alg. 10.31 */
 102 bernsteinSplit(coeffList, c, d,e) :=
 103   block([p : length(coeffList)-1,
 104          b : make_array('any),
 105          bern_ce,
 106          bern_ed,
 107          alpha, beta,
 108          i,j],
 109
 110   b : make_array('any,p+1),
 111   for i : 0 thru p do
 112     b[i] : coeffList[i+1],
 113
 114   bern_ce : [b[0]],
 115   bern_ed : [b[p]],
 116   alpha : (d-e)/(d-c),
 117   beta : (e-c)/(d-c),
 118   for i : 1 thru p do
 119     (
 120     for j : 0 thru p-i do
 121        b[j] : alpha * b[j] + beta * b[j+1],
 122     bern_ce : endcons(b[0],bern_ce),
 123
 124     bern_ed : cons(b[p-i],bern_ed)
 125
 126     ), /* end for */
 127   return([bern_ce,bern_ed])
 128   ); /* end proc */
 129
 130
 131 /* Algorithm 10.34 */
 132 specialBernsteinSplit(coeffList,c,d) :=
 133  block([p : length(coeffList)-1,
 134         bern_first,
 135         bern_second,
 136         b : make_array('any),
 137         i,j],
 138
 139   b : make_array('any,p+1),
 140   for i : 0 thru p do
 141     b[i] : coeffList[i+1],
 142
 143   bern_first : [2^p * b[0]],
 144   bern_second : [2^p * b[p]],
 145
 146   for i : 1 thru p do
 147     (
 148     for j : 0 thru p-i do
 149         b[j] : b[j] + b[j+1],
 150
 151    bern_first : endcons(2^(p-i) * b[0],bern_first),
 152
 153    bern_second : cons(2^(p-i)*b[p-i],bern_second)
 154
 155     ), /* end for */
 156
 157   return([bern_first,bern_second])
 158   ); /* end proc */
 159
 160
 161
 162 /* Square free part of a polynomial */
 163 squareFree(pol,var) := gcdFreePartWithZ(pol,diff(pol,var),var)[2];
 164
 165 /* ---------------------------------------------- */
 166 /* Auxiliary routines for polynomial manipulation */
 167
 168 translated(pol,var,c) := subst(var-c,var,pol);
 169
 170 reciprocal(pol,var,p) := var^p * subst(1/var,var,pol);
 171
 172 /* contraction */
 173 /* (defined for l different from zero) */
 174 contracted(pol,var,l) := subst(var*l,var,pol);
 175
 176 /* auxiliary coefficients as defined in page 340 */
 177 /* (in the book notation: a_(c,k))               */
 178 specialCoeff(l,c,k) :=
 179   (-2^(l+k) + c*2^(l+1))/(2^k);
 180
 181
 182
 183 /* Generic base logarithm */
 184 logInBase(val,base) :=
 185    log(val)/log(base);
 186
 187 discreteLogInBaseTwo(val) :=
 188    block([i],
 189    i : 0,
 190    while 2^i < val do
 191      i : i + 1,
 192    return(i)
 193    ); /* end proc */
 194
 195 bernsteinCoeffList(pol,var,c,d) :=
 196   convert2Bernstein(pol,degree(pol,var),var,c,d);
 197
 198 /* Bernstein coefficients */
 199 /* computed by using Proposition 10.24 */
 200 convert2Bernstein(p,polDeg,var,c,d) :=
 201   block([pol,auxPol,transMinusC, conDMinusC, recP, transMinusOne,res,i],
 202
 203 pol:p,
 204
 205     if polDeg<degree(pol,var) then ( print("bad degree, d must be more than degree(P)"),
 206     return(false))
 207
 208     else
 209
 210     pol : expandIf(pol),
 211     polDeg : max(degree(pol,var),polDeg),
 212
 213     transMinusC : expand(translated(pol,var,-c)),
 214     conDMinusC : expand(contracted(transMinusC,var,d-c)),
 215
 216     recP : expand(reciprocal(conDMinusC,var,polDeg)),
 217
 218     transMinusOne : expand(translated(recP,var,-1)),
 219
 220     auxPol : translated(reciprocal(
 221              contracted(translated(
 222              pol,var,-c),var,d-c),var,polDeg),var,-1),
 223     auxPol : ratsimp(auxPol,var),
 224
 225     res : makelist(ratcoeff(auxPol,var,polDeg-i)/
 226                    binomial(polDeg,i),i,0,polDeg),
 227     return(res)
 228     ); /* end proc */
 229
 230
 231 /* For a given squarefree polynonial SFP of degree p  */
 232 /* Bernstein coefficients of the polynomial           */
 233 /* 2^(k*p)*SFP in the Bernstein basis                 */
 234 /* (specialCoeff(l,c,k),specialCoeff(l,c+1,k)))           */
 235 bernsteinSpecialCoeffList(pol,var,leftEnd,rightEnd) :=
 236    bernsteinCoeffList(pol,var,leftEnd,rightEnd);
 237
 238
 239
 240
 241
 242 /* --------------------------------------------------- */
 243 /* Real root isolation */
 244
 245 isListLess(lhs,rhs) :=
 246   if length(lhs) = 1 then
 247     if length(rhs) = 1 then
 248       lhs[1]<rhs[1]
 249     else
 250       lhs[1]<=rhs[1]
 251   else
 252     if length(rhs) = 1 then
 253       lhs[2]<=rhs[1]
 254     else
 255       lhs[1]<rhs[1];
 256
 257
 258 /* Order Isolating List */
 259 isListSort(isList) :=
 260   sort(isList,isListLess);
 261
 262
 263 isListWithSignsLess(lhs,rhs) :=
 264   isListLess(first(lhs),first(rhs));
 265
 266 /* Order Isolating list with signs */
 267 isListWithSignsSort(isListWithSigns) :=
 268   sort(isListWithSigns,isListWithSignsLess);
 269
 270
 271
 272 /* Real root isolation in an archimedean real closed field */
 273 /* the output is a list of elements of the following form: */
 274 /* single point : [p] */
 275 /* interval : [[a,b],[c,k]], where c, k are the coeffs in a_{c,k} */
 276
 277 /* Aliases for "deCasteljauIsolateRoots..." */
 278
 279 /* Aliases for "deCasteljauIsolateRootsGeneral" */
 280
 281 deCasteljauIsolateRootsWithARCF(pol,var) := /* ARCF stands for Archimedean Real Closed Field */
 282   deCasteljauIsolateRootsGenStruct(pol,var,false);
 283
 284 deCasteljauIsolateRootsWithZ(pol,var) :=
 285   deCasteljauIsolateRootsGenStruct(pol,var,true);
 286
 287 /* Aliases for "deCasteljauIsolateRootsGenStructVerbose" */
 288
 289 deCasteljauIsolateRootsGenStruct(pol,var,integer_flag) :=
 290   deCasteljauFindRootsGenStruct(pol,var,inf,integer_flag);
 291
 292 deCasteljauIsolateRootsBetweenWithARCF(pol,var,a,b) :=
 293   deCasteljauFindRootsGenIntGenStructCert(pol,var,inf,a,b,false,false)[1];
 294
 295 deCasteljauIsolateRootsBetweenWithZ(pol,var,a,b) :=
 296   deCasteljauFindRootsGenIntGenStructCert(pol,var,inf,a,b,true,false)[1];
 297
 298 deCasteljauIsolateRootsBetweenWithZCert(pol,var,a,b) :=
 299   deCasteljauFindRootsGenIntGenStructCert(pol,var,inf,a,b,true,true);
 300
 301
 302 /* Aliases for "deCasteljauFindRootsGeneral..." */
 303
 304 /* Aliases for "deCasteljauFindRootsGeneral" */
 305
 306 deCasteljauFindRootsWithARCF(pol,var,threshold) :=
 307   deCasteljauFindRootsGenStruct(pol,var,threshold,false);
 308
 309 deCasteljauFindRootsWithZ(pol,var,threshold) :=
 310   deCasteljauFindRootsGenStruct(pol,var,threshold,true);
 311
 312 /* Aliases for "deCasteljauFindRootsGenStruct" */
 313
 314 deCasteljauFindRootsBetweenWithZ(pol,var,threshold,a,b) :=
 315   deCasteljauFindRootsGenIntGenStructCert(pol,var,threshold,a,b,true,false)[1];
 316
 317 deCasteljauFindRootsWithZ(pol,var,threshold) :=
 318   deCasteljauFindRootsGenStruct(pol,var,threshold,true);
 319
 320 deCasteljauFindRootsGenStruct(pol,var,threshold,integer_flag) :=
 321   deCasteljauFindRootsGenIntGenStructCert(pol,var,threshold,false,0,integer_flag,false)[1];
 322
 323 deCasteljauFindRootsGenStructCert(pol,var,threshold,integer_flag) :=
 324   deCasteljauFindRootsGenIntGenStructCert(pol,var,threshold,false,0,integer_flag,true)[1];
 325
 326 /* definition de deCasteljauFindRootsGenIntGenStructCert */
 327
 328 deCasteljauFindRootsGenIntGenStructCert(pol,var,threshold,a,b,integer_flag,modifCert) :=
 329   block([leftEnd,rightEnd,midPoint,n,sqPol,bCoeff,pos,deg,
 330          resList,item,sgnCh,splitRes,c,k,count,gcdTest,normFact,searchLeftEnd, searchRightEnd,mod,sg,certificate,P,resCert,bernIntSet,bsplit,lhsSplit, rhsSplit,bernInt],
 331
 332
 333   P : expandIf(pol),
 334   deg : degree(P,var),
 335
 336
 337   if a = false then
 338      (
 339      n : discreteLogInBaseTwo(cauchyRootUpperBound(P,var)),
 340      searchLeftEnd : -(2^n),
 341      searchRightEnd : 2^n
 342      ) /* end if */
 343   else
 344      (
 345      searchLeftEnd : a,
 346      searchRightEnd : b
 347      ), /* end else */
 348
 349   sqPol : expand(factor(squareFree(P,var))),
 350
 351   resList : [],
 352   resCert : [],
 353   if modifCert then (
 354   bCoeff : bernsteinCoeffList(P,var,searchLeftEnd,searchRightEnd)
 355
 356      ) /* end if */
 357  else (
 358        bCoeff : bernsteinCoeffList(sqPol,var,searchLeftEnd,searchRightEnd)
 359              ), /* end else */
 360   sg: sgn(bCoeff[1]),
 361   normFact : 0,
 362
 363 /* normFact is the multiplicative coefficient which makes the berstein coefficients integers*/
 364
 365   if integer_flag then
 366      (
 367      normFact : apply(lcm,map(denom,bCoeff)),
 368      bCoeff : normFact*bCoeff
 369      ), /* end if */
 370
 371
 372   bernIntSet : [[[searchLeftEnd,searchRightEnd],normFact,bCoeff]],
 373   certificate : true,
 374   while not(bernIntSet=[]) and certificate do
 375     (
 376     item : first(bernIntSet),
 377     bernIntSet : rest(bernIntSet),
 378
 379     bernInt : first(item),
 380     bCoeff : third(item),
 381
 382     leftEnd : first(bernInt),
 383     rightEnd : second(bernInt),
 384
 385   if modifCert then
 386    (
 387     if first(bCoeff)=0 then
 388       (
 389       certificate:false,
 390       resCert:[false,[[leftEnd],sqPol],true,false]
 391       )
 392     else
 393       (
 394       if last(bCoeff)=0 then
 395        (
 396        certificate:false,
 397        resCert:[false,[[rightEnd],sqPol],true,false]
 398        ) /* end if */
 399      ) /* end else */
 400     ), /* end if */
 401     sgnCh : signChanges(bCoeff),
 402    if (sgnCh = 1) and (rightEnd-leftEnd)<threshold then
 403
 404         (if modifCert then      /* Loop which is useful for the function in the certificate file */
 405          ((certificate : false),
 406          resCert : [false,[[leftEnd,rightEnd],sqPol],true,false]
 407                ) /* end if */
 408          else
 409           resList : endcons([leftEnd,rightEnd],resList)
 410             ) /* end if */
 411
 412     else
 413        (
 414
 415        midPoint : (leftEnd+rightEnd)/2,
 416        if (sgnCh>=1) then
 417           (
 418           if subst(midPoint,var,pol)=0 then
 419              (
 420              resList : endcons([midPoint],resList)
 421                           ), /* end if */
 422
 423          bsplit : bernsteinSplitWithZ(bCoeff,leftEnd,rightEnd,midPoint),
 424          lhsSplit : first(bsplit),
 425
 426         rhsSplit : second(bsplit),
 427
 428         bernIntSet : endcons([[leftEnd,midPoint],
 429                              normFact,lhsSplit],bernIntSet),
 430         bernIntSet : endcons([[midPoint,rightEnd],
 431                              normFact,rhsSplit],bernIntSet),
 432
 433         if rightEnd-midPoint<1 then
 434
 435          ( normFact : normFact*2^deg ) /* end if */
 436
 437
 438
 439                    ) /* end if */
 440       else ( if modifCert then
 441                resCert : endcons([[leftEnd,rightEnd],normFact,bCoeff],resCert)
 442
 443               )  /* end else */
 444        ) /* end else */
 445     ), /* end while */
 446   if certificate then (
 447         if modifCert then
 448             if isListCertSort(resCert)[1][3][1]>0 then
 449                   return([positive,compressCertificate(isListCertSort(resCert)),0,false,true])
 450             else
 451                   return([negative,compressCertificate(isListCertSort(resCert)),0,false,true])
 452         else
 453         return([isListSort(resList),sqPol,false,false])
 454          ) /* end if */
 455   else
 456   return(resCert)
 457   );
 458  /* end proc */
 459
 460
 461 deCasteljauNumberOfRoots(pol,var) :=
 462   length(deCasteljauIsolateRoots(pol,var));
 463
 464 deCasteljauNumberOfRootsBetween(pol,var,a,b) :=
 465   length(deCasteljauIsolateRootsBetween(pol,var,a,b));
 466
 467 /* Refine an open interval describing a unique real root */
 468 refineRoot(rootInterval,pol,var) :=
 469   block([leftEnd,rightEnd,midPoint,n,sqPol,bCoeffList,pos,
 470          resList,item,sgnCh,splitRes,c,k,count,gcdTest,
 471          lhs,rhs],
 472   leftEnd : rootInterval[1],
 473   rightEnd : rootInterval[2],
 474   midPoint : (leftEnd+rightEnd)/2,
 475   if subst(midPoint,var,pol)=0 then
 476     return([midPoint])
 477   else
 478     (
 479     resList : [],
 480     sqPol : expand(factor(squareFree(pol,var))),
 481
 482
 483     bCoeffList : bernsteinCoeffList(sqPol,var,leftEnd,rightEnd),
 484
 485     splitRes : specialBernsteinSplit(bCoeffList,leftEnd,rightEnd),
 486     lhs : splitRes[1],
 487     rhs : splitRes[2],
 488     if signChanges(lhs)=1 then
 489       return([leftEnd,midPoint])
 490     else
 491       return([midPoint,rightEnd])
 492     ) /* end else */
 493 ); /* end proc */
 494
 495
 496
 497 /* Finds an intermediate point (already in the right order) */
 498 intermidiatePoint(lhs,rhs,poly,var) :=
 499   if length(lhs)=1 then
 500     if length(rhs) = 1 then
 501       (lhs[1]+rhs[1])/2
 502     else
 503        if not(lhs[1]=rhs[1]) then
 504          (lhs[1]+rhs[1])/2
 505        else
 506          intermidiatePoint(lhs,refineRoot(rhs,poly,var),poly,var)
 507   else
 508     if length(rhs)=1 then
 509       if not(lhs[2]=rhs[1]) then
 510         (lhs[2]+rhs[1])/2
 511       else
 512         intermidiatePoint(refineRoot(lhs,poly,var),rhs,poly,var)
 513     else
 514       if not(lhs[2]=rhs[1]) then
 515         (lhs[2]+rhs[1])/2
 516       else
 517         lhs[2];
 518
 519
 520 /* -------------------------------------------------------- */
 521 /* Sign or real roots */
 522
 523 /* Sign of real roots */
 524 /* in an archimedean real closed field */
 525
 526 deCasteljauRootsSign(isolListForP,p,q,var) :=
 527   block([sfP,sfQ,g,commonPtList,nonCommonPtList,commonInList,nonCommonInList,
 528          candList,nCom,sfPCoeffList,gCoeffList,
 529          bound,i,l,c,k,
 530          leftEnd,rightEnd,midPt,item,bernG, bernSfP, bernQ,
 531          sgnM,sgnL,sgnR, bernSplitSfP, bernSplitG],
 532
 533   p : expandIf(p),
 534   q : expandIf(q),
 535   sfP : expand(factor(squareFree(p,var))),
 536   g : gcd(sfP,q),
 537
 538   commonPtList : [],  /* In the book's notation : N(P) */
 539   commonInList : [],
 540   nonCommonPtList : [],
 541   nonCommonInList : [],
 542   candList : [], /* In the book's notation : Pos */
 543   nCom : [],
 544   for i : 1 thru length(isolListForP) do
 545     (
 546
 547     if length(isolListForP[i])=1 then
 548        (
 549        sgnM : sgn(subst(isolListForP[i][1],var,q)),
 550        if sgnM = 0 then
 551           commonPtList : endcons([isolListForP[i],sgnM],commonPtList)
 552        else
 553           nonCommonPtList : endcons([isolListForP[i],sgnM],nonCommonPtList)
 554        ) /* end if */
 555     else
 556        (
 557
 558        leftEnd : isolListForP[i][1],
 559        rightEnd : isolListForP[i][2],
 560
 561        candList : endcons([[bernsteinCoeffList(sfP,var,leftEnd,rightEnd),
 562                          bernsteinCoeffList(g,var,leftEnd,rightEnd)],
 563                          isolListForP[i]],
 564                          candList)
 565        ) /* end else */
 566     ), /* end for */
 567
 568   while not(candList=[]) do
 569     (
 570     item : first(candList),
 571
 572     candList : rest(candList),
 573     bernG : item[1][2],
 574     bernSfP : item[1][1],
 575     leftEnd : item[2][1],
 576     rightEnd : item[2][2],
 577
 578     if signChanges(bernG) = 1 then
 579       commonInList : endcons([item[2],0],commonInList)
 580     else
 581       if signChanges(bernG) = 0 then
 582         nCom : endcons([bernSfP,[leftEnd,rightEnd]],nCom)
 583       else
 584         (
 585         midPt : (leftEnd+rightEnd)/2,
 586         sgnL : subst(leftEnd,var,sqP),
 587         sgnM : subst(midPt,var,sqP),
 588         sgnR : subst(rightEnd,var,sqP),
 589
 590         if sgnM = 0 then
 591            (
 592            if sgn(subst(midPt,var,q)) = 0 then
 593              commonPtList : endcons([[midPt],0],commonPtList)
 594            else
 595              nonCommonPtList : endcons([[midPt],sgn(subst(midPt,var,q))],nonPtCommonList)
 596            ) /* end if */
 597         else
 598            (
 599            bernSplitSfP : specialBernsteinSplit(bernSfP,leftEnd,rightEnd),
 600            bernSplitG : specialBernsteinSplit(bernG,leftEnd,rightEnd),
 601            if sgnL*sgnM<0 then
 602               candList : endcons([[bernSplitSfP[1],bernSplitG[1]],
 603                                [leftEnd,midPt]
 604                               ],candList)
 605
 606            else
 607               candList : endcons([[bernSplitSfP[2],bernSplitG[2]],
 608                                 [midPt,rightEnd]
 609                               ],candList)
 610            ) /* end else */
 611
 612         ) /* end else */
 613      ), /* end while */
 614
 615 /* second part */
 616   candList : nCom,
 617
 618   while not(candList=[]) do
 619     (
 620     item : first(candList),
 621     candList : rest(candList),
 622     leftEnd : item[2][1],
 623     rightEnd : item[2][2],
 624
 625     bernSfP : item[1],
 626
 627     bernQ : bernsteinCoeffList(q,var,leftEnd,rightEnd),
 628
 629     if signChanges(bernQ)=0 then
 630       (
 631
 632       i : 1,
 633       while(bernQ[i]=0) do
 634         i : i+1,
 635
 636
 637       nonCommonInList : endcons([[leftEnd,rightEnd],
 638                                sgn(bernQ[i])],nonCommonInList)
 639       ) /* end if */
 640     else
 641       (
 642       midPt : (leftEnd+rightEnd)/2,
 643       sgnM : sgn(subst(midPt,var,sfP)),
 644
 645       i : 1,
 646       while(bernSfP[i]=0) do
 647         i : i+1,
 648       sgnL : sgn(bernSfP[i]),
 649
 650       i : length(bernSfP),
 651       while(bernSfP[i]=0) do
 652         i : i - 1,
 653       sgnR : sgn(bernSfP[i]),
 654
 655
 656       if sgnM=0 then
 657          nonCommonPtList : endcons([[midPt],sgn(subst(midPt,var,q))],nonCommonPtList)
 658       else
 659          (
 660          if sgnL*sgnM<0 then
 661             candList : endcons([bernsteinCoeffList(sfP,var,leftEnd,midPt),
 662                              [leftEnd,midPt]
 663                             ],candList)
 664          else
 665             candList : endcons([bernsteinCoeffList(sfP,var,midPt,rightEnd),
 666                              [midPt,rightEnd]
 667                             ],candList)
 668          ) /* end else */
 669       ) /* end ele */
 670     ), /* end if */
 671
 672   return([[nonCommonPtList,nonCommonInList],[commonPtList,commonInList]])
 673   ); /* end proc */
 674
 675
 676 deCasteljauSignsAtRoots(isolListForP,p,q,var) :=
 677   isListWithSignsSort(lambda([l],append(l[1][1],l[1][2],
 678                     l[2][1],l[2][2]))(
 679          deCasteljauRootsSign(isolListForP,p,q,var)));
 680
 681
 682 /* --------------------------------------- */
 683 /* Sign Comparison of real roots */
 684 /* in the Z ring */
 685
 686
 687 /* Rewrite a set of isolating points with signs as intervals and removes the sign */
 688 makeAsIn(ptList) :=
 689   if length(ptList)=0 then
 690     []
 691   else
 692     cons([ptList[1][1][1],ptList[1][1][1]],makeAsIn(rest(ptList)));
 693
 694 /* Rewrite dummy intervals as points */
 695 makeAsPt(inList) :=
 696   if length(inList)=0 then
 697    []
 698   else
 699    if inList[1][1]=inList[1][2] then
 700      cons([inList[1][1]],makeAsPt(rest(inList)))
 701    else
 702      cons(inList[1],makeAsPt(rest(inList)));
 703
 704
 705 /* Removes the sign */
 706 removeSign(inList) :=
 707   if length(inList)=0 then
 708     []
 709   else
 710     cons(first(first(inList)),removeSign(rest(inList)));
 711
 712 /* Interval Intersection */
 713 intersectIn(inList1,inList2) :=
 714   block([i,len1,len2, leftEnd1,rightEnd1,leftEnd2,rightEnd2,res],
 715
 716   len1 : length(inList1),
 717   len2 : length(inList2),
 718
 719
 720   if len1=0 then
 721     return(inList2)
 722   else
 723     if len2=0 then
 724       return(intList1),
 725
 726   res : [],
 727   i : 1,
 728
 729   while (i<=len1)  do
 730     (
 731     leftEnd1 : inList1[i][1],
 732     leftEnd2 : inList2[i][1],
 733     rightEnd1 : inList1[i][2],
 734     rightEnd2 : inList2[i][2],
 735
 736     if (leftEnd1<=leftEnd2) and
 737        (leftEnd2<=rightEnd1) then
 738        (
 739        res : endcons([leftEnd2,min(rightEnd1,rightEnd2)],res),
 740        i : i + 1
 741        )   /* end if */
 742     else
 743       if (leftEnd2<=leftEnd1) and
 744          (leftEnd1<=rightEnd2) then
 745          (
 746          res : endcons([leftEnd1,min(rightEnd1,rightEnd2)],res),
 747          i : i + 1
 748          ) /* end if */
 749       else
 750          (
 751          print("impossible configuration"),
 752          return([])
 753          ) /* end else */
 754     ), /* end while */
 755     return(res)
 756
 757
 758   );  /* end proc */
 759
 760
 761
 762 /* Evaluation of the signs of two polynomials at their roots */
 763 /* in an archimedean real closed field */
 764
 765 deCasteljauEvaluateSignsAtRoots(p,q,var):=
 766   block([l,isolP,isolQ,
 767          signComPtP,signComInP,signComPtQ,signComInQ,
 768          signNComPtP,signNComInP,signNComPtQ,signNComInQ,
 769          signP,signQ,signComP,signNComP,signComQ,signNComQ,
 770          signComAsInP,signComAsInQ,comIn],
 771
 772   p : expandIf(p),
 773   q : expandIf(q),
 774
 775   isolP : deCasteljauIsolateRoots(p,var),
 776
 777   isolQ : deCasteljauIsolateRoots(q,var),
 778   signP : deCasteljauRootsSign(isolP,p,q,var),
 779   signComP : second(signP),
 780
 781   signComPtP : first(signComP),
 782   signComInP : second(signComP),
 783   signNComP : first(signP),
 784   signNComP : sort(append(first(signNComP),second(signNComP))),
 785
 786   signQ : deCasteljauRootsSign(isolQ,q,p,var),
 787   signComQ : second(signQ),
 788
 789   signComPtQ : signComQ[1],
 790   signComInQ : signComQ[2],
 791   signNComQ : first(signQ),
 792   signNComQ : sort(append(first(signNComQ),second(signNComQ))),
 793
 794   signComAsInP : sort(append(makeAsIn(signComPtP),removeSign(signComInP))),
 795   signComAsInQ : sort(append(makeAsIn(signComPtQ),removeSign(signComInQ))),
 796
 797 comIn : makeAsPt(intersectIn(signComAsInP,signComAsInQ)),
 798 return([comIn,signNComP,signNComQ])
 799   ); /* end proc */
 800
 801
 802 /* Real root sample points */
 803 /* for an archimedean real closed field */
 804 realSamplePointsARCF(polList) :=
 805   block([i,j,listSize,compRes],
 806   listSize : length(polList),
 807   for i : 1 thru listSize do
 808     for j : 1 thru listSize do
 809       compRes : realRootComparisonARCF(p,q,var),
 810   return([0])
 811   );
 812
 813
 814