//************************************************************************
//
//  Purpose:
//    Compute elliptic division polynomials via SLPs.
//    Compute also the partial derivatives wrt X, A, B, also via SLPs.
//  Author:
//    Pierrick Gaudry
//
//  The formulae are \overline{f_m} in Blake-Seroussi-Smart, page 40-41
//
//***********************************************************************
//
//  Prototypes:
//
//    SLPDivisionPolynomial(A, B, X, n);
//    SLPDiffXDivisionPolynomial(A, B, X, n);
//    SLPDiffADivisionPolynomial(A, B, X, n);
//    SLPDiffBDivisionPolynomial(A, B, X, n);
//
//***********************************************************************



procedure SLPDivisionPolynomialRemember(A, B, X, n, ~TablePol, ~TableIndices)
  ind := Position(TableIndices, n);
  if ind ne 0 then
    return ;
  end if;
  
  if IsOdd(n) then
    m := n div 2;  // n = 2*m + 1
    SLPDivisionPolynomialRemember(A, B, X, m+2, ~TablePol, ~TableIndices);
    SLPDivisionPolynomialRemember(A, B, X, m, ~TablePol, ~TableIndices);
    SLPDivisionPolynomialRemember(A, B, X, m-1, ~TablePol, ~TableIndices);
    SLPDivisionPolynomialRemember(A, B, X, m+1, ~TablePol, ~TableIndices);
    ind := Position(TableIndices, m+2); Pmp2 := TablePol[ind];
    ind := Position(TableIndices, m); Pm := TablePol[ind];
    ind := Position(TableIndices, m-1); Pmm1 := TablePol[ind];
    ind := Position(TableIndices, m+1); Pmp1 := TablePol[ind];
    if IsOdd(m) then
      assert m ge 3;
      Psi := Pmp2*Pm^3 - 16*(X^3 + A*X + B)^2*Pmm1*Pmp1^3;
    else // m is even
      assert m ge 2;  
      Psi := 16*(X^3 + A*X + B)^2*Pmp2*Pm^3 - Pmm1*Pmp1^3;
    end if;
  else // n is even
    m := n div 2;  // n = 2*m
    assert m gt 2;
    SLPDivisionPolynomialRemember(A, B, X, m+2, ~TablePol, ~TableIndices);
    SLPDivisionPolynomialRemember(A, B, X, m-2, ~TablePol, ~TableIndices);
    SLPDivisionPolynomialRemember(A, B, X, m, ~TablePol, ~TableIndices);
    SLPDivisionPolynomialRemember(A, B, X, m-1, ~TablePol, ~TableIndices);
    SLPDivisionPolynomialRemember(A, B, X, m+1, ~TablePol, ~TableIndices);
    ind := Position(TableIndices, m+2); Pmp2 := TablePol[ind];
    ind := Position(TableIndices, m-2); Pmm2 := TablePol[ind];
    ind := Position(TableIndices, m); Pm := TablePol[ind];
    ind := Position(TableIndices, m-1); Pmm1 := TablePol[ind];
    ind := Position(TableIndices, m+1); Pmp1 := TablePol[ind];
    
    Psi := Pm*(Pmp2*Pmm1^2 - Pmm2*Pmp1^2);
  end if;
  Append(~TablePol, Psi);
  Append(~TableIndices, n);
end procedure;

function SLPDivisionPolynomial(A, B, X, n)
  myring := Parent(X+A - (X+A));
  
  TableIndices := [0,1,2,3,4];
  TablePol := [ myring | 0, 1, 1, 3*X^4 + 6*A*X^2 + 12*B*X - A^2,
    2*X^6 + 10*A*X^4 + 40*B*X^3 - 10*A^2*X^2 -8*A*B*X - 16*B^2 - 2*A^3 ];
  
  ind := Position(TableIndices, n);
  if ind ne 0 then 
    return TablePol[ind];
  end if;

  SLPDivisionPolynomialRemember(A, B, X, n, ~TablePol, ~TableIndices);
  ind := Position(TableIndices, n);
  assert ind ne 0;
  return TablePol[ind];
end function;


procedure SLPDiffXDivisionPolynomialRemember(A, B, X, n, ~TablePol,
  ~TableIndices, ~TableDiffXPol, ~TableDiffXIndices)
  ind := Position(TableDiffXIndices, n);
  if ind ne 0 then
    return ;
  end if;
 
  if IsOdd(n) then
    m := n div 2;  // n = 2*m + 1
    SLPDivisionPolynomialRemember(A, B, X, m+2, ~TablePol, ~TableIndices);
    SLPDivisionPolynomialRemember(A, B, X, m, ~TablePol, ~TableIndices);
    SLPDivisionPolynomialRemember(A, B, X, m-1, ~TablePol, ~TableIndices);
    SLPDivisionPolynomialRemember(A, B, X, m+1, ~TablePol, ~TableIndices);
    ind := Position(TableIndices, m+2); Pmp2 := TablePol[ind];
    ind := Position(TableIndices, m); Pm := TablePol[ind];
    ind := Position(TableIndices, m-1); Pmm1 := TablePol[ind];
    ind := Position(TableIndices, m+1); Pmp1 := TablePol[ind];
    SLPDiffXDivisionPolynomialRemember(A, B, X, m+2, ~TablePol, ~TableIndices,
      ~TableDiffXPol, ~TableDiffXIndices);
    SLPDiffXDivisionPolynomialRemember(A, B, X, m, ~TablePol, ~TableIndices,
      ~TableDiffXPol, ~TableDiffXIndices);
    SLPDiffXDivisionPolynomialRemember(A, B, X, m-1, ~TablePol, ~TableIndices,
      ~TableDiffXPol, ~TableDiffXIndices);
    SLPDiffXDivisionPolynomialRemember(A, B, X, m+1, ~TablePol, ~TableIndices,
      ~TableDiffXPol, ~TableDiffXIndices);
    ind := Position(TableDiffXIndices, m+2); DPmp2 := TableDiffXPol[ind];
    ind := Position(TableDiffXIndices, m); DPm := TableDiffXPol[ind];
    ind := Position(TableDiffXIndices, m-1); DPmm1 := TableDiffXPol[ind];
    ind := Position(TableDiffXIndices, m+1); DPmp1 := TableDiffXPol[ind];
    if IsOdd(m) then
      assert m ge 3;
      DPsi := DPmp2*Pm^3 + 3*Pmp2*DPm*Pm^2
	- 2*16*(3*X^2 + A)*(X^3 + A*X + B)*Pmm1*Pmp1^3
	- 16*(X^3 + A*X + B)^2*(DPmm1*Pmp1^3 + 3*Pmm1*DPmp1*Pmp1^2);
    else // m is even
      assert m ge 2;
      DPsi := 2*16*(3*X^2 + A)*(X^3 + A*X + B)*Pmp2*Pm^3 +
	16*(X^3 + A*X + B)^2*(DPmp2*Pm^3 + 3*Pmp2*DPm*Pm^2)
	- (DPmm1*Pmp1^3 + 3*Pmm1*DPmp1*Pmp1^2);
    end if;
  else // n is even
    m := n div 2;  // n = 2*m
    assert m gt 2;
    SLPDivisionPolynomialRemember(A, B, X, m+2, ~TablePol, ~TableIndices);
    SLPDivisionPolynomialRemember(A, B, X, m-2, ~TablePol, ~TableIndices);
    SLPDivisionPolynomialRemember(A, B, X, m, ~TablePol, ~TableIndices);
    SLPDivisionPolynomialRemember(A, B, X, m-1, ~TablePol, ~TableIndices);
    SLPDivisionPolynomialRemember(A, B, X, m+1, ~TablePol, ~TableIndices);
    ind := Position(TableIndices, m+2); Pmp2 := TablePol[ind];
    ind := Position(TableIndices, m-2); Pmm2 := TablePol[ind];
    ind := Position(TableIndices, m); Pm := TablePol[ind];
    ind := Position(TableIndices, m-1); Pmm1 := TablePol[ind];
    ind := Position(TableIndices, m+1); Pmp1 := TablePol[ind];
    SLPDiffXDivisionPolynomialRemember(A, B, X, m+2, ~TablePol, ~TableIndices,
      ~TableDiffXPol, ~TableDiffXIndices);
    SLPDiffXDivisionPolynomialRemember(A, B, X, m-2, ~TablePol, ~TableIndices,
      ~TableDiffXPol, ~TableDiffXIndices);
    SLPDiffXDivisionPolynomialRemember(A, B, X, m, ~TablePol, ~TableIndices,
      ~TableDiffXPol, ~TableDiffXIndices);
    SLPDiffXDivisionPolynomialRemember(A, B, X, m-1, ~TablePol, ~TableIndices,
      ~TableDiffXPol, ~TableDiffXIndices);
    SLPDiffXDivisionPolynomialRemember(A, B, X, m+1, ~TablePol, ~TableIndices,
      ~TableDiffXPol, ~TableDiffXIndices);
    ind := Position(TableDiffXIndices, m+2); DPmp2 := TableDiffXPol[ind];
    ind := Position(TableDiffXIndices, m-2); DPmm2 := TableDiffXPol[ind];
    ind := Position(TableDiffXIndices, m); DPm := TableDiffXPol[ind];
    ind := Position(TableDiffXIndices, m-1); DPmm1 := TableDiffXPol[ind];
    ind := Position(TableDiffXIndices, m+1); DPmp1 := TableDiffXPol[ind];

    DPsi := DPm*(Pmp2*Pmm1^2 - Pmm2*Pmp1^2)+
      Pm*(
	DPmp2*Pmm1^2 + 2*Pmp2*DPmm1*Pmm1
	  - (DPmm2*Pmp1^2 + 2*Pmm2*DPmp1*Pmp1)
	    );
  end if;
  Append(~TableDiffXPol, DPsi);
  Append(~TableDiffXIndices, n);
end procedure;


function SLPDiffXDivisionPolynomial(A, B, X, n)
  myring := Parent(X+A - (X+A));
  
  TableDiffXIndices := [0,1,2,3,4];
  TableDiffXPol := [ myring | 0, 0, 0, 12*X^3 + 12*A*X + 12*B,
    12*X^5 + 40*A*X^3 + 120*B*X^2 - 20*A^2*X -8*A*B ];
  
  TableIndices := [0,1,2,3,4];
  TablePol := [ myring | 0, 1, 1, 3*X^4 + 6*A*X^2 + 12*B*X - A^2,
    2*X^6 + 10*A*X^4 + 40*B*X^3 - 10*A^2*X^2 -8*A*B*X - 16*B^2 - 2*A^3 ];
    
  ind := Position(TableDiffXIndices, n);
  if ind ne 0 then 
    return TableDiffXPol[ind];
  end if;

  SLPDiffXDivisionPolynomialRemember(A, B, X, n, ~TablePol, ~TableIndices,
    ~TableDiffXPol, ~TableDiffXIndices);
  ind := Position(TableDiffXIndices, n);
  assert ind ne 0;
  return TableDiffXPol[ind];
end function;




procedure SLPDiffADivisionPolynomialRemember(A, B, X, n, ~TablePol,
  ~TableIndices, ~TableDiffAPol, ~TableDiffAIndices)
  ind := Position(TableDiffAIndices, n);
  if ind ne 0 then
    return ;
  end if;
 
  if IsOdd(n) then
    m := n div 2;  // n = 2*m + 1
    SLPDivisionPolynomialRemember(A, B, X, m+2, ~TablePol, ~TableIndices);
    SLPDivisionPolynomialRemember(A, B, X, m, ~TablePol, ~TableIndices);
    SLPDivisionPolynomialRemember(A, B, X, m-1, ~TablePol, ~TableIndices);
    SLPDivisionPolynomialRemember(A, B, X, m+1, ~TablePol, ~TableIndices);
    ind := Position(TableIndices, m+2); Pmp2 := TablePol[ind];
    ind := Position(TableIndices, m); Pm := TablePol[ind];
    ind := Position(TableIndices, m-1); Pmm1 := TablePol[ind];
    ind := Position(TableIndices, m+1); Pmp1 := TablePol[ind];
    SLPDiffADivisionPolynomialRemember(A, B, X, m+2, ~TablePol, ~TableIndices,
      ~TableDiffAPol, ~TableDiffAIndices);
    SLPDiffADivisionPolynomialRemember(A, B, X, m, ~TablePol, ~TableIndices,
      ~TableDiffAPol, ~TableDiffAIndices);
    SLPDiffADivisionPolynomialRemember(A, B, X, m-1, ~TablePol, ~TableIndices,
      ~TableDiffAPol, ~TableDiffAIndices);
    SLPDiffADivisionPolynomialRemember(A, B, X, m+1, ~TablePol, ~TableIndices,
      ~TableDiffAPol, ~TableDiffAIndices);
    ind := Position(TableDiffAIndices, m+2); DPmp2 := TableDiffAPol[ind];
    ind := Position(TableDiffAIndices, m); DPm := TableDiffAPol[ind];
    ind := Position(TableDiffAIndices, m-1); DPmm1 := TableDiffAPol[ind];
    ind := Position(TableDiffAIndices, m+1); DPmp1 := TableDiffAPol[ind];
    if IsOdd(m) then
      assert m ge 3;
      DPsi := DPmp2*Pm^3 + 3*Pmp2*DPm*Pm^2
	- 32*X*(X^3 + A*X + B)*Pmm1*Pmp1^3
	- 16*(X^3 + A*X + B)^2*(DPmm1*Pmp1^3 + 3*Pmm1*DPmp1*Pmp1^2);
    else // m is even
      assert m ge 2;
      DPsi := 32*X*(X^3 + A*X + B)*Pmp2*Pm^3 +
	16*(X^3 + A*X + B)^2*(DPmp2*Pm^3 + 3*Pmp2*DPm*Pm^2)
	- (DPmm1*Pmp1^3 + 3*Pmm1*DPmp1*Pmp1^2);
    end if;
  else // n is even
    m := n div 2;  // n = 2*m
    assert m gt 2;
    SLPDivisionPolynomialRemember(A, B, X, m+2, ~TablePol, ~TableIndices);
    SLPDivisionPolynomialRemember(A, B, X, m-2, ~TablePol, ~TableIndices);
    SLPDivisionPolynomialRemember(A, B, X, m, ~TablePol, ~TableIndices);
    SLPDivisionPolynomialRemember(A, B, X, m-1, ~TablePol, ~TableIndices);
    SLPDivisionPolynomialRemember(A, B, X, m+1, ~TablePol, ~TableIndices);
    ind := Position(TableIndices, m+2); Pmp2 := TablePol[ind];
    ind := Position(TableIndices, m-2); Pmm2 := TablePol[ind];
    ind := Position(TableIndices, m); Pm := TablePol[ind];
    ind := Position(TableIndices, m-1); Pmm1 := TablePol[ind];
    ind := Position(TableIndices, m+1); Pmp1 := TablePol[ind];
    SLPDiffADivisionPolynomialRemember(A, B, X, m+2, ~TablePol, ~TableIndices,
      ~TableDiffAPol, ~TableDiffAIndices);
    SLPDiffADivisionPolynomialRemember(A, B, X, m-2, ~TablePol, ~TableIndices,
      ~TableDiffAPol, ~TableDiffAIndices);
    SLPDiffADivisionPolynomialRemember(A, B, X, m, ~TablePol, ~TableIndices,
      ~TableDiffAPol, ~TableDiffAIndices);
    SLPDiffADivisionPolynomialRemember(A, B, X, m-1, ~TablePol, ~TableIndices,
      ~TableDiffAPol, ~TableDiffAIndices);
    SLPDiffADivisionPolynomialRemember(A, B, X, m+1, ~TablePol, ~TableIndices,
      ~TableDiffAPol, ~TableDiffAIndices);
    ind := Position(TableDiffAIndices, m+2); DPmp2 := TableDiffAPol[ind];
    ind := Position(TableDiffAIndices, m-2); DPmm2 := TableDiffAPol[ind];
    ind := Position(TableDiffAIndices, m); DPm := TableDiffAPol[ind];
    ind := Position(TableDiffAIndices, m-1); DPmm1 := TableDiffAPol[ind];
    ind := Position(TableDiffAIndices, m+1); DPmp1 := TableDiffAPol[ind];

    DPsi := DPm*(Pmp2*Pmm1^2 - Pmm2*Pmp1^2)+
      Pm*(
	DPmp2*Pmm1^2 + 2*Pmp2*DPmm1*Pmm1
	  - (DPmm2*Pmp1^2 + 2*Pmm2*DPmp1*Pmp1)
	    );
  end if;
  Append(~TableDiffAPol, DPsi);
  Append(~TableDiffAIndices, n);
end procedure;


function SLPDiffADivisionPolynomial(A, B, X, n)
  myring := Parent(X+A - (X+A));
  
  TableDiffAIndices := [0,1,2,3,4];
  TableDiffAPol := [ myring | 0, 0, 0, 6*X^2 - 2*A, 
    10*X^4 - 20*A*X^2 -8*B*X - 6*A^2 ];
  
  TableIndices := [0,1,2,3,4];
  TablePol := [ myring | 0, 1, 1, 3*X^4 + 6*A*X^2 + 12*B*X - A^2,
    2*X^6 + 10*A*X^4 + 40*B*X^3 - 10*A^2*X^2 -8*A*B*X - 16*B^2 - 2*A^3 ];
    
  ind := Position(TableDiffAIndices, n);
  if ind ne 0 then 
    return TableDiffAPol[ind];
  end if;

  SLPDiffADivisionPolynomialRemember(A, B, X, n, ~TablePol, ~TableIndices,
    ~TableDiffAPol, ~TableDiffAIndices);
  ind := Position(TableDiffAIndices, n);
  assert ind ne 0;
  return TableDiffAPol[ind];
end function;


procedure SLPDiffBDivisionPolynomialRemember(A, B, X, n, ~TablePol,
  ~TableIndices, ~TableDiffBPol, ~TableDiffBIndices)
  ind := Position(TableDiffBIndices, n);
  if ind ne 0 then
    return ;
  end if;
 
  if IsOdd(n) then
    m := n div 2;  // n = 2*m + 1
    SLPDivisionPolynomialRemember(A, B, X, m+2, ~TablePol, ~TableIndices);
    SLPDivisionPolynomialRemember(A, B, X, m, ~TablePol, ~TableIndices);
    SLPDivisionPolynomialRemember(A, B, X, m-1, ~TablePol, ~TableIndices);
    SLPDivisionPolynomialRemember(A, B, X, m+1, ~TablePol, ~TableIndices);
    ind := Position(TableIndices, m+2); Pmp2 := TablePol[ind];
    ind := Position(TableIndices, m); Pm := TablePol[ind];
    ind := Position(TableIndices, m-1); Pmm1 := TablePol[ind];
    ind := Position(TableIndices, m+1); Pmp1 := TablePol[ind];
    SLPDiffBDivisionPolynomialRemember(A, B, X, m+2, ~TablePol, ~TableIndices,
      ~TableDiffBPol, ~TableDiffBIndices);
    SLPDiffBDivisionPolynomialRemember(A, B, X, m, ~TablePol, ~TableIndices,
      ~TableDiffBPol, ~TableDiffBIndices);
    SLPDiffBDivisionPolynomialRemember(A, B, X, m-1, ~TablePol, ~TableIndices,
      ~TableDiffBPol, ~TableDiffBIndices);
    SLPDiffBDivisionPolynomialRemember(A, B, X, m+1, ~TablePol, ~TableIndices,
      ~TableDiffBPol, ~TableDiffBIndices);
    ind := Position(TableDiffBIndices, m+2); DPmp2 := TableDiffBPol[ind];
    ind := Position(TableDiffBIndices, m); DPm := TableDiffBPol[ind];
    ind := Position(TableDiffBIndices, m-1); DPmm1 := TableDiffBPol[ind];
    ind := Position(TableDiffBIndices, m+1); DPmp1 := TableDiffBPol[ind];
    if IsOdd(m) then
      assert m ge 3;
      DPsi := DPmp2*Pm^3 + 3*Pmp2*DPm*Pm^2
	- 32*(X^3 + A*X + B)*Pmm1*Pmp1^3
	- 16*(X^3 + A*X + B)^2*(DPmm1*Pmp1^3 + 3*Pmm1*DPmp1*Pmp1^2);
    else // m is even
      assert m ge 2;
      DPsi := 32*(X^3 + A*X + B)*Pmp2*Pm^3 +
	16*(X^3 + A*X + B)^2*(DPmp2*Pm^3 + 3*Pmp2*DPm*Pm^2)
	- (DPmm1*Pmp1^3 + 3*Pmm1*DPmp1*Pmp1^2);
    end if;
  else // n is even
    m := n div 2;  // n = 2*m
    assert m gt 2;
    SLPDivisionPolynomialRemember(A, B, X, m+2, ~TablePol, ~TableIndices);
    SLPDivisionPolynomialRemember(A, B, X, m-2, ~TablePol, ~TableIndices);
    SLPDivisionPolynomialRemember(A, B, X, m, ~TablePol, ~TableIndices);
    SLPDivisionPolynomialRemember(A, B, X, m-1, ~TablePol, ~TableIndices);
    SLPDivisionPolynomialRemember(A, B, X, m+1, ~TablePol, ~TableIndices);
    ind := Position(TableIndices, m+2); Pmp2 := TablePol[ind];
    ind := Position(TableIndices, m-2); Pmm2 := TablePol[ind];
    ind := Position(TableIndices, m); Pm := TablePol[ind];
    ind := Position(TableIndices, m-1); Pmm1 := TablePol[ind];
    ind := Position(TableIndices, m+1); Pmp1 := TablePol[ind];
    SLPDiffBDivisionPolynomialRemember(A, B, X, m+2, ~TablePol, ~TableIndices,
      ~TableDiffBPol, ~TableDiffBIndices);
    SLPDiffBDivisionPolynomialRemember(A, B, X, m-2, ~TablePol, ~TableIndices,
      ~TableDiffBPol, ~TableDiffBIndices);
    SLPDiffBDivisionPolynomialRemember(A, B, X, m, ~TablePol, ~TableIndices,
      ~TableDiffBPol, ~TableDiffBIndices);
    SLPDiffBDivisionPolynomialRemember(A, B, X, m-1, ~TablePol, ~TableIndices,
      ~TableDiffBPol, ~TableDiffBIndices);
    SLPDiffBDivisionPolynomialRemember(A, B, X, m+1, ~TablePol, ~TableIndices,
      ~TableDiffBPol, ~TableDiffBIndices);
    ind := Position(TableDiffBIndices, m+2); DPmp2 := TableDiffBPol[ind];
    ind := Position(TableDiffBIndices, m-2); DPmm2 := TableDiffBPol[ind];
    ind := Position(TableDiffBIndices, m); DPm := TableDiffBPol[ind];
    ind := Position(TableDiffBIndices, m-1); DPmm1 := TableDiffBPol[ind];
    ind := Position(TableDiffBIndices, m+1); DPmp1 := TableDiffBPol[ind];

    DPsi := DPm*(Pmp2*Pmm1^2 - Pmm2*Pmp1^2)+
      Pm*(
	DPmp2*Pmm1^2 + 2*Pmp2*DPmm1*Pmm1
	  - (DPmm2*Pmp1^2 + 2*Pmm2*DPmp1*Pmp1)
	    );
  end if;
  Append(~TableDiffBPol, DPsi);
  Append(~TableDiffBIndices, n);
end procedure;


function SLPDiffBDivisionPolynomial(A, B, X, n)
  myring := Parent(X+A - (X+A));
  
  TableDiffBIndices := [0,1,2,3,4];
  TableDiffBPol := [ myring | 0, 0, 0, 12*X, 
    40*X^3 - 8*A*X - 32*B ];
    
  TableIndices := [0,1,2,3,4];
  TablePol := [ myring | 0, 1, 1, 3*X^4 + 6*A*X^2 + 12*B*X - A^2,
    2*X^6 + 10*A*X^4 + 40*B*X^3 - 10*A^2*X^2 -8*A*B*X - 16*B^2 - 2*A^3 ];
    
  ind := Position(TableDiffBIndices, n);
  if ind ne 0 then 
    return TableDiffBPol[ind];
  end if;

  SLPDiffBDivisionPolynomialRemember(A, B, X, n, ~TablePol, ~TableIndices,
    ~TableDiffBPol, ~TableDiffBIndices);
  ind := Position(TableDiffBIndices, n);
  assert ind ne 0;
  return TableDiffBPol[ind];
end function;