##############################################################################
### INPUT: a trivariate rational function f in s,t,x.
### OUTPUT: a linear differential operator L(x,Dx) annihilating Diag(f).
### NOTES on the algorithm: it is based on Remark 3 in [Lipshitz89]: let
### F(s,t,x) = 1/s/t * f(s,t/s,x/t). If P(x;Dx,Ds,Dt) annihilates F, then
### P(x;Dx,0,0) kills Diag(f).
### The idea is to determine such a P by undetermined coefficients: what makes
### things work is the observation that P(F)=0 amounts to a linear homogeneous
### system in the undetermined coefficients of P, having nb_eqs = O(N^3)
### equations and nb_unk = O(N^4) unknowns.
### In our approach (which is more in Fasenmyer's style), the notable
### difference is that the equation P(F)=0 is also equivalent  to a
### homogeneous linear system in the (undetermined) coefficients c_{ijk}(x) of
### P = sum(c_{ijk}*Dx^i*Ds^j*Dt^k). The system has nb_eqs = O(N^2) equations
### and nb_unk = O(N^3) unknowns, these unknowns being polynomials in x,
### instead of scalars.
### The algorithm works by increasing the value of N until nb_unk > nb_eqs.
### This value of N is typically smaller than the one in Lipshitz's approach.
### On our particular example of 3D Rooks, we have:
### nb_eqs:=(13*N+14)*(N+1)/2
### nb_unk:=binomial(N+3,3);
##############################################################################

Fasenmyer3rat:=proc(f, delay := 0, $)
local F, Q, L, N, NB_eqs, NB_unk, del, NumDer, j, i, k, eqs, sys,
unk, nb_eqs, nb_unk, x0, sol0, sol, Eq, Lstu, pol_eqs, pol_unk, tmp, Qx, Qs, Qt, zob;
  F:=normal(1/s/t*subs({x=s,s=t/s,t=x/t},f));
  Q:=denom(F);
  Qx := diff(Q, x) ;
  Qs := diff(Q, s) ;
  Qt := diff(Q, t) ;
  L:=0;
  NB_eqs:=[]:
  NB_unk:=NB_eqs:
  del := delay ;

  NumDer[0, 0, 0] := numer(F) ;

  # Invariant: NumDer[i, j, k] / Q^(N+1) = diff(F, [x$i, s$j, t$k]).
  for N while L = 0 or del > 0 do
    if L <> 0 then
      del := del - 1
    end if ;
  print(N);
  print("setting up the linear system");
  print("derivatives");

  tmp := NumDer[0, N-1, 0] ;
  NumDer[0, N, 0] := expand( diff(tmp, s) * Q - N * tmp * Qs ) ;
  tmp := NumDer[0, 0, N-1] ;
  NumDer[0, 0, N] := expand( diff(tmp, t) * Q - N * tmp * Qt ) ;

  for i to N do
    tmp := NumDer[0, i-1, N-i] ;
    NumDer[0, i, N-i]:=expand( diff(tmp, s) * Q - N * tmp * Qs )
  od ;

  for i from 1 to N do
    for j from 0 to N-i do
      tmp := NumDer[i-1, j, N-i-j] ;
      NumDer[i, j, N-i-j] := expand( diff(tmp, x) * Q - N * tmp * Qx )
    od
  od ;

  NumDer[0, 0, 0] := normal(NumDer[0, 0, 0] * Q) ;
  for j from 0 to N-2 do
    NumDer[0, j+1, 0] := normal(NumDer[0, j+1, 0] * Q) ;
    NumDer[0, 0, j+1] := normal(NumDer[0, 0, j+1] * Q) ;
  od ;
  for i to N do
    for j to N-i-1 do
      NumDer[0, i, j] := normal(NumDer[0, i, j] * Q)
    od
  od ;
  for i from 1 to N do
    for j from 0 to N-i do
      for k from 0 to N-i-j-1 do
        NumDer[i, j, k] := normal(NumDer[i, j, k] * Q)
      od
    od
  od ;

  if L <> 0 and del > 0 then
    next
  end if ;

  print("add equations");

  eqs:=normal(add(add(add(c[i, j, k]*NumDer[i, j, k],
    k=0..N-j-i),j=0..N-i),i=0..N)):

  print("extract equations");

  eqs := expand(eqs) ;
  sys := [coeffs(eqs, {s, t})] ;

  unk := subs(x = NULL, indets(sys)) ;

	nb_eqs:=nops(sys); NB_eqs:=[op(NB_eqs),nb_eqs]:
	nb_unk:=nops(unk); NB_unk:=[op(NB_unk),nb_unk]:

  print("number of equations and unkowns", nb_eqs, nb_unk);
  print("degree in x", max(seq(degree(eq,x), eq=eqs)));

    print("solving the linear system");
    x0:=nextprime(rand(2^10)(1));
    print("at a random point x =", x0);
    sol0 := SolveTools:-Linear(eval(sys, x = x0), unk) ;

    if map2(op, 2, sol0) <> {0} and L * del = 0 then
      print("it seems that there is a solution");
      sol := SolveTools:-Linear(sys, unk) ;

      print("computing annihilating operator");
      Eq:=0;

      for i from 0 to N do
        for j from 0 to N-i do
          for k from 0 to N-i-j do
            Eq:= normal(Eq+c[i, j, k]*Dx^i*Ds^j*Dt^k);
          od
        od
      od;

      Lstu:=subs(sol,numer(Eq));
      if Lstu<>0 then
        L:=subs({Ds=0,Dt=0},Lstu);
      end if;

    end if;

  if N = 6 then
  pol_eqs:=factor(CurveFitting[PolynomialInterpolation]([$1..nops(NB_eqs)], NB_eqs,Z));
  pol_unk:=factor(CurveFitting[PolynomialInterpolation]([$1..nops(NB_unk)], NB_unk,Z));
  print("number of equations",pol_eqs);
  print("number of unknowns", pol_unk);
  end if

od ;

  return collect(numer(L),Dx), numer(Lstu), pol_eqs, pol_unk;
end proc;

####

MinimalFasenmyer3rat := proc(f, delay := 0, $)
local tt, FR, L, ind, L1, L2, G;
  tt:=time(): FR:=Fasenmyer3rat(f, delay); print("time Fasenmyer3rat", time()-tt);
  L:=FR[1];
  ind := indets(L) minus {x, Dx} ;
  L1:=subs([seq(s=rand(100)(1),s=ind)],L):
  L2:=subs([seq(s=rand(100)(1),s=ind)],L):
  tt:=time(): G:=DEtools[GCRD](L1,L2,[Dx,x]); print("time GCRD", time()-tt);
return G;
end proc;


#### simplest example
#f:=1/(1-x-s-t):
#MinimalFasenmyer3rat(f);

#### 3D rook
f:=1/(1 - s/(1-s) - t/(1-t) - x/(1-x) );
MinimalFasenmyer3rat(f);