share/contrib/sarag/rootCounting.mac

   1 /* ---------------------------------------------------------------- */
   2 /* SARAG - Root Counting                                            */
   3 /* by Fabrizio Caruso  modified by Alexandre Le meur and Marie-Françoise Roy              */
   4
   5
   6 /* Verbosity levels */
   7 NONVERBOSE : 0;
   8 SUMMARY : 1;
   9 NORMAL : 2;
  10 VERY : 3;
  11 EXTRA : 4;
  12
  13 padWithZeros(seq,size) :=
  14   append(seq,makelist(0,i,1,size-length(seq)));
  15
  16 /* signed pseudo-remainder of p divised by q*/
  17 sPsRem(p,q,var):=
  18 if degree(p,var)<degree(q,var) then p
  19 else (if remainder(degree(p,var)-degree(q,var),2)=1 then
  20 remainder(leadCoeff(q,var)**(degree(p,var)-degree(q,var)+1)*p,q,var)
  21 else remainder(leadCoeff(q,var)**(degree(p,var)-degree(q,var)+2)*p,q,var)
  22 ) /* end else */
  23 ;
  24
  25 /* Signed Remainder Sequence */
  26 sRem(a,b,var) :=
  27   block([aa,bb,rOld,rNew,quo,res,seq,srsRes],
  28
  29   aa : expandIf(a),
  30   bb : expandIf(b),
  31
  32   rOld : aa,
  33   rNew : bb,
  34   srsRes : [aa],
  35
  36   while not(rNew = 0) do
  37     (
  38     srsRes : endcons(rNew,srsRes),
  39
  40     res : remainder(rOld,rNew,var),
  41
  42     rOld : rNew,
  43     rNew : -res
  44     ), /* end while */
  45
  46   return(srsRes)
  47   ); /* end proc */
  48
  49
  50
  51 /* Extended Signed Remainder Sequence */
  52 sRemExt(a,b,var) :=
  53   block([rOld,rNew,uOld,uNew,vOld,vNew,rRes, uRes, vRes,quo,rList,uList,vList,rAux,uAux,vAux],
  54
  55   a : expandIf(a),
  56   b : expandIf(b),
  57
  58   rOld : a,
  59   rNew : b,
  60   uOld : 1, vOld : 0,
  61   uNew : 0, vNew : 1,
  62   rList : [a],
  63   uList : [uOld],
  64   vList : [vOld],
  65
  66   while not(rNew = 0 ) do
  67      (
  68      rList : endcons(rNew,rList),
  69      uList : endcons(uNew,uList),
  70      vList : endcons(vNew,vList),
  71
  72      rRes : divide(rOld,rNew,var),
  73
  74      quo : first(rRes),
  75
  76      rAux : rNew,
  77      uAux : uNew,
  78      vAux : vNew,
  79
  80      uNew : ratexpand(-uOld + quo*uNew),
  81      vNew : ratexpand(-vOld + quo*vNew),
  82
  83      rOld : rAux,
  84      uOld : uAux,
  85      vOld : vAux,
  86      rNew : -second(rRes)
  87
  88      ),/* end while */
  89   return([rList,endcons(uNew,uList),endcons(vNew,vList)])
  90   );/* end proc */
  91
  92
  93 /* It counts the sign change of a determinant after */
  94 /* i consecutive row changes */
  95 revRowsCount(i) := (-1)^(i*(i-1)/2);
  96
  97 /* ------------------------------------------------------- */
  98 /* Part concerning the computation of signed subresultants */
  99
 100
 101 sSubResPol(p,q,var) :=
 102   block([pp,qq,degA,degB,i,j,k,h,lcA,lcB,delta,
 103          SR,s,sOld,res],
 104
 105     pp : expandIf(p),
 106     qq : expandIf(q),
 107
 108     degA : degree(pp,var),
 109     degB : degree(qq,var),
 110     if degA<=degB then
 111       (
 112       print("sSubRes) Error: inconsistent degrees"),
 113       return(false)
 114       ), /* end if */
 115
 116     SR : make_array('any,degA+1),
 117     s : make_array('any,degA+1),
 118     sOld : make_array('any,degA+1),
 119
 120
 121     lcA : leadCoeff(pp,var),
 122     lcB : leadCoeff(qq,var),
 123
 124     SR[degA] : pp,
 125   s[degA] : 1,
 126
 127     sOld[degA] : 1,
 128
 129     SR[degA-1] : qq,
 130     sOld[degA-1] : lcB,
 131
 132     i : degA + 1,
 133     j : degA,
 134
 135     while j>=1 and not(SR[j-1]=0) do
 136       (
 137
 138       k : degree(SR[j-1],var),
 139
 140       if k = j-1 then
 141         (
 142         s[j-1] : sOld[j-1],
 143
 144         if k>=1 then
 145            SR[k-1] : NORM_ALGORITHM(-remainder(s[j-1]^2*SR[i-1],SR[j-1],var)/
 146                          (s[j]*sOld[i-1]))
 147         ), /* end if */
 148       if k < j-1 then
 149         (
 150         /* Computation of s[k] */
 151         for delta : 1 thru min(j-k-1,j-1) do
 152            sOld[j-delta-1] : (-1)^delta * (sOld[j-1]*sOld[j-delta])/s[j],
 153         s[k] : sOld[k],
 154         s[j-1] : 0,
 155         SR[k] : NORM_ALGORITHM(s[k]*SR[j-1]/sOld[j-1]),
 156
 157         for h : k+1 thru j-2 do
 158           (
 159           SR[h] : 0,
 160           s[h] : 0
 161           ), /* end for */
 162
 163         /* Computation of SR[k-1] */
 164         /* ratexpand or expand */
 165         if k>= 1 then
 166           SR[k-1] : NORM_ALGORITHM(-remainder(sOld[j-1]*s[k]*SR[i-1],
 167                                              SR[j-1],var)/
 168                          (s[j]*sOld[i-1]))
 169         ), /* end if */
 170
 171       if k>= 1 then
 172          sOld[k-1] : leadCoeff(SR[k-1],var),
 173       i : j,
 174       j : k
 175
 176       ), /*end while */
 177
 178
 179 res : padWithZeros(makelist(SR[degA-i],i,0,degA-j),degA+1),
 180
 181      return(res)
 182  ); /* end proc */
 183
 184
 185
 186 sSubResCoeff(p,q,var) :=
 187 block([s,pp,qq],
 188 pp:expand(p),
 189 qq:expand(q),
 190 s :  sSubResPol(pp,qq,var),
 191 return(
 192 makelist(coeff(s[j],var,degree(pp,var)+1-j),j,1,degree(pp,var)+1)
 193 )
 194 );
 195
 196
 197
 198 /* Extended Signed Subresultant Sequence      */
 199 /* Algorithm 8.75                             */
 200 /* (A variation of this algorithm is used in  */
 201 /* Algorithm 9.47 and Algorithm 10.17         */
 202
 203
 204 /* Extended Signed Subresultant Sequence*/
 205
 206 sSubResExt(a,b,var) :=
 207   sSubResExtVerbose(a,b,var,DEFAULT_VERBOSITY);
 208
 209 sSubResExtVerbose(a,b,var,verbosity) :=
 210    block([q,degA,degB,i,j,k,lcA,lcB,delta,
 211           SR:make_array('any),s:make_array('any),sOld:make_array('any),
 212           u:make_array('any),v:make_array('any)],
 213
 214    a : expandIf(a),
 215    b : expandIf(b),
 216
 217    degA : degree(a,var),
 218    degB : degree(b,var),
 219    if degA<=degB then
 220       (
 221       print("sSubResExt) Error: inconsistent degrees"),
 222       return(false)
 223       ),
 224
 225
 226    SR : make_array('any,degA+1),
 227    s : make_array('any,degA+1),
 228    sOld : make_array('any,degA+1),
 229    u : make_array('any,degA+2),
 230    v : make_array('any,degA+2),
 231
 232
 233    lcA : leadCoeff(a,var),
 234    lcB : leadCoeff(b,var),
 235
 236    SR[degA] : a,
 237
 238
 239       s[degA] : 1,
 240
 241    sOld[degA] : s[degA],
 242    SR[degA-1] : b,
 243
 244    sOld[degA-1] : lcB,
 245
 246    u[degA+1] : 1,
 247    v[degA+1] : 0,
 248    u[degA-1+1] : 0,
 249    v[degA-1+1] : 1,
 250
 251    i : degA + 1,
 252    j : degA,
 253
 254    while j>=1 and not(SR[j-1]=0) do
 255      (
 256      k : degree(SR[j-1],var),
 257      if verbosity>=EXTRA then
 258         (
 259         print("sSubResExt) j : ", j),
 260         print("sSubResExt) k : ", k)
 261         ),
 262      if k = j-1 then
 263        (
 264
 265        s[j-1] : sOld[j-1],
 266        q : quotient(s[j-1]^2 * SR[i-1], SR[j-1],var),
 267        if verbosity >= VERY then
 268          print("sSubResExt) q : ", q),
 269
 270        if k>=1 then
 271          SR[k-1] : NORM_ALGORITHM((-s[j-1]^2 *
 272                                  SR[i-1] + q * SR[j-1])/(s[j]*sOld[i-1])),
 273
 274
 275        u[k-1+1] : NORM_ALGORITHM((-s[j-1]^2 *
 276                                 u[i-1+1] + q * u[j-1+1])/(s[j]*sOld[i-1])),
 277        v[k-1+1] : NORM_ALGORITHM((-s[j-1]^2 *
 278                                 v[i-1+1] + q * v[j-1+1])/(s[j]*sOld[i-1]))
 279
 280        ),
 281
 282      if k<j-1 then
 283        (
 284        s[j-1] : 0,
 285        /*  Computation of s[k] */
 286        for delta : 1 thru min(j-k-1,j-1) do
 287           sOld[j-delta-1] : (-1)^delta * (sOld[j-1]*sOld[j-delta])/s[j],
 288
 289        if k>=0 then
 290          (
 291          s[k] : sOld[k],
 292
 293          SR[k] : NORM_ALGORITHM(s[k]*SR[j-1]/sOld[j-1])
 294          ),
 295        u[k+1] : NORM_ALGORITHM(s[k]*u[j-1+1]/sOld[j-1]),
 296        v[k+1] : NORM_ALGORITHM(s[k]*v[j-1+1]/sOld[j-1]),
 297
 298
 299        for h : k+1 thru j-2 do
 300           (
 301           SR[h] : 0,
 302           u[h+1] : 0,
 303           v[h+1] : 0,
 304           s[h] : 0
 305           ),
 306
 307        /* Computation of SR[k-1], u[k-1], v[k-1] */
 308        q : quotient(s[k]*sOld[j-1]*SR[i-1],SR[j-1],var),
 309        if k>=1 then
 310          SR[k-1] : NORM_ALGORITHM((-sOld[j-1] * s[k] * SR[i-1] +
 311                            q * SR[j-1])/(s[j]*sOld[i-1])),
 312
 313        u[k-1+1] : NORM_ALGORITHM((-sOld[j-1] * s[k] * u[i-1+1] +
 314                            q * u[j-1+1])/(s[j]*sOld[i-1])),
 315        v[k-1+1] : NORM_ALGORITHM((-sOld[j-1] * s[k] * v[i-1+1] +
 316                           q * v[j-1+1])/(s[j]*sOld[i-1]))
 317
 318        ),
 319
 320      if k>=1 then
 321         sOld[k-1] : leadCoeff(SR[k-1],var),
 322      i : j,
 323      j : k,
 324
 325      if verbosity>= NORMAL and k>= 1 then
 326         print("sSubResExt) SR[", k-1,"] : ", SR[k-1])
 327      ),
 328 for h : 0 thru j-2 do
 329    (
 330    SR[h] : 0,
 331    s[h] : 0,
 332    u[h+1] : 0,
 333    v[h+1] : 0
 334    ),
 335
 336
 337   s[degA]:lcA,
 338
 339 return(
 340        [padWithZeros(makelist(SR[degA-i],i,0,degA-j),degA+1),
 341         makelist(u[degA-i+1],i,0,degA-j+1),
 342         makelist(v[degA-i+1],i,0,degA-j+1)])
 343 );
 344
 345
 346 /* last non zero element of a sequence */
 347
 348 lastNonZero(seq) :=
 349 block([i],
 350   for i:1 thru length(seq) do
 351     if not(seq[length(seq)-i+1] = 0) then
 352       return(seq[length(seq)-i+1])
 353 );  /* end proc */
 354
 355 /* index of last non zero element of a sequence */
 356 lastNonZeroIndex(seq) :=
 357 block([i],
 358   for i:1 thru length(seq) do
 359     if not(seq[length(seq)-i+1] = 0) then
 360       return(length(seq)-i+1)
 361 );  /* end proc */
 362
 363 /* index of first non zero element of a sequence */
 364 firstNonZeroIndex(seq) :=
 365 block([i],
 366   for i:1 thru length(seq) do
 367     if not(seq[i] = 0) then
 368       return(i)
 369 );  /* end proc */
 370
 371 /* Gcd and Gcd-Free part computation */
 372 /* It outputs : */
 373 /* gcd(p,q) and p/gcd(p,q) */
 374
 375 gcdFreePart(p,q,var) :=
 376 block([r],
 377 r:sSubResExt(ratexpand(p),sPsRem(ratexpand(q),ratexpand(p),var),var),
 378 return ([lastNonZero(r[1]),lastNonZero(r[3])])
 379 ); /* end proc */
 380
 381
 382 gcdFreePartWithZ(p,q,var) :=
 383 block([r,u,v,c],
 384 r:sSubResExt(p,sPsRem(ratexpand(q),ratexpand(p),var),var),
 385 c:leadCoeff(expand(p),var),
 386 v:lastNonZero(r[3]),
 387 return ([ratexpand(c*lastNonZero(r[1])/leadCoeff(lastNonZero(r[1]),var)),ratexpand(c*v/leadCoeff(v,var))])
 388 ); /*  end proc */
 389
 390
 391
 392 sSubResCoeffLast(p,q,var) :=
 393   last(sSubResCoeff(p,q,var));
 394
 395 sylvesterResultant(p,q,var) :=
 396   block([aux,pp ,qq],
 397     pp : expandIf(p),
 398     qq : expandIf(q),
 399     if degree(pp,var)>degree(qq,var) then
 400       return(epsilon(degree(pp,var))*sSubResCoeffLast(pp,qq,var))
 401     else
 402       if degree(pp,var)<degree(qq,var) then
 403         return((-1)^(degree(pp,var)*degree(pp,var))*sylvesterResultant(qq,pp,var))
 404       else
 405         (
 406         aux : expand(leadCoeff(pp,var)*qq-leadCoeff(qq,var)*pp),
 407
 408         return(expand(epsilon(degree(pp,var))*sSubResCoeffLast(pp,aux,var)/leadCoeff(pp,var)^degree(aux,var)))
 409         ) /* end else */
 410       ); /* end proc */
 411
 412
 413
 414 subDiscr(p,var) :=
 415   block([pcRes,pRes,cRes],
 416     pcRes : sSubRes(p,diff(p,var),var),
 417     pRes : first(pcRes/leadCoeff(expand(p),var)),
 418     pRes : cons(first(pRes),expand(rest(pRes)/leadCoeff(expand(p),var))),
 419     cRes : second(pcRes),
 420     cRes :  cons(first(cRes),expand(rest(cRes)/leadCoeff(expand(p),var))),
 421     return([pRes,cRes])
 422     ); /* end proc */
 423
 424 subDiscrCoeff(p,var):=
 425   block([cRes],
 426   cRes : sSubResCoeff(p,diff(p,var),var),
 427   cRes : cons(first(cRes)/leadCoeff(expand(p),var),expand(rest(cRes)/leadCoeff(expand(p),var))),
 428   return(cRes)
 429   ); /* end proc */
 430
 431 discriminant(p,var):=
 432 subDiscrCoeff(p,var)[degree(expand(p),var)+1];
 433
 434
 435 /* --------------------------------------------- */
 436 /* Part concerning the signed remainder sequence */
 437
 438 /* In the following S(P,Q) is the signed remainder sequence of P and Q */
 439
 440 /* Number of sign changes (notation "V") */
 441 /* of sequence containing non-zero elements */
 442
 443
 444 nonZeroSignChanges(seq) :=
 445 block([variations,i],
 446  variations:0,
 447  if length(seq) > 1 then
 448    for i:1 thru length(seq)-1 do
 449      if seq[i]*seq[i+1] < 0 then
 450        variations:variations+1,
 451 return(variations)
 452 ); /* end proc */
 453
 454 signChanges(seq) :=
 455   nonZeroSignChanges(trimZeros(seq));
 456
 457
 458 /* Number of sign changes of a polynomial at a certain value */
 459 /* Notation "V(P;a)" */
 460 /* We use proposition 2.4 for the signs at infinities */
 461 signChangesAt(seq,var,a) :=
 462   if a = inf then
 463      signChanges(makelist(leadCoeff(seq[i],var),i,1,length(seq)))
 464   else
 465      if a = -inf then
 466        signChanges(makelist((-1)^(degree(seq[i],var))*
 467                             leadCoeff(seq[i],var),i,1,length(seq)))
 468      else
 469        signChanges(makelist(subst(a,var,seq[i]),i,1,length(seq)));
 470 /* end proc */
 471
 472
 473 /* Difference between the no. of changes a two points */
 474 /* Notation "V(P;a,b)" */
 475 signChangesDiff(seq,var,a,b) :=
 476   signChangesAt(seq,var,a) - signChangesAt(seq,var,b);
 477
 478 /* Sign changes of the i-th coefficients of a polynomial */
 479 /* Notation "V(P)" (P is a polynomial) */
 480 signChangesPolyCoeff(poly,var) :=
 481   signChanges(makelist(coeff(poly,var,i),i,0,degree(poly,var)));
 482
 483 /* Sequence of the derivatives of a polynomial */
 484 /* Notation "Der(P)" */
 485 derSeq(f,var) :=
 486   makelist(diff(f,var,i),i,0,degree(f,var));
 487
 488 /* Cauchy Index */
 489 /* Theorem 2.52 */
 490 /* V(S(P,Q);a,b) = Ind(Q/P;a,b) */
 491
 492 sRemCauchyIndexBetween(num,den,var,a,b) :=
 493   signChangesDiff(sRem(den,num,var),var,a,b);
 494
 495 sRemCauchyIndex(num,den,var) :=
 496   signChangesDiff(sRem(den,num,var),
 497                       var,-inf,+inf);
 498
 499 /* Tarski query computed by Sylvester's theorem's formula */
 500 /* Theorem 2.55 */
 501 /* V(S(P,P'Q);a,b) = SQ(Q,P;a,b) */
 502
 503 sRemTarskiQueryBetween(q,p,var,a,b) :=
 504   signChangesDiff(sRem(p,diff(p,var)*q,var),var,a,b);
 505
 506 sRemTarskiQuery(q,p,var) :=
 507   signChangesDiff(sRem(p,diff(p,var)*q,var),
 508                        var,-inf,+inf);
 509
 510
 511 /* Number of roots counted by Tarski's theorem's formula */
 512 /* Theorem 2.56 */
 513 /* V(S(P,P');a,b) */
 514 sRemNumberOfRootsBetween(p,var,a,b) :=
 515   sRemTarskiQueryBetween(1,p,var,a,b);
 516
 517 sRemNumberOfRoots(p,var) :=
 518   sRemTarskiQuery(1,p,var);
 519
 520 /* Tarski sequence */
 521 /* defined as: S(P,P') */
 522 sturmSequence(p,var) :=
 523    sRem(p,diff(p,var),var);
 524
 525
 526 /* ------------------------------------------------- */
 527 /* Part concerning signed subresultants coefficients */
 528 /* and the Cauchy Index in all R(Subsection 9.1.2)   */
 529
 530
 531 /* Notation taken from Remark 4.36 and used by Notation 9.4 */
 532 epsilon(x) :=
 533   (-1)^(x*(x-1)/2);
 534
 535 /* It counts and removes the trailing zeros out of a sequence */
 536
 537 trimZeroCount(seq,count) :=
 538 block([list,count,i],
 539 list:[],
 540 count:0,
 541 if length(seq)>0 then
 542   for i:1 thru length(seq) do
 543      if seq[i] = 0 then
 544 count:count+1
 545      else
 546 list:append(list,[seq[i]]),
 547 return([list,count])
 548 ); /* end proc */
 549
 550 /* Interface function */
 551 trimZero(seq) :=
 552   trimZeroCount(seq,0);
 553
 554
 555
 556 /* Notation 9.4 ("D(s)") */
 557 /* When there are no zeros this gives the difference between */
 558 /* the sign permanencies and the number of sign changes */
 559
 560 genPermMVar(seq) :=
 561 block([res,count,i],
 562   res:0, i:1,
 563    ( while not(i>length(seq)-1) do
 564 if not(seq[i]=0) then
 565 (count:1,
 566   (while (seq[i+count]=0 and not(i+count>length(seq)-1)) do count:count+1), /* end while */
 567     ( if  remainder(count,2)=1 then
 568 res:res+(-1)^((count)*(count-1)/2)*sgn(seq[i]*seq[i+count])),
 569 i:i+count) /* end if */
 570 else i:i+1), /* end while */
 571 return(res)
 572 ); /* end proc */
 573
 574
 575 /* Cauchy Index by signed subresultant */
 576 /* Algorithm 9.27 */
 577
 578 sSubResCauchyIndex(q,p,var) :=
 579 genPermMVar(sSubResCoeff(p,sPsRem(ratexpand(q),ratexpand(p),var),var));
 580
 581 /* Tarski Query by signed subresultant */
 582 /* Algorithm 9.28 */
 583 sSubResTarskiQuery(q,p,var) :=
 584 genPermMVar(sSubResCoeff(p,sPsRem(ratexpand(diff(p,var)*q),p,var),var));
 585
 586 sSubResNumberOfRoots(p,var) :=
 587 sSubResTarskiQuery(1,p,var);
 588
 589
 590
 591 /* ----------------------------------------------  */
 592 /* Part concerning the Bezoutian (subsection 9.1.3) */
 593
 594 /* Bezoutian (Notation 9.14) */
 595 bez(p,q,var,x,y) :=
 596    (subst(y,var,q) * subst(x,var,p) -
 597     subst(x,var,q) * subst(y,var,p))/(x-y); /* end proc */
 598
 599
 600
 601 /* ----------------------------------------------- */
 602 /* Part concerning the Cauchy Index on an interval */
 603 /* (Subsection 9.1.5)                              */
 604
 605 /* Modified number of sign changes */
 606 /* Notation 9.29, ("W(s)") */
 607 /* Modified Sign Changes: */
 608 /* counting as two sign changes the groups +,0,0,+ and -,0,0,- */
 609
 610
 611 modifiedSignChanges(seq) :=
 612 block([variations,i],
 613  variations:signChanges(seq),
 614  if length(seq) > 1 then
 615    for i:1 thru length(seq)-3 do
 616 if (seq[i+1]=0 and seq[i+2]=0 and seq[i]*seq[i+3] > 0) then
 617        variations:variations+2,
 618 return(variations)
 619 ); /* end proc */
 620
 621
 622 /* Remove identically zero elements from a list */
 623
 624
 625 removeZeros(seq) :=
 626 block([list,i],
 627 list:[],
 628 if length(seq)>0 then
 629   for i:1 thru length(seq) do
 630      if not (seq[i] = 0) then
 631 list:append(list,[seq[i]]),
 632 return(list)
 633 ); /* end proc */
 634
 635
 636 /* Modified Sign Changes at a certain value */
 637 /* Notation "W(P;a)" */
 638 modifiedSignChangesAt(seq,var,value) :=
 639    if value = inf then
 640       signChanges(makelist(leadCoeff(seq[i],var),i,1,length(seq)))
 641    else
 642      if value = -inf then
 643        signChanges(makelist((-1)^(degree(seq[i],var))*
 644                             leadCoeff(seq[i],var),i,1,length(seq)))
 645      else
 646        modifiedSignChanges(subst(value,var,removeZeros(seq))); /* end proc */
 647
 648
 649 /* Modified Sign Changes Difference */
 650 /* Notation "W(P;a,b)" */
 651 modifiedSignChangesDiff(seq,var,a,b) :=
 652    modifiedSignChangesAt(seq,var,a)-
 653    modifiedSignChangesAt(seq,var,b);
 654
 655
 656
 657 /* Cauchy Index in an interval computed by subresultants */
 658 /* Theorem 9.30 */
 659 /* W(SR(P,Q);a,b) = Ind(Q/P;a,b) */
 660
 661 sSubResCauchyIndexBetween(num,den,var,a,b) :=
 662   modifiedSignChangesDiff(sSubResPol(den,sPsRem(num,den,var),var),var,a,b);
 663
 664 /* Tarski Query by subresultants */
 665 /* Corallary 9.32 */
 666 /* W(SR(P,Remainder(P,Q));a,b) = SQ(Q,P;a,b) */
 667 sSubResTarskiQueryBetween(q,p,var,a,b) :=
 668   modifiedSignChangesDiff(sSubResPol(p,
 669                               sPsRem(ratexpand(diff(p,var)*q),p,var),var),var,a,b);
 670
 671 /* Number of roots on an interval by signed subresultants */
 672 /* Consequence of Corollary 9.32 */
 673 sSubResNumberOfRootsBetween(p,var,a,b) :=
 674   sSubResCauchyIndexBetween(diff(p,var),p,var,a,b);
 675
 676
 677 /* ------------------------------- */
 678 /* Part concerning Hankel Matrices */
 679
 680
 681 /* It computes the signature of a Hankel quadratic form */
 682 /* Algorithm 9.47 */
 683
 684
 685 hankelSignature(seq) :=
 686   block([m,n,sLen,trimRes,trimmed,count,c,p,q,t,i,SR,j,defective],
 687    sLen : length(seq),
 688    if evenp(sLen) then
 689      (
 690      print("hankelSignature) Number of elements must be odd"),
 691      return(false)
 692      ), /* end if */
 693    n : (sLen+1)/2,
 694
 695    trimRes : trimZero(seq),
 696
 697    if length(trimRes) = 0 then
 698       return(0)
 699    else
 700        c:firstNonZeroIndex(seq),
 701       (if c>=n then
 702          (if oddp(c) then
 703             return(sgn(seq[c]))
 704          else
 705             return(0)) /* end if */
 706         else m:lastNonZeroIndex(seq),
 707          (if m <= n then
 708            (if oddp(2*n-m) then
 709               return(sgn(seq[m]))
 710            else
 711               return(0)) /* end if */
 712          else
 713            (
 714            p : t^m,
 715            q : sum(seq[i]*t^(m-i),i,1,m),
 716
 717            SR : sSubResPol(p,q,t),
 718
 719            /* check of defectiveness */
 720            j : m-n,
 721            defective : true,
 722            while defective do
 723               if degree(SR[m-j+1],t) = j then
 724                  defective : false
 725               else
 726                  j : j-1,
 727
 728            return(genPermMVar(
 729                    makelist(sSubResCoeff(p,q,t)[i],
 730                      i,1,m-j+1)))
 731
 732            ) /* end else */
 733          ) /* end if */
 734
 735       ) /* end if */
 736
 737 ); /* end proc */
 738
 739
 740 /* It build a Hankel matrix out of an odd-numbered sequence */
 741 hankelMatrix(seq) :=
 742   block([n,len],
 743   len : length(seq),
 744   n : (len+1)/2,
 745   return(makelist(
 746            makelist(
 747                   seq[i+j+1],
 748                j,0,n-1),
 749            i,0,n-1)
 750          )
 751   ); /* end proc */
 752
 753
 754
 755 /* ----------------------------------------------------- */
 756 /* Part concerning complex roots with negative real part */
 757
 758
 759 evenPart(pol,var) :=
 760   sum(ratcoeff(pol,var,i*2)*var^i,
 761       i,0,floor(degree(pol,var)/2));
 762
 763 oddPart(pol,var) :=
 764   sum(ratcoeff(pol,var,i*2+1)*var^i,
 765       i,0,ceiling(degree(pol,var)/2));
 766
 767 lienardChipartConditions(pol,var) :=
 768   block([F,G,SR,degP,m],
 769   F : evenPart(pol,var),
 770
 771   G : oddPart(pol,var),
 772
 773   degP : degree(pol,var),
 774
 775   if evenp(degP) then
 776     (
 777     m = p/2,
 778     SR : sSubResCoeff(F,G,var),
 779
 780     return(setify(append(poly2list(pol,var),SR)))
 781     ) /* end if */
 782   else
 783     (
 784     m = (p-1)/2,
 785     SR : sSubResCoeff(var*G,F,var),
 786
 787     return(setify(append(poly2list(pol,var),SR)))
 788     ) /* end else */
 789 ); /* end proc */
 790
 791
 792
 793
 794 /* It computes the difference between the number of roots of P */
 795 /* with positive real parts and the number of roots with negative real parts */
 796 /* Theorem 9.48 (notation "n(P)") */
 797 posNegDiff(a,var) :=
 798   block([p,degP,m,f,g,i],
 799   p : expandIf(a),
 800   degP : degree(p,var),
 801
 802   if evenp(degP) then
 803      (
 804      m : degP/2,
 805      f : sum(coeff(p,var,i*2)* var^i,i,0,m),
 806      g : sum(coeff(p,var,i*2+1)*var^i,i,0,m-1),
 807
 808
 809      return(-sSubResCauchyIndex(g,f,var)+
 810              sSubResCauchyIndexBetween(sPsRem(var*g,f,x),f,var))
 811      ) /* end if */
 812   else
 813      (
 814      m : (degP-1)/2,
 815      f : sum(coeff(p,var,i*2)* var^i,i,0,m),
 816      g : sum(coeff(p,var,2*i+1)* var^i,i,0,m),
 817
 818
 819      return(-sRemCauchyIndexBetween(f,NORM_ALGORITHM(var*g),
 820             var,-inf,+inf)+
 821              sRemCauchyIndexBetween(f,g,
 822             var,-inf,+inf))
 823      ) /* end else */
 824   ); /* end proc */
 825
 826   /* New name for Horner evaluation */
 827
 828 hornerEval(pol,var,value) :=saragHorner(pol,var,value) ;
 829
 830 saragHorner(pol,var,value) :=
 831   block([res,i,degP],
 832
 833   degP : degree(pol,var),
 834
 835   res : leadCoeff(pol,var),
 836
 837   for i : 1 thru degP do
 838      res : value*res+ratcoeff(pol,var,degP-i),
 839
 840   return(res)
 841   ); /* end proc */