share/contrib/sarag/signDetermination.mac

   1 /* ---------------------------------------------------------------- */
   2 /* SARAG - Sign Determination                                       */
   3 /* by Fabrizio Caruso   modified by Alexandre Le Meur and Marie-Françoise Roy */
   4
   5
   6 /* ----------------------------------------- */
   7 /* Backward compatibility issues since 5.9.2 */
   8
   9 /* Auxiliary floor function */
  10 /* unsed */
  11 /*
  12 comp_floor(x) :=
  13   fix(x);
  14 */
  15
  16 /* Auxiliary ceiling function */
  17 comp_ceiling(x) :=
  18   if x=float(fix(x)) then
  19      fix(x)
  20   else
  21      fix(x)+1;
  22
  23 /* Uncenter the modulo representation */
  24 /* It works for both odd and even positive integer m */
  25 comp_uncenter(x,m) :=
  26   if(x >= 0) then
  27      x
  28   else
  29      m + x;
  30
  31
  32 /* ----------------------------------------------------------- */
  33 /* Realization of sign conditions                              */
  34
  35 ceilFract(num,den) :=
  36   block([i],
  37    i : 1,
  38    while(i*den<num) do
  39      i: i+1,
  40    return(i)
  41    ); /* end proc */
  42
  43 /* Auxiliary block matrix computation */
  44 /* defined for blocks of same size    */
  45 /* */
  46 /* nRows,nCols : number of rows/cols of the result */
  47 /* nBlockRows, nBlockCols : number of rows/cols in the blocks */
  48 /* inBlockRowIndex, inBlockColIndex : row/col in the block */
  49 blockMatrix(descr) :=
  50   block([nRows,rCols,nBlockRows,i,j,rowRes,res,
  51          blockRow,blockCol,inBlockRows, inBlockCols,
  52          inBlockRowIndex,inBlockColIndex, nCols],
  53      nBlockRows : length(descr),
  54      nBlockCols : length(first(descr)),
  55      inBlockRows : length(first(first(descr))),
  56      nRows : nBlockRows*inBlockRows,
  57      inBlockCols : length(first(first(first(descr)))),
  58      nCols : nBlockCols*inBlockCols,
  59      res : [],
  60
  61      for i : 1 thru nRows do
  62        (
  63        rowRes : [],
  64        for j : 1 thru nCols do
  65          (
  66          blockRow : ceilFract(i,inBlockRows),
  67          blockCol : ceilFract(j,inBlockCols),
  68
  69          if inBlockRows = 1 then
  70             inBlockRowIndex : 1
  71          else
  72          /* after 5.9.2 */
  73            /*
  74             inBlockRowIndex : mod(i-1,inBlockRows)+1,
  75            */
  76          /* before 5.9.2 but also forwards-compatible */
  77             inBlockRowIndex : comp_uncenter(mod(i-1,inBlockRows),inBlockRows)+1,
  78          if inBlockCols = 1 then
  79             inBlockColIndex : 1
  80          else
  81          /* after 5.9.2 */
  82             /*
  83             inBlockColIndex : mod(j-1,inBlockCols)+1,
  84             */
  85          /* before 5.9.2 but also forwards-compatible */
  86             inBlockColIndex : comp_uncenter(mod(j-1,inBlockCols),inBlockRows)+1,
  87          rowRes : endcons(
  88                   descr[blockRow][blockCol][inBlockRowIndex][inBlockColIndex],
  89                   rowRes)
  90          ), /* end for */
  91        res : endcons(rowRes,res)
  92        ), /* end for */
  93      return(res)
  94      ); /* end proc */
  95
  96
  97
  98 /* Matrix tensor product */
  99 /* Notation as in Notation 2.63 */
 100 tensorProduct(m,mPrime) :=
 101   block([i,j,nBlockRows,nBlockCols,bm],
 102     nBlockRows : length(mPrime),
 103     nBlockCols : length(first(mPrime)),
 104     bm : makelist(makelist(mPrime[i][j]*m,j,1,nBlockCols),i,1,nBlockRows),
 105     return(blockMatrix(bm))
 106     ); /* end proc */
 107
 108
 109
 110 /* Matrix of signs for the case Sign(P,Z) = {-1,1} */
 111 matrixOfSignsForOnePolyNonZero :
 112   [[1,1],[1,-1]];
 113
 114 /* Matrix of signs for the case Sign(P,Z) = {0,1}  */
 115 matrixOfSignsForOnePolyNonNegative :
 116   [[1,1],[0,1]];
 117
 118 /* Matrix of signs for the case Sign(P,Z) = {0,-1} */
 119 matrixOfSignsForOnePolyNonPositive :
 120   [[1,1],[0,-1]];
 121
 122
 123 /* Matrix of signs for one polynomial */
 124 matrixOfSignsForOnePolyAllSigns :
 125   [[1,1,1],[0,1,-1],[0,1,1]];
 126
 127 /* Matrix of signs of P^A */
 128 naiveMatrixOfSigns(nPols) :=
 129   if nPols = 1 then
 130     matrixOfSignsForOnePolyAllSigns
 131   else
 132     tensorProduct(naiveMatrixOfSigns(nPols-1),
 133                   matrixOfSignsForOnePolyAllSigns);
 134
 135
 136
 137 /* It produces the product coefficient terms  */
 138 /* i.e. the terms of the coefficients of the products */
 139 composeVector(lhs,rhs) :=
 140   block([lhsLen,rhsLen,i,j,prodRes,listRes],
 141     lhsLen : length(lhs),
 142     rhsLen : length(rhs),
 143     listRes : [],
 144     for i : 1 thru rhsLen do
 145       (
 146       for j : 1 thru lhsLen do
 147          (
 148          prodRes : lhs[j]*rhs[i],
 149          listRes : endcons(prodRes,listRes)
 150          ) /* end for */
 151       ), /* end for */
 152     return(listRes)
 153     ); /* end proc */
 154
 155
 156
 157 /* Auxiliary function */
 158 tarskiQueryComposition(polList) :=
 159   block([i,j,res,prodRes,card],
 160    if length(polList) = 1 then
 161      [1,polList[1],polList[1]^2]
 162    else
 163      return(composeVector([1,first(polList),first(polList)^2],
 164                             tarskiQueryComposition(rest(polList)))
 165                           )
 166    ); /* end proc */
 167
 168
 169 /* Tarski Query Vector */
 170 /* In the book: SQ(P^a,Z) */
 171 /* tarskiQueryAlg can be for instance : sSubResTarskiQuery(Q,P,var) */
 172 naiveTarskiQueryVector(polList,ptSet,tarskiQueryAlg,var) :=
 173   block([sQConv, res],
 174
 175     sQConv : tarskiQueryComposition(reverse(polList)),
 176
 177     res : map(lambda([y],tarskiQueryAlg(expand(y),expand(ptSet),var)),
 178          sQConv),
 179
 180     return(res)
 181     ); /* end proc */
 182
 183
 184
 185 /* Naive Sign Determination */
 186
 187 naiveSignDetermination(polList,ptSet,sqAlg,var) :=
 188   naiveSignConv(solveSys(naiveMatrixOfSigns(length(polList)),
 189                 naiveTarskiQueryVector(polList,ptSet,sqAlg,var),LINEAR_SOLVER));
 190
 191
 192
 193 collectSimilarAux(item,checkList,newCheckList,res) :=
 194  if checkList=[] then
 195    [res,newCheckList]
 196  else
 197    if rest(item)=rest(first(checkList)) then
 198      collectSimilarAux(item,rest(checkList),
 199                        newCheckList,endcons(first(checkList),res))
 200    else
 201      collectSimilarAux(item,rest(checkList),
 202                        endcons(first(checkList),newCheckList),res);
 203
 204 collectSimilar(item,checkList) :=
 205   collectSimilarAux(item,checkList,[],[]);
 206
 207 extractExt(extSignDet) :=
 208   block([signList,listLen,item,simRes,simList,ext2,ext3],
 209    signList : extSignDet,
 210    ext2 : [],
 211    ext3 : [],
 212    while not(signList=[]) do
 213       (
 214
 215       listLen : length(signList),
 216       item : first(signList),
 217       signList : rest(signList),
 218       simRes : collectSimilar(item,signList),
 219       simList : simRes[1],
 220       signList : simRes[2],
 221
 222       if length(simList)>=1 then
 223         ext2 : endcons(rest(item),ext2),
 224       if length(simList)=2 then
 225         ext3 : endcons(rest(item),ext3)
 226
 227      ), /* end while */
 228    return([ext2,ext3])
 229    ); /* end proc */
 230
 231 signs2exps(signCondList) := signCondListConv(signCondList);
 232
 233 exps2signs(exps) := signCondListBackConv(exps);
 234
 235 signConv(singleSign) :=
 236   if singleSign = -1 then
 237     2
 238   else
 239     singleSign;
 240
 241 signCondConv(signList) :=
 242   map(signConv,signList);
 243
 244 signCondListConv(signCondList) :=
 245   map(signCondConv,signCondList);
 246
 247 signBackConv(singleSign) :=
 248   if singleSign = 2 then
 249     -1
 250   else
 251     singleSign;
 252
 253 signCondBackConv(signList) :=
 254   map(signBackConv,signList);
 255
 256 signCondListBackConv(signCondList) :=
 257   map(signCondBackConv,signCondList);
 258
 259 genExp(base,expo) :=
 260   if base=0 and expo = 0 then
 261     1
 262   else
 263     base^expo;
 264
 265 signAt(signCond,adaExp) :=
 266   product(genExp(signCond[i],adaExp[i]),i,1,length(signCond));
 267
 268 matrixOfSigns(adaExpList,signCondList) :=
 269  block([i,j,nSigns,nExp,newRow,res],
 270    nSigns : length(signCondList),
 271    nExp : length(adaExpList),
 272    signCondList : signCondListBackConv(
 273                       sort(signCondListConv(signCondList))),
 274    adaExpList : sort(adaExpList),
 275
 276    res : [],
 277    for i : 1 thru nExp do
 278      (
 279      newRow : makelist(signAt(signCondList[j],adaExpList[i]),j,1,nSigns),
 280
 281      res : endcons(newRow,res)
 282      ), /* end for */
 283    return(res)
 284    ); /* end proc */
 285
 286 tarskiQueryVector(adaExpList,polList,var,sqAlg,ptSet) :=
 287   block([i,j,y],
 288   adaExpList : sort(adaExpList),
 289   return(
 290            map(lambda([y],sqAlg(expand(y),expand(ptSet),var)),
 291            makelist(
 292            product(polList[i]^adaExpList[j][i],
 293            i,1,length(polList)),j,1,length(adaExpList))))
 294   ); /* end proc */
 295
 296 extendByOneExponent(newExp,oldExpList) :=
 297   if oldExpList = [] then
 298     []
 299   else
 300     sort(makelist(cons(newExp,oldExpList[i]),i,1,length(oldExpList)));
 301
 302
 303 extendExponents(newExpList,oldExpList) :=
 304   block([i,j,res],
 305   res : [],
 306   if length(newExpList) = 0 then
 307     return(oldExpList),
 308   if length(oldExpList) = 0 then
 309     return(makelist([newExpList[i]],i,1,length(newExpList))),
 310   for i : 1 thru length(newExpList) do
 311     for j : 1 thru length(oldExpList) do
 312       res : cons(cons(newExpList[i],oldExpList[j]),res),
 313   return(sort(res))
 314   ); /* end proc */
 315
 316
 317 extendSigns(newSgnList,oldSgnList) :=
 318   exps2signs(extendExponents(signCondConv(newSgnList),signs2exps(oldSgnList)));
 319
 320
 321 extendAdaptedFamily(extSignDet,ada) :=
 322   block([ext2,ext3,r1,r2,r3,m2,m3,ada2,ada3,extRes,trN,i,j,indList,trM],
 323   if extSignDet=[] then
 324      return(ada),
 325   if length(extSignDet[1])=1 then
 326      if length(extSignDet)= 1 then
 327        return([[0]])
 328      else
 329         if length(extSignDet)=2 then
 330           return([[0],[1]])
 331         else
 332           return([[0],[1],[2]]),
 333   extRes : extractExt(extSignDet),
 334
 335
 336   ext2 : extRes[1],
 337   ext3 : extRes[2],
 338
 339   r2 : length(ext2),
 340   r3 : length(ext3),
 341   m2 : matrixOfSigns(ada,ext2),
 342   trM : gaussElim(m2),
 343
 344   i : 1,
 345   indList : [],
 346   while length(indList) < r2 /* and i<= length(trM[1]) */ do
 347     (
 348     if not(member(i-1,trM[3])) then
 349       indList : endcons(i,indList),
 350     i : i +1
 351     ), /* end while */
 352
 353   ada2 : makelist(ada[indList[j]],j,1,length(indList)), /* it was r2 before */
 354   m3 : matrixOfSigns(ada,ext3),
 355   trM : gaussElim(m3),
 356
 357   i : 1,
 358   indList : [],
 359   while length(indList)<r3 /* and i<=length(trM[1]) */ do /* it was r2 before */
 360     (
 361     if not(member(i-1,trM[3])) then
 362       indList : endcons(i,indList),
 363     i : i +1
 364     ), /* end while */
 365
 366  ada3 : makelist(ada[indList[j]],j,1,length(indList)), /* it was r3 before */
 367
 368
 369   return(sort(append(extendByOneExponent(0,ada),
 370                      extendByOneExponent(1,ada2),
 371                      extendByOneExponent(2,ada3))))
 372
 373   ); /* end proc */
 374
 375
 376 smartSignDetermination(polList,ptSet,sqAlg,var) :=
 377   block([sqOne,sqPexp1,sqPexp2,cVect,i,j,
 378          nPols,adaOnePoly,sgnsOnePoly,sgns,oldSgns,allSgns,
 379          nSigns,solOne,mtx,mtxOnePoly,ada,
 380          extMtx,extSysSol,sqVect,solExt,allSigns],
 381
 382    sqOne : sqAlg(1,ptSet,var),
 383    if sqOne = 0 then
 384      return([]),
 385    nPols : length(polList),
 386    ada : [],
 387    sgns : [], oldSgns : [],
 388    for i : nPols thru 1 step -1 do
 389       (
 390       sqPexp1 : sqAlg(polList[i],ptSet,var),
 391
 392       sqPexp2 : sqAlg(expand(polList[i]^2),ptSet,var),
 393       solOne : solveSys(matrixOfSignsForOnePolyAllSigns,
 394                [sqOne,sqPexp1,sqPexp2],LINEAR_SOLVER),
 395
 396       if solOne[1]>0 then
 397          sgnsOnePoly : [0]
 398       else
 399          sgnsOnePoly : [],
 400       if solOne[2]>0 then
 401          sgnsOnePoly : endcons(1,sgnsOnePoly),
 402       if solOne[3]>0 then
 403          sgnsOnePoly : endcons(-1,sgnsOnePoly),
 404       nSigns : length(sgnsOnePoly),
 405       adaOnePoly : makelist(j,j,0,nSigns-1),
 406
 407       mtxOnePoly : matrixOfSigns(makelist([adaOnePoly[j]],
 408                                           j,1,length(adaOnePoly)),
 409                                  makelist([sgnsOnePoly[j]],
 410                                           j,1,length(sgnsOnePoly))),
 411
 412       if i < nPols then
 413          (
 414          sqVect : tarskiQueryVector(extendExponents(adaOnePoly,ada),
 415                                    makelist(polList[j],j,i,nPols),
 416                                    var,sqAlg,ptSet),
 417
 418          extMtx : tensorProduct(mtx,mtxOnePoly),
 419
 420
 421          solExt : solveSys(extMtx,sqVect,LINEAR_SOLVER),
 422
 423          allSigns : extendSigns(sgnsOnePoly,sgns),
 424          sgns : [],
 425          for j : 1 thru length(allSigns) do
 426            if not(solExt[j] = 0) then
 427               sgns : endcons(allSigns[j],sgns),
 428          ada : extendAdaptedFamily(sgns,ada),
 429
 430
 431          mtx : matrixOfSigns(ada,sgns)
 432
 433          ) /* end if */
 434       else
 435          (
 436          sgns : makelist([sgnsOnePoly[j]],j,1,length(sgnsOnePoly)),
 437          mtx : mtxOnePoly,
 438          ada : makelist([adaOnePoly[j]],j,1,length(adaOnePoly))
 439          ) /* end else */
 440       ), /* end for */
 441    return(sgns)
 442    ); /* end proc */
 443
 444
 445 log3aux(n,c) :=
 446   if 3^c = n then
 447     c
 448   else
 449     if c>n then
 450       false
 451     else
 452       log3aux(n,c+1);
 453
 454 log3(n) := log3aux(n,0);
 455
 456 naiveSignConv(naiveSignList) :=
 457   block([res,i,j,nPols,tsign,tpol,item,propag,copyItem,prInd],
 458    nPols : log3(length(naiveSignList)),
 459    res : [],
 460    tpol : nPols,
 461    tsign : 0,
 462    item : makelist(0,i,1,nPols),
 463    for i : 1 thru length(naiveSignList) do
 464      (
 465      if tsign = 3 then
 466        (
 467        tsign : 0,
 468        propag : false,
 469        prInd : nPols-1,
 470        while not(propag) do
 471          (
 472          if item[prInd] < 2 then
 473            (
 474            item[prInd] : item[prInd]+1,
 475            propag : true
 476            ) /* end if */
 477          else
 478            (
 479            item[prInd] : 0,
 480            prInd : prInd-1
 481            ) /* end else */
 482          ) /* end while */
 483      ), /* end if */
 484      if naiveSignList[i] > 0 then
 485         (
 486         copyItem : copylist(item),
 487
 488         res : cons(copyItem,res)
 489         ), /* end if */
 490      tsign : tsign+1,
 491
 492
 493      item[nPols] : tsign
 494      ), /* end for */
 495    return(exps2signs(res))
 496    ); /* end proc */
 497
 498
 499
 500 adaptedFamily(signDet) :=
 501   block([partSignDet,i,partAda],
 502   numPols : length(signDet[1]),
 503
 504   for i : numPols thru 1 step - 1 do
 505     (
 506     partSignDet : listify(setify(map(lambda([x],rest(x,i-1)),signDet))),
 507
 508     if i = numPols then
 509       partAda : makelist([i],i,0,length(partSignDet)-1)
 510     else
 511       partAda : extendAdaptedFamily(partSignDet,partAda)
 512
 513     ), /* end for */
 514   return(partAda)
 515   );
 516
 517 /* -------------------------------------------------- */
 518 /* Thom Encoding related-functions                    */
 519
 520
 521 /* Thom Encoding sorting */
 522
 523 thomLess(lhs,rhs) :=
 524   thomLessAux(rest(lhs,-1),rest(rhs,-1),last(lhs));
 525
 526
 527 thomLessAux(lhs,rhs,lastSgn) :=
 528   if lhs = [] and rhs = [] then
 529      false
 530   else
 531     if last(lhs)=last(rhs) then
 532       thomLessAux(rest(lhs,-1),rest(rhs,-1),last(lhs))
 533     else
 534       if lastSgn=1 then
 535          if last(lhs)<last(rhs) then
 536             true
 537          else
 538             false
 539       else
 540          if last(lhs)<last(rhs) then
 541             false
 542          else
 543             true;
 544
 545 thomSort(thomList) :=
 546   sort(thomList,thomLess);
 547
 548
 549 /* Thom Encoding */
 550 thomEncoding(pol,var) :=
 551    thomSort(THOM_SIGN_DET_ALG(rest(derSeq(pol,var)),pol,
 552                               sSubResTarskiQuery,var));
 553
 554
 555
 556 extThomEncoding(p,q,var) :=
 557   block([prepro],
 558     prepro : append(derSeq(q,var),rest(derSeq(p,var))),
 559
 560   return(THOM_SIGN_DET_ALG(prepro,p,sSubResTarskiQuery,var))
 561           ); /* end proc */
 562
 563
 564
 565 /* Merge two thom-sorted encodings and remove duplicates */
 566 thomMerge(lhs,rhs) :=
 567   mergeSorted(lhs,rhs,thomLess);
 568
 569
 570 extThomMerge(lhs,degLhs,rhs,degRhs) :=
 571   block([i,j,k,lhsLen,rhsLen,res,lhsComm,rhsComm,
 572          lhsSign,lhsSignAtRhs,rhsSign,rhsSignAtLhs,numOfSigns],
 573   lhsLen : length(lhs),
 574   rhsLen : length(rhs),
 575   numOfSigns : degLhs + degRhs+1,
 576   i : 1,
 577   j : 1,
 578   res : [],
 579   while (i<= lhsLen) and (j<= rhsLen) do
 580     (
 581     lhsSign : makelist(lhs[i][k],k,degRhs+2,numOfSigns),
 582     lhsComm : lhs[i][1],
 583     lhsSignAtRhs : makelist(lhs[i][k],k,2,degRhs+1),
 584
 585     rhsSign : makelist(rhs[j][k],k,degLhs+2,numOfSigns),
 586     rhsComm : rhs[j][1],
 587     rhsSignAtLhs : makelist(rhs[j][k],k,2,degLhs+1),
 588
 589     if lhsComm=0 then
 590        (
 591         while(not(rhs[j][1]=0)) do
 592           (
 593           res : endcons([2,rhsSign],res),
 594           j : j + 1,
 595           rhsSign : makelist(rhs[j][k],k,degLhs+2,numOfSigns)
 596           ), /* end while */
 597         res : endcons([0,lhsSign,rhsSign],res),
 598         i : i + 1,
 599         j : j + 1
 600         ) /* end if */
 601     else
 602       if rhsComm=0 then
 603         (
 604         while(not(lhs[i][1]=0)) do
 605           (
 606           res : endcons([1,lhsSign],res),
 607           i : i + 1,
 608           lhsSign : makelist(lhs[i][k],k,degRhs+2,numOfSigns)
 609           ), /* end while */
 610         res : endcons([0,lhsSign,rhsSign],res),
 611         i : i + 1,
 612         j : j + 1
 613         ) /* end if */
 614
 615       else
 616        if not(lhsSign=rhsSignAtLhs) then
 617          if thomLess(lhsSign,rhsSignAtLhs) then
 618            (
 619
 620            res : endcons([1,lhsSign],res),
 621            i : i + 1
 622            ) /* end if */
 623          else
 624            (
 625            res : endcons([2,rhsSign],res),
 626            j : j + 1
 627            ) /* end else */
 628        else
 629          if thomLess(rhsSign,lhsSignAtRhs) then
 630            (
 631
 632            res : endcons([2,rhsSign],res),
 633            j : j + 1
 634            ) /* end if */
 635          else
 636            (
 637
 638            res : endcons([1,lhsSign],res),
 639            i : i + 1
 640            ) /* end else */
 641
 642      ), /* end while */
 643
 644   if i<= lhsLen then
 645     res : append(res, makelist([1,
 646                          makelist(lhs[k][j],j,degRhs+2,numOfSigns)],
 647                       k,i,lhsLen))
 648   else
 649      if j <= rhsLen then
 650         res : append(res, makelist([2,
 651                              makelist(rhs[k][j],j,degLhs+2,numOfSigns)],
 652                           k,j,rhsLen)),
 653 return(res)
 654  ); /* end proc */
 655
 656
 657 /* It verifies whether a root of p is a common root */
 658 /* with the polynomial used for the extended Thom encoding */
 659 thomComm(pol,var,extThomEnc) :=
 660   if(extThomEnc[degree(pol,var)+1]=0) then
 661     true
 662   else
 663     false;
 664
 665
 666 /* Root Comparison */
 667 thomCompare(p,q,var) :=
 668   block([pSign,qSign, exP, exQ],
 669
 670   exP: expandIf(p),
 671   exQ: expandIf(q),
 672
 673   pSign : thomSort(extThomEncoding(exP,exQ,var)),
 674
 675   qSign : thomSort(extThomEncoding(exQ,exP,var)),
 676
 677   return(extThomMerge(pSign,degree(exP,var),qSign,degree(exQ,var)))
 678
 679   ); /* end proc */
 680
 681 thomComparePrettyPrint(thomCompareOutput) :=
 682     map(lambda([x],
 683         if x[1]=1 then
 684            ["1st",x[2]]
 685         else
 686            if x[1]=2 then
 687              ["2nd",x[2]]
 688            else
 689              ["Com",x[2]]),thomCompareOutput);
 690
 691
 692 /* For every Thom encoding of P it computes the signs of q */
 693 /*currently buggy*/
 694 thomSignsAtRoots(thomP,p,q,var) :=
 695   block([adaQAux,qSigns,adaQ,adaDerPmP,
 696          taQVect,mtxQ,mtxDerPmP,mtxExt,solExt,allSigns,
 697          res,i],
 698     qSigns : smartSignDetermination([q],p,sSubResTarskiQuery,var),
 699
 700     adaQ : adaptedFamily(qSignsAux),
 701     adaQAux : makelist(first(adaQ[i]),i,1,length(adaQ)),
 702
 703     adaDerPmP : adaptedFamily(thomP), /* Ada(Der(P)\{P}) */
 704
 705     taQVect : tarskiQueryVector(extendExponents(adaQAux,adaDerPmP),
 706                           cons(q,makelist(diff(p,var,i),i,1,degree(p,var))),
 707                           var,sSubResTarskiQuery,p),
 708
 709     mtxQ : matrixOfSigns(adaQ,qSigns),
 710
 711
 712     mtxDerPmP : matrixOfSigns(adaDerPmP,thomP),
 713
 714
 715     mtxExt : tensorProduct(mtxQ,mtxDerPmP),
 716
 717     solExt : solveSys(mtxExt,taQVect,LINEAR_SOLVER),
 718
 719
 720     allSigns : extendSigns(qSigns,thomP),
 721
 722     res : [],
 723     for i : 1 thru length(allSigns) do
 724        if not(solExt[i]=0) then
 725          res : endcons([rest(allSigns[i]),first(allSigns[i])],res),
 726
 727
 728     return(res)
 729     ); /* end proc */
 730
 731
 732 /* It computes the thom Encodings of p with */
 733 /* with the sign of q at each Thom Encoding */
 734 thomEncodingWithSigns(p,q,var) :=
 735   block([signDetRes],
 736
 737     signDetRes : smartSignDetermination(append(rest(derSeq(p,var)),[q]),p,
 738                                    sSubResTarskiQuery,
 739                                    var),
 740
 741     return(map(lambda([x],[rest(x,-1),last(x)]),
 742               signDetRes))
 743
 744     ); /* end proc */