// Fasenmyer.
F<n,k>:=FunctionField(RationalField(),2);

DecalageN:=hom<F->F|n+1,k>;
DecalageK:=hom<F->F|n,k+1>;
DecalageInverseK:=hom<F->F|n,k-1>;

function Facteur(Rn,Rk,i,j)
  if j gt 0 then return Rk*DecalageK($$(Rn,Rk,i,j-1)); end if;
  if i gt 0 then return Rn*DecalageN($$(Rn,1,i-1,0)); end if;
  return 1;
end function;

function Fasenmyer(Rn,Rk,S)
  C:=[F!0:i in [1..#S]];
  commonDenominator:=Denominator(F!1);
  for i:=1 to #S do
    C[i]:=Facteur(Rn,Rk,S[i][1],S[i][2]);
    commonDenominator:=Lcm(commonDenominator,Denominator(C[i]));
  end for;
  for i:=1 to #S do
    C[i]*:=commonDenominator;
  end for;
  C:=[Numerator(C[i]):i in [1..#C]];
  commonDivisor:=C[1];
  for i:=2 to #S do
    commonDivisor:=Gcd(commonDivisor,C[i]);
  end for;
  for i:=1 to #S do
    C[i]/:=commonDivisor;
  end for;

/*
  C:=[Facteur(Rn,Rk,s[1],s[2]):s in S];
  commonDenominator:=Lcm([Denominator(c):c in C]);
  C:=[Numerator(commonDenominator*c):c in C];
  commonDivisor:=Gcd(C);
  C:=[Numerator(c/commonDenominator):c in C];
*/

/*
  maximumDegree:=Maximum([Degree(c,Numerator(k)):c in C]);
  order:=Maximum([s[1]:s in S]);
*/

  maximumDegree:=0;
  order:=0;
  for i:=1 to #S do
    maximumDegree:=Maximum(maximumDegree,Degree(C[i],Numerator(k)));
    order:=Maximum(order,S[i][1]);
  end for;

  C:=Matrix(#S,maximumDegree+1,
    [Coefficient(C[i],Numerator(k),j):j in [0..maximumDegree],i in [1..#S]]);
  N:=NullSpace(C);
  B:=Basis(N);
  R:=Matrix(#B,order+1,[F!0:i in [1..(order+1)*#B]]);
  for i:=1 to #B do
    for j:=1 to #S do
      R[i,order+1-S[j][1]]+:=B[i][j];
    end for;
  end for;
  R:=EchelonForm(R);
  return R;
end function;

// Binome.
Fasenmyer((n+1)/(n+1-k), (n-k)/(k+1), [<i,j>:i in [0..1],j in [0..1]]);

// Dixon.
Fasenmyer(
  ((n+1)/(n+1-k))^3, -((n-k)/(k+1))^3, [<i,j>:i in [0..3],j in [0..3]]
);


// Bof.
// Gosper d'origine.
function FormeDePetkovsek(Rk)
  num:=Numerator(Rk);
  den:=Denominator(Rk);
  p:=Gcd(Numerator(DecalageInverseK(num)),den);
  return Rk*p/DecalageK(p);
end function;

FormeDePetkovsek(-(4*k+5)*(2*k-1)/2/(4*k+1)/(k+2));


// Abramov.
// Prend en entree [a_0,...,a_d].
function Abramov(listeDesCoeff)
  
end function;