Changeset cc07066 in git


Ignore:
Timestamp:
Aug 31, 2020, 11:57:27 AM (4 years ago)
Author:
Janko Boehm <bohem@…>
Branches:
(u'spielwiese', '5b153614cbc72bfa198d75b1e9e33dab2645d9fe')
Children:
fdea544a15d2eb16f952ba088de11e5efa5494d9
Parents:
6924d452049dc8bfbb851672c4b3ceac05081b6c
Message:
update
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • Singular/LIB/pfd.lib

    r6924d4 rcc07066  
    55LIBRARY: pfd.lib Multivariate Partial Fraction Decomposition
    66
    7 AUTHOR: Marcel Wittmann, e-mail: mwittman@rhrk.uni-kl.de
     7AUTHOR: Marcel Wittmann, e-mail: mwittman@mathematik.uni-kl.de
    88
    99OVERVIEW:
     
    1212\"smaller\" numerators and denominators.
    1313This can be used to shorten the IBP reduction coeffcients of multi-loop Feynman
    14 integrals as shown in [J. Böhm, M. Wittmann, Z. Wu, Y. Xu, Y. Zhang: 'IBP
     14integrals as shown in [J. Boehm, M. Wittmann, Z. Wu, Y. Xu, Y. Zhang: 'IBP
    1515reduction coefficients made simple'] (preprint 2020). For this application,
    1616we also provide a procedure that applies the algorithm to all entries of a
     
    3030  getStringpfd();         turn a decomposition gotten from @code{pfd} into one
    3131                          string
    32   getStringpfd_indexed(); like @code{getStringpfd}, but writing the denominator
    33                           factors just as @code{q1},@code{q2},...
     32  getStringpfd_indexed(); like @code{getStringpfd}, but writes the denominator
     33                          factors just as @code{q1}, @code{q2}, ...
    3434  readInputTXT();         read a matrix of rational functions from a txt-file
    3535  pfdMat();               apply @code{pfd} to a matrix of rational functions
     
    5050static proc mod_init()
    5151{
     52  printlevel = 2;
     53  // export some static functions so they can be run in parallel using parallelWaitAll:
    5254  if (!defined(Tasks)) {LIB "tasks.lib";}
    53   // export some static functions so they can be run in parallel using parallelWaitAll
    5455  exportto(Tasks,pfdWrap);
    5556  importfrom(Tasks,pfdWrap);
    5657  exportto(Pfd,pfdWrap);
    57 
    5858  exportto(Tasks,testEntry);
    5959  importfrom(Tasks,testEntry);
     
    6666          pfd(arguments[, parallelize]);     arguments list, parallelize int
    6767RETURN:   a partial fraction decomposition of f/g as a list @code{l} where
    68           @code{l[1]} is an ideal generate by irreducible polynomials and
     68          @code{l[1]} is an ideal generated by irreducible polynomials and
    6969          @code{l[2]} is a list of fractions.
    7070          Each fraction is represented by a list of
    7171       @* 1) the numerator polynomial
    72        @* 2) an intvec giving the indices @code{i} for which @code{l[1][i]}
    73              occurs as a factor in the denominator
     72       @* 2) an intvec of indices @code{i} for which @code{l[1][i]} occurs
     73             as a factor in the denominator
    7474       @* 3) an intvec containing the exponents of those irreducible factors.
    75        @* Setting debug to a positive integer measures runtimes and
    76           creates a log file (default: debug=0).
     75       @* Setting @code{debug} to a positive integer measures runtimes and
     76          creates a log file (default: @code{debug=0}).
    7777       @* The denominator g can also be given in factorized form as a list of
    7878          an ideal of irreducible non constant polynomials and an intvec of
    7979          exponents. This can save time since the first step in the algorithm is
    80           to factorize g. (g=list(ideal(0),intvec(0:0)) represents a denominator
    81           of 1.)
     80          to factorize g. (A list of the zero-ideal and an empty intvec
     81          represents a denominator of 1.)
    8282       @* If instead of f and g, the input is a single list (or even a list of
    83           lists) containing elements of the form list(f,g[,debug]) (f,g[,debug]
    84           as above), the algorithm is applied to all entries in parallel (using
    85           @ref{parallel_lib}), if @code{parallelize=1} (default) and in sequence
    86           if @code{parallelize=0}. A list of the results is returned.
     83          lists) containing elements of the form @code{list(f,g[,debug])}
     84          (@code{f,g,debug} as above), the algorithm is applied to all entries
     85          in parallel (using @ref{parallel_lib}), if @code{parallelize=1}
     86          (default) and in sequence if @code{parallelize=0}. A list (or list of
     87          lists) of the results is returned.
    8788NOTE:     The result depends on the monomial ordering. For \"small\" results
    88           use dp.
     89          use @code{dp}.
    8990EXAMPLE:  example pfd; shows an example
    9091SEE ALSO: checkpfd, evaluatepfd, displaypfd, displaypfd_long, pfdMat"
     
    141142
    142143  int debug=0; link l=":w ";
    143   if(size(#)>2)
     144  if(size(#)>3) {ERROR("wrong number of arguments, expected 2 or 3");}
     145  if(size(#)==3)
    144146  {
    145147    debug=#[3];
     
    160162    if(debug)
    161163      {fprintf(l,"%ntotal: 0 ms (numerator was 0)",0); close(l);}
    162     if(voice<3)
    163       {displaypfd(dec);}
     164    if(voice<=printlevel) {displaypfd(dec);}
    164165    return(dec);
    165166  }
     
    171172      if(debug)
    172173        {fprintf(l,"%ntotal: 0 ms (denominator was constant)",0); close(l);}
    173       if(voice<3)
    174         {displaypfd(dec);}
     174      if(voice<=printlevel) {displaypfd(dec);}
    175175      return(dec);
    176176    }
    177177
    178     // (1) factorization of the denominator:
     178    // (1) factorization of the denominator ////////////////////////////////////
    179179    if(debug) {int t1 = rtimer; write(l,"factorizing ");}
    180180    list factor = factorize(g);
    181181    number lcoeff;
    182     for(i=2; i<=size(factor[1]); i++) // making polynomials unique
     182    for(i=2; i<=size(factor[1]); i++)
    183183    {
    184184      lcoeff = leadcoef(factor[1][i]);
     
    203203      if(debug)
    204204        {fprintf(l,"%ntotal: 0 ms (denominator was constant)",0); close(l);}
    205       if(voice<3)
    206         {displaypfd(dec);}
     205      if(voice<=printlevel) {displaypfd(dec);}
    207206      return(dec);
    208207    }
     
    221220  }
    222221  else
    223     {ERROR("wrong type for second argument, expected poly or list(ideal,intvec)");}}
    224 
    225   // (2) Nullstellensatz decomposition
     222  {ERROR("wrong type for second argument, expected poly or list(ideal,intvec)");}}
     223
     224  // (2) Nullstellensatz decomposition /////////////////////////////////////////
    226225  if(debug)
    227226    {int t2 = rtimer; write(l,"Nullstellensatz decomposition ");}
     
    248247    {t2 = rtimer-t2; fprintf(l,"done! (%s ms, %s terms)", t2, size(dec));}
    249248
    250   // (3) short numerator decomposition
     249  // (3) short numerator decomposition /////////////////////////////////////////
    251250  if(debug)
    252251    {int t3 = rtimer; write(l,"short numerator decompositions "); counter=0;}
     
    277276    {t3 = rtimer-t3; fprintf(l,"done! (%s ms, %s terms)", t3, size(dec));}
    278277
    279   // (4) algebraic dependence decomposition
     278  // (4) algebraic dependence decomposition ////////////////////////////////////
    280279  if(debug)
    281280    {int t4 = rtimer; write(l,"algebraic dependence decomposition "); counter=0;}
     
    305304    {t4 = rtimer-t4; fprintf(l,"done! (%s ms, %s terms)", t4, size(dec));}
    306305
    307   // (5) numerator decomposition
     306  // (5) numerator decomposition ///////////////////////////////////////////////
    308307  if(debug)
    309308    {int t5 = rtimer; write(l,"numerator decompositions "); counter=0;}
     
    335334  dec = list(q,dec);
    336335  if(debug) {fprintf(l,"%ntotal: %s ms", t1+t2+t3+t4+t5); close(l);}
    337   if(voice<3) {displaypfd(dec);}
     336  if(voice<=printlevel) {displaypfd(dec);}
    338337  return(dec);
    339338}
     
    366365  // a more complicated example
    367366  ring S = 0,(s12,s15,s23,s34,s45),dp;
    368   poly f = (7*s12^4*s15^2 + 11*s12^3*s15^3 + 4*s12^2*s15^4 - 10*s12^4*s15*s23
     367  poly f = 7*s12^4*s15^2 + 11*s12^3*s15^3 + 4*s12^2*s15^4 - 10*s12^4*s15*s23
    369368    - 14*s12^3*s15^2*s23 - 4*s12^2*s15^3*s23 + 3*s12^4*s23^2 + 3*s12^3*s15*s23^2
    370369    + 13*s12^4*s15*s34 + 12*s12^3*s15^2*s34 + 2*s12^2*s15^3*s34
     
    402401    - 8*s12*s15*s34*s45^3 - 10*s15^2*s34*s45^3 + 9*s12*s23*s34*s45^3
    403402    + s12*s34^2*s45^3 + 8*s15*s34^2*s45^3 - 9*s23*s34^2*s45^3 - s34^3*s45^3
    404     - s15^2*s45^4 + s15*s34*s45^4);
     403    - s15^2*s45^4 + s15*s34*s45^4;
    405404  poly g = 4*s12*s15*(s12 + s15 - s34)*(s15 - s23 - s34)*(s12 + s23 - s45)
    406405                                 *(s12 - s34 - s45)*(s12 + s15 - s34 - s45)*s45;
     
    424423
    425424  if(m==0) // denominator is 1
    426   {
    427     return(list(l)); //do nothing, return input
    428   }
     425    {return(list(l));} // do nothing, return input
    429426
    430427  ideal qe = q[indices];
    431428  for(int i=1; i<=m; i++)
    432   {
    433     qe[i] = qe[i]^e[i];
    434   }
     429    {qe[i] = qe[i]^e[i];}
    435430  matrix T;
    436431  ideal qe_std = liftstd(qe,T);
     
    438433  if(deg(qe_std) == 0)
    439434  {
    440     T = T/qe_std[1]; // now 1 = T[1,1]*qe[1] + ... + T[m,1]*qe[m] is a Nullstellensatz certificate
     435    T = T/qe_std[1];
     436      // now 1 = T[1,1]*qe[1] + ... + T[m,1]*qe[m] is a Nullstellensatz certificate
    441437    list dec;
    442438    poly h;
     
    445441      h = T[i,1];
    446442      if(h != 0)
    447       {
    448         dec[size(dec)+1] = list(f*h,delete(indices,i),delete(e,i));
    449       }
     443        {dec[size(dec)+1] = list(f*h,delete(indices,i),delete(e,i));}
    450444    }
    451445    return(dec);
     
    453447  else
    454448  {
    455     return(list(l)); //do nothing, return input
     449    return(list(l)); // do nothing, return input
    456450  }
    457451}
     
    472466    if(debug)
    473467      {fprintf(ll,"      shortNumeratorDecompStep: %s ms (m=%s, e=%s) "
    474                   +"--> constant denominator", rtimer-tt, m, e);}
    475     return(list(list(),list(l))); //do nothing, return input
    476   }
    477 
    478   ideal q_denom = q[indices]; //factors in the denominator
     468                    +"--> constant denominator", rtimer-tt, m, e);}
     469    return(list(list(),list(l))); // do nothing, return input
     470  }
     471
     472  ideal q_denom = q[indices]; // factors occuring in the denominator
    479473  matrix T;
    480474  ideal q_std = liftstd(q_denom,T);
     
    488482                  +"--> remainder is nonzero", rtimer-tt, m, e);}
    489483    return(list(list(),list(l))); // if there is a rest, decomposing further would
    490                                   // not help in the next step (alg. depend. decomposition)
    491   }
     484  }                               // not help in the next step (alg. depend. decomposition)
    492485
    493486  matrix a = divrem[1]/divrem[3][1,1]; // now f = r + a[1,1]*q_std[1] + ... +a[m,1]*q_std[m]
    494   a = T*a; // lift coefficients     ==>    now f = r + a[1,1]*q[1]     + ... +a[m,1]*q[m]
     487  a = T*a;   // lift coefficients   ==>   now f = r + a[1,1]*q[1]     + ... +a[m,1]*q[m]
    495488
    496489  // reduce w.r.t. groebner basis of syz(q) to make the numerators "smaller":
     
    520513    }
    521514    dec[size(dec)+1] = fraction;
    522     }
     515  }
    523516
    524517  if(debug)
     
    544537    if(debug)
    545538      {fprintf(ll,"      algDependDecompStep: %s ms (m=%s, e=%s)", rtimer-tt, m, e);}
    546     return(list(l)); //do nothing, return input
     539    return(list(l)); // do nothing, return input
    547540  }
    548541
    549542  if(m<=d)
    550543  {
    551     if(size(syz(module(transpose(jacob(ideal(q[indices]))))))==0) //jacobian criterion
     544    if(size(syz(module(transpose(jacob(ideal(q[indices]))))))==0) // jacobian criterion
    552545    {
    553546      if(debug)
    554547        {fprintf(ll,"      algDependDecompStep: %s ms (m=%s, e=%s) "
    555548                    +"--> alg. indep.", rtimer-tt, m, e);}
    556       return(list(l)); //do nothing, return input
     549      return(list(l)); // do nothing, return input
    557550    }
    558551  }
     
    567560  ideal I;
    568561  for(i=1; i<=m; i++)
    569   {
    570     I[i] = y(i)-q[indices[i]];//^e[i];
    571   }
     562    {I[i] = y(i)-q[indices[i]];}
    572563
    573564  ideal annihilatingPolys = eliminate(I,intvec(1..d));
     
    575566  poly g = annihilatingPolys[1];
    576567
    577   poly tail = g[size(g)]; // term of lowest dp-order (thus lowest degree)
     568  poly tail = g[size(g)];    // term of lowest dp-order (thus lowest degree)
    578569  number tcoeff = leadcoef(tail);
    579570  intvec texpon = leadexp(tail);
     
    597588    for(i=1; i<=m; i++)
    598589    {
    599       //pow = e[i]*(expon[i]-texpon[i]-1);
    600590      pow = expon[i]-texpon[i]-e[i];
    601591      if(pow>=0)
     
    634624      {fprintf(ll,"      numeratorDecompStep: %s ms (m=%s, e=%s) "
    635625                  +"--> constant denominator", rtimer-tt, m, e);}
    636     return(list(list(),list(l))); //do nothing, return input
    637   }
    638 
    639   ideal q_denom = q[indices]; //factors in the denominator
     626    return(list(list(),list(l))); // do nothing, return input
     627  }
     628
     629  ideal q_denom = q[indices]; // factors in the denominator
    640630  matrix T;
    641631  ideal q_std = liftstd(q_denom,T);
     
    643633  matrix a = divrem[1]/divrem[3][1,1];
    644634  poly r = divrem[2][1]/divrem[3][1,1]; // now f = r + a[1,1]*q_std[1] + ... +a[m,1]*q_std[m]
    645   a = T*a; // lift coefficients     ==>    now f = r + a[1,1]*q[1]     + ... +a[m,1]*q[m]
     635  a = T*a;   // lift coefficients    ==>   now f = r + a[1,1]*q[1]     + ... +a[m,1]*q[m]
    646636
    647637  // reduce w.r.t. groebner basis of syz(q) to make the numerators "smaller"
     
    652642  list fraction,dec,rest;
    653643  if(r!=0)
    654   {
    655     rest[1] = list(r,indices,e);
    656   }
     644    {rest[1] = list(r,indices,e);}
    657645  for(i=1; i<=m; i++)
    658646  {
     
    684672
    685673static proc mergepfd(list dec1, list dec2)
    686 // assumes dec1 is already sorted w.r.t. dp in the denominator exponents
    687 {
     674{
     675  // Note: assumes dec1 is already sorted w.r.t. dp in the denominator exponents
    688676  int n1=size(dec1);
    689677  int n2=size(dec2);
     
    717705    }
    718706    if(is_dp_smaller(dec1[a][2], dec1[a][3], entry[2], entry[3]))
    719     {         // numerator exponent vector of dec1[a] is dp-smaller than entry
    720       dec1=insert(dec1,entry,a); n1++;
    721     }
     707      {dec1=insert(dec1,entry,a); n1++;}
    722708    else{if(entry[2]==dec1[a][2] && entry[3]==dec1[a][3])
    723709    {
    724       dec1[a][1] = dec1[a][1] + entry[1]; //same denominator: add numerators
     710      dec1[a][1] = dec1[a][1] + entry[1];   //same denominator: add numerators
    725711      if(dec1[a][1]==0)
    726       {
    727         dec1 = delete(dec1,a); n1--;
    728       }
     712        {dec1 = delete(dec1,a); n1--;}
    729713    }
    730714    else
    731     {
    732       dec1=insert(dec1,entry,a-1); n1++;
    733     }}
     715      {dec1=insert(dec1,entry,a-1); n1++;}}
    734716  }
    735717  return(dec1);
     
    758740"USAGE:   checkpfd(list(f,g),dec[,N,C]);   f,g poly, dec list, N,C int
    759741RETURN:   0 or 1
    760 PURPOSE:  tests for (mathematical) equality of f/g and a partial fraction
     742PURPOSE:  test for (mathematical) equality of f/g and a partial fraction
    761743          decomposition dec. The list dec has to have the same structure as the
    762744          output of @ref{pfd}.
    763745       @* The denominator g can also be given in factorized form as a list of
    764746          an ideal of irreducible non constant polynomials and an intvec of
    765           exponents. (g=list(ideal(0),intvec(0:0)) represents a denominator
    766           of 1.)
     747          exponents. This can save time since the first step in the algorithm is
     748          to factorize g. (a list of the zero-ideal and an empty intvec
     749          represents a denominator of 1.)
    767750       @* By default the test is done (exactly) by bringing all terms of the
    768751          decomposition on the same denominator and comparing to f/g.
    769        @* If additional parameters N [, C] are given and if N>0, a probabilistic
    770           method is chosen: evaluation at N random points with coordinates
    771           between -C and C. This may be faster for big polynomials.
     752       @* If additional parameters N [, C] are given and if @code{N>0}, a
     753          probabilistic method is chosen: evaluation at N random points with
     754          coordinates between -C and C. This may be faster for big polynomials.
    772755SEE ALSO: pfd
    773756EXAMPLE:  example checkpfd; shows an example"
    774757{
    775   poly f=fraction[1];
    776   def g=fraction[2];
     758  poly f = fraction[1];
     759  def g = fraction[2];
    777760  if(size(#)>0)
    778761  {
     
    793776
    794777        if(typeof(g)=="poly")
    795         {
    796           val1 = number(substitute(g,vars,values));
    797         }
    798         else{if(typeof(g)=="list") //denominator given in factorized form
     778          {val1 = number(substitute(g,vars,values));}
     779        else{if(typeof(g)=="list") // denominator given in factorized form
    799780        {
    800781          val1 = number(1);
    801782          for(j=1; j<=size(g[1]); j++)
    802           {
    803             val1 = val1 * number(substitute(g[1][j]^g[2][j],vars,values));
    804           }
     783            {val1 = val1 * number(substitute(g[1][j]^g[2][j],vars,values));}
    805784        }
    806785        else {ERROR("wrong type for second argument, expected poly or list");}}
     
    811790        if(div_by_0) {continue;}
    812791        if(val1 != val2)
    813         {
    814           return(0);
    815         }
     792          {return(0);}
    816793      }
    817794      return(1);
     
    848825    num = dec[i][1];
    849826    for(j=1; j<=m; j++)
    850     {
    851       num = num * q[j]^(e[j]);
    852     }
     827      {num = num * q[j]^(e[j]);}
    853828    sum_of_numerators = sum_of_numerators + num;
    854829  }
    855 
    856830  // the decomposition is now equal to sum_of_numerators/product(q[i]^e_max[i]) (i from 1 to imax)
    857831  // now: check if this equals f/g:
    858 
    859832  if(typeof(g)=="poly")
    860833  {
     
    864837    num = f/fact[1][1];
    865838  }
    866   else{if(typeof(g)=="list") //denominator given in factorized form
     839  else{if(typeof(g)=="list")    //denominator given in factorized form
    867840  {
    868841    ideal q_g = g[1];
     
    877850  for(i=1; i<=m_g; i++)
    878851  {
    879 
    880852    j=0;
    881853    for(k=1; k<=m; k++)
     
    884856      if(c*q_g[i]==q[k]) {j=k; break;}
    885857    }
    886 
    887858    if(j==0)
    888     {
    889       sum_of_numerators = sum_of_numerators*q_g[i]^e_g[i];
    890     }
     859      {sum_of_numerators = sum_of_numerators*q_g[i]^e_g[i];}
    891860    else
    892861    {
    893       num = num*(c^e_g[i]); //fix lead coefficient
    894 
     862      num = num*(c^e_g[i]);    //fix lead coefficient
    895863      expon = e_g[i]-e_max[j];
    896864      if(expon>0)
    897       {
    898         sum_of_numerators = sum_of_numerators*q[j]^expon;
    899       }
     865        {sum_of_numerators = sum_of_numerators*q[j]^expon;}
    900866      else{if(expon<0)
    901       {
    902         num = num*q[j]^(-expon);
    903       }}
     867        {num = num*q[j]^(-expon);}}
    904868    }
    905869  }
     
    933897"USAGE:   evaluatepfd(dec, values[, mode]);   dec list, values ideal, mode int
    934898RETURN:   the number gotten by substituting the numbers generating the ideal
    935           values for the variables in the partial fraction decomposition dec.
    936           The list dec has to have the same structure as the output of @ref{pfd}.
    937        @* mode = 1: raise Error in case of division by 0 (default)
    938        @* mode = 2: return a second integer which is 1 if the denominator
     899          @code{values} for the variables in the partial fraction decomposition
     900          @code{dec}. The list @code{dec} has to have the same structure as the
     901          output of @ref{pfd}.
     902       @* @code{mode=1}: raise Error in case of division by 0 (default)
     903       @* @code{mode=2}: return a second integer which is 1 if the denominator
    939904          becomes 0, and 0 otherwise.
    940905SEE ALSO: pfd
     
    944909  if(size(#)>0) {mode = #[1];}
    945910
    946   ideal vars = maxideal(1); // ideal generate by ring variables
     911  ideal vars = maxideal(1); // ideal generated by ring variables
    947912  number val=0;
    948913  number denom;
     
    978943
    979944  displaypfd_long(dec);
     945  // evaluation at x=2, y=1:
    980946  ideal values = 2,1;
    981   evaluate(dec,values);
    982   // compare: f(2,1)/g(2,1) = (2+2*1)/(2^2-1^1) = 4/7
     947  evaluatepfd(dec,values);
     948
     949  // compare: f(2,1)/g(2,1) = (2+2*1)/(2^2-1^1) = 4/3
    983950}
    984951
     
    987954  int n=size(l);
    988955  for(int i=1; i<=n; i++)
    989   {
    990     if(entry==l[i]) {return(i);}
    991   }
     956    {if(entry==l[i]) {return(i);}}
    992957  return(0);
    993958}
     
    995960proc displaypfd(list dec)
    996961"USAGE:   displaypfd(dec);   dec list
    997 PURPOSE:  print a partial fraction decomposition dec in a readable way.
    998           The list dec has to have the same structure as the output of @ref{pfd}.
     962PURPOSE:  print a partial fraction decomposition @code{dec} in a readable way.
     963          The list @code{dec} has to have the same structure as the output of
     964          @ref{pfd}.
    999965SEE ALSO: pfd, displaypfd_long, getStringpfd, getStringpfd_indexed
    1000966EXAMPLE:  example displaypfd; shows an example"
     
    1032998    printf("q%s = %s",i,q[i]);
    1033999  }
    1034   printf("(%s terms)%n", size(dec), 0);
     1000  if(size(dec)==1) {printf("(%s terms)%n", 1, 0);}
     1001  else {printf("(%s terms)%n", size(dec), 0);}
    10351002}
    10361003example
     
    10571024  print("  "+getStringFraction(dec[1],q));
    10581025  for(int i=2; i<=size(dec); i++)
    1059   {
    1060     print("+ "+getStringFraction(dec[i],q));
    1061   }
    1062   printf("(%s terms)%n", size(dec), 0);
     1026    {print("+ "+getStringFraction(dec[i],q));}
     1027  if(size(dec)==1) {printf("(%s terms)%n", 1, 0);}
     1028  else {printf("(%s terms)%n", size(dec), 0);}
    10631029}
    10641030example
     
    10771043proc getStringpfd(list dec)
    10781044"USAGE:   getStringpfd(dec);   dec list
    1079 PURPOSE:  turn a partial fraction decomposition dec into one string. The list
    1080           dec has to have the same structure as the output of @ref{pfd}.
     1045PURPOSE:  turn a partial fraction decomposition @code{dec} into one string. The
     1046          list @code{dec} has to have the same structure as the output of
     1047          @ref{pfd}.
    10811048SEE ALSO: pfd, getStringpfd_indexed, displaypfd, displaypfd_long
    10821049EXAMPLE:  example getStringpfd; shows an example"
     
    10861053  string s = getStringFraction(dec[1],q);
    10871054  for(int i=2; i<=size(dec); i++)
    1088   {
    1089     s = s+" + "+getStringFraction(dec[i],q);
    1090   }
     1055    {s = s+" + "+getStringFraction(dec[i],q);}
    10911056  return(s);
    10921057}
     
    11071072proc getStringpfd_indexed(list dec)
    11081073"USAGE:   getStringpfd_indexed(dec);   dec list
    1109 PURPOSE:  turn a partial fraction decomposition dec into one string, writing the
    1110           denominator factors just as q1,q2,... . The list dec has to have the
    1111           same structure as the output of @ref{pfd}.
     1074PURPOSE:  turn a partial fraction decomposition @code{dec} into one string,
     1075          writing the denominator factors just as @code{q1},@code{q2},... . The
     1076          list @code{dec} has to have the same structure as the output of
     1077          @ref{pfd}.
    11121078SEE ALSO: pfd, getStringpfd, displaypfd, displaypfd_long
    11131079EXAMPLE:  example getStringpfd_indexed; shows an example"
     
    11161082  string s = getStringFraction_indexed(dec[1]);
    11171083  for(int i=2; i<=size(dec); i++)
    1118   {
    1119     s = s+" + "+getStringFraction_indexed(dec[i]);
    1120   }
     1084    {s = s+" + "+getStringFraction_indexed(dec[i]);}
    11211085  return(s);
    11221086}
     
    11781142PURPOSE:  read matrix of rational functions from a txt-file and turn it into a
    11791143          matrix (i.e. a list of lists) of pairs of polynomials (numerators and
    1180           denominators). The string file is the [directory +] name of the file,
    1181           r is the ring chosen for the numerator/denominator polynomials.
    1182        @* The input file should be a list of lists using the brackets \"{\",
    1183           \"}\" and \",\" as seperation, e.g.
    1184        @* \"{{(x+y)/(x^2-x*y), -(x^2*y+1)/(y), x^2}, {(x+1)/y, y/x, 0}}\".
     1144          denominators). The string @code{file} should be the [directory +] name
     1145          of the file in the form \"@code{<path-to-file>/<filename>.txt}\".
     1146       @* The input file should be a list of lists separated by the characters
     1147          \"{\", \"}\" and \",\". Example:
     1148       @* \"{{(x+y)/(x^2-x*y), -(x^2*y+1)/(y), x^2}, {(x+1)/y, y/x, 0}}\"
     1149       @* Each rational function has to be an expression of the form \"a\",
     1150          \"(a)/(b)\", \"(b)^(-n)\" or \"(a)*(b)^(-n)\" where \"a\",\"b\" stand
     1151          for polynomials (i.e. strings, that can be parsed as a polynomial with
     1152          the @code{execute} command) and \"n\" stands for a positive integer. A
     1153          minus sign \"-\" followed by such an expression is also allowed.
     1154       @* IMPORTANT: The strings \"a\",\"b\" must NOT contain the symbol \"/\".
     1155          (So in case the coefficient field is the rationals, all denominators
     1156          in the coefficients of numerator and denominator polynomials should be
     1157          cleared.)
    11851158       @* The file should contain less than 2^31 characters (filesize < 2 GB).
    1186           For bigger files the matrix should be split row-wise in multiple
    1187           matrices and saved in different files. A list of the filenames
    1188           (in the right order) can then be given as first argument instead.
     1159          For bigger files the matrix should be split row-wise into multiple
     1160          matrices and saved in different files (each smaller than 2 GB). A list
     1161          of the filenames (in the right order) can then be given as first
     1162          argument instead.
    11891163       @* Also the basering has to match the variable names used in the
    11901164          input file(s).
    1191        @* mode = 1: save result to an ssi file of the same name (default)
    1192        @* mode = 2: return result
    1193        @* mode = 3: both save to ssi file and return result
     1165       @* @code{mode=1} (default): save result to an ssi-file of the same name
     1166       @* @code{mode=2}: return result
     1167       @* @code{mode=3}: save to ssi-file AND return result
    11941168SEE ALSO: pfdMat"
    11951169{
     
    11981172  if(!defined(basering))
    11991173    {ERROR("no basering defined!");}
    1200   int left__,right__,pos1__,pos2__,tmp__,i__,j__,t__,tt__;
     1174  int left__,right__,pos1__,mid__,pos2__,tmp__,i__,j__,k__,t__,tt__,depth__;
    12011175  int mode__=1;
    12021176  if(size(#)>0) {mode__=#[1];}
     
    12071181    for(i__=1;i__<=n;i__++)
    12081182    {
    1209       printf("  file %s of %s:",i__,n);
     1183      dbprint(sprintf("  file %s of %s:",i__,n));
    12101184      if(find(file__[i__],".txt")==0) {ERROR("wrong file type, expected txt");}
     1185      printlevel = printlevel+1;
    12111186      mat__ = mat__ + readInputTXT(file__[i__],2);
     1187      printlevel = printlevel-1;
    12121188    }
    12131189
     
    12151191
    12161192    string filename__ = file__[1][1,find(file__[1],".txt")-1];
    1217     print("  saving to file "+filename__+".ssi "); t__ = rtimer;
     1193    dbprint("  saving to file "+filename__+".ssi "); t__ = rtimer;
    12181194    write("ssi:w "+filename__+".ssi", mat__);
    1219     printf("  done! (%s ms)", rtimer-t__);
     1195    dbprint(sprintf("  done! (%s ms)", rtimer-t__));
    12201196
    12211197    if(mode__==3) {return(mat__);}
     
    12251201  if(find(file__,".txt")==0) {ERROR("wrong file type, expected txt");}
    12261202
    1227   print("  reading file "); t__=rtimer;
     1203  dbprint("  reading file "); t__=rtimer;
    12281204  string data__ = read(":r "+file__);
    1229   printf("  done! (%s ms)", rtimer-t__);
    1230 
    1231 
    1232   print("  processing input "); t__ = rtimer;
     1205  dbprint(sprintf("  done! (%s ms)", rtimer-t__));
     1206
     1207  dbprint("  processing input "); t__ = rtimer;
    12331208  left__ = find(data__,"{");
    12341209  right__ = find(data__,"}");
     
    12361211  if(left__<tmp__ && tmp__<right__) {left__ = tmp__;}
    12371212
    1238   int finished__;
     1213  int finished__,n__;
    12391214  poly p__,q__;
    12401215  list row__,mat__;
     1216  string s__,ss__;
    12411217  i__=0;
    12421218  while(1)
     
    12511227    {
    12521228      j__++;
     1229      s__ = "";
    12531230
    12541231      pos1__ = pos2__+1;
    1255       pos2__ = find(data__,"/",pos1__);
    1256       tmp__ = find(data__,",",pos1__);
    1257       if(pos2__==0 || pos2__>right__) //no denominator
    1258       {
    1259         if(tmp__==0 || tmp__>right__) {pos2__ = right__; finished__ = 1;}
    1260         else {pos2__ = tmp__;}
    1261         execute("p__=" + fixBrackets(data__[pos1__,pos2__-pos1__]));
    1262         q__=1;
    1263       }
    1264       else{if(tmp__>0 && tmp__<pos2__) // no denominator
    1265       {
    1266         pos2__ = tmp__;
    1267         execute("p__=" + fixBrackets(data__[pos1__,pos2__-pos1__]));
    1268         q__=1;
    1269       }
    1270       else
    1271       {
    1272         execute("p__=" + fixBrackets(data__[pos1__,pos2__-pos1__]));
    1273 
    1274         pos1__ = pos2__+1;
    1275         pos2__ = tmp__;
    1276         if(pos2__==0 || pos2__>right__)
     1232      pos2__ = find(data__,",",pos1__);
     1233      if(pos2__==0 || pos2__>right__)    // end of row
     1234      {
     1235        finished__ = 1;
     1236        pos2__ = right__;
     1237      }
     1238      s__ = data__[pos1__,pos2__-pos1__];
     1239      mid__ = find(s__,"/");
     1240      if(mid__==0)
     1241      {
     1242        tmp__ = find(s__,"^(-");
     1243        if(tmp__==0)    //no denominator
    12771244        {
    1278           pos2__=right__;
    1279           finished__=1;
     1245          execute("p__=" + s__);
     1246          q__=1;
     1247          row__[j__] = list(p__,q__);
     1248          continue;
    12801249        }
    1281         execute("q__=" + fixBrackets(data__[pos1__,pos2__-pos1__]));
    1282       }}
     1250        else    // denominator is given by a negative exponent
     1251        {
     1252          if(find(s__,"^(-",tmp__+1)>0)
     1253            {ERROR(sprintf("invalid syntax in (%s,%s)-th entry:"
     1254                              +" more than one negative exponent",i__,j__));}
     1255          for(k__=tmp__+3; s__[k__]!=")"; k__++) {}
     1256          execute("n__=" + s__[tmp__+3,k__-tmp__-3]);
     1257          while(k__<size(s__))
     1258          {
     1259            k__++;
     1260            if(s__[k__]!=" ")
     1261            {ERROR(sprintf("invalid syntax in (%s,%s)-th entry",i__,j__));}
     1262          }
     1263          s__ = s__[1,tmp__-1];
     1264          depth__=0;
     1265          for(k__=tmp__-1; 1; k__--)
     1266          {
     1267            ss__ = s__[k__];
     1268            if(ss__ == ")") {depth__++;}
     1269            else{if(ss__ == "(") {depth__--;}}
     1270            if(depth__==0) {break;}
     1271          }
     1272          if(k__>1)
     1273          {
     1274            while(1)
     1275            {
     1276              if(k__==0)
     1277                {ERROR(sprintf("invalid syntax in (%s,%s)-th entry",i__,j__));}
     1278              k__--;
     1279              ss__ = s__[k__];
     1280              if(ss__=="*" || ss__==" " || ss__=="-") {break;}
     1281            }
     1282          }
     1283          s__ = s__[1,k__] + "1/" + s__[k__+1,size(s__)-k__] + "^" + string(n__);
     1284          mid__ = k__+2;    // position of the character "/"
     1285        }
     1286      }
     1287      if(find(s__,"/",mid__+1)>0)
     1288        {ERROR(sprintf("invalid syntax in (%s,%s)-th entry:%n"
     1289        +"no '/' allowed in the string representing the polynomials",i__,j__,0));}
     1290
     1291      execute("p__=" + fixBrackets(s__[1,mid__-1]));
     1292      execute("q__=" + fixBrackets(s__[mid__+1,size(s__)-mid__]));
    12831293
    12841294      row__[j__] = list(p__,q__);
    12851295    }
    1286 
    1287     mat__[i__] = row__; // append row to matrix
    1288     printf("    row %s done! (%s ms)",i__,rtimer-tt__);
     1296    mat__[i__] = row__;    // append row to matrix
     1297    dbprint(sprintf("    row %s done! (%s ms)",i__,rtimer-tt__));
    12891298
    12901299    left__ = find(data__,"{",right__);
     
    12921301    right__ = find(data__,"}",left__);
    12931302  }
    1294   printf("  done! (%s ms)", rtimer-t__);
     1303  dbprint(sprintf("  done! (%s ms)", rtimer-t__));
    12951304
    12961305  if(mode__==2) {return(mat__);}
    12971306
    12981307  string filename__ = file__[1,find(file__,".txt")-1];
    1299   print("  saving to file "+filename__+".ssi "); t__ = rtimer;
     1308  dbprint("  saving to file "+filename__+".ssi "); t__ = rtimer;
    13001309  write("ssi:w "+filename__+".ssi", mat__);
    1301   printf("  done! (%s ms)", rtimer-t__);
     1310  dbprint(sprintf("  done! (%s ms)", rtimer-t__));
    13021311
    13031312  if(mode__==3) {return(mat__);}
    13041313}
    1305 /* new version:
    1306 {
    1307   system("--ticks-per-sec",1000);
    1308   if(!defined(basering))
    1309     {ERROR("no basering defined!");}
    1310   int left__,right__,pos1__,pos2__,tmp__,i__,j__,k__,d__,t__,tt__;
    1311   int mode__=1;
    1312   if(size(#)>0) {mode__=#[1];}
    1313   if(typeof(file__)=="list") // list of filenames given --> apply function to each
    1314   {                          // file and concatenate the resulting matrices
    1315     int n = size(file__);
    1316     list mat__ = list();
    1317     for(i__=1;i__<=n;i__++)
    1318     {
    1319       printf("  file %s of %s:",i__,n);
    1320       if(find(file__[i__],".txt")==0) {ERROR("wrong file type, expected txt");}
    1321       mat__ = mat__ + readInputTXT(file__[i__],2);
    1322     }
    1323 
    1324     if(mode__==2) {return(mat__);}
    1325 
    1326     string filename__ = file__[1][1,find(file__[1],".txt")-1];
    1327     print("  saving to "+filename__+".ssi "); t__ = rtimer;
    1328     write("ssi:w "+filename__+".ssi", mat__);
    1329     printf("  done! (%s ms)", rtimer-t__);
    1330 
    1331     if(mode__==3) {return(mat__);}
    1332   }
    1333   if(typeof(file__)!="string")
    1334     {ERROR("wrong type for first argument (expected string or list)");}
    1335   if(find(file__,".txt")==0) {ERROR("wrong file type, expected txt");}
    1336 
    1337   print("  reading input matrix as string from "+file__); t__=rtimer;
    1338   string data__ = read(":r "+file__);
    1339   printf("  done! (%s ms)", rtimer-t__);
    1340 
    1341 
    1342   print("  processing input "); t__ = rtimer;
    1343   left__ = find(data__,"{");
    1344   right__ = find(data__,"}");
    1345   tmp__ = find(data__,"{",left__+1);
    1346   if(left__<tmp__ && tmp__<right__) {left__ = tmp__;}
    1347 
    1348   int finished__;
    1349   poly p__,q__;
    1350   int depth__;
    1351   list mat__;
    1352   string s__,s1__,s2__;
    1353   for(i__=1;1;i__++)
    1354   {
    1355     tt__ = rtimer;
    1356     mat__[i__] = list();
    1357     pos2__ = left__;
    1358     finished__ = 0;
    1359     for(j__=1; not finished__; j__++)
    1360     {
    1361       pos1__ = pos2__+1;
    1362       pos2__ = find(data__,",",pos1__);
    1363       if(pos2__==0 || pos2__>right__) {finished__=1; pos2__=right__;}
    1364       s__ = data__[pos1__,pos2__-pos1__];
    1365       for(k__=1; s__[k__]==" "; k__++) {}
    1366       s__ = s__[k__,size(s__)-k__+1]; // remove spaces
    1367 
    1368       tmp__=find(s__,"/");
    1369       //printf("s: %s",s__);
    1370       if(tmp__>0)
    1371       {
    1372         depth__=0;
    1373         d__=0;
    1374         for(k__=tmp__+1;k__<=size(s__);k__++)
    1375         {
    1376           if(s__[k__]=="(") {depth__++; k__++; continue;}
    1377           if(s__[k__]==")") {depth__--; k__++; continue;}
    1378           if(s__[k__]=="/" && depth__<=d__) {tmp__=k__; d__=depth__;}
    1379         }
    1380         s1__ = s__[1,tmp__-1];
    1381         s2__ = s__[tmp__+1,size(s__)-tmp__];
    1382         depth__=0;
    1383         for(k__=size(s1__); k__>1; k__--)
    1384         {
    1385           if(s1__[k__]==")") {depth__++; k__--; continue;}
    1386           if(s1__[k__]=="(") {depth__--; k__--; continue;}
    1387           if(depth__==0 && (s1__[k__]=="+" || s1__[k__]=="-")) {tmp__=0; break}
    1388         }
    1389         depth__=0;
    1390         for(k__=1; k__<=size(s2__); k__++)
    1391         {
    1392           if(s2__[k__]=="(") {depth__++; k__++; continue;}
    1393           if(s2__[k__]==")") {depth__--; k__++; continue;}
    1394           if(depth__==0 && (s2__[k__]=="+" || s2__[k__]=="-" || s2__[k__]=="*")) {tmp__=0; break}
    1395         }
    1396         //printf("s1: %s,%ns2: %s",s1__,s2__);
    1397       }
    1398       if(tmp__==0) // no denominator
    1399       {
    1400         execute("p__=" + s__);
    1401         q__=1;
    1402       }
    1403       else // s1__ = numerator, s2__ = denominator
    1404       {
    1405         execute("p__=" + fixBrackets(s1__));
    1406         execute("q__=" + fixBrackets(s2__));
    1407         if(deg(q__)==0) {p__=p__/q__; q__=1;}
    1408       }
    1409       ;
    1410 
    1411       mat__[i__][j__] = list(p__,q__);
    1412     }
    1413     printf("    row %s done! (%s ms)",i__,rtimer-tt__);
    1414 
    1415     left__ = find(data__,"{",right__);
    1416     if(left__==0) {break;}
    1417     right__ = find(data__,"}",left__);
    1418   }
    1419   printf("  done! (%s ms)", rtimer-t__);
    1420 
    1421   if(mode__==2) {return(mat__);}
    1422 
    1423   string filename__ = file__[1,find(file__,".txt")-1];
    1424   print("  saving to "+filename__+".ssi "); t__ = rtimer;
    1425   write("ssi:w "+filename__+".ssi", mat__);
    1426   printf("  done! (%s ms)", rtimer-t__);
    1427 
    1428   if(mode__==3) {return(mat__);}
    1429 } */
    1430 
    14311314
    14321315static proc fixBrackets(string data)
     
    14891372          file string, dotest,ignore_nonlin,output_mode,parallelize int
    14901373PURPOSE:  apply @code{pfd} to all entries of a matrix of rational functions
    1491           saved in a file. The string file is the [directory +] name of the
    1492           file.
     1374          saved in a txt-file. The string @code{file} should be the
     1375          [directory +] name of the file.
    14931376       @* The input file can either be a txt-file or an ssi-file created with
    14941377          @code{readInputTXT}. In case of a txt-file, the base ring has to match
    14951378          and the matrix has to be in the same format specified in
    1496           @ref{readInputTXT}.
    1497           Also, txt-files that are bigger than 2 GB should be split as described
    1498           for @code{readInputTXT} and a list of the filenames can be given as
    1499           first argument instead.
    1500        @* The result is written in two txt-files: once with the denominators
    1501           written out and once with indexed denominator factors q1,q2,... .
    1502           The polynomials q1,q2,... are saved in a separate txt-file.
     1379          @ref{readInputTXT}. Also, txt-files that are bigger than 2 GB should
     1380          be split as described for @code{readInputTXT} and a list of the
     1381          filenames can be given as first argument instead.
     1382       @* The result is saved in multiple txt- (and ssi-) files (see below)
     1383          within the directory of the input file.
    15031384       @* Also a logfile is created, which protocols the memory used and the
    1504           runtimes of for each matrix entry in real-time.
     1385          runtimes of @code{pfd} for each matrix entry in real-time.
    15051386
    15061387       @* There are also 4 optional arguments:
    15071388       @* If @code{dotest} is nonzero, test the results with checkpfd:
    1508             dotest<0: exact test (may be slow),
    1509             dotest>0: do this amount of probabilistic tests for each entry
    1510                       (see @ref{checkpfd}).
    1511        @* (default: @code{dotest=-1})
    1512        @* If @code{ignore_nonlin} is nonzero, for each denominator, the
    1513           nonlinear denominator factors are saved into a seperate file and
    1514           removed before applying @code{pfd}. The nonlinear factors thus do not
    1515           appear in the output files. (So the input matrix is equal to the
    1516           output matrix divided element-wise by the matrix containing the
    1517           nonlinear factors.)
    1518           (default: @code{ignore_nonlin=1})
    1519        @* If @code{output_mode} is nonzero, the input and result are also
    1520           saved to ssi-files containing matrices (as list of lists) of the input
    1521           rational functions (as lists of numerator and denominator poly) and
    1522           decompositions respectively. For the ssi-file containing the result,
    1523           the first entry is an ideal @code{q} of denominator factors and the
    1524           second entry is a matrix (as list of lists) containing the
    1525           decompositions, each of which is a list of terms, where a term is
    1526           represented as in the result of @ref{pfd} by a list @code{l} containing
    1527        @*   1) the numerator polynomial
    1528        @*   2) an intvec giving the indices @code{i} for which @code{q[i]}
    1529                occurs as a factor in the denominator
    1530        @*   3) an intvec containing the exponents of those irreducible factors.
    1531        @* Also, the results of @code{pfd} will be saved directly in ssi-files
    1532           named pfd_results_i_j where i,j are the matrix indices in the current
    1533           directory. (default: @code{output_mode=0})
    1534        @* If @code{parallelize} is nonzero, the decompositions are calculated in
    1535           parallel using @ref{parallel_lib}. (default: @code{parallelize=1})
     1389       @* @code{dotest<0} (default): exact test (may be slow),
     1390       @* @code{dotest>0}: do this amount of probabilistic tests for each entry
     1391                           (see @ref{checkpfd}).
     1392
     1393       @* If @code{ignore_nonlin} is nonzero (default), for each denominator,
     1394          the nonlinear factors in the factorization are removed before applying
     1395          @code{pfd} (and added back in in the output files).
     1396
     1397       @* If @code{parallelize} is nonzero (default), the decompositions are
     1398          calculated in parallel using @ref{parallel_lib}.
     1399
     1400       @* The parameter @code{output_mode} controls the output files created:
     1401       @* @code{output_mode=1} (default): The result consists of two files:
     1402          @code{<filename>_pfd_indexed.txt} contains the matrix of all
     1403          decompositions (as list of lists separated by the characters \"{\",
     1404          \"}\" and \",\") where all the denominators are written in factorized
     1405          form depending on irreducible factors @code{q1}, @code{q2}, ... .
     1406          The file @code{<filename>_denominator_factors.txt} lists all the
     1407          polynomials @code{q1}, @code{q2}, ... .
     1408       @* @code{output_mode=2}: Additionally to mode 1, the file
     1409          @code{<filename>_pfd.txt} is created, which also contains the matrix
     1410          of decompositions but the factors in the denominators are written out.
     1411       @* @code{output_mode=3}: Additionally to mode 2, the result and some
     1412          intermediate results are saved as SINGULAR objects in ssi-files:
     1413       @* @code{<filename>.ssi}: contains the result of @code{readInputTXT} in
     1414          case a txt-file was given as input.
     1415       @* @code{<filename>_factorized_denominators.ssi}: like the first file,
     1416          but the denominators are saved in factorized form, that is as a list
     1417          of an ideal of irreducible non constant polynomials and an intvec of
     1418          exponents.
     1419       @* @code{<filename>_linear_part.ssi} (only if @code{ignore_nonlin} is
     1420          nonzero): like the previous file, but all the irreducible denominator
     1421          factors are removed
     1422       @* @code{<filename>_non_linear_factors.ssi} (only if @code{ignore_nonlin}
     1423          is nonzero): a list of an ideal @code{p} generated by irreducible
     1424          polynomials and a matrix (list of lists) of the nonlinear denominator
     1425          factors of each entry of the input matrix. These are represented as
     1426          lists of an intvec of indices @code{i} for which @code{p[i]} occurs
     1427          as a (nonlinear) factor in the denominator and an intvec containing
     1428          the exponents of those factors.
     1429       @* @code{<filename>_pfd.ssi}: a list, where the first entry is an ideal
     1430          @code{q} of denominator factors and the second entry is a matrix (as
     1431          list of lists) containing the decompositions, each of which is a list
     1432          of terms, where a term is represented as in the result of @ref{pfd}
     1433          by a list containing
     1434       @* 1) the numerator polynomial
     1435       @* 2) an intvec of indices @code{i} for which @code{q[i]} occurs
     1436             as a factor in the denominator
     1437       @* 3) an intvec containing the exponents of those irreducible factors.
     1438       @* IMPORTANT: If @code{ignore_nonlin} is nonzero, this file contains the
     1439          decompositions of the entries of the matrix in
     1440          @code{<filename>_linear_part.ssi}. Thus the nonlinear factors, are
     1441          NOT contained in this file.
     1442       @* @code{output_mode=4}: Additionally to mode 3, the direct output of
     1443          each call of @code{pfd} is saved in separate ssi-files called
     1444          @code{pfd_results_i_j.ssi} where i,j are the matrix indices. This
     1445          creates a lot of files, but may be useful in case the algorithm does
     1446          not terminate in time for every matrix entry. Other than the files
     1447          created in mode 1-3, these files are saved in the current directory,
     1448          rather than the directory of the input file.
    15361449SEE ALSO: readInputTXT, pfd, checkpfd, checkpfdMat"
    15371450{
    15381451  system("--ticks-per-sec",1000);
    15391452  short = 0;
    1540   int dotest,ignore_nonlin,output_mode,parallelize = -1,1,0,1;
     1453  int dotest,ignore_nonlin,output_mode,parallelize = -1,1,1,1;
    15411454  if(size(#)>0) {dotest = #[1];}
    15421455  if(size(#)>1) {ignore_nonlin = #[2];}
     
    15471460  list arguments,results;
    15481461
    1549   print(newline+"reading data "); int t0=rtimer;
     1462  dbprint(newline+"reading data "); int t0=rtimer;
    15501463
    15511464  if(typeof(infile)=="list")
    15521465  {
     1466    printlevel = printlevel+1;
    15531467    if(output_mode>2) {list mat = readInputTXT(infile,3);}
    15541468    else {list mat = readInputTXT(infile,2);}
     1469    printlevel = printlevel-1;
    15551470    int pos=find(infile[1],".txt");
    15561471    string filename = infile[1][1,pos-1];
     
    15611476    if(pos!=0)
    15621477    {
     1478      printlevel = printlevel+1;
    15631479      if(output_mode>2) {list mat = readInputTXT(infile,3);}
    15641480      else {list mat = readInputTXT(infile,2);}
     1481      printlevel = printlevel-1;
    15651482    }
    15661483    else
    15671484    {
    15681485      pos=find(infile,".ssi");
    1569       if(pos!=0)
    1570       {
    1571         list mat = read("ssi:r "+infile);
    1572       }
     1486      if(pos!=0) {list mat = read("ssi:r "+infile);}
    15731487      else {ERROR("invalid file type, expected ssi or txt");}
    15741488    }
     
    15791493  int n = size(mat);
    15801494  int m = size(mat[1]);
    1581   printf("done! (%s ms)",rtimer-t0);
    1582 
    1583   if(typeof(mat[1][1][2])!="list") // denominators are not yet factorized
    1584   {
    1585     print("factorizing the denominators "); t0=rtimer;
     1495  dbprint(sprintf("done! (%s ms)",rtimer-t0));
     1496
     1497  if(typeof(mat[1][1][2])!="list")    // denominators are not yet factorized
     1498  {
     1499    dbprint("factorizing the denominators "); t0=rtimer;
     1500    printlevel = printlevel+1;
    15861501    mat = FactDenom(mat);
     1502    printlevel = printlevel-1;
    15871503    if(output_mode>2)
    15881504      {write("ssi:w "+filename+"_factorized_denominators.ssi",mat);}
    1589     printf("done! (%s ms)",rtimer-t0);
     1505    dbprint(sprintf("done! (%s ms)",rtimer-t0));
    15901506  }
    15911507
    15921508  if(ignore_nonlin)
    15931509  {
    1594     print("removing nonlinear denominator factors before pfd is applied");
    1595     //       +" (%s_non_linear_factors[_indexed].txt):", filename);
    1596     list denom_factors;
     1510    dbprint("removing nonlinear denominator factors before pfd is applied");
     1511    list nonlin_denom_factors;
    15971512    ideal p;
    1598     mat, denom_factors, p = removeNonlinearFactors(mat, filename);
     1513    printlevel = printlevel+1;
     1514    mat, nonlin_denom_factors, p = removeNonlinearFactors(mat, filename);
     1515    printlevel = printlevel-1;
    15991516    if(output_mode>2)
    16001517    {
    1601       print("saving nonlinear factors to "+filename+"_non_linear_factors.ssi "); t0 = rtimer;
    1602       write("ssi:w "+filename+"_non_linear_factors.ssi", denom_factors);
    1603       printf("done! (%s ms)",rtimer-t0);
    1604       print("saving input matrix without the nonlinear factors to "+filename+"_linear_part.ssi "); t0 = rtimer;
     1518      dbprint("saving nonlinear factors to "+filename+"_non_linear_factors.ssi ");
     1519      t0 = rtimer;
     1520      write("ssi:w "+filename+"_non_linear_factors.ssi", list(p,nonlin_denom_factors));
     1521      dbprint(sprintf("done! (%s ms)",rtimer-t0));
     1522      dbprint("saving input matrix without the nonlinear factors to "
     1523              +filename+"_linear_part.ssi ");
     1524      t0 = rtimer;
    16051525      write("ssi:w "+filename+"_linear_part.ssi", mat);
    1606       printf("done! (%s ms)",rtimer-t0);
    1607     }
    1608     //ideal p = saveNonlinearFactorsTXT(denom_factors, filename);
    1609     //filename = filename + "_linear_part";
     1526      dbprint(sprintf("done! (%s ms)",rtimer-t0));
     1527    }
    16101528  }
    16111529
    16121530  if(parallelize)
    16131531  {
    1614     print("creating tasks "); t0=rtimer;
     1532    dbprint("creating tasks "); t0=rtimer;
    16151533    write(":w "+filename+"_pfdMat_logfile.txt","finished matrix entries with runtimes "
    16161534          +"(calculated in parallel on "+string(getcores())+" cores):");
     
    16241542      }
    16251543    }
    1626     printf("done! (%s ms)",rtimer-t0);
    1627 
    1628     print("applying pfd to each matrix entry "); t0 = rtimer;
     1544    dbprint(sprintf("done! (%s ms)",rtimer-t0));
     1545
     1546    dbprint("applying pfd to each matrix entry "); t0 = rtimer;
    16291547    results = parallelWaitAll("pfdWrap",arguments);
    16301548    arguments = list();
    1631     printf("done! (%s ms)",rtimer-t0);
     1549    dbprint(sprintf("done! (%s ms)",rtimer-t0));
    16321550
    16331551    write(logfile,"decomposition: "+string(rtimer-t0)+" ms and "+string(memory(2))
    16341552                   +" Byte Memory max. (after calling pfd on each matrix entry)");
    16351553
    1636     print("writing results in matrix shape "); t0 = rtimer;
     1554    dbprint("writing results in matrix shape "); t0 = rtimer;
    16371555    list dec_mat;
    16381556    for(i=1; i<=n; i++)
     
    16461564    }
    16471565    results = list();
    1648     printf("done! (%s ms)",rtimer-t0);
     1566    dbprint(sprintf("done! (%s ms)",rtimer-t0));
    16491567  }
    16501568  else
    16511569  {
    1652     print("applying pfd to each matrix entry "); t0=rtimer;
     1570    dbprint("applying pfd to each matrix entry "); t0=rtimer;
    16531571    write(":w "+filename+"_pfdMat_logfile.txt",
    16541572          "finished matrix entries with runtimes (no parallelization):");
     
    16631581      }
    16641582    }
    1665     printf("done! (%s ms)",rtimer-t0);
     1583    dbprint(sprintf("done! (%s ms)",rtimer-t0));
    16661584    write(logfile,"decomposition: "+string(rtimer-t0)+" ms and "+string(memory(2))
    16671585          +" Byte Memory max. (after calling pfd on each matrix entry)");
    16681586 }
    16691587
    1670   print("making one single list of denominator factors "); t0 = rtimer;
     1588  dbprint("making one single list of denominator factors "); t0 = rtimer;
    16711589  ideal q,new_q;
    16721590  intvec dict;
     
    16911609      {
    16921610        if(size(dec_mat[i][j][k][2])>0)
    1693         {
    1694           dec_mat[i][j][k][2] = intvec(dict[dec_mat[i][j][k][2]]);
    1695         }
    1696       }
    1697     }
    1698     printf("  row %s complete!",i);
    1699   }
    1700   printf("done! (%s ms)",rtimer-t0);
     1611          {dec_mat[i][j][k][2] = intvec(dict[dec_mat[i][j][k][2]]);}
     1612      }
     1613    }
     1614    dbprint(sprintf("  row %s complete!",i));
     1615  }
     1616  dbprint(sprintf("done! (%s ms)",rtimer-t0));
    17011617
    17021618  if(output_mode>2)
    17031619  {
    1704     print("saving result to "+filename+"_pfd_only_linear_factors.ssi "); t0 = rtimer;
    1705     write("ssi:w "+filename+"_pfd_only_linear_factors.ssi", list(q,dec_mat));
    1706     printf("done! (%s ms)",rtimer-t0);
     1620    dbprint("saving result to "+filename+"_pfd.ssi "); t0 = rtimer;
     1621    write("ssi:w "+filename+"_pfd.ssi", list(q,dec_mat));
     1622    dbprint(sprintf("done! (%s ms)",rtimer-t0));
    17071623  }
    17081624
     
    17141630    {
    17151631      for(j=1; j<=m; j++)
    1716         {
    1717           denom_factors[i][j][1] = denom_factors[i][j][1]+ind; // adjust indices
    1718         }
     1632      {
     1633        nonlin_denom_factors[i][j][1] = nonlin_denom_factors[i][j][1]+ind; // adjust indices
     1634      }
    17191635    }
    17201636  }
    17211637
    17221638  if(ignore_nonlin)
    1723     {printf("creating readable .txt-files (including the nonlinear factors again)");}
     1639    {dbprint(sprintf("creating readable .txt-files (including the nonlinear factors again)"));}
    17241640  else
    1725     {printf("creating readable .txt-files ");}
     1641    {dbprint(sprintf("creating readable .txt-files "));}
    17261642  t0 = rtimer;
    1727   print(" indexed ("+filename+"_pfd_indexed.txt):");
     1643  dbprint(" indexed ("+filename+"_pfd_indexed.txt):");
     1644  printlevel = printlevel+1;
    17281645  if(ignore_nonlin)
    1729     {saveResultTXT_indexed(dec_mat, filename+"_pfd_indexed", denom_factors);}
     1646    {saveResultTXT_indexed(dec_mat, filename+"_pfd_indexed", nonlin_denom_factors);}
    17301647  else {saveResultTXT_indexed(dec_mat, filename+"_pfd_indexed");}
     1648  printlevel = printlevel-1;
    17311649  for(i=1; i<=size(q); i++) {fprintf(qfile, "q%s = %s;", i, q[i]);}
    17321650  if(output_mode>1)
    17331651  {
    1734     print(" denominators written out ("+filename+"_pfd.txt):");
     1652    dbprint(" denominators written out ("+filename+"_pfd.txt):");
     1653    printlevel = printlevel+1;
    17351654    if(ignore_nonlin)
    1736       {saveResultTXT(dec_mat, q, filename+"_pfd", denom_factors);}
     1655      {saveResultTXT(dec_mat, q, filename+"_pfd", nonlin_denom_factors);}
    17371656    else {saveResultTXT(dec_mat, q, filename+"_pfd");}
    1738   }
    1739   printf("done! (%s ms)",rtimer-t0);
     1657    printlevel = printlevel-1;
     1658  }
     1659  dbprint(sprintf("done! (%s ms)",rtimer-t0));
    17401660
    17411661  if(dotest)
    17421662  {
    17431663    if(dotest<0)
    1744       {print("checking for correctness (exact test) ");}
     1664      {dbprint("checking for correctness (exact test) ");}
    17451665    else
    1746       {printf("checking for correctness (%s random evaluations per entry) ",dotest);}
     1666      {dbprint(sprintf("checking for correctness (%s random evaluations per entry) ",dotest));}
    17471667    t0 = rtimer;
    17481668    if(parallelize)
     
    17671687      }
    17681688    }
    1769 
    1770     printf("%s out of %s = %sx%s decompositions are correct! (%s ms)%n",
    1771             sum(results),n*m,n,m,rtimer-t0,0);
     1689    dbprint(sprintf("%s out of %s = %sx%s decompositions are correct! (%s ms)%n",
     1690            sum(results),n*m,n,m,rtimer-t0,0));
    17721691    write(logfile,"checking for correctness: "+string(rtimer-t0)+" ms and "
    17731692          +string(memory(2))+" Byte Memory max. (at the end of pfdMat), "
     
    18021721  {
    18031722    t = rtimer;
    1804     for(i=1; i<=n; i++) // create argument list
     1723    for(i=1; i<=n; i++)   // create argument list
    18051724    {
    18061725      for(j=1; j<=m; j++)
    1807       {
    1808         arguments[m*(i-1)+j] = list(denom[i][j]);
    1809       }
     1726        {arguments[m*(i-1)+j] = list(denom[i][j]);}
    18101727    }
    18111728
     
    18131730    arguments = list();
    18141731
    1815     for(i=1; i<=n; i++) // update q, mat, denom
     1732    for(i=1; i<=n; i++)   // update q, mat, denom
    18161733    {
    18171734      for(j=1; j<=m; j++)
     
    18821799      }
    18831800    }
    1884 
    18851801    timeout = timeout*2;
    1886 
    1887     printf("  %s out of %s denominators factorized completely (%s ms)",counter,n*m,rtimer-t);
     1802    dbprint(sprintf("  %s out of %s denominators factorized completely (%s ms)",
     1803                                                         counter,n*m,rtimer-t));
    18881804  }
    18891805  return(mat);
     
    19241840    k = size(s);
    19251841    s[k-1] = "}"; s[k] = ","; s = s+" ";
    1926     printf("  row %s done! (%s ms)",i,rtimer-t);
     1842    dbprint(sprintf("  row %s done! (%s ms)",i,rtimer-t));
    19271843  }
    19281844  k = size(s);
     
    19651881    k = size(s);
    19661882    s[k-1] = "}"; s[k] = ","; s = s+" ";
    1967     printf("  row %s done! (%s ms)",i,rtimer-t);
     1883    dbprint(sprintf("  row %s done! (%s ms)",i,rtimer-t));
    19681884  }
    19691885  k = size(s);
     
    19781894  int i,j,k,t0,ind;
    19791895  ideal p;
    1980   list denom_factors;
     1896  list nonlin_denom_factors;
    19811897  intvec factors,exponents;
    19821898  list fac;
     
    19841900  {
    19851901    t0 = rtimer;
    1986     denom_factors[i] = list();
     1902    nonlin_denom_factors[i] = list();
    19871903    for(j=1; j<=m; j++)
    19881904    {
     
    20001916            for(ind=1;1;ind++)
    20011917            {
    2002               //print("   "+string(ind));
    20031918              if(p[ind]==fac[1][k]) {break;}
    20041919              if(ind==size(p))
     
    20181933      }
    20191934      fractions[i][j][2] = fac;
    2020 
    2021       denom_factors[i][j] = list(factors,exponents);
    2022     }
    2023     printf("  row %s done! (%s ms)",i,rtimer-t0);
    2024   }
    2025   return(fractions, denom_factors, p);
     1935      nonlin_denom_factors[i][j] = list(factors,exponents);
     1936    }
     1937    dbprint(sprintf("  row %s done! (%s ms)", i, rtimer-t0));
     1938  }
     1939  return(fractions, nonlin_denom_factors, p);
    20261940}
    20271941
     
    20291943"USAGE:   checkpfdMat(input, output, denomFactors[, N, parallelize]);
    20301944          input,output,denomFactors string, N,parallelize int
    2031           checkpfdMat(input, output, denomFactors, nonlinFactors[, N, parallelize]);
    2032           input,output,denomFactors,nonlinFactors string, N,parallelize int
    20331945PURPOSE:  test the output files of @code{pfdMat} for correctness. Input and
    2034           output (indexed) txt-file names are given as strings. If nonlinear
    2035           denominator factors were ignored in the calculation (as described in
    2036           @ref{pfdMat}), the txt-file containing the nonlinear factors must be
    2037           given as fourth argument @code{nonlinFactors}. The output and
    2038           nonlinear Factors should be indexed and @code{denomFactors} is the
    2039           file containing the denominator factors @code{q1},@code{q2},... .
     1946          output (indexed) txt-files have to be given as strings in the form
     1947          \"@code{<path-to-file>/<filename>.txt}\". The output should be indexed
     1948          (that is the output file ending in @code{..._pfd_indexed.txt}) and
     1949          @code{denomFactors} has to be the file containing the denominator
     1950          factors @code{q1}, @code{q2}, ... (the txt-file ending in
     1951          @code{..._denominator_factors.txt}).
    20401952       @* As for @code{readInputTXT} and @code{pfdMat}, the basering has to
    20411953          match the variable names used in the input file, which has to be in
     
    20431955          than 2 GB have to be split as described for @code{readInputTXT} and a
    20441956          list of filenames can be given as first argument instead.
    2045        @* If an integer N is given, the test is done probabilistically by N
    2046           random substitutions for each entry of the matrix. If N is omitted,
    2047           the fractions in the decomposition will be expanded symbolically and
    2048           compared to the input (may be slower).
    2049        @* If @code{parallelize} is nonzero, the tests are run in parallel using
    2050           @ref{parallel_lib}. (default: @code{parallelize=1})
     1957       @* If a positive integer N is given, the test is done probabilistically by
     1958          evaluation at N random points for each entry of the matrix. If N is
     1959          nonpositive (default), the fractions in the decompositions will be
     1960          expanded symbolically and compared to the input (may be slower).
     1961       @* If @code{parallelize} is nonzero (default), the tests are run in
     1962          parallel using @ref{parallel_lib}.
    20511963       @* The result is printed and as in @code{pfdMat} a logfile is created
    20521964          showing the results for each matrix entry.
     
    20711983  if(size(#)>2) {ERROR("too many arguments");}
    20721984
    2073   print(newline+"reading input file:");
     1985  dbprint(newline+"reading input file:");
     1986  printlevel = printlevel+1;
    20741987  list frac = readInputTXT(input,2);
     1988  printlevel = printlevel-1;
    20751989  if(typeof(input)=="string") {string filename=input;}
    20761990  if(typeof(input)=="list")   {string filename=input[1];}
    20771991  filename = filename[1,find(filename,".txt")-1];
    20781992
    2079   print("factorizing the denominators "); t0=rtimer;
     1993  dbprint("factorizing the denominators "); t0=rtimer;
     1994  printlevel = printlevel+1;
    20801995  frac = FactDenom(frac);
    2081   printf("done! (%s ms)",rtimer-t0);
    2082 
    2083   print("reading output files:"); t0=rtimer;
    2084   print(" reading list of denominator factors from "+qfile);
     1996  printlevel = printlevel-1;
     1997  dbprint(sprintf("done! (%s ms)",rtimer-t0));
     1998
     1999  dbprint("reading output files:"); t0=rtimer;
     2000  dbprint(" reading list of denominator factors from "+qfile);
    20852001  ideal q;
    20862002  q = readQfileTXT(qfile);
    2087   //printf("denominator factors: %s",q);
    2088   print(" done!");
    2089 
    2090   print(" reading (indexed) output decompositions ");
     2003  dbprint(" done!");
     2004
     2005  dbprint(" reading (indexed) output decompositions ");
    20912006  list dec,nonlin;
     2007  printlevel = printlevel+1;
    20922008  dec,nonlin = readOutputTXT(output);
     2009  printlevel = printlevel-1;
    20932010
    20942011  if(parallelize)
    2095     {printf("done! (%s ms)%n%ncreating tasks",rtimer-t0,0); t0=rtimer;}
     2012    {dbprint(sprintf("done! (%s ms)%n%ncreating tasks",rtimer-t0,0)); t0=rtimer;}
    20962013  else
    20972014  {
    2098     printf("done! (%s ms)%n%n",rtimer-t0,0);
     2015    dbprint(sprintf("done! (%s ms)%n",rtimer-t0,0));
    20992016    if(N<=0)
    2100       {print("checking for correctness (exact test) ");}
     2017      {dbprint("checking for correctness (exact test) ");}
    21012018    else
    2102       {printf("checking for correctness (%s random evaluations per entry) ",N);}
     2019      {dbprint(sprintf("checking for correctness (%s random evaluations per entry) ",N));}
    21032020    t0=rtimer;
    21042021  }
    21052022
    2106   fprintf(":w "+filename+"_checkpfd_logfile.txt","Input file (matrix of rational functions): %s"
    2107           +"%nOutput file (decompositions): %s%nlist of all denominator factors: %s"
    2108           +"%n%nResults of checkpfdMat:",input,output,qfile,0);
    2109   link logfile = ":a "+filename+"_checkpfd_logfile.txt";
     2023  fprintf(":w "+filename+"_checkpfdMat_logfile.txt","Input file (matrix of rational functions):"
     2024          +" %s%nOutput file (decompositions): %s%nlist of all denominator factors:"
     2025          +" %s%n%nResults of checkpfdMat:",input,output,qfile,0);
     2026  link logfile = ":a "+filename+"_checkpfdMat_logfile.txt";
    21102027
    21112028  int n=size(frac);
     
    21352052  if(parallelize)
    21362053  {
    2137     printf("done! (%s ms)%n%n",rtimer-t0,0);
     2054    dbprint(sprintf("done! (%s ms)%n",rtimer-t0,0));
    21382055    if(N<=0)
    2139       {print("checking for correctness (exact test) ");}
     2056      {dbprint("checking for correctness (exact test) ");}
    21402057    else
    2141       {printf("checking for correctness (%s random evaluations per entry) ",N);}
     2058      {dbprint(sprintf("checking for correctness (%s random evaluations per entry) ",N));}
    21422059    t0=rtimer;
    21432060    list results = parallelWaitAll("testEntry",arguments);
    21442061  }
    2145   else {printf("done! (%s ms)",rtimer-t0);}
    2146   printf("%s out of %s = %sx%s decompositions are correct! (%s ms)%n",
    2147           sum(results),n*m,n,m,rtimer-t0,0);
     2062  else {dbprint(sprintf("done! (%s ms)",rtimer-t0));}
     2063  dbprint(sprintf("%s out of %s = %sx%s decompositions are correct! (%s ms)%n",
     2064          sum(results),n*m,n,m,rtimer-t0,0));
    21482065  fprintf(logfile,"%s out of %s = %sx%s decompositions are correct! (%s ms)%n",
    21492066          sum(results),n*m,n,m,rtimer-t0,0);
     
    21532070{
    21542071  system("--ticks-per-sec",1000);
    2155   print("  reading matrix of decompositions from file "+filename__);
     2072  dbprint("  reading matrix of decompositions from file "+filename__);
    21562073  int t__=rtimer;
    21572074  int tt__;
    21582075  string data__ = read(":r "+filename__);
    2159   printf("  done! (%s ms)", rtimer-t__);
    2160 
    2161   print("  processing input "); t__ = rtimer;
     2076  dbprint(sprintf("  done! (%s ms)", rtimer-t__));
     2077
     2078  dbprint("  processing input "); t__ = rtimer;
    21622079  list mat__;
    21632080  list nonlin__=list();
     
    21902107      mat__[i__][j__] = list();
    21912108
    2192 
    21932109      factors__ = intvec(0:0);
    21942110      exponents__ = intvec(0:0);
     
    21962112      if(l__>0)
    21972113      {
    2198         ss__ = s__[1,l__-1];  //ss__ contains the nonlinear factors
     2114        ss__ = s__[1,l__-1];    //ss__ contains the nonlinear factors
    21992115        s__ = s__[l__+3,size(s__)-l__-2];
    22002116        for(p1__=1;s__[p1__]!="(";p1__++) {}
     
    22122128          p2__ = find(ss__,"^",p1__);
    22132129          tmp__ = find(ss__,"*",p1__);
    2214           if((p2__>tmp__ && tmp__>0) || (p2__==0)) //exponent is 1
     2130          if((p2__>tmp__ && tmp__>0) || (p2__==0))    //exponent is 1
    22152131          {
    22162132            execute("factors__[l__]="+ss__[p1__,tmp__-p1__]);
     
    22382154        if((s__[k__]=="+" && depth__==0) || k__==max__)
    22392155        {
    2240           if(tmp2__==0) // no denominator
     2156          if(tmp2__==0)    // no denominator
    22412157          {
    22422158            execute("numerator__="+s__[tmp__,k__-tmp__]);
     
    22492165          p1__ = find(ss__,"(");
    22502166          p2__ = find(ss__,")");
    2251           ss__ = ss__[p1__+1,p2__-p1__-1]; // now ss__ is only the denominator
     2167          ss__ = ss__[p1__+1,p2__-p1__-1];   // now ss__ is only the denominator
    22522168
    22532169          ss__ = ss__+"*";
     
    22622178            p2__ = find(ss__,"^",p1__);
    22632179            tmp__ = find(ss__,"*",p1__);
    2264             if((p2__>tmp__ && tmp__>0) || (p2__==0)) //exponent is 1
     2180            if((p2__>tmp__ && tmp__>0) || (p2__==0))    //exponent is 1
    22652181            {
    22662182              execute("factors__[l__]="+ss__[p1__,tmp__-p1__]);
     
    22792195      }
    22802196    }
    2281     printf("    row %s done! (%s ms)",i__,rtimer-tt__);
    2282   }
    2283   printf("  done! (%s ms)", rtimer-t__);
     2197    dbprint(sprintf("    row %s done! (%s ms)",i__,rtimer-tt__));
     2198  }
     2199  dbprint(sprintf("  done! (%s ms)", rtimer-t__));
    22842200
    22852201  return(mat__,nonlin__);
     
    22912207  ideal q__;
    22922208  int pos1__,pos2__=1,1;
    2293 
    22942209  while(1)
    22952210  {
  • Tst/Short/pfd.tst

    r6924d4 rcc07066  
    66example pfd;
    77
    8 /* pfd(poly,poly) */
     8// pfd(poly,poly) //////////////////////////////////////////////////////////////
    99ring r1 = 0, x(1..5), dp;
    1010poly f = -2*x(1)*x(3)+3*x(3)*x(4)+x(2)*x(5)+x(3)*x(5)-x(4)*x(5);
     
    3232kill dec;
    3333
     34
     35// pfd(poly,list) (denominator = ideal of factors & intvec of exponents) ///////
    3436ring r3 = 5, (x,y,z), dp;
    35 poly f = 2*x^3-x^2*y+x*y^2+y^3-2*x^2*z-2*x*y*z-y^2*z+2*x*z^2+y*z^2-2*z^3;
    36 poly g = (x^2+3*y)^2*(x-2*y^2)*(x+y+1)^3*(2*x+3*y+4)*(x^2-x*y+y^2);
     37poly f = x+y+z+1;
     38list g = list(ideal((x^2+y^2+z^2),(x+y^2),(y+z^2),(z+x^2)), intvec(2,1,1,1));
     39list dec = pfd(f,g);
     40displaypfd(dec);
     41checkpfd(list(f,g), dec);
     42checkpfd(list(f,g), dec, 10);
     43kill dec;
     44
     45// different ordering, same polynomials:
     46ring r4 = 5, (x,y,z), lp;
     47poly f = fetch(r3,f);
     48list g = fetch(r3,g);
    3749list dec = pfd(f,g);
    3850displaypfd(dec);
     
    4153kill dec;
    4254
    43 ring r4 = 5, (x,y,z), lp;
    44 poly f = fetch(r3,f);
    45 poly g = fetch(r3,g);
    46 list dec = pfd(f,g);
    47 displaypfd(dec);
    48 checkpfd(list(f,g), dec);
    49 checkpfd(list(f,g), dec,10);
    50 kill dec;
    5155
    52 /* pfd(list) */
     56// pfd(list) ///////////////////////////////////////////////////////////////////
    5357ring s1 = 0, (x,y,z), dp;
    5458poly f1 = x*y+y*z+z*x-x-y-z+1;
     
    6872kill dec;
    6973
    70 /* pfd(matrix) */
     74// pfd(matrix) /////////////////////////////////////////////////////////////////
    7175ring s2 = 3, (x,y,z), dp;
    7276poly f11 = (x+y+z+1)^3;
     
    8589displaypfd(dec[2][1]);
    8690displaypfd(dec[2][2]);
    87 checkpfd(list(f11,g11), dec[1][1]);
    88 checkpfd(list(f12,g12), dec[1][2]);
    89 checkpfd(list(f21,g21), dec[2][1]);
    90 checkpfd(list(f22,g22), dec[2][2]);
     91checkpfd(list(f1,g1), dec[1][1]);
     92checkpfd(list(f2,g2), dec[1][2]);
     93checkpfd(list(f3,g3), dec[2][1]);
     94checkpfd(list(f4,g4), dec[2][2]);
    9195kill dec;
    9296
Note: See TracChangeset for help on using the changeset viewer.