\\ Elliptic division polynomials
\\  magma Author:
\\    Pierrick Gaudry
\\ ported to gp/pari by unknown
\\ usage:
\\ gp
\\ \r DP-pub.gp
{
Position(li,el)=
local(i,l);
for(i=1,length(li),if (el == li[i],return(i)));
return(0);
}

{
SLPDivisionPolynomialRemember(A, B, X, n, XXXTablePol, XXXTableIndices)=
local(ind,Psi,m,Pmm1,Pmm2,Pmp1,Pmp2,Pm);
  ind = Position(TableIndices, n);
\\print("remember ",ind,TableIndices);
  if(ind != 0 , return(X));
  
  if( (n%2)!= 0,
    m = floor(n / 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((m%2!=0) , 
      \\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;
    )
  , \\ else n is even
    m = floor(n / 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);
  );
  listput(TablePol, Psi);
  listput(TableIndices, n);
}

{
SLPDivisionPolynomial(A, B, X, n)=
local(ind);
\\  myring = Parent(X+A - (X+A));
  
\\  TableIndices = [0,1,2,3,4];
if(type(TablePol)=="t_LIST",listkill(TablePol));
if(type(TableIndices)=="t_LIST",listkill(TableIndices));

TableIndices=listcreate(MAXLI);
TablePol=listcreate(MAXLI);

	listput(TableIndices,0);

	listput(TableIndices,0);
	listput(TableIndices,1);
	listput(TableIndices,2);
	listput(TableIndices,3);
	listput(TableIndices,4);

	listput(TablePol,0);

	listput(TablePol,0);
	listput(TablePol,1);
	listput(TablePol,1);
	listput(TablePol,3*X^4 + 6*A*X^2 + 12*B*X - A^2);
	listput(TablePol,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 );
\\  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 != 0 , return(TablePol[ind]));

  SLPDivisionPolynomialRemember(A, B, X, n, TablePol, TableIndices);
  ind = Position(TableIndices, n);
  \\assert ind != 0;
  return(TablePol[ind]);
}

\\ n-th chebyshev polynomial of 1st kind
\\ author: www.jjj.de
{
cheby1f(n,x)=
return( ([1,x]*[0, -1;1,2*x]^n)[1]);
}


MAXLI=20000; \\XXX lame
\\global(TableIndices,TablePol);
TableIndices=listcreate(MAXLI);
TablePol=listcreate(MAXLI);
po1=SLPDivisionPolynomial(1,42,x,17)
tt1=SLPDivisionPolynomial(1,42,Mod(2,113),1005)