////////////////////////////////////////////////////////////////////
version="version pfd.lib 4.1.3.2 Aug_2020 ";
category="??";
info="
LIBRARY: pfd.lib Multivariate Partial Fraction Decomposition

AUTHOR: Marcel Wittmann, e-mail: mwittman@mathematik.uni-kl.de

OVERVIEW:
This Library implements an algorithm based on the work of E. K. Leinartas to
write rational functions in multiple variables as a sum of functions with
\"smaller\" numerators and denominators.
This can be used to shorten the IBP reduction coefficients of multi-loop Feynman
integrals. For this application,
we also provide a procedure that applies the algorithm to all entries of a
matrix of rational functions given as one (possibly very big) txt-file.
If you use the library pfd.lib, please cite the corresponding paper
[J. Boehm, M. Wittmann, Z. Wu, Y. Xu, Y. Zhang: 'IBP reduction coefficients 
made simple'] (preprint 2020).

KEYWORDS: partial fraction; decomposition; Leinartas

PROCEDURES:
  pfd();                  calculate a partial fraction decomposition
                          of a rational function
  checkpfd();             test if a decomposition is equal to a rational
                          function given by numerator/denominator polynomials
  evaluatepfd();          substitute values in a partial fraction
                          decomposition gotten from @code{pfd}
  displaypfd();           print a decomposition gotten as output of @code{pfd}
  displaypfd_long();      like @code{display}, but denominators are written out
  getStringpfd();         turn a decomposition gotten from @code{pfd} into one
                          string
  getStringpfd_indexed(); like @code{getStringpfd}, but writes the denominator
                          factors just as @code{q1}, @code{q2}, ...
  readInputTXT();         read a matrix of rational functions from a txt-file
  pfdMat();               apply @code{pfd} to a matrix of rational functions
                          in parallel (using @ref{parallel_lib}) and save result
                          as easy-to-read txt-files.
  checkpfdMat();          test output files of @code{pfdMat} for correctness
";
/////////////////////////////////////////////////////////////////////////////

LIB "random.lib";
LIB "ring.lib";
LIB "parallel.lib";
LIB "elim.lib";
LIB "polylib.lib";

/////////////////////////////////////////////////////////////////////////////

static proc mod_init()
{
  printlevel = 2;
  // export some static functions so they can be run in parallel using parallelWaitAll:
  if (!defined(Tasks)) {LIB "tasks.lib";}
  exportto(Tasks,pfdWrap);
  importfrom(Tasks,pfdWrap);
  exportto(Pfd,pfdWrap);
  exportto(Tasks,testEntry);
  importfrom(Tasks,testEntry);
  exportto(Pfd,testEntry);
}

proc pfd(list #)
"USAGE:   pfd(f,g[,debug]);   f,g poly, debug int
          pfd(f,g[,debug]);   f poly, g list, debug int
          pfd(arguments[, parallelize]);     arguments list, parallelize int
RETURN:   a partial fraction decomposition of f/g as a list @code{l} where
          @code{l[1]} is an ideal generated by irreducible polynomials and
          @code{l[2]} is a list of fractions.
          Each fraction is represented by a list of
       @* 1) the numerator polynomial
       @* 2) an intvec of indices @code{i} for which @code{l[1][i]} occurs
             as a factor in the denominator
       @* 3) an intvec containing the exponents of those irreducible factors.
       @* Setting @code{debug} to a positive integer measures runtimes and
          creates a log file (default: @code{debug=0}).
       @* The denominator g can also be given in factorized form as a list of
          an ideal of irreducible non constant polynomials and an intvec of
          exponents. This can save time since the first step in the algorithm is
          to factorize g. (A list of the zero-ideal and an empty intvec
          represents a denominator of 1.)
       @* If instead of f and g, the input is a single list (or even a list of
          lists) containing elements of the form @code{list(f,g[,debug])}
          (@code{f,g,debug} as above), the algorithm is applied to all entries
          in parallel (using @ref{parallel_lib}), if @code{parallelize=1}
          (default) and in sequence if @code{parallelize=0}. A list (or list of
          lists) of the results is returned.
NOTE:     The result depends on the monomial ordering. For \"small\" results
          use @code{dp}.
SEE ALSO: checkpfd, evaluatepfd, displaypfd, displaypfd_long, pfdMat
EXAMPLE:  example pfd; shows an example"
{
  short = 0;
  int i;
  if(typeof(#[1])=="list")
  {
    if(size(#)>1)
    {
      if(typeof(#[2])=="list")
        {list arguments = #;    int parallelize = 1;}
      else{if(typeof(#[2])=="int")
        {list arguments = #[1]; int parallelize = #[2];}
      else {ERROR("wrong type for second argument, expected int");}}
    }
    if(parallelize)
    {
      for(i=1; i<=size(arguments); i++)
      {
        if(typeof(arguments[i][1])=="list") //input is list of lists
          {arguments[i] = list(arguments[i],1);}
      }
      return(parallelWaitAll("pfd",arguments));
    }
    else
    {
      list results;
      for(i=1; i<=size(arguments); i++)
      {
        if(typeof(arguments[i][1])=="list") //input is list of lists
          {results[i] = pfd(arguments[i],0);}
        else
        {
          if(size(arguments[i])==2)
            {results[i] = pfd(arguments[i][1],arguments[i][2]);}
          else{if(size(arguments[i])==3)
            {results[i] = pfd(arguments[i][1],arguments[i][2],arguments[i][3]);}
          else {ERROR("wrong number of arguments, expected 2 or 3");}}
        }
      }
      return(results);
    }
  }
  poly f = #[1];
  if(typeof(#[2])=="list")
  {
    list g=#[2];
  }
  else
  {
    poly g=#[2];
  }

  int debug=0; link l=":w ";
  if(size(#)>3) {ERROR("wrong number of arguments, expected 2 or 3");}
  if(size(#)==3)
  {
    debug=#[3];
    l=":a "+string(debug)+"_log_"+datetime()+".txt";
    system("--ticks-per-sec",1000);
  }
  if(debug)
  {
    fprintf(l,"debug: %s", debug);
    fprintf(l,"size(string(f)) = %s, size(string(g)) = %s %n",
               size(string(f)), size(string(g)), 0);
    int counter,tt,ttt;
  }

  if(f==0)
  {
    list dec = list(ideal(),list(list(poly(0),intvec(0:0),intvec(0:0))));
    if(debug)
      {fprintf(l,"%ntotal: 0 ms (numerator was 0)",0); close(l);}
    if(voice<=printlevel) {displaypfd(dec);}
    return(dec);
  }
  if(typeof(g)=="poly")
  {
    if(deg(g)==0)
    {
      list dec = list(ideal(),list(list(f/g,intvec(0:0),intvec(0:0))));
      if(debug)
        {fprintf(l,"%ntotal: 0 ms (denominator was constant)",0); close(l);}
      if(voice<=printlevel) {displaypfd(dec);}
      return(dec);
    }

    // (1) factorization of the denominator ////////////////////////////////////
    if(debug) {int t1 = rtimer; write(l,"factorizing ");}
    list factor = factorize(g);
    number lcoeff;
    for(i=2; i<=size(factor[1]); i++)
    {
      lcoeff = leadcoef(factor[1][i]);
      factor[1][i] = factor[1][i]/lcoeff;
      factor[1][1] = factor[1][1]*(lcoeff^factor[2][i]); // polynomial is monic (thus unique)
      lcoeff = content(factor[1][i]);
      factor[1][i] = factor[1][i]/lcoeff;
      factor[1][1] = factor[1][1]*(lcoeff^factor[2][i]); // polynomial has nice coefficients
    }
    ideal q = factor[1];
    f = f/q[1];
    q=delete(q,1);
    intvec e = factor[2]; e=delete(e,1);
    int m = size(q);
    if(debug) {t1 = rtimer-t1; fprintf(l,"done! (%s ms)", t1);}
  }
  else{if(typeof(g)=="list")
  {
    if(size(g[1])==0)
    {
      list dec = list(ideal(),list(list(f,intvec(0:0),intvec(0:0))));
      if(debug)
        {fprintf(l,"%ntotal: 0 ms (denominator was constant)",0); close(l);}
      if(voice<=printlevel) {displaypfd(dec);}
      return(dec);
    }

    // denominator is already factorized
    for(i=1;i<=size(g[1]);i++)
    {
      if(size(factorize(g[1][i])[1])>2)
        {ERROR("factors should be irreducible");}
    }

    ideal q = g[1];
    intvec e = g[2];
    int m = size(q);
    if(debug) {int t1=0;}
  }
  else
  {ERROR("wrong type for second argument, expected poly or list(ideal,intvec)");}}

  // (2) Nullstellensatz decomposition /////////////////////////////////////////
  if(debug)
    {int t2 = rtimer; write(l,"Nullstellensatz decomposition ");}
  list terms = list(list(poly(f),1..m,e));
  list dec,newterms,result;
  int imax;
  while(size(terms)>0)
  {
    if(debug) {ttt = rtimer; counter++;}
    imax = size(terms);
    newterms = list();
    for(i=1; i<=imax; i++)
    {
      result = NSSdecompStep(terms[i],q);
      if(size(result)>1) {newterms = mergepfd(newterms, result);}
      else                    {dec = mergepfd(dec,      result);}
    }
    terms = newterms;
    if(debug)
      {fprintf(l,"  %s: %s ms, %s terms, %s unfinished",
       counter, rtimer-ttt, size(terms)+size(dec), size(terms));}
  }
  if(debug)
    {t2 = rtimer-t2; fprintf(l,"done! (%s ms, %s terms)", t2, size(dec));}

  // (3) short numerator decomposition /////////////////////////////////////////
  if(debug)
    {int t3 = rtimer; write(l,"short numerator decompositions "); counter=0;}
  terms = dec;
  dec = list();
  while(size(terms)>0)
  {
    if(debug) {tt = rtimer; counter++;}
    imax = size(terms);
    newterms = list();
    for(i=1; i<=imax; i++)
    {
      result = shortNumeratorDecompStep(terms[i],q,debug,l);

      if(debug) {ttt = rtimer;}
      newterms = mergepfd(newterms, result[1]);
      dec = mergepfd(dec, result[2]);
      if(debug)
        {fprintf(l,"      merging: %s ms, %s new terms",
         rtimer-ttt, size(result[1]));}
    }
    terms = newterms;
    if(debug)
      {fprintf(l,"  %s: %s ms, %s terms, %s unfinished",
       counter, rtimer-tt, size(terms)+size(dec), size(terms));}
  }
  if(debug)
    {t3 = rtimer-t3; fprintf(l,"done! (%s ms, %s terms)", t3, size(dec));}

  // (4) algebraic dependence decomposition ////////////////////////////////////
  if(debug)
    {int t4 = rtimer; write(l,"algebraic dependence decomposition "); counter=0;}
  terms = dec;
  dec = list();
  while(size(terms)>0)
  {
    if(debug) {tt = rtimer; counter++;}
    imax = size(terms);
    newterms = list();
    for(i=1; i<=imax; i++)
    {
      result = algDependDecompStep(terms[i],q,debug,l);

      if(debug) {ttt = rtimer;}
      if(size(result)>1) {newterms = mergepfd(newterms, result);}
      else                    {dec = mergepfd(dec,      result);}
      if(debug)
        {fprintf(l,"      merging: %s ms, %s new terms", rtimer-ttt, size(result));}
    }
    terms = newterms;
    if(debug)
      {fprintf(l,"  %s: %s ms, %s terms, %s unfinished",
       counter, rtimer-tt, size(terms)+size(dec),size(terms));}
  }
  if(debug)
    {t4 = rtimer-t4; fprintf(l,"done! (%s ms, %s terms)", t4, size(dec));}

  // (5) numerator decomposition ///////////////////////////////////////////////
  if(debug)
    {int t5 = rtimer; write(l,"numerator decompositions "); counter=0;}
  terms = dec;
  dec = list();
  while(size(terms)>0)
  {
    if(debug) {tt = rtimer; counter++;}
    imax = size(terms);
    newterms = list();
    for(i=1; i<=imax; i++)
    {
      result = numeratorDecompStep(terms[i],q,debug,l);

      if(debug) {ttt = rtimer;}
      newterms = mergepfd(newterms, result[1]);
      dec = mergepfd(dec, result[2]);
      if(debug)
        {fprintf(l,"      merging: %s ms, %s new terms", rtimer-ttt, size(result[1]));}
    }
    terms = newterms;
    if(debug)
      {fprintf(l,"  %s: %s ms, %s terms, %s unfinished",
       counter, rtimer-tt, size(terms)+size(dec), size(terms));}
  }
  if(debug)
    {t5 = rtimer-t5; fprintf(l,"done! (%s ms, %s terms)", t5, size(dec));}

  dec = list(q,dec);
  if(debug) {fprintf(l,"%ntotal: %s ms", t1+t2+t3+t4+t5); close(l);}
  if(voice<=printlevel) {displaypfd(dec);}
  return(dec);
}
example
{
  "EXAMPLE:";
  echo=voice;
  ring R = 0,(x,y),dp;
  poly f = x^3+3*x^2*y+2*y^2-x^2+4*x*y;
  poly g = x^2*y*(x-1)*(x-y)^2;

  list dec = pfd(f,g);

  displaypfd_long(dec);   // display result
  checkpfd(list(f,g),dec);   // check for equality to f/g

  // calculate decompositions of a 2x2 matrix of rational functions at once:
  list arguments = list(list(f, g),          list(1, f)     ),
                        list(list(x*y, y+1), list(1, x^2-y^2));

  dec = pfd(arguments);

  // the result has the same shape as the
  // input (2x2 matrix as list of lists):
  displaypfd_long(dec[1][1]);
  displaypfd_long(dec[1][2]);
  displaypfd_long(dec[2][1]);
  displaypfd_long(dec[2][2]);

  // a more complicated example
  ring S = 0,(s12,s15,s23,s34,s45),dp;
  poly f = 7*s12^4*s15^2 + 11*s12^3*s15^3 + 4*s12^2*s15^4 - 10*s12^4*s15*s23
    - 14*s12^3*s15^2*s23 - 4*s12^2*s15^3*s23 + 3*s12^4*s23^2 + 3*s12^3*s15*s23^2
    + 13*s12^4*s15*s34 + 12*s12^3*s15^2*s34 + 2*s12^2*s15^3*s34
    - 5*s12^4*s23*s34 + 33*s12^3*s15*s23*s34 + 49*s12^2*s15^2*s23*s34
    + 17*s12*s15^3*s23*s34 - 17*s12^3*s23^2*s34 - 19*s12^2*s15*s23^2*s34
    - 5*s12*s15^2*s23^2*s34 - 24*s12^3*s15*s34^2 - 15*s12^2*s15^2*s34^2
    + 2*s12*s15^3*s34^2 + 15*s12^3*s23*s34^2 - 34*s12^2*s15*s23*s34^2
    - 31*s12*s15^2*s23*s34^2 + 2*s15^3*s23*s34^2 + 33*s12^2*s23^2*s34^2
    + 29*s12*s15*s23^2*s34^2 + 5*s15^2*s23^2*s34^2 + 9*s12^2*s15*s34^3
    - 4*s12*s15^2*s34^3 - 15*s12^2*s23*s34^3 + 9*s12*s15*s23*s34^3
    - 4*s15^2*s23*s34^3 - 27*s12*s23^2*s34^3 - 13*s15*s23^2*s34^3
    + 2*s12*s15*s34^4 + 5*s12*s23*s34^4 + 2*s15*s23*s34^4 + 8*s23^2*s34^4
    - 6*s12^3*s15^2*s45 - 9*s12^2*s15^3*s45 - 2*s12*s15^4*s45
    + 30*s12^3*s15*s23*s45 + 56*s12^2*s15^2*s23*s45 + 24*s12*s15^3*s23*s45
    - 12*s12^3*s23^2*s45 - 23*s12^2*s15*s23^2*s45 - 10*s12*s15^2*s23^2*s45
    - 30*s12^3*s15*s34*s45 - 32*s12^2*s15^2*s34*s45 - 6*s12*s15^3*s34*s45
    + 7*s12^3*s23*s34*s45 - 86*s12^2*s15*s23*s34*s45 - 104*s12*s15^2*s23*s34*s45
    - 15*s15^3*s23*s34*s45 + 41*s12^2*s23^2*s34*s45 + 51*s12*s15*s23^2*s34*s45
    + 10*s15^2*s23^2*s34*s45 - 5*s12^3*s34^2*s45 + 33*s12^2*s15*s34^2*s45
    + 14*s12*s15^2*s34^2*s45 - 2*s15^3*s34^2*s45 - 21*s12^2*s23*s34^2*s45
    + 62*s12*s15*s23*s34^2*s45 + 28*s15^2*s23*s34^2*s45 - 46*s12*s23^2*s34^2*s45
    - 28*s15*s23^2*s34^2*s45 + 10*s12^2*s34^3*s45 - s12*s15*s34^3*s45
    + 4*s15^2*s34^3*s45 + 21*s12*s23*s34^3*s45 - 6*s15*s23*s34^3*s45
    + 17*s23^2*s34^3*s45 - 5*s12*s34^4*s45 - 2*s15*s34^4*s45 - 7*s23*s34^4*s45
    - 6*s12^2*s15^2*s45^2 - 5*s12*s15^3*s45^2 - 2*s15^4*s45^2
    - 28*s12^2*s15*s23*s45^2 - 42*s12*s15^2*s23*s45^2 - 10*s15^3*s23*s45^2
    + 9*s12^2*s23^2*s45^2 + 10*s12*s15*s23^2*s45^2 + 24*s12^2*s15*s34*s45^2
    + 36*s12*s15^2*s34*s45^2 + 10*s15^3*s34*s45^2 - 11*s12^2*s23*s34*s45^2
    + 31*s12*s15*s23*s34*s45^2 + 25*s15^2*s23*s34*s45^2
    - 18*s12*s23^2*s34*s45^2 - 10*s15*s23^2*s34*s45^2 + 4*s12^2*s34^2*s45^2
    - 29*s12*s15*s34^2*s45^2 - 17*s15^2*s34^2*s45^2 + 27*s12*s23*s34^2*s45^2
    + 2*s15*s23*s34^2*s45^2 + 9*s23^2*s34^2*s45^2 - 3*s12*s34^3*s45^2
    + 10*s15*s34^3*s45^2 - 16*s23*s34^3*s45^2 - s34^4*s45^2 + 6*s12*s15^2*s45^3
    + 3*s15^3*s45^3 + 8*s12*s15*s23*s45^3 + 10*s15^2*s23*s45^3
    - 8*s12*s15*s34*s45^3 - 10*s15^2*s34*s45^3 + 9*s12*s23*s34*s45^3
    + s12*s34^2*s45^3 + 8*s15*s34^2*s45^3 - 9*s23*s34^2*s45^3 - s34^3*s45^3
    - s15^2*s45^4 + s15*s34*s45^4;
  poly g = 4*s12*s15*(s12 + s15 - s34)*(s15 - s23 - s34)*(s12 + s23 - s45)
                                 *(s12 - s34 - s45)*(s12 + s15 - s34 - s45)*s45;

  list dec = pfd(f,g);

  displaypfd(dec);
  checkpfd(list(f,g),dec);

  // size comparison:
  size(string(f)) + size(string(g));
  size(getStringpfd(dec));
}

static proc NSSdecompStep(list l, ideal q)
{
  poly f=l[1];
  intvec indices=l[2];
  intvec e=l[3];
  int m = size(indices);

  if(m==0) // denominator is 1
    {return(list(l));} // do nothing, return input

  ideal qe = q[indices];
  for(int i=1; i<=m; i++)
    {qe[i] = qe[i]^e[i];}
  matrix T;
  ideal qe_std = liftstd(qe,T);

  if(deg(qe_std) == 0)
  {
    T = T/qe_std[1];
      // now 1 = T[1,1]*qe[1] + ... + T[m,1]*qe[m] is a Nullstellensatz certificate
    list dec;
    poly h;
    for(i=1; i<=m; i++)
    {
      h = T[i,1];
      if(h != 0)
        {dec[size(dec)+1] = list(f*h,delete(indices,i),delete(e,i));}
    }
    return(dec);
  }
  else
  {
    return(list(l)); // do nothing, return input
  }
}

static proc shortNumeratorDecompStep(list l, ideal q, list #)
{
  int debug=0;
  if(size(#)>0) {debug=#[1]; link ll=#[2];}
  if(debug) {system("--ticks-per-sec",1000); int tt=rtimer;}

  poly f=l[1];
  intvec indices=l[2];
  intvec e=l[3];
  int m = size(indices);

  if(m==0) // denominator is 1
  {
    if(debug)
      {fprintf(ll,"      shortNumeratorDecompStep: %s ms (m=%s, e=%s) "
                    +"--> constant denominator", rtimer-tt, m, e);}
    return(list(list(),list(l))); // do nothing, return input
  }

  ideal q_denom = q[indices]; // factors occurring in the denominator
  matrix T;
  ideal q_std = liftstd(q_denom,T);
  list divrem = division(f,q_std);
  poly r = divrem[2][1]/divrem[3][1,1];

  if(r!=0)
  {
    if(debug)
      {fprintf(ll,"      shortNumeratorDecompStep: %s ms (m=%s, e=%s) "
                  +"--> remainder is nonzero", rtimer-tt, m, e);}
    return(list(list(),list(l))); // if there is a rest, decomposing further would
  }                               // not help in the next step (alg. depend. decomposition)

  matrix a = divrem[1]/divrem[3][1,1]; // now f = r + a[1,1]*q_std[1] + ... +a[m,1]*q_std[m]
  a = T*a;   // lift coefficients   ==>   now f = r + a[1,1]*q[1]     + ... +a[m,1]*q[m]

  // reduce w.r.t. groebner basis of syz(q) to make the numerators "smaller":
  vector v;
  for(int i=1; i<=m; i++) {v = v + gen(i)*a[i,1];}
  v = reduce(v, std(syz(q_denom)));

  list fraction,dec;
  for(i=1; i<=m; i++)
  {
    if(v[i] == 0)
    {
      i++;
      continue;
    }
    fraction[1] = v[i];
    if(e[i]==1)
      {
        fraction[2] = delete(indices,i);
        fraction[3] = delete(e,i);
      }
    else
    {
      fraction[2] = indices;
      fraction[3] = e;
      fraction[3][i] = fraction[3][i] - 1;
    }
    dec[size(dec)+1] = fraction;
  }

  if(debug)
    {fprintf(ll,"      shortNumeratorDecompStep: %s ms (m=%s, e=%s, deg(v)=%s, size(v)=%s)",
     rtimer-tt, m, e, deg(v), size(v));}
  return(list(dec,list()));
}

static proc algDependDecompStep(list l, ideal q, list #)
{
  int debug=0;
  if(size(#)>0) {debug=#[1]; link ll=#[2];}
  if(debug) {system("--ticks-per-sec",1000); int tt=rtimer;}
  def br = basering;
  int d = nvars(br);
  intvec indices=l[2];
  int m = size(indices);
  intvec e=l[3];
  int i;

  if(m==0) // denominator is 1
  {
    if(debug)
      {fprintf(ll,"      algDependDecompStep: %s ms (m=%s, e=%s)", rtimer-tt, m, e);}
    return(list(l)); // do nothing, return input
  }

  if(m<=d)
  {
    if(size(syz(module(transpose(jacob(ideal(q[indices]))))))==0) // jacobian criterion
    {
      if(debug)
        {fprintf(ll,"      algDependDecompStep: %s ms (m=%s, e=%s) "
                    +"--> alg. indep.", rtimer-tt, m, e);}
      return(list(l)); // do nothing, return input
    }
  }

  def R = changeord(list(list("dp",m+d)),extendring(m, "y(", "dp", 1, changevar("x()",br)));
  setring(R);

  list l = fetch(br,l);
  ideal q = fetch(br,q);
  poly f=l[1];

  ideal I;
  for(i=1; i<=m; i++)
    {I[i] = y(i)-q[indices[i]];}

  ideal annihilatingPolys = eliminate(I,intvec(1..d));

  poly g = annihilatingPolys[1];

  poly tail = g[size(g)];    // term of lowest dp-order (thus lowest degree)
  number tcoeff = leadcoef(tail);
  intvec texpon = leadexp(tail);
  texpon = texpon[(d+1)..(d+m)];
  g = g-tail;
  poly term;
  number coeff;
  intvec expon;
  list fraction,dec;
  int pow;
  int jmax = size(g);
  for(int j=1; j<=jmax; j++)
  {
    term = g[j];
    coeff = leadcoef(term);
    expon = leadexp(term);
    expon = expon[(d+1)..(d+m)];
    fraction[1] = -f*coeff/tcoeff;
    fraction[2] = intvec(0:0);
    fraction[3] = intvec(0:0);
    for(i=1; i<=m; i++)
    {
      pow = expon[i]-texpon[i]-e[i];
      if(pow>=0)
      {
        fraction[1] = fraction[1]*q[indices[i]]^pow;
      }
      else
      {
        fraction[2][size(fraction[2])+1] = indices[i];
        fraction[3][size(fraction[3])+1] = -pow;
      }
    }
    dec[size(dec)+1] = fraction;
  }
  setring(br);
  list dec = fetch(R,dec);
  if(debug)
    {fprintf(ll,"      algDependDecompStep: %s ms (m=%s, e=%s)", rtimer-tt, m, e);}
  return(dec);
}

static proc numeratorDecompStep(list l, ideal q, list #)
{
  int debug=0;
  if(size(#)>0) {debug=#[1]; link ll=#[2];}
  if(debug) {system("--ticks-per-sec",1000); int tt=rtimer;}

  poly f=l[1];
  intvec indices=l[2];
  intvec e=l[3];
  int m = size(indices);

  if(m==0) // denominator is 1
  {
    if(debug)
      {fprintf(ll,"      numeratorDecompStep: %s ms (m=%s, e=%s) "
                  +"--> constant denominator", rtimer-tt, m, e);}
    return(list(list(),list(l))); // do nothing, return input
  }

  ideal q_denom = q[indices]; // factors in the denominator
  matrix T;
  ideal q_std = liftstd(q_denom,T);
  list divrem = division(f,q_std);
  matrix a = divrem[1]/divrem[3][1,1];
  poly r = divrem[2][1]/divrem[3][1,1]; // now f = r + a[1,1]*q_std[1] + ... +a[m,1]*q_std[m]
  a = T*a;   // lift coefficients    ==>   now f = r + a[1,1]*q[1]     + ... +a[m,1]*q[m]

  // reduce w.r.t. groebner basis of syz(q) to make the numerators "smaller"
  vector v;
  for(int i=1; i<=m; i++) {v = v + gen(i)*a[i,1];}
  v = reduce(v, std(syz(q_denom)));

  list fraction,dec,rest;
  if(r!=0)
    {rest[1] = list(r,indices,e);}
  for(i=1; i<=m; i++)
  {
    if(v[i] == 0)
    {
      i++;
      continue;
    }
    fraction[1] = v[i];
    if(e[i]==1)
      {
        fraction[2] = delete(indices,i);
        fraction[3] = delete(e,i);
      }
    else
    {
      fraction[2] = indices;
      fraction[3] = e;
      fraction[3][i] = fraction[3][i] - 1;
    }
    dec[size(dec)+1] = fraction;
  }

  if(debug)
    {fprintf(ll,"      numeratorDecompStep: %s ms (m=%s, e=%s, deg(v)=%s, size(v)=%s)",
     rtimer-tt, m, e, deg(v), size(v));}
  return(list(dec,rest));
}

static proc mergepfd(list dec1, list dec2)
{
  // Note: assumes dec1 is already sorted w.r.t. dp in the denominator exponents
  int n1=size(dec1);
  int n2=size(dec2);
  if(n2==0) {return(dec1);}
  int i;
  int a,b,m;
  list entry;
  for(i=1; i<=n2; i++)
  {
    entry = dec2[i];
    if(n1==0)
    {
      dec1=list(entry); n1++;
      i++; continue;
    }

    a=1;
    b=n1;
    m = (a+b) div 2;
    while(b>a)
    {
      if(is_dp_smaller(dec1[m][2], dec1[m][3], entry[2], entry[3]))
      {
        a=m+1;
      }
      else
      {
        b=m;
      }
      m = (a+b) div 2;
    }
    if(is_dp_smaller(dec1[a][2], dec1[a][3], entry[2], entry[3]))
      {dec1=insert(dec1,entry,a); n1++;}
    else{if(entry[2]==dec1[a][2] && entry[3]==dec1[a][3])
    {
      dec1[a][1] = dec1[a][1] + entry[1];   //same denominator: add numerators
      if(dec1[a][1]==0)
        {dec1 = delete(dec1,a); n1--;}
    }
    else
      {dec1=insert(dec1,entry,a-1); n1++;}}
  }
  return(dec1);
}

static proc is_dp_smaller(intvec indices1, intvec e1, intvec indices2, intvec e2)
{
  if(size(e2)==0) {return(0);}
  if(size(e1)==0) {return(1);}
  int s1,s2 = sum(e1),sum(e2);
  if(s1<s2) {return(1);}
  if(s1>s2) {return(0);}
  int n1,n2 = size(indices1),size(indices2);
  int imax = min(n1,n2);
  for(int i=0; i<imax; i++)
  {
    if(indices1[n1-i]>indices2[n2-i]) {return(1);}
    if(indices1[n1-i]<indices2[n2-i]) {return(0);}
    if(e1[n1-i]>e2[n2-i]) {return(1);}
    if(e1[n1-i]<e2[n2-i]) {return(0);}
  }
  return(0);
}

proc checkpfd(list fraction, list dec, list #)
"USAGE:   checkpfd(list(f,g),dec[,N,C]);   f,g poly, dec list, N,C int
RETURN:   0 or 1
PURPOSE:  test for (mathematical) equality of f/g and a partial fraction
          decomposition dec. The list dec has to have the same structure as the
          output of @ref{pfd}.
       @* The denominator g can also be given in factorized form as a list of
          an ideal of irreducible non constant polynomials and an intvec of
          exponents. This can save time since the first step in the algorithm is
          to factorize g. (a list of the zero-ideal and an empty intvec
          represents a denominator of 1.)
       @* By default the test is done (exactly) by bringing all terms of the
          decomposition on the same denominator and comparing to f/g.
       @* If additional parameters N [, C] are given and if @code{N>0}, a
          probabilistic method is chosen: evaluation at N random points with
          coordinates between -C and C. This may be faster for big polynomials.
SEE ALSO: pfd
EXAMPLE:  example checkpfd; shows an example"
{
  poly f = fraction[1];
  def g = fraction[2];
  if(size(#)>0)
  {
    if(#[1]>0)
    {
      int N = #[1]; // number of random tests
      int max_val=16;
      if(size(#)>1) {max_val = #[2];}
      ideal values;
      ideal vars = maxideal(1);
      int d=nvars(basering);
      number val1,val2;
      int div_by_0;
      int i,j;
      for(i=1; i<=N; i++)
      {
        values = ideal(random(max_val,1,d));

        if(typeof(g)=="poly")
          {val1 = number(substitute(g,vars,values));}
        else{if(typeof(g)=="list") // denominator given in factorized form
        {
          val1 = number(1);
          for(j=1; j<=size(g[1]); j++)
            {val1 = val1 * number(substitute(g[1][j]^g[2][j],vars,values));}
        }
        else {ERROR("wrong type for second argument, expected poly or list");}}

        if(val1==0) {continue;}
        val1 = number(substitute(f,vars,values))/val1;
        val2, div_by_0 = evaluatepfd(dec,values,2);
        if(div_by_0) {continue;}
        if(val1 != val2)
          {return(0);}
      }
      return(1);
    }
  }
  ideal q = dec[1];
  dec = dec[2];

  int m = size(q);
  int jmax,j,ind,k;
  intvec e_max; e_max[m]=0;
  int imax = size(dec);
  for(int i=1; i<=imax; i++)
  {
    jmax = size(dec[i][2]);
    for(j=1; j<=jmax; j++)
    {
      ind = dec[i][2][j];
      e_max[ind] = max(e_max[ind],dec[i][3][j]);
    }
  }
  poly num;
  poly sum_of_numerators = 0;
  intvec e;
  for(i=1; i<=imax; i++)
  {
    e = e_max;
    jmax = size(dec[i][2]);
    for(j=1; j<=jmax; j++)
    {
      ind = dec[i][2][j];
      e[ind] = e[ind]-dec[i][3][j];
    }
    num = dec[i][1];
    for(j=1; j<=m; j++)
      {num = num * q[j]^(e[j]);}
    sum_of_numerators = sum_of_numerators + num;
  }
  // the decomposition is now equal to sum_of_numerators/product(q[i]^e_max[i]) (i from 1 to imax)
  // now: check if this equals f/g:
  if(typeof(g)=="poly")
  {
    list fact = factorize(g);
    ideal q_g = delete(fact[1],1);
    intvec e_g = delete(fact[2],1);
    num = f/fact[1][1];
  }
  else{if(typeof(g)=="list")    //denominator given in factorized form
  {
    ideal q_g = g[1];
    intvec e_g = g[2];
    num = f;
  }
  else {ERROR("wrong type for second argument, expected poly or list");}}

  int m_g = size(q_g);
  int expon;
  number c;
  for(i=1; i<=m_g; i++)
  {
    j=0;
    for(k=1; k<=m; k++)
    {
      c = leadcoef(q[k])/leadcoef(q_g[i]);
      if(c*q_g[i]==q[k]) {j=k; break;}
    }
    if(j==0)
      {sum_of_numerators = sum_of_numerators*q_g[i]^e_g[i];}
    else
    {
      num = num*(c^e_g[i]);    //fix lead coefficient
      expon = e_g[i]-e_max[j];
      if(expon>0)
        {sum_of_numerators = sum_of_numerators*q[j]^expon;}
      else{if(expon<0)
        {num = num*q[j]^(-expon);}}
    }
  }
  return(sum_of_numerators==num);
}
example
{
  "EXAMPLE:";
  echo=voice;
  ring R = 0,(x,y),dp;
  poly f = x^3+3*x^2*y+2*y^2-x^2+4*x*y;
  poly g = x^2*y*(x-1)*(x-y)^2;

  // partial fraction decomposition of f/g:
  list dec = pfd(f,g);
  // some other decomposition (not equal to f/g):
  list wrong_dec = pfd(f+1,g);

  displaypfd_long(dec);
  list fraction = f,g;

  // exact test:
  checkpfd(fraction,dec);
  checkpfd(fraction,wrong_dec);
  // probabilistic test (evaluation at 10 random points):
  checkpfd(fraction,dec,10);
  checkpfd(fraction,wrong_dec,10);
}

proc evaluatepfd(list dec, ideal values, list #)
"USAGE:   evaluatepfd(dec, values[, mode]);   dec list, values ideal, mode int
RETURN:   the number gotten by substituting the numbers generating the ideal
          @code{values} for the variables in the partial fraction decomposition
          @code{dec}. The list @code{dec} has to have the same structure as the
          output of @ref{pfd}.
       @* @code{mode=1}: raise Error in case of division by 0 (default)
       @* @code{mode=2}: return a second integer which is 1 if the denominator
          becomes 0, and 0 otherwise.
SEE ALSO: pfd
EXAMPLE:  example evaluatepfd; shows an example"
{
  int mode = 1;
  if(size(#)>0) {mode = #[1];}

  ideal vars = maxideal(1); // ideal generated by ring variables
  number val=0;
  number denom;
  int m,i,j;
  for(i=1; i<=size(dec[2]); i++)
  {
    denom = 1;
    m = size(dec[2][i][2]);
    for(j=1; j<=m; j++)
    {
      denom = denom * (number(substitute(dec[1][dec[2][i][2][j]],vars,values)))^dec[2][i][3][j];
      if(denom == 0)
      {
        if(mode==1) {ERROR("division by 0");}
        else {return(0,1);}
      }
    }
    val = val + number(substitute(dec[2][i][1],vars,values))/denom;
  }
  if(mode==1) {return(val);}
  else {return(val,0);}
}
example
{
  "EXAMPLE:";
  echo=voice;
  ring R = 0,(x,y),dp;
  poly f = x+2*y;
  poly g = x^2-y^2;

  // partial fraction decomposition of f/g:
  list dec = pfd(f,g);

  displaypfd_long(dec);
  // evaluation at x=2, y=1:
  ideal values = 2,1;
  evaluatepfd(dec,values);

  // compare: f(2,1)/g(2,1) = (2+2*1)/(2^2-1^1) = 4/3
}

static proc find_entry(def l, def entry)
{
  int n=size(l);
  for(int i=1; i<=n; i++)
    {if(entry==l[i]) {return(i);}}
  return(0);
}

proc displaypfd(list dec)
"USAGE:   displaypfd(dec);   dec list
PURPOSE:  print a partial fraction decomposition @code{dec} in a readable way.
          The list @code{dec} has to have the same structure as the output of
          @ref{pfd}.
SEE ALSO: pfd, displaypfd_long, getStringpfd, getStringpfd_indexed
EXAMPLE:  example displaypfd; shows an example"
{
  ideal q = dec[1];
  dec = dec[2];
  int imax=size(dec);
  int jmax,j;
  string s1,s2;
  for(int i=1; i<=imax; i++)
  {
    s1 = "(" + string(dec[i][1]);
    if(i>1) {s1 = "+ " + s1;}
    else {s1 = "  " + s1;}
    if(size(dec[i][2])>0)
    {
      s2 = ") / (";
      jmax = size(dec[i][2]);
      for(j=1; j<=jmax; j++)
      {
        s2 = s2 + "q" + string(dec[i][2][j]);
        if(dec[i][3][j] != 1) {s2 = s2 + "^" + string(dec[i][3][j]);}
        if(j<jmax) {s2 = s2 + "*";}
      }
      s2 = s2 + ")";
    }
    else {s2 = ")";}

    if(size(s1)+size(s2)>192) {s1=s1[1..(192-size(s2))]; s1 = s1 + "... ";}
    print(s1+s2);
  }
  print("where");
  for(i=1; i<=size(q); i++)
  {
    printf("q%s = %s",i,q[i]);
  }
  if(size(dec)==1) {printf("(%s terms)%n", 1, 0);}
  else {printf("(%s terms)%n", size(dec), 0);}
}
example
{
  "EXAMPLE:";
  echo=voice;
  ring R = 0,(x,y),dp;
  poly f = x^3+3*x^2*y+2*y^2-x^2+4*x*y;
  poly g = x^2*y*(x-1)*(x-y)^2;

  list dec = pfd(f,g);

  displaypfd(dec);
}

proc displaypfd_long(list dec)
"USAGE:   displaypfd_long(dec);   dec list
PURPOSE:  like @ref{displaypfd}, but denominators are written out, not indexed.
SEE ALSO: pfd, displaypfd, getStringpfd, getStringpfd_indexed
EXAMPLE:  example displaypfd_long; shows an example"
{
  ideal q = dec[1];
  dec = dec[2];
  print("  "+getStringFraction(dec[1],q));
  for(int i=2; i<=size(dec); i++)
    {print("+ "+getStringFraction(dec[i],q));}
  if(size(dec)==1) {printf("(%s terms)%n", 1, 0);}
  else {printf("(%s terms)%n", size(dec), 0);}
}
example
{
  "EXAMPLE:";
  echo=voice;
  ring R = 0,(x,y),dp;
  poly f = x^3+3*x^2*y+2*y^2-x^2+4*x*y;
  poly g = x^2*y*(x-1)*(x-y)^2;

  list dec = pfd(f,g);

  displaypfd_long(dec);
}

proc getStringpfd(list dec)
"USAGE:   getStringpfd(dec);   dec list
PURPOSE:  turn a partial fraction decomposition @code{dec} into one string. The
          list @code{dec} has to have the same structure as the output of
          @ref{pfd}.
SEE ALSO: pfd, getStringpfd_indexed, displaypfd, displaypfd_long
EXAMPLE:  example getStringpfd; shows an example"
{
  ideal q = dec[1];
  dec = dec[2];
  string s = getStringFraction(dec[1],q);
  for(int i=2; i<=size(dec); i++)
    {s = s+" + "+getStringFraction(dec[i],q);}
  return(s);
}
example
{
  "EXAMPLE:";
  echo=voice;
  ring R = 0,(x,y),dp;
  poly f = x^3+3*x^2*y+2*y^2-x^2+4*x*y;
  poly g = x^2*y*(x-1)*(x-y)^2;

  list dec = pfd(f,g);

  displaypfd_long(dec);
  getStringpfd(dec);
}

proc getStringpfd_indexed(list dec)
"USAGE:   getStringpfd_indexed(dec);   dec list
PURPOSE:  turn a partial fraction decomposition @code{dec} into one string,
          writing the denominator factors just as @code{q1},@code{q2},... . The
          list @code{dec} has to have the same structure as the output of
          @ref{pfd}.
SEE ALSO: pfd, getStringpfd, displaypfd, displaypfd_long
EXAMPLE:  example getStringpfd_indexed; shows an example"
{
  if(typeof(dec[1])=="ideal") {dec = dec[2];}
  string s = getStringFraction_indexed(dec[1]);
  for(int i=2; i<=size(dec); i++)
    {s = s+" + "+getStringFraction_indexed(dec[i]);}
  return(s);
}
example
{
  "EXAMPLE:";
  echo=voice;
  ring R = 0,(x,y),dp;
  poly f = x^3+3*x^2*y+2*y^2-x^2+4*x*y;
  poly g = x^2*y*(x-1)*(x-y)^2;

  list dec = pfd(f,g);

  displaypfd(dec);
  getStringpfd_indexed(dec);
}

static proc getStringFraction(list fraction, ideal q)
{
  int jmax,j;
  string s = "(" + string(fraction[1]);
  if(size(fraction[2])>0)
  {
    s = s + ")/(";
    jmax = size(fraction[2]);
    for(j=1; j<=jmax; j++)
    {
      s = s + "(" + string(q[fraction[2][j]]) + ")";
      if(fraction[3][j] != 1) {s = s + "^" + string(fraction[3][j]);}
      if(j<jmax) {s = s + "*";}
    }
  }
  s = s + ")";
  return(s);
}

static proc getStringFraction_indexed(list fraction)
{
  int jmax,j;
  string s = "(" + string(fraction[1]);
  if(size(fraction[2])>0)
  {
    s = s + ")/(";
    jmax = size(fraction[2]);
    for(j=1; j<=jmax; j++)
    {
      s = s + "q" + string(fraction[2][j]);
      if(fraction[3][j] != 1) {s = s + "^" + string(fraction[3][j]);}
      if(j<jmax) {s = s + "*";}
    }
  }
  s = s + ")";
  return(s);
}

proc readInputTXT(def file__, list #)
"USAGE:   readInputTXT(file[, mode]), file string, mode int
          readInputTXT(filelist[, mode]), filelist list, mode int
PURPOSE:  read matrix of rational functions from a txt-file and turn it into a
          matrix (i.e. a list of lists) of pairs of polynomials (numerators and
          denominators). The string @code{file} should be the [directory +] name
          of the file in the form \"@code{<path-to-file>/<filename>.txt}\".
       @* The input file should be a list of lists separated by the characters
          \"{\", \"}\" and \",\". Example:
       @* \"{{(x+y)/(x^2-x*y), -(x^2*y+1)/(y), x^2}, {(x+1)/y, y/x, 0}}\"
       @* Each rational function has to be an expression of the form \"a\",
          \"(a)/(b)\", \"(b)^(-n)\" or \"(a)*(b)^(-n)\" where \"a\",\"b\" stand
          for polynomials (i.e. strings, that can be parsed as a polynomial with
          the @code{execute} command) and \"n\" stands for a positive integer. A
          minus sign \"-\" followed by such an expression is also allowed.
       @* IMPORTANT: The strings \"a\",\"b\" must NOT contain the symbol \"/\".
          (So in case the coefficient field is the rationals, all denominators
          in the coefficients of numerator and denominator polynomials should be
          cleared.)
       @* The file should contain less than 2^31 characters (filesize < 2 GB).
          For bigger files the matrix should be split row-wise into multiple
          matrices and saved in different files (each smaller than 2 GB). A list
          of the filenames (in the right order) can then be given as first
          argument instead.
       @* Also the basering has to match the variable names used in the
          input file(s).
       @* @code{mode=1} (default): save result to an ssi-file of the same name
       @* @code{mode=2}: return result
       @* @code{mode=3}: save to ssi-file AND return result
SEE ALSO: pfdMat
"
{
  system("--ticks-per-sec",1000);
  short = 0;
  if(!defined(basering))
    {ERROR("no basering defined!");}
  int left__,right__,pos1__,mid__,pos2__,tmp__,i__,j__,k__,t__,tt__,depth__;
  int mode__=1;
  if(size(#)>0) {mode__=#[1];}
  if(typeof(file__)=="list") // list of filenames given --> apply function to each
  {                          // file and concatenate the resulting matrices
    int n = size(file__);
    list mat__ = list();
    for(i__=1;i__<=n;i__++)
    {
      dbprint(sprintf("  file %s of %s:",i__,n));
      if(find(file__[i__],".txt")==0) {ERROR("wrong file type, expected txt");}
      printlevel = printlevel+1;
      mat__ = mat__ + readInputTXT(file__[i__],2);
      printlevel = printlevel-1;
    }

    if(mode__==2) {return(mat__);}

    string filename__ = file__[1][1,find(file__[1],".txt")-1];
    dbprint("  saving to file "+filename__+".ssi "); t__ = rtimer;
    write("ssi:w "+filename__+".ssi", mat__);
    dbprint(sprintf("  done! (%s ms)", rtimer-t__));

    if(mode__==3) {return(mat__);}
  }
  if(typeof(file__)!="string")
    {ERROR("wrong type for first argument (expected string or list)");}
  if(find(file__,".txt")==0) {ERROR("wrong file type, expected txt");}

  dbprint("  reading file "); t__=rtimer;
  string data__ = read(":r "+file__);
  dbprint(sprintf("  done! (%s ms)", rtimer-t__));

  dbprint("  processing input "); t__ = rtimer;
  left__ = find(data__,"{");
  right__ = find(data__,"}");
  tmp__ = find(data__,"{",left__+1);
  if(left__<tmp__ && tmp__<right__) {left__ = tmp__;}

  int finished__,n__;
  poly p__,q__;
  list row__,mat__;
  string s__,ss__;
  i__=0;
  while(1)
  {
    i__++;
    tt__ = rtimer;
    row__ = list();
    pos2__ = left__;
    finished__ = 0;
    j__=0;
    while(not finished__)
    {
      j__++;
      s__ = "";

      pos1__ = pos2__+1;
      pos2__ = find(data__,",",pos1__);
      if(pos2__==0 || pos2__>right__)    // end of row
      {
        finished__ = 1;
        pos2__ = right__;
      }
      s__ = data__[pos1__,pos2__-pos1__];
      mid__ = find(s__,"/");
      if(mid__==0)
      {
        tmp__ = find(s__,"^(-");
        if(tmp__==0)    //no denominator
        {
          execute("p__=" + s__);
          q__=1;
          row__[j__] = list(p__,q__);
          continue;
        }
        else    // denominator is given by a negative exponent
        {
          if(find(s__,"^(-",tmp__+1)>0)
            {ERROR(sprintf("invalid syntax in (%s,%s)-th entry:"
                              +" more than one negative exponent",i__,j__));}
          for(k__=tmp__+3; s__[k__]!=")"; k__++) {}
          execute("n__=" + s__[tmp__+3,k__-tmp__-3]);
          while(k__<size(s__))
          {
            k__++;
            if(s__[k__]!=" ")
            {ERROR(sprintf("invalid syntax in (%s,%s)-th entry",i__,j__));}
          }
          s__ = s__[1,tmp__-1];
          depth__=0;
          for(k__=tmp__-1; 1; k__--)
          {
            ss__ = s__[k__];
            if(ss__ == ")") {depth__++;}
            else{if(ss__ == "(") {depth__--;}}
            if(depth__==0) {break;}
          }
          if(k__>1)
          {
            while(1)
            {
              if(k__==0)
                {ERROR(sprintf("invalid syntax in (%s,%s)-th entry",i__,j__));}
              k__--;
              ss__ = s__[k__];
              if(ss__=="*" || ss__==" " || ss__=="-") {break;}
            }
          }
          s__ = s__[1,k__] + "1/" + s__[k__+1,size(s__)-k__] + "^" + string(n__);
          mid__ = k__+2;    // position of the character "/"
        }
      }
      if(find(s__,"/",mid__+1)>0)
        {ERROR(sprintf("invalid syntax in (%s,%s)-th entry:%n"
        +"no '/' allowed in the string representing the polynomials",i__,j__,0));}

      execute("p__=" + fixBrackets(s__[1,mid__-1]));
      execute("q__=" + fixBrackets(s__[mid__+1,size(s__)-mid__]));

      row__[j__] = list(p__,q__);
    }
    mat__[i__] = row__;    // append row to matrix
    dbprint(sprintf("    row %s done! (%s ms)",i__,rtimer-tt__));

    left__ = find(data__,"{",right__);
    if(left__==0) {break;}
    right__ = find(data__,"}",left__);
  }
  dbprint(sprintf("  done! (%s ms)", rtimer-t__));

  if(mode__==2) {return(mat__);}

  string filename__ = file__[1,find(file__,".txt")-1];
  dbprint("  saving to file "+filename__+".ssi "); t__ = rtimer;
  write("ssi:w "+filename__+".ssi", mat__);
  dbprint(sprintf("  done! (%s ms)", rtimer-t__));

  if(mode__==3) {return(mat__);}
}

static proc fixBrackets(string data)
{
  int pos=0;
  int left_brackets =0;
  int right_brackets=0;
  int n=size(data);
  while(pos<n)
  {
    pos = find(data,"(",pos+1);
    if(pos==0) {break;}
    left_brackets++;
  }
  pos=0;
  while(pos<n)
  {
    pos = find(data,")",pos+1);
    if(pos==0) {break;}
    right_brackets++;
  }
  int difference = left_brackets-right_brackets;
  if(difference>0)
  {
    for(int i=1; i<=difference; i++) {data = data+")";}
  }
  if(difference<0)
  {
    for(int i=1; i<=-difference; i++) {data = "("+data;}
  }
  return(data);
}

static proc pfdWrap(poly f, def g, int i, int j, link logfile, int output_mode)
{
  system("--ticks-per-sec",1000);
  if(output_mode>3)
    {write("ssi:w pfd_results_"+string(i)+"_"+string(j)
           +".ssi","task started, but not finished yet");}
  int t0 = rtimer;
  list result = pfd(f,g);
  fprintf(logfile,"_[%s,%s]: %s ms",i,j,rtimer-t0);
  if(output_mode>3)
    {write("ssi:w pfd_results_"+string(i)+"_"+string(j)+".ssi",result);}
  return(result);
}

static proc testEntry(int i, int j, list fraction, list dec, link logfile, int N);
{
  int t=rtimer;
  system("--ticks-per-sec",1000);
  int result = checkpfd(fraction,dec,N);
  if(result==1) {fprintf(logfile, " _[%s,%s]: correct (%s ms)",i,j,rtimer-t);}
  else          {fprintf(logfile, " _[%s,%s]: WRONG! (%s ms)", i,j,rtimer-t);}
  return(result);
}

proc pfdMat(def infile, list #)
"USAGE:   pfdMat(file[, dotest, ignore_nonlin, output_mode, parallelize]);
          file string, dotest,ignore_nonlin,output_mode,parallelize int
PURPOSE:  apply @code{pfd} to all entries of a matrix of rational functions
          saved in a txt-file. The string @code{file} should be the
          [directory +] name of the file.
       @* The input file can either be a txt-file or an ssi-file created with
          @code{readInputTXT}. In case of a txt-file, the base ring has to match
          and the matrix has to be in the same format specified in
          @ref{readInputTXT}. Also, txt-files that are bigger than 2 GB should
          be split as described for @code{readInputTXT} and a list of the
          filenames can be given as first argument instead.
       @* The result is saved in multiple txt- (and ssi-) files (see below)
          within the directory of the input file.
       @* Also a logfile is created, which protocols the memory used and the
          runtimes of @code{pfd} for each matrix entry in real-time.

       @* There are also 4 optional arguments:
       @* If @code{dotest} is nonzero, test the results with checkpfd:
       @* @code{dotest<0} (default): exact test (may be slow),
       @* @code{dotest>0}: do this amount of probabilistic tests for each entry
                           (see @ref{checkpfd}).

       @* If @code{ignore_nonlin} is nonzero (default), for each denominator,
          the nonlinear factors in the factorization are removed before applying
          @code{pfd} (and added back in in the output files).

       @* If @code{parallelize} is nonzero (default), the decompositions are
          calculated in parallel using @ref{parallel_lib}.

       @* The parameter @code{output_mode} controls the output files created:
       @* @code{output_mode=1} (default): The result consists of two files:
          @code{<filename>_pfd_indexed.txt} contains the matrix of all
          decompositions (as list of lists separated by the characters \"{\",
          \"}\" and \",\") where all the denominators are written in factorized
          form depending on irreducible factors @code{q1}, @code{q2}, ... .
          The file @code{<filename>_denominator_factors.txt} lists all the
          polynomials @code{q1}, @code{q2}, ... .
       @* @code{output_mode=2}: Additionally to mode 1, the file
          @code{<filename>_pfd.txt} is created, which also contains the matrix
          of decompositions but the factors in the denominators are written out.
       @* @code{output_mode=3}: Additionally to mode 2, the result and some
          intermediate results are saved as SINGULAR objects in ssi-files:
       @* @code{<filename>.ssi}: contains the result of @code{readInputTXT} in
          case a txt-file was given as input.
       @* @code{<filename>_factorized_denominators.ssi}: like the first file,
          but the denominators are saved in factorized form, that is as a list
          of an ideal of irreducible non constant polynomials and an intvec of
          exponents.
       @* @code{<filename>_linear_part.ssi} (only if @code{ignore_nonlin} is
          nonzero): like the previous file, but all the irreducible denominator
          factors are removed
       @* @code{<filename>_non_linear_factors.ssi} (only if @code{ignore_nonlin}
          is nonzero): a list of an ideal @code{p} generated by irreducible
          polynomials and a matrix (list of lists) of the nonlinear denominator
          factors of each entry of the input matrix. These are represented as
          lists of an intvec of indices @code{i} for which @code{p[i]} occurs
          as a (nonlinear) factor in the denominator and an intvec containing
          the exponents of those factors.
       @* @code{<filename>_pfd.ssi}: a list, where the first entry is an ideal
          @code{q} of denominator factors and the second entry is a matrix (as
          list of lists) containing the decompositions, each of which is a list
          of terms, where a term is represented as in the result of @ref{pfd}
          by a list containing
       @* 1) the numerator polynomial
       @* 2) an intvec of indices @code{i} for which @code{q[i]} occurs
             as a factor in the denominator
       @* 3) an intvec containing the exponents of those irreducible factors.
       @* IMPORTANT: If @code{ignore_nonlin} is nonzero, this file contains the
          decompositions of the entries of the matrix in
          @code{<filename>_linear_part.ssi}. Thus the nonlinear factors, are
          NOT contained in this file.
       @* @code{output_mode=4}: Additionally to mode 3, the direct output of
          each call of @code{pfd} is saved in separate ssi-files called
          @code{pfd_results_i_j.ssi} where i,j are the matrix indices. This
          creates a lot of files, but may be useful in case the algorithm does
          not terminate in time for every matrix entry. Other than the files
          created in mode 1-3, these files are saved in the current directory,
          rather than the directory of the input file.
SEE ALSO: readInputTXT, pfd, checkpfd, checkpfdMat
"
{
  system("--ticks-per-sec",1000);
  short = 0;
  int dotest,ignore_nonlin,output_mode,parallelize = -1,1,1,1;
  if(size(#)>0) {dotest = #[1];}
  if(size(#)>1) {ignore_nonlin = #[2];}
  if(size(#)>2) {output_mode = #[3];}
  if(size(#)>3) {parallelize = #[4];}

  int i,j,k,l,ind;
  list arguments,results;

  dbprint(newline+"reading data "); int t0=rtimer;

  if(typeof(infile)=="list")
  {
    printlevel = printlevel+1;
    if(output_mode>2) {list mat = readInputTXT(infile,3);}
    else {list mat = readInputTXT(infile,2);}
    printlevel = printlevel-1;
    int pos=find(infile[1],".txt");
    string filename = infile[1][1,pos-1];
  }
  else
  {
    int pos=find(infile,".txt");
    if(pos!=0)
    {
      printlevel = printlevel+1;
      if(output_mode>2) {list mat = readInputTXT(infile,3);}
      else {list mat = readInputTXT(infile,2);}
      printlevel = printlevel-1;
    }
    else
    {
      pos=find(infile,".ssi");
      if(pos!=0) {list mat = read("ssi:r "+infile);}
      else {ERROR("invalid file type, expected ssi or txt");}
    }
    string filename = infile[1,pos-1];
  }

  link qfile = ":w "+filename+"_denominator_factors.txt";
  int n = size(mat);
  int m = size(mat[1]);
  dbprint(sprintf("done! (%s ms)",rtimer-t0));

  if(typeof(mat[1][1][2])!="list")    // denominators are not yet factorized
  {
    dbprint("factorizing the denominators "); t0=rtimer;
    printlevel = printlevel+1;
    mat = FactDenom(mat);
    printlevel = printlevel-1;
    if(output_mode>2)
      {write("ssi:w "+filename+"_factorized_denominators.ssi",mat);}
    dbprint(sprintf("done! (%s ms)",rtimer-t0));
  }

  if(ignore_nonlin)
  {
    dbprint("removing nonlinear denominator factors before pfd is applied");
    list nonlin_denom_factors;
    ideal p;
    printlevel = printlevel+1;
    mat, nonlin_denom_factors, p = removeNonlinearFactors(mat, filename);
    printlevel = printlevel-1;
    if(output_mode>2)
    {
      dbprint("saving nonlinear factors to "+filename+"_non_linear_factors.ssi ");
      t0 = rtimer;
      write("ssi:w "+filename+"_non_linear_factors.ssi", list(p,nonlin_denom_factors));
      dbprint(sprintf("done! (%s ms)",rtimer-t0));
      dbprint("saving input matrix without the nonlinear factors to "
              +filename+"_linear_part.ssi ");
      t0 = rtimer;
      write("ssi:w "+filename+"_linear_part.ssi", mat);
      dbprint(sprintf("done! (%s ms)",rtimer-t0));
    }
  }

  if(parallelize)
  {
    dbprint("creating tasks "); t0=rtimer;
    write(":w "+filename+"_pfdMat_logfile.txt","finished matrix entries with runtimes "
          +"(calculated in parallel on "+string(getcores())+" cores):");
    link logfile = ":a "+filename+"_pfdMat_logfile.txt";
    for(i=1; i<=n; i++)
    {
      for(j=1; j<=m; j++)
      {
        ind = m*(i-1)+j;
        arguments[ind] = list(mat[i][j][1],mat[i][j][2],i,j,logfile,output_mode);
      }
    }
    dbprint(sprintf("done! (%s ms)",rtimer-t0));

    dbprint("applying pfd to each matrix entry "); t0 = rtimer;
    results = parallelWaitAll("pfdWrap",arguments);
    arguments = list();
    dbprint(sprintf("done! (%s ms)",rtimer-t0));

    write(logfile,"decomposition: "+string(rtimer-t0)+" ms and "+string(memory(2))
                   +" Byte Memory max. (after calling pfd on each matrix entry)");

    dbprint("writing results in matrix shape "); t0 = rtimer;
    list dec_mat;
    for(i=1; i<=n; i++)
    {
      dec_mat[i] = list();
      for(j=1; j<=m; j++)
      {
        ind = m*(i-1)+j;
        dec_mat[i][j] = results[ind];
      }
    }
    results = list();
    dbprint(sprintf("done! (%s ms)",rtimer-t0));
  }
  else
  {
    dbprint("applying pfd to each matrix entry "); t0=rtimer;
    write(":w "+filename+"_pfdMat_logfile.txt",
          "finished matrix entries with runtimes (no parallelization):");
    link logfile = ":a "+filename+"_pfdMat_logfile.txt";
    list dec_mat;
    for(i=1; i<=n; i++)
    {
      dec_mat[i] = list();
      for(j=1; j<=m; j++)
      {
        dec_mat[i][j] = pfdWrap(mat[i][j][1],mat[i][j][2],i,j,logfile,output_mode);
      }
    }
    dbprint(sprintf("done! (%s ms)",rtimer-t0));
    write(logfile,"decomposition: "+string(rtimer-t0)+" ms and "+string(memory(2))
          +" Byte Memory max. (after calling pfd on each matrix entry)");
 }

  dbprint("making one single list of denominator factors "); t0 = rtimer;
  ideal q,new_q;
  intvec dict;
  for(i=1; i<=n; i++)
  {
    for(j=1; j<=m; j++)
    {
      new_q = dec_mat[i][j][1];
      dec_mat[i][j] = dec_mat[i][j][2];
      dict = 0:0;
      for(k=1; k<=size(new_q); k++)
      {
        ind = find_entry(q,new_q[k]);
        if(ind==0)
        {
          ind = size(q)+1;
          q[ind] = new_q[k];
        }
        dict[k] = ind;
      }
      for(k=1; k<=size(dec_mat[i][j]); k++)
      {
        if(size(dec_mat[i][j][k][2])>0)
          {dec_mat[i][j][k][2] = intvec(dict[dec_mat[i][j][k][2]]);}
      }
    }
    dbprint(sprintf("  row %s complete!",i));
  }
  dbprint(sprintf("done! (%s ms)",rtimer-t0));

  if(output_mode>2)
  {
    dbprint("saving result to "+filename+"_pfd.ssi "); t0 = rtimer;
    write("ssi:w "+filename+"_pfd.ssi", list(q,dec_mat));
    dbprint(sprintf("done! (%s ms)",rtimer-t0));
  }

  if(ignore_nonlin)
  {
    ind = size(q);
    for(i=1; i<=size(p); i++) {q[ind+i]=p[i];} // add nonlin. polynomials to q
    for(i=1; i<=n; i++)
    {
      for(j=1; j<=m; j++)
      {
        nonlin_denom_factors[i][j][1] = nonlin_denom_factors[i][j][1]+ind; // adjust indices
      }
    }
  }

  if(ignore_nonlin)
    {dbprint(sprintf("creating readable .txt-files (including the nonlinear factors again)"));}
  else
    {dbprint(sprintf("creating readable .txt-files "));}
  t0 = rtimer;
  dbprint(" indexed ("+filename+"_pfd_indexed.txt):");
  printlevel = printlevel+1;
  if(ignore_nonlin)
    {saveResultTXT_indexed(dec_mat, filename+"_pfd_indexed", nonlin_denom_factors);}
  else {saveResultTXT_indexed(dec_mat, filename+"_pfd_indexed");}
  printlevel = printlevel-1;
  for(i=1; i<=size(q); i++) {fprintf(qfile, "q%s = %s;", i, q[i]);}
  if(output_mode>1)
  {
    dbprint(" denominators written out ("+filename+"_pfd.txt):");
    printlevel = printlevel+1;
    if(ignore_nonlin)
      {saveResultTXT(dec_mat, q, filename+"_pfd", nonlin_denom_factors);}
    else {saveResultTXT(dec_mat, q, filename+"_pfd");}
    printlevel = printlevel-1;
  }
  dbprint(sprintf("done! (%s ms)",rtimer-t0));

  if(dotest)
  {
    if(dotest<0)
      {dbprint("checking for correctness (exact test) ");}
    else
      {dbprint(sprintf("checking for correctness (%s random evaluations per entry) ",dotest));}
    t0 = rtimer;
    if(parallelize)
    {
      for(i=1; i<=n; i++)
      {
        for(j=1; j<=m; j++)
        {
          arguments[m*(i-1)+j]=list(i,j,mat[i][j],list(q,dec_mat[i][j]),logfile,dotest);
        }
      }
      results = parallelWaitAll("testEntry",arguments);
    }
    else
    {
      for(i=1; i<=n; i++)
      {
        for(j=1; j<=m; j++)
        {
          results[m*(i-1)+j]=testEntry(i,j,mat[i][j],list(q,dec_mat[i][j]),logfile,dotest);
        }
      }
    }
    dbprint(sprintf("%s out of %s = %sx%s decompositions are correct! (%s ms)%n",
            sum(results),n*m,n,m,rtimer-t0,0));
    write(logfile,"checking for correctness: "+string(rtimer-t0)+" ms and "
          +string(memory(2))+" Byte Memory max. (at the end of pfdMat), "
          +string(sum(results))+" correct out of "+string(n*m));
  }
}

static proc FactDenom(list mat)
{
  system("--ticks-per-sec",1000);
  int i,j,k,ind,t,counter;
  int n = size(mat);
  int m = size(mat[1]);
  list denom = list();
  for(i=1; i<=n; i++)
  {
    denom[i] = list();
    for(j=1; j<=m; j++)
    {
      denom[i][j] = mat[i][j][2];
      mat[i][j][2] = list(ideal(),intvec(0:0));
    }
  }
  int expon;
  list fact;
  number lcoeff;
  int timeout = 60000;
  int finished = 0;
  list arguments,results;
  ideal q;
  while(!finished)
  {
    t = rtimer;
    for(i=1; i<=n; i++)   // create argument list
    {
      for(j=1; j<=m; j++)
        {arguments[m*(i-1)+j] = list(denom[i][j]);}
    }

    results = parallelWaitAll("factorize",arguments,timeout);
    arguments = list();

    for(i=1; i<=n; i++)   // update q, mat, denom
    {
      for(j=1; j<=m; j++)
      {
        ind = m*(i-1)+j;
        if(typeof(results[ind]) != "none")
        {
          fact = results[ind];
          for(k=2; k<=size(fact[1]); k++)
          {
            lcoeff = leadcoef(fact[1][k]);
            fact[1][k] = fact[1][k]/lcoeff;
            fact[1][1] = fact[1][1]*(lcoeff^fact[2][k]); // polynomial is monic (thus unique)
            lcoeff = content(fact[1][k]);
            fact[1][k] = fact[1][k]/lcoeff;
            fact[1][1] = fact[1][1]*(lcoeff^fact[2][k]); // polynomial has nice coefficients

            // add a new factor to q:
            ind = find_entry(q,fact[1][k]);
            if(ind==0) {q[size(q)+1] = fact[1][k];}
            // complete the factorization of the i,j-th denominator:
            ind = find_entry(mat[i][j][2][1],fact[1][k]);
            if(ind==0)
            {
              ind = size(mat[i][j][2][1])+1;
              mat[i][j][2][1][ind] = fact[1][k];
              mat[i][j][2][2][ind] = fact[2][k];
            }
            else
            {
              mat[i][j][2][2][ind] = mat[i][j][2][2][ind] + fact[2][k];
            }
          }
          mat[i][j][1] = mat[i][j][1]/fact[1][1];
          denom[i][j] = 1;
        }
      }
    }
    results = list();

    finished = 1;
    counter = n*m;
    for(i=1; i<=n; i++) // factorize by any known factors (from q)
    {
      for(j=1; j<=m; j++)
      {
        for(k=1; k<=size(q); k++)
        {
          expon=0;
          while(reduce(denom[i][j],q[k])==0)
          {
            denom[i][j] = denom[i][j]/q[k];
            expon++;
          }
          if(expon>0)
          {
            ind = size(mat[i][j][2][1])+1;
            mat[i][j][2][1][ind] = q[k];
            mat[i][j][2][2][ind] = expon;
          }
        }
        if(deg(denom[i][j])==0)
        {
          mat[i][j][1] = mat[i][j][1]/denom[i][j];
          denom[i][j] = poly(1);
        }
        if(denom[i][j]!=poly(1)) {finished = 0; counter--;}
      }
    }
    timeout = timeout*2;
    dbprint(sprintf("  %s out of %s denominators factorized completely (%s ms)",
                                                         counter,n*m,rtimer-t));
  }
  return(mat);
}

static proc saveResultTXT_indexed(list dec, string filename, list #)
{
  // expect dec = list(list(dec11, dec12, ...), list(dec21, dec22, ...), ...),
  // where dec11, dec12, ... are decompositions of form list(list(poly, intvec, intvec), ...)
  system("--ticks-per-sec",1000);
  list nonlinFactors = list();
  if(size(#)>0) {nonlinFactors = #;}
  link file = ":w "+filename+".txt";

  int i,j;
  int t;
  string s="{";
  int n=size(dec);
  int m=size(dec[1]);
  int k;
  for(i=1; i<=n; i++)
  {
    t = rtimer;
    s = s+"{";
    for(j=1; j<=m; j++)
    {
      if(size(nonlinFactors)>0)
      {
        if(size(nonlinFactors[i][j][1])>0)
        {
          s = s + getStringFraction_indexed(list(poly(1))+nonlinFactors[i][j])
                + " * (" + getStringpfd_indexed(dec[i][j]) + "), ";
          j++; continue;
        }
      }
      s = s + getStringpfd_indexed(dec[i][j]) + ", ";
    }
    k = size(s);
    s[k-1] = "}"; s[k] = ","; s = s+" ";
    dbprint(sprintf("  row %s done! (%s ms)",i,rtimer-t));
  }
  k = size(s);
  s[k-1] = "}"; s[k] = " ";
  write(file,s);
}

static proc saveResultTXT(list dec, ideal q, string filename, list #)
{
  // expect dec = list(list(dec11, dec12, ...), list(dec21, dec22, ...), ...),
  // where dec11, dec12, ... are decompositions of form list(list(poly, intvec, intvec), ...)
  system("--ticks-per-sec",1000);
  list nonlinFactors = list();
  if(size(#)>0) {nonlinFactors = #;}
  link file = ":w "+filename+".txt";

  int i,j;
  int t;
  string s="{";
  int n=size(dec);
  int m=size(dec[1]);
  int k;
  for(i=1; i<=n; i++)
  {
    t = rtimer;
    s = s+"{";
    for(j=1; j<=m; j++)
    {
      if(size(nonlinFactors)>0)
      {
        if(size(nonlinFactors[i][j][1])>0)
        {
          s = s + getStringFraction(list(poly(1))+nonlinFactors[i][j],q)
                + " * (" + getStringpfd(list(q,dec[i][j])) + "), ";
          j++; continue;
        }
      }
      s = s + getStringpfd(list(q,dec[i][j])) + ", ";
    }
    k = size(s);
    s[k-1] = "}"; s[k] = ","; s = s+" ";
    dbprint(sprintf("  row %s done! (%s ms)",i,rtimer-t));
  }
  k = size(s);
  s[k-1] = "}"; s[k] = " ";
  write(file,s);
}

static proc removeNonlinearFactors(list fractions, string filename)
{
  int n=size(fractions);
  int m=size(fractions[1]);
  int i,j,k,t0,ind;
  ideal p;
  list nonlin_denom_factors;
  intvec factors,exponents;
  list fac;
  for(i=1; i<=n; i++)
  {
    t0 = rtimer;
    nonlin_denom_factors[i] = list();
    for(j=1; j<=m; j++)
    {
      fac = fractions[i][j][2];
      factors = 0:0;
      exponents = 0:0;
      for(k=1; k<=size(fac[1]); k++)
      {
        if(deg(poly(fac[1][k]))>1)
        {
          // add the nonlin. factor fac[1][k] to p if necessary:
          if(size(p)==0) {ind = 1; p[ind]=fac[1][k];}
          else
          {
            for(ind=1;1;ind++)
            {
              if(p[ind]==fac[1][k]) {break;}
              if(ind==size(p))
              {
                ind++;
                p[ind]=fac[1][k];
                break;
              }
            }
          }
          factors[size(factors)+1] = ind;
          exponents[size(exponents)+1] = fac[2][k];
          fac[1] = delete(fac[1],k);
          fac[2] = delete(fac[2],k);
          continue;
        }
      }
      fractions[i][j][2] = fac;
      nonlin_denom_factors[i][j] = list(factors,exponents);
    }
    dbprint(sprintf("  row %s done! (%s ms)", i, rtimer-t0));
  }
  return(fractions, nonlin_denom_factors, p);
}

proc checkpfdMat(def input, string output, string qfile, list #)
"USAGE:   checkpfdMat(input, output, denomFactors[, N, parallelize]);
          input,output,denomFactors string, N,parallelize int
PURPOSE:  test the output files of @code{pfdMat} for correctness. Input and
          output (indexed) txt-files have to be given as strings in the form
          \"@code{<path-to-file>/<filename>.txt}\". The output should be indexed
          (that is the output file ending in @code{..._pfd_indexed.txt}) and
          @code{denomFactors} has to be the file containing the denominator
          factors @code{q1}, @code{q2}, ... (the txt-file ending in
          @code{..._denominator_factors.txt}).
       @* As for @code{readInputTXT} and @code{pfdMat}, the basering has to
          match the variable names used in the input file, which has to be in
          the same format specified in @ref{readInputTXT}. Also, files bigger
          than 2 GB have to be split as described for @code{readInputTXT} and a
          list of filenames can be given as first argument instead.
       @* If a positive integer N is given, the test is done probabilistically by
          evaluation at N random points for each entry of the matrix. If N is
          nonpositive (default), the fractions in the decompositions will be
          expanded symbolically and compared to the input (may be slower).
       @* If @code{parallelize} is nonzero (default), the tests are run in
          parallel using @ref{parallel_lib}.
       @* The result is printed and as in @code{pfdMat} a logfile is created
          showing the results for each matrix entry.
SEE ALSO: readInputTXT, pfd, checkpfd, pfdMat
"
{
  system("--ticks-per-sec",1000);
  short = 0;
  int t0;
  int N=0;
  int parallelize=1;
  if(size(#)==1)
  {
    if(typeof(#[1])=="int") {N=#[1];}
    else {ERROR("invalid argument type: "+typeof(#[1]));}
  }
  if(size(#)==2)
  {
    if(typeof(#[1])=="int" && typeof(#[2])=="int")
      {N=#[1]; parallelize=#[2];}
    else {ERROR("invalid argument types: "+typeof(#[1])+", "+typeof(#[2]));}
  }
  if(size(#)>2) {ERROR("too many arguments");}

  dbprint(newline+"reading input file:");
  printlevel = printlevel+1;
  list frac = readInputTXT(input,2);
  printlevel = printlevel-1;
  if(typeof(input)=="string") {string filename=input;}
  if(typeof(input)=="list")   {string filename=input[1];}
  filename = filename[1,find(filename,".txt")-1];

  dbprint("factorizing the denominators "); t0=rtimer;
  printlevel = printlevel+1;
  frac = FactDenom(frac);
  printlevel = printlevel-1;
  dbprint(sprintf("done! (%s ms)",rtimer-t0));

  dbprint("reading output files:"); t0=rtimer;
  dbprint(" reading list of denominator factors from "+qfile);
  ideal q;
  q = readQfileTXT(qfile);
  dbprint(" done!");

  dbprint(" reading (indexed) output decompositions ");
  list dec,nonlin;
  printlevel = printlevel+1;
  dec,nonlin = readOutputTXT(output);
  printlevel = printlevel-1;

  if(parallelize)
    {dbprint(sprintf("done! (%s ms)%n%ncreating tasks",rtimer-t0,0)); t0=rtimer;}
  else
  {
    dbprint(sprintf("done! (%s ms)%n",rtimer-t0,0));
    if(N<=0)
      {dbprint("checking for correctness (exact test) ");}
    else
      {dbprint(sprintf("checking for correctness (%s random evaluations per entry) ",N));}
    t0=rtimer;
  }

  fprintf(":w "+filename+"_checkpfdMat_logfile.txt","Input file (matrix of rational functions):"
          +" %s%nOutput file (decompositions): %s%nlist of all denominator factors:"
          +" %s%n%nResults of checkpfdMat:",input,output,qfile,0);
  link logfile = ":a "+filename+"_checkpfdMat_logfile.txt";

  int n=size(frac);
  int m=size(frac[1]);
  int i,j,k,ind;
  list arguments;
  for(i=1;i<=n;i++)
  {
    for(j=1;j<=m;j++)
    {
      for(k=1;k<=size(nonlin[i][j][1]);k++)
      {
        ind = find_entry(frac[i][j][2][1],q[nonlin[i][j][1][k]]);
        if(ind==0) {ERROR("nonlinear factors are wrong");}
        if(frac[i][j][2][2][ind]!=nonlin[i][j][2][k])
          {ERROR("nonlinear factors are wrong");}
        frac[i][j][2][1] = delete(frac[i][j][2][1],ind);
        frac[i][j][2][2] = delete(frac[i][j][2][2],ind);
      }
      if(parallelize)
        {arguments[(i-1)*m+j] = list(i,j,frac[i][j],list(q,dec[i][j]),logfile,N);}
      else
        {results[(i-1)*m+j] = testEntry(i,j,frac[i][j],list(q,dec[i][j]),logfile,N);}
    }
  }

  if(parallelize)
  {
    dbprint(sprintf("done! (%s ms)%n",rtimer-t0,0));
    if(N<=0)
      {dbprint("checking for correctness (exact test) ");}
    else
      {dbprint(sprintf("checking for correctness (%s random evaluations per entry) ",N));}
    t0=rtimer;
    list results = parallelWaitAll("testEntry",arguments);
  }
  else {dbprint(sprintf("done! (%s ms)",rtimer-t0));}
  dbprint(sprintf("%s out of %s = %sx%s decompositions are correct! (%s ms)%n",
          sum(results),n*m,n,m,rtimer-t0,0));
  fprintf(logfile,"%s out of %s = %sx%s decompositions are correct! (%s ms)%n",
          sum(results),n*m,n,m,rtimer-t0,0);
}

static proc readOutputTXT(string filename__)
{
  system("--ticks-per-sec",1000);
  dbprint("  reading matrix of decompositions from file "+filename__);
  int t__=rtimer;
  int tt__;
  string data__ = read(":r "+filename__);
  dbprint(sprintf("  done! (%s ms)", rtimer-t__));

  dbprint("  processing input "); t__ = rtimer;
  list mat__;
  list nonlin__=list();
  int nonlin_factors_seperated__=0;
  if(find(data__," * ")>0) {nonlin_factors_seperated__=1;}
  int left__,right__=0,0;
  left__ = find(data__,"{");
  int pos1__,pos2__=0,0;
  int p1__,p2__;
  int tmp__,tmp2__,depth__;
  int i__,j__,k__,l__,max__;
  intvec factors__,exponents__;
  poly numerator__;
  string s__,ss__;
  for(i__=1;1;i__++)
  {
    tt__ = rtimer;
    left__ = find(data__,"{",left__+1);
    if(left__==0) {break;}
    right__ = find(data__,"}",left__);
    mat__[i__]=list();
    nonlin__[i__]=list();
    pos2__=left__;
    for(j__=1;pos2__<right__&&pos2__>0;j__++)
    {
      pos1__ = pos2__+1;
      pos2__ = find(data__,",",pos1__);
      if(pos2__==0||pos2__>right__) {s__ = data__[pos1__,right__-pos1__];}
      else {s__ = data__[pos1__,pos2__-pos1__];}
      mat__[i__][j__] = list();

      factors__ = intvec(0:0);
      exponents__ = intvec(0:0);
      l__ = find(s__," * ");
      if(l__>0)
      {
        ss__ = s__[1,l__-1];    //ss__ contains the nonlinear factors
        s__ = s__[l__+3,size(s__)-l__-2];
        for(p1__=1;s__[p1__]!="(";p1__++) {}
        for(p2__=size(s__);s__[p2__]!=")";p2__--) {}
        s__ = s__[p1__+1,p2__-p1__-1];
        l__ = find(ss__,"q");
        ss__ = ss__[l__,find(ss__,")",l__)-l__];
        ss__ = ss__+"*";
        p1__=0;
        for(l__=1;1;l__++)
        {
          p1__ = find(ss__,"q",p1__+1);
          if(p1__==0) {break;}
          p1__++;
          p2__ = find(ss__,"^",p1__);
          tmp__ = find(ss__,"*",p1__);
          if((p2__>tmp__ && tmp__>0) || (p2__==0))    //exponent is 1
          {
            execute("factors__[l__]="+ss__[p1__,tmp__-p1__]);
            exponents__[l__] = 1;
          }
          else
          {
            execute("factors__[l__]="+ss__[p1__,p2__-p1__]);
            execute("exponents__[l__]="+ss__[p2__+1,tmp__-p2__-1]);
          }
        }
      }
      nonlin__[i__][j__] = list(factors__,exponents__);

      depth__ = 0;
      s__=s__+" ";
      max__ = size(s__);
      tmp__ = 1;
      tmp2__ = 0;
      for(k__=1;k__<=max__;k__++)
      {
        if(s__[k__]=="(") {depth__++;k__++;continue;}
        if(s__[k__]==")") {depth__--;k__++;continue;}
        if(s__[k__]=="/" && depth__==0) {tmp2__ = k__;}
        if((s__[k__]=="+" && depth__==0) || k__==max__)
        {
          if(tmp2__==0)    // no denominator
          {
            execute("numerator__="+s__[tmp__,k__-tmp__]);
            mat__[i__][j__][size(mat__[i__][j__])+1] = list(numerator__,intvec(0:0),intvec(0:0));
            tmp__ = k__+1;
            k__++; continue;
          }
          execute("numerator__="+s__[tmp__,tmp2__-tmp__]);
          ss__ = s__[tmp2__+1,k__-tmp2__-1];
          p1__ = find(ss__,"(");
          p2__ = find(ss__,")");
          ss__ = ss__[p1__+1,p2__-p1__-1];   // now ss__ is only the denominator

          ss__ = ss__+"*";
          factors__ = intvec(0:0);
          exponents__ = intvec(0:0);
          p1__=0;
          for(l__=1;1;l__++)
          {
            p1__ = find(ss__,"q",p1__+1);
            if(p1__==0) {break;}
            p1__++;
            p2__ = find(ss__,"^",p1__);
            tmp__ = find(ss__,"*",p1__);
            if((p2__>tmp__ && tmp__>0) || (p2__==0))    //exponent is 1
            {
              execute("factors__[l__]="+ss__[p1__,tmp__-p1__]);
              exponents__[l__] = 1;
            }
            else
            {
              execute("factors__[l__]="+ss__[p1__,p2__-p1__]);
              execute("exponents__[l__]="+ss__[p2__+1,tmp__-p2__-1]);
            }
          }
          mat__[i__][j__][size(mat__[i__][j__])+1] = list(numerator__,factors__,exponents__);
          tmp__ = k__+1;
          tmp2__ = 0;
        }
      }
    }
    dbprint(sprintf("    row %s done! (%s ms)",i__,rtimer-tt__));
  }
  dbprint(sprintf("  done! (%s ms)", rtimer-t__));

  return(mat__,nonlin__);
}

static proc readQfileTXT(string filename__)
{
  string data__ = read(":r "+filename__);
  ideal q__;
  int pos1__,pos2__=1,1;
  while(1)
  {
    pos1__=find(data__,"=",pos2__);
    if(pos1__==0) {break;}
    pos1__++;
    pos2__=find(data__,";",pos1__);
    execute("q__[size(q__)+1]="+data__[pos1__,pos2__-pos1__]);
  }
  return(q__);
}