share/contrib/sarag/sarag_linear_algebra.mac

   1 /* ------------------------------------------------------------------ */
   2 /* SARAG - Linear Algebra                                             */
   3 /* by Fabrizio Caruso   modified by Alexandre Le Meur and Marie-Françoise Roy                           */
   4
   5
   6 /* It counts the sign change of a determinant after */
   7 /* i consecutive row changes */
   8 revRowsCount(i) := (-1)^(i*(i-1)/2);
   9
  10 exchange(m,i,j) :=
  11   block([aux,size],
  12
  13   size : length(m),
  14
  15   aux : m[i],
  16   m[i] : m[j],
  17   m[j] : aux,
  18
  19   return(m)
  20   ); /* end proc */
  21
  22 exchangeCol(m,i,j) :=
  23    block([aux,size,k],
  24    size : first(third(arrayinfo(m)))+1,
  25    for k : 0 thru size-1 do
  26       (
  27       aux : m[k,i],
  28       m[k,i] : m[k,j],
  29       m[k,j] : aux
  30       ), /* end for */
  31    return(m)
  32    ); /* end proc */
  33
  34
  35
  36 /* Fix columns exchanges */
  37 fixColumns(m,colEx) :=
  38   block([i],
  39    colEx : reverse(colEx),
  40    for i : 1 thru colEx do
  41      exchangeCols(m,colEx[i][1],colEx[i][2]),
  42    return(m)
  43    );  /* end proc */
  44
  45
  46 removeElements(lst,blackList) :=
  47   removeElementsAux(lst,blackList,[],1);
  48
  49 removeElementsAux(lst,blackList,res,indx) :=
  50    if blackList=[] then
  51      append(res,lst)
  52    else
  53      if first(blackList)=indx then
  54        removeElementsAux(rest(lst,1),rest(blackList,1),res,indx+1)
  55      else
  56        if indx<first(blackList) then
  57           removeElementsAux(rest(lst,1),blackList,
  58                             endcons(first(lst),res),indx+1)
  59        else
  60           removeElementsAux(rest(lst,1),rest(blackList),
  61                             endcons(first(lst),res),indx+1);
  62
  63
  64
  65 gaussTriangularize(m) :=
  66   lambda([lst],
  67          [removeElements(array2list(first(lst)),
  68                          third(lst)),
  69           second(lst)])(gaussElimArray(m));
  70
  71 gaussElim(m) :=
  72   lambda([lst],
  73          [array2list(first(lst)),second(lst),third(lst)])(gaussElimArray(m));
  74
  75 /* It performs Gaussian elimiation on a matrix */
  76 /* by columns exchange */
  77 gaussElimArray(m) :=
  78   block([
  79          nRows : length(m),
  80          nCols : length(m[1]),
  81          g : make_array('any),
  82          k,i,j,r,flag,colExList,
  83          zeroRows,offset],
  84
  85   offset : 0, /* it compensates for zero rows */
  86
  87   zeroRows : [],
  88
  89   g : make_array('any,nRows,nCols),
  90
  91   colExList : [],
  92
  93   /* Initialization */
  94   for i : 0 thru nRows-1 do
  95      for j : 0 thru nCols-1 do
  96         g[i,j] : m[i+1][j+1],
  97
  98
  99
 100   /* Main loop */
 101   k : 0, /* index of the row of the pivot */
 102   r : 0, /* index of the column of the pivot */
 103
 104   while r<= nCols-1 and k <= nRows-1 do
 105     (
 106
 107     /* Test for zero row */
 108     flag : true,
 109     j : r+1,
 110     while flag and j <= nCols do
 111       if not(g[k+1-1,j-1]=0) then
 112          flag : false
 113       else
 114          j : j+1,
 115     if flag then
 116       (
 117       zeroRows : endcons(k,zeroRows),
 118       offset : offset + 1
 119       ), /* end if */
 120
 121     if not(j=r+1) and not(flag) then
 122       (
 123       colExList : endcons([r+1-1,j-1],colExList),
 124       exchangeCol(g,r+1-1,j-1)
 125       ),  /* end if */
 126     if not(flag) then
 127       (
 128       for i : k+2 thru nRows do
 129         (
 130         for j : r+2 thru nCols do
 131           (
 132           g[i-1,j-1] :
 133                 g[i-1,j-1] - g[i-1,r+1-1]*
 134                              g[k+1-1,j-1]/g[k+1-1,r+1-1]
 135           ), /* end for */
 136
 137         g[i-1,r+1-1] : 0
 138
 139         ), /* end for */
 140       r : r + 1
 141       ), /* end if */
 142     k : k +1
 143     ), /* end main while loop */
 144
 145 /* add zero rows below the upper square part to zeroRows*/
 146
 147   if k<=nRows-1 then
 148      zeroRows : append(zeroRows, makelist(i,i,k,nRows-1)),
 149
 150   return([g,colExList,zeroRows])
 151   ); /* end proc */
 152
 153
 154 gaussElimMod(m,gmod) :=
 155   gaussElimWithDivisionMod(m,gmod);
 156
 157 /*
 158   lambda([lst],
 159          [array2list(first(lst)),second(lst),third(lst)])(gaussElimModArray(m,gmod));
 160 */
 161
 162 gaussElimWithDivisionMod(m,gmod) :=
 163   lambda([lst],
 164          [array2list(first(lst)),second(lst),third(lst)])(gaussElimWithDivisionModArray(m,gmod));
 165
 166 /* It performs Gaussian elimiation on a matrix */
 167 /* by columns exchange */
 168 gaussElimWithDivisionModArray(m,gmod) :=
 169   block([
 170          nRows : length(m),
 171          nCols : length(m[1]),
 172          g : make_array('any),
 173          k,i,j,r,flag,colExList,
 174          zeroRows,offset,modInv],
 175
 176   offset : 0, /* it compensates for zero rows */
 177
 178   modulus : gmod, /* it sets the modular computation */
 179
 180   zeroRows : [],
 181
 182   g : make_array('any,nRows,nCols),
 183
 184   colExList : [],
 185
 186   /* Initialization */
 187   for i : 0 thru nRows-1 do
 188      for j : 0 thru nCols-1 do
 189         g[i,j] : m[i+1][j+1],
 190
 191
 192   /* Main loop */
 193   k : 0, /* index of the row of the pivot */
 194   r : 0, /* index of the column of the pivot */
 195
 196   while r<= nCols-1 and k <= nRows-1 do
 197     (
 198
 199     /* Test for zero row */
 200     flag : true,
 201     j : r+1,
 202     while flag and j <= nCols do
 203       if not(g[k+1-1,j-1]=0) then
 204          flag : false
 205       else
 206          j : j+1,
 207     if flag then
 208       (
 209       zeroRows : endcons(k,zeroRows),
 210       offset : offset + 1
 211       ), /* end if */
 212
 213     if not(j=r+1) and not(flag) then
 214       (
 215       colExList : endcons([r+1-1,j-1],colExList),
 216       exchangeCol(g,r+1-1,j-1)
 217       ),    /* end if */
 218     if not(flag) then
 219       (
 220       for i : k+2 thru nRows do
 221         (
 222         /* it computes the inverse for the column */
 223         modInv : rat(1/g[k+1-1,r+1-1]),
 224
 225         for j : r+2 thru nCols do
 226           (
 227
 228           g[i-1,j-1] :
 229                 ratexpand(g[i-1,j-1] - g[i-1,r+1-1]*g[k+1-1,j-1]*modInv)
 230
 231           ), /* end for */
 232
 233         g[i-1,r+1-1] : 0
 234
 235         ), /* end for */
 236       r : r + 1
 237       ), /* end if */
 238     k : k +1
 239     ), /* end main while loop */
 240
 241 /* add zero rows below the upper square part to zeroRows*/
 242
 243   if k<=nRows-1 then
 244      zeroRows : append(zeroRows, makelist(i,i,k,nRows-1)),
 245
 246   modulus : false, /* it resets the computation over the rationals */
 247   return([g,colExList,zeroRows])
 248   ); /* end proc */
 249
 250
 251 gaussElimWithoutDivisionMod(m,gmod) :=
 252   lambda([lst],
 253          [array2list(first(lst)),
 254           second(lst),third(lst)])(gaussElimWithoutDivisionModArray(m,gmod));
 255
 256 gaussElimWithoutDivisionModArray(m,gmod) :=
 257   block([
 258          nRows : length(m),
 259          nCols : length(m[1]),
 260          g : make_array('any),
 261          k,i,j,r,flag,colExList,
 262          zeroRows,offset],
 263
 264   offset : 0, /* it compensates for zero rows */
 265
 266   modulus : gmod, /* it sets the modular computation */
 267
 268   zeroRows : [],
 269
 270   g : make_array('any,nRows,nCols),
 271
 272   colExList : [],
 273
 274   /* Initialization */
 275   for i : 0 thru nRows-1 do
 276      for j : 0 thru nCols-1 do
 277         g[i,j] : m[i+1][j+1],
 278
 279
 280   /* Main loop */
 281   k : 0, /* index of the row of the pivot */
 282   r : 0, /* index of the column of the pivot */
 283
 284   while r<= nCols-1 and k <= nRows-1 do
 285     (
 286
 287     /* Test for zero row */
 288     flag : true,
 289     j : r+1,
 290     while flag and j <= nCols do
 291       if not(g[k+1-1,j-1]=0) then
 292          flag : false
 293       else
 294          j : j+1,
 295     if flag then
 296       (
 297       zeroRows : endcons(k,zeroRows),
 298       offset : offset + 1
 299       ), /* end if */
 300
 301     if not(j=r+1) and not(flag) then
 302       (
 303       colExList : endcons([r+1-1,j-1],colExList),
 304       exchangeCol(g,r+1-1,j-1)
 305       ),    /* end if */
 306     if not(flag) then
 307       (
 308       for i : k+2 thru nRows do
 309         (
 310         /* it computes the inverse for the column */
 311 /*
 312         modInv : rat(1/g[k+1-1,r+1-1]),
 313 */
 314         for j : r+2 thru nCols do
 315           (
 316
 317           g[i-1,j-1] :
 318                 ratexpand(g[k+1-1,r+1-1]*g[i-1,j-1] - g[i-1,r+1-1]*g[k+1-1,j-1])
 319
 320           ), /* end for */
 321
 322         g[i-1,r+1-1] : 0
 323
 324         ), /* end for */
 325       r : r + 1
 326       ), /* end if */
 327     k : k +1
 328     ), /* end main while loop */
 329
 330 /* add zero rows below the upper square part to zeroRows*/
 331
 332   if k<=nRows-1 then
 333      zeroRows : append(zeroRows, makelist(i,i,k,nRows-1)),
 334
 335   modulus : false, /* it resets the computation over the rationals */
 336   return([g,colExList,zeroRows])
 337   ); /* end proc */
 338
 339 echelonMod(m,gmod) :=
 340   echelonWithDivisionMod(m,gmod);
 341
 342
 343 echelonWithDivisionMod(m,gmod) :=
 344   lambda([lst],
 345          [array2list(first(lst)),second(lst),third(lst)])(echelonWithDivisionModArray(m,gmod));
 346
 347 /* It performs Gaussian elimiation on a matrix */
 348 /* by columns exchange */
 349 echelonWithDivisionModArray(m,gmod) :=
 350   block([
 351          nRows : length(m),
 352          nCols : length(m[1]),
 353          g : make_array('any),
 354          k,i,j,r,flag,colExList,
 355          zeroRows,offset,modInv],
 356
 357   offset : 0, /* it compensates for zero rows */
 358
 359   modulus : gmod, /* it sets the modular computation */
 360
 361   zeroRows : [],
 362
 363   g : make_array('any,nRows,nCols),
 364
 365   colExList : [],
 366   /* Initialization */
 367   for i : 0 thru nRows-1 do
 368      for j : 0 thru nCols-1 do
 369         g[i,j] : m[i+1][j+1],
 370
 371
 372   /* Main loop */
 373   k : 0, /* index of the row of the pivot */
 374   r : 0, /* index of the column of the pivot */
 375
 376   while r<= nCols-1 and k <= nRows-1 do
 377     (
 378
 379     /* Test for zero row */
 380     flag : true,
 381     j : r+1,
 382     while flag and j <= nCols do
 383       (
 384       if not(g[k+1-1,j-1]=0) then
 385          flag : false
 386       else
 387          j : j+1
 388       ), /* end while */
 389
 390     if flag then
 391       (
 392       zeroRows : endcons(k,zeroRows),
 393       offset : offset + 1
 394       ), /* end if */
 395
 396     if not(j=r+1) and not(flag) then
 397       (
 398       colExList : endcons([r+1-1,j-1],colExList),
 399       exchangeCol(g,r+1-1,j-1)
 400       ),  /* end if */
 401     if not(flag) then
 402       (
 403
 404       for i : 1 thru k do
 405         (
 406         /* it computes the inverse for the column */
 407         modInv : rat(1/g[k+1-1,r+1-1]),
 408
 409
 410         for j : r+2 thru nCols do
 411           (
 412
 413           g[i-1,j-1] :
 414                 ratexpand(g[i-1,j-1] - g[i-1,r+1-1]*g[k+1-1,j-1]*modInv)
 415
 416           ), /* end for */
 417
 418         g[i-1,r+1-1] : 0
 419
 420         ), /* end for */
 421
 422
 423       for i : k+2 thru nRows do
 424         (
 425         /* it computes the inverse for the column */
 426         modInv : rat(1/g[k+1-1,r+1-1]),
 427
 428
 429         for j : r+2 thru nCols do
 430           (
 431
 432           g[i-1,j-1] :
 433                 ratexpand(g[i-1,j-1] - g[i-1,r+1-1]*g[k+1-1,j-1]*modInv)
 434
 435           ), /* end for */
 436
 437         g[i-1,r+1-1] : 0
 438
 439         ), /* end for */
 440       r : r + 1
 441       ), /* end if */
 442     k : k +1
 443     ), /* end main while loop */
 444
 445 /* add zero rows below the upper square part to zeroRows*/
 446
 447   if k<=nRows-1 then
 448      zeroRows : append(zeroRows, makelist(i,i,k,nRows-1)),
 449
 450   modulus : false, /* it resets the computation over the rationals */
 451   return([g,colExList,zeroRows])
 452   ); /* end proc */
 453
 454
 455 /* It computes the determinant by Gauss elimination */
 456 /* Algorithm 8.26                                   */
 457 gaussDet(m) :=
 458   block([
 459          g : make_array('any),
 460          k,i,j,flag,colEx,size],
 461
 462   size : length(m),
 463   g : list2array(m),
 464   colEx : 0,
 465
 466   /* Initialization */
 467   for i : 0 thru size-1 do
 468      for j : 0 thru size-1 do
 469         g[i,j] : m[i+1][j+1],
 470
 471   /* Main loop */
 472   for k : 0 thru size-2 do
 473     (
 474     /* Test for zero row */
 475     flag : true,
 476     j : k+1,
 477     while flag and j <= size do
 478       if not(g[k+1-1,j-1]=0) then
 479          flag : false
 480       else
 481          j : j+1,
 482     if flag then
 483       return(0),
 484
 485     /* column exchange check */
 486     if not(j=k+1) then
 487       (
 488       /* column exchange required */
 489       colEx : colEx+1,
 490       exchangeCol(g,k+1-1,j-1)
 491       ), /* end if */
 492     /* elimination */
 493     for i : k+2 thru size do
 494       (
 495       for j : k+2 thru size do
 496         g[i-1,j-1] : g[i-1,j-1] - g[i-1,k+1-1]*g[k+1-1,j-1]/g[k+1-1,k+1-1],
 497       g[i-1,k+1-1] : 0
 498       ) /* end for */
 499     ), /* end for */
 500   return(expand((-1)^colEx * product(g[i,i],i,0,size-1)))
 501   ); /* end proc */
 502
 503
 504
 505 /* It computes the determinant by Dogdson-Jordan-Bareiss method */
 506 /* Algorithm 8.30                                               */
 507 bareissDet(m) :=
 508   block([
 509          g : make_array('any),
 510          k,i,j,flag,colEx,oldPivot, size],
 511   size : length(m),
 512   g : list2array(m),
 513   colEx : 0,
 514   /* Initialization */
 515   for i : 0 thru size-1 do
 516      for j : 0 thru size-1 do
 517         g[i,j] : m[i+1][j+1],
 518   oldPivot : 1,
 519   /* Main loop */
 520   for k : 0 thru size-2 do
 521     (
 522     /* Test for zero row */
 523     flag : true,
 524     j : k+1,
 525     while flag and j <= size do
 526       if not(g[k+1-1,j-1]=0) then
 527          flag : false
 528       else
 529          j : j+1,
 530     if flag then
 531       return(0),
 532     /* column exchange check */
 533     if not(j=k+1) then
 534       (
 535       /* column exchange required */
 536       colEx : colEx+1,
 537       exchangeCol(g,k+1-1,j-1)
 538       ), /* end if */
 539     /* elimination */
 540     for i : k+2 thru size do /* column loop */
 541       (
 542       for j : k+2 thru size do /* row loop */
 543         (
 544         if k= 0 then
 545           oldPivot : 1
 546         else
 547           oldPivot : g[k-1,k-1],
 548         g[i-1,j-1] : (g[k+1-1,k+1-1]*g[i-1,j-1]-g[i-1,k+1-1]*g[k+1-1,j-1])/oldPivot
 549         ), /* end for */
 550
 551       g[i-1,k+1-1] : 0
 552       ) /* end for */
 553     ), /* end if */
 554   return(expand((-1)^colEx * g[size-1,size-1]))
 555   ); /* end proc */
 556
 557
 558 princSubMat(m,ord) :=
 559   block([res:make_array('any,ord,ord),i,j],
 560
 561   for i : 0 thru ord-1 do
 562     for j : 0 thru ord-1 do
 563       res[i,j] : m[i,j],
 564   return(res)
 565   ); /* end proc */
 566
 567 extractCol(m,j,size) :=
 568   block([res:make_array('any,size),k],
 569
 570   for k : 0 thru size-1 do
 571     res[k] : m[k,j],
 572
 573   return(res)
 574   ); /* end proc */
 575
 576 extractRow(m,i,size) :=
 577   block([res:make_array('any,size),k],
 578
 579   for k : 0 thru size-1 do
 580     res[k] : m[i,k],
 581
 582   return(res)
 583   ); /* end proc */
 584
 585
 586 /* Linear algebra support routines */
 587 solveSys(mat,vec,solver) :=
 588   block([sol,unk,i,j,sys,nRows,nCols,newEq],
 589   sys : [],
 590   nRows : length(mat),
 591   nCols : length(mat[1]),
 592
 593   for i : 1 thru nRows do
 594     (
 595     newEq : sum(mat[i][j]*unk[j],j,1,nCols)-vec[i]=0,
 596     sys : endcons(newEq,sys)
 597     ), /* end for */
 598   sol : solver(sys,makelist(unk[i],i,1,nCols)),
 599
 600   if sol = [] then
 601     return([])
 602   else
 603     return(map(second,sol))
 604   ); /* end proc */
 605
 606
 607
 608 newtonFromPoly(pol,var,len) :=
 609   array2singleList(newtonArrFromPoly(pol,var,len));
 610
 611 newtonArrFromPoly(pol,var,len) :=
 612   block([newtonRes,degPol,i,j],
 613     degPol : degree(pol,var),
 614     newtonRes : make_array('any,len+1),
 615
 616     newtonRes[0] : degPol,
 617
 618     for i : 1 thru len do
 619       (
 620
 621       newtonRes[i] : (degPol-i)*ratcoeff(pol,var,degPol-i)-
 622                      sum(ratcoeff(pol,var,degPol-j)*
 623                      newtonRes[i-j],j,1,
 624                      min(i,degPol))
 625       ), /* end for */
 626     return(newtonRes)
 627     ); /* end proc */
 628
 629
 630 polyFromNewton(ns,var) :=
 631   list2poly(getCoeffFromNewton(ns),var);
 632
 633 getCoeffFromNewton(ns) :=
 634   array2singleList(getCoeffFromNewtonArray(singleList2array(ns)));
 635
 636 getCoeffFromNewtonArray(ns) :=
 637   block([i,j,degRes : arrayLength(ns)-1,resCoeff:make_array('any)],
 638
 639    resCoeff : make_array('any,degRes+1),
 640
 641    resCoeff[degRes]:1,
 642
 643    for i : 1 thru degRes do
 644      resCoeff[degRes-i] : -1/i*sum(ns[j+1-1]*resCoeff[degRes+j-i],j,1,i),
 645
 646
 647    return(resCoeff)
 648    ); /* end proc */
 649
 650
 651 getCoeffFromNewtonList(ns) :=
 652   array2singleList(getCoeffFromNewton(ns));
 653
 654 matrixProd(a,b) :=
 655   array2list(matrixProdArray(list2array(a),list2array(b)));
 656
 657 matrixProdArray(a,b) :=
 658   block([i,j,k,n,m,l,res:mame_array('any)],
 659
 660   n : first(third(arrayinfo(a)))+1,
 661
 662   m : second(third(arrayinfo(a)))+1,
 663
 664  if not(first(third(arrayinfo(b)))+1=m) then
 665     (
 666     print("matrixProd) incompatible matrices"),
 667     print("matrixProd) a : ", a),
 668     print("matrixProd) b : ", b),
 669     return(false)
 670     ) /* end if */
 671   else
 672     (
 673   l : second(third(arrayinfo(b)))+1,
 674
 675   res : make_array('any,n,l),
 676
 677   for i : 1 thru n do
 678     for k : 1 thru l do
 679
 680        res[i-1,k-1] : sum(a[i-1,j-1]*b[j-1,k-1],j,1,m),
 681
 682
 683   return(res)
 684     ) /* end else */
 685   ); /* end proc */
 686
 687 matrixTrace(A) :=
 688   block([ALen,i,res],
 689     ALen : numOfCols(A),
 690
 691     res : 0,
 692     for i : 1 thru ALen do
 693       res : res + A[i-1,i-1],
 694     return(res)
 695     ); /* end proc */
 696
 697 discreteSquareRoot(n) :=
 698  block([i,res],
 699  res : 1,
 700  /* after 5.9.2 */
 701  /*
 702  for i : 1 thru ceiling(bfloat(n/2)) do
 703    if i*i >= n then
 704      return(i)
 705  */
 706  /* before 5.9.2 but also forwards-compatible*/
 707  for i : 1 thru comp_ceiling(bfloat(n/2)) do
 708    if i*i >= n then
 709      return(i)
 710  ); /* end proc */
 711
 712
 713 babyGiantCharPol(A,var) :=
 714   babyGiantCharPolAux(list2array(A),var);
 715
 716 babyGiantCharPolAux(A,var) :=
 717   list2poly(array2singleList(getCoeffFromNewtonArray(
 718          getNewtonFromMatrix(A))),
 719         var);
 720
 721
 722 matrixMinus(A,B) :=
 723   makelist(
 724     makelist(A[i][j]-B[i][j],j,1,length(A[1])),
 725     i,1,length(A));
 726
 727 diagonal(elem,len) :=
 728   makelist(
 729     makelist(if(i=j) then elem else 0,j,1,len),
 730     i,1,len);
 731
 732 gaussCharPol(A,var) :=
 733   expandIf(gaussDet(matrixMinus(A,diagonal(var,length(A)))));
 734
 735 bareissCharPol(A,var) :=
 736   expandIf(bareissDet(matrixMinus(A,diagonal(var,length(A)))));
 737
 738 /* Input : a bidimensional array */
 739 /* Output : the Newton sums of the corresponding characteristic polynomial */
 740 getNewtonFromMatrix(A) :=
 741    block([i,j,k,r,n,dsr,
 742           B:make_array('any),G:make_array('any),ss:make_array('any)],
 743
 744    n : numOfRows(A),
 745    ss : make_array('any,n +1 /* +1 */),
 746
 747    dsr : discreteSquareRoot(n),
 748    if dsr*dsr = n then
 749      r : dsr+1
 750    else
 751      r : dsr,
 752    B : make_array('any,r-1),
 753    for i : 1 thru r-1 do
 754      B[i-1] : make_array('any,n,r),
 755    /* baby step */
 756    B[1-1] : A,
 757    ss[0] : n,
 758    ss[1] : matrixTrace(A),
 759    /* baby loop */
 760    for i : 1 thru r-2 do
 761      (
 762      B[i+1 -1] : matrixProdArray(A,B[i -1]),
 763      ss[i+1] : matrixTrace(B[i+1 -1])
 764
 765      ),  /* end for */
 766    G : make_array('any,r-1),
 767    for j : 1 thru r-1 do
 768      G[j-1] : make_array('any,n,r),
 769    /* giant step */
 770    G[1 -1] : matrixProdArray(A,B[r-1 -1]),
 771    ss[r] : matrixTrace(G[1 -1]),
 772
 773    for j : 1 thru r-2 do
 774      (
 775      if (j+1)*r <= n then
 776        (
 777        G[j+1 -1] : matrixProdArray(G[1-1],G[j -1]),
 778        ss[(j+1)*r] : matrixTrace(G[j+1 -1])
 779        ) /* end if */
 780      ),   /* end for */
 781    /* Newton's sums */
 782    for i : 1 thru r-1 do
 783       for j : 1 thru r-1 do
 784          (
 785          if j*r+i <= n then
 786            ss[j*r+i] : matrixTrace(matrixProdArray(B[i -1],G[j -1]))
 787          ), /* end for */
 788
 789 return(ss)
 790 ); /* end proc */
 791
 792 companion(pol,var) :=
 793   array2list(companion(pol,var));
 794
 795 companionArray(pol,var) :=
 796   block(
 797     [m : make_array('any),n, i,j],
 798
 799     n : degree(pol,var),
 800     m: make_array('any,n,n),
 801
 802     for j : 0 thru n - 2 do
 803       m[0,j] : 0,
 804     m[0,n-1] : -coeff(pol,var,0),
 805
 806
 807     for i : 1 thru n-1 do
 808       (
 809       for j : 0 thru i-2 do
 810         m[i,j] : 0,
 811       m[i,i-1] : 1,
 812       for j : i thru n-2 do
 813         m[i,j] : 0,
 814       m[i,n-1] : -coeff(pol,var,i)
 815       ), /* end for */
 816     return(m)
 817     ); /* end proc */
 818
 819
 820 /* This should be changed with the algorithm shown in Proposition 4.7 */
 821 companionPoly2Newton(pol,var) :=
 822   array2singleList(getNewtonFromMatrix(companionArray(pol,var)));
 823
 824 newton2Poly(newtonLst) :=
 825   array2singleList(newtonArr2Poly(singleList2array(newtonLst)));
 826
 827 newtonArr2Poly(newtonArr) :=
 828    block([res,i,j,degP],
 829    degP : arrayLength(newtonArr)-1,
 830
 831    res : make_array('any,degP+1),
 832    res[degP] : 1,
 833    for i : 1 thru degP do
 834      (
 835
 836      for j : 1 thru i do
 837        res[degP-i] : sum(res[degP-i+j]*newtonArr[j],j,1,i)*(-1/i)
 838
 839      ), /* end for */
 840    return(res)
 841    ); /* end proc */
 842
 843 /* Descartes Signature */
 844 descartesSignature(mtx) :=
 845   block([charP,x,revCoeffList,charPDeg,i],
 846     charP : babyGiantCharPol(mtx,x), /* DEBUGGING */
 847     charPDeg : length(mtx),
 848     revCoeffList :  reverse(poly2list(charP,x)),
 849
 850     return(signChanges(revCoeffList)-
 851            signChanges(
 852              makelist((-1)^(charPDeg-i)*revCoeffList[i+1],i,0,
 853                                                         charPDeg)))
 854
 855     ); /* end proc */
 856