reduce-historical/r38/packages/algint/nagell.red

module nagell;

% Author: James H. Davenport.

fluid '(!*tra !*trmin intvar);

exports lutz!-nagell;

symbolic procedure lutz!-nagell(divisor);
begin
  scalar ans,places,mults,save!*tra;
  for each u in divisor do <<
    places:=(car u).places;
    mults :=(cdr u).mults >>;
  ans:=lutz!-nagell!-2(places,mults);
  if ans eq 'infinite
     then return 'provably!-impossible;
  save!*tra:=!*tra;
  if !*trmin
    then !*tra:=nil;
  ans:=coates!-multiple(places,mults,ans);
  !*tra:=save!*tra;
  return ans
  end;


symbolic procedure lutz!-nagell!-2(places,mults);
begin
  scalar wst,x,y,equation,point,a;
  wst:=weierstrass!-form getsqrtsfromplaces places;
  x:=car wst;
  y:=cadr wst;
  equation:=caddr wst;
  equation:=!*q2f !*multsq(equation,equation);
  equation:=makemainvar(equation,intvar);
  if ldeg equation = 3
    then equation:=red equation
    else interr "Equation not of correct form";
  if mvar equation eq intvar
    then if ldeg equation = 1
      then <<
        a:=(lc equation) ./ 1;
        equation:=red equation >>
      else interr "Equation should not have a x**2 term"
    else a:=nil ./ 1;
  equation:= a . (equation ./ 1);
  places:=for each u in places collect
        wst!-convert(u,x,y);
  point:=elliptic!-sum(places,mults,equation);
  a:=lutz!-nagell!-bound(point,equation);
  if !*tra or !*trmin then <<
    princ "Point actually is of order ";
    printc a >>;
  return a
  end;


symbolic procedure wst!-convert(place,x,y);
begin
  x:=subzero(xsubstitutesq(x,place),intvar);
  y:=subzero(xsubstitutesq(y,place),intvar);
  return x.y
  end;


symbolic procedure elliptic!-sum(places,mults,equation);
begin
  scalar point;
  point:=elliptic!-multiply(car places,car mults,equation);
  places:=cdr places;
  mults:=cdr mults;
  while places do <<
    point:=elliptic!-add(point,
             elliptic!-multiply(car places,car mults,
                        equation),
                         equation);
    places:=cdr places;
    mults:=cdr mults >>;
  return point
  end;


symbolic procedure elliptic!-multiply(point,n,equation);
if n < 0
  then elliptic!-multiply( (car point) . (negsq cdr point),
                           -n,
                           equation)
  else if n = 0
    then interr "N=0 in elliptic!-multiply"
    else if n = 1
      then point
      else begin
        scalar q,r;
        q:=divide(n,2);
        r:=cdr q;
        q:=car q;
    q:=elliptic!-multiply(elliptic!-add(point,point,equation),q,
                        equation);
        if r = 0
          then return q
      else return elliptic!-add(point,q,equation)
        end;


symbolic procedure elliptic!-add(p1,p2,equation);
begin
  scalar x1,x2,y1,y2,x3,y3,inf,a,b,lhs,rhs;
  a:=car equation;
  b:=cdr equation;
  inf:=!*kk2q 'infinite;
  x1:=car p1;
  y1:=cdr p1;
  x2:=car p2;
  y2:=cdr p2;
  if x1 = x2
    then if y1 = y2
      then <<
    % this is the doubling case.
    x3:= multsq(4 ./ 1,
		!*addsq(b,!*multsq(x1,!*addsq(a, !*exptsq(x1,2)))));
    if null numr x3 then return inf . inf;
    % We doubled a point and got infinity. 
    x3:=!*multsq(!*addsq(!*addsq(!*multsq(a,a),
		     !*exptsq(x1,4)),
		 !*addsq(multsq(-8 ./ 1,!*multsq(x1,b)),
		     !*multsq(!*multsq(x1,x1),
					      multsq(-2 ./ 1,a)))),
	     !*invsq x3);
    y3:=!*addsq(y1,!*multsq(!*multsq(!*addsq(x3,negsq x1),
		     !*addsq(a,multsq(3 ./ 1,
			     !*multsq(x1,x1)))),
		 !*invsq multsq(2 ./ 1,
						y1))) >>
      else x3:=(y3:=inf)
    else if x1 = inf
      then <<
	x3:=x2;
	y3:=y2 >>
      else if x2 = inf
	then <<
	  x3:=x1;
	  y3:=y1 >>
	else <<
      x3:=!*multsq(!*addsq(!*multsq(a,!*addsq(x1,x2)),
		   !*addsq(multsq(2 ./ 1,b),
		       !*addsq(!*multsq(!*multsq(x1,x2),
			    !*addsq(x1,x2)),
					       multsq(-2 ./ 1,
			    !*multsq(y1,y2))))),
	       !*invsq !*exptsq(!*addsq(x1,negsq x2),2));
      y3:=!*multsq(!*addsq(!*multsq(!*addsq(y2,negsq y1),x3),
		   !*addsq(!*multsq(x2,y1),
		       !*multsq(x1,negsq y2))),
	       !*invsq !*addsq(x1,negsq x2)) >>;
  if x3 = inf
    then return x3.y3;
  lhs:=!*multsq(y3,y3);
  rhs:=!*addsq(b,!*multsq(x3,!*addsq(a,!*multsq(x3,x3))));
  if numr !*addsq(lhs,negsq rhs) % We can't just compare them
                  % since they're algebraic numbers.
                  % JHD Jan 14th. 1987.
    then <<
      prin2t "Point defined by X and Y as follows:";
      printsq x3;
      printsq y3;
      prin2t "on the curve defined by A and B as follows:";
      printsq a;
      printsq b;
      prin2t "gives a consistency check between:";
      printsq lhs;
      printsq rhs;
      interr "Consistency check failed in elliptic!-add" >>;
  return x3.y3
  end;

symbolic procedure infinitep u;
   kernp u and (mvar numr u eq 'infinite);

symbolic procedure lutz!-nagell!-bound(point,equation);
begin
  scalar x,y,a,b,lutz!-alist,n,point2,p,l,ans;
    % THE LUTZ!-ALIST is an association list of elements of the form
    % [X-value].([Y-value].[value of N for this point])
    % See thesis, chapter 7, algorithm LUTZ!-NAGELL, step [1].
  x:=car point;
  y:=cdr point;
  if !*tra or !*trmin then <<
    printc "Point to have torsion investigated is";
    printsq x;
    printsq y >>;
  a:=car equation;
  b:=cdr equation;
  if denr y neq 1 then <<
    l:=denr y;
    % we can in fact make l an item whose cube is > denr y.
    y:=!*multsq(y,!*exptf(l,3) ./ 1);
    x:=!*multsq(x,!*exptf(l,2) ./ 1);
    a:=!*multsq(a,!*exptf(l,4) ./ 1);
    b:=!*multsq(b,!*exptf(l,6) ./ 1) >>;
  if denr x neq 1 then <<
    l:=denr x;
    % we can in fact make l an item whose square is > denr x.
    y:=!*multsq(y,!*exptf(l,3) ./ 1);
    x:=!*multsq(x,!*exptf(l,2) ./ 1);
    a:=!*multsq(a,!*exptf(l,4) ./ 1);
    b:=!*multsq(b,!*exptf(l,6) ./ 1) >>;
  % we now have integral co-ordinates for x,y.
  lutz!-alist:=list (x . (y . 0));
  if (x neq car point) and (!*tra or !*trmin) then <<
    printc "Point made integral as ";
    printsq x;
    printsq y;
    printc "on the curve with coefficients";
    printsq a;
    printsq b >>;
  point:=x.y;
  equation:=a.b;
  n:=0;
loop:
  n:=n+1;
  point2:=elliptic!-multiply(x.y,2,equation);
  x:=car point2;
  y:=cdr point2;
  if infinitep x
    then return 2**n;
  if denr x neq 1
    then go to special!-denr;
  if a:=assoc(x,lutz!-alist)
    then if y = cadr a
      then return (ans:=lutz!-reduce(point,equation,2**n-2**(cddr a)))
      else if null numr !*addsq(y,cadr a)
    then return (ans:=lutz!-reduce(point,equation,2**n+2**(cddr a)))
    else interr "Cannot have 3 points here";
  lutz!-alist:=(x.(y.n)).lutz!-alist;
  if ans
    then return ans;
  go to loop;
special!-denr:
  p:=denr x;
  if not primep p
    then return 'infinite;
  n:=1;
  n:=1;
loop2:
  point:=elliptic!-multiply(point,p,equation);
  n:=n*p;
  if infinitep car point
    then return n;
  if quotf(p,denr car point)
    then go to loop2;
  return 'infinite
  end;


symbolic procedure lutz!-reduce(point,equation,power);
begin
  scalar n;
  if !*tra or !*trmin then <<
    princ "Point is of order dividing ";
    printc power >>;
  n:=1;
  while evenp power do <<
    power:=power/2;
    n:=n*2;
    point:=elliptic!-add(point,point,equation) >>;
    % we know that all the powers of 2 must appear in the answer.
  if power = 1
    then return n;
  if primep power
    then return n*power;
  return n*lutz!-reduce2(point,equation,power,3)
  end;



symbolic procedure lutz!-reduce2(point,equation,power,prime);
if power = 1
  then if infinitep car point
    then 1
    else nil
  else if infinitep car point
    then power
    else begin
      scalar n,prime2,u,ans;
      n:=0;
      while cdr divide(power,prime)=0 do <<
        n:=n+1;
        power:=power/prime >>;
      prime2:=nextprime prime;
      for i:=0:n do <<
    u:=lutz!-reduce2(point,equation,power,prime2);
        if u
          then <<
              ans:=u*prime**i;
              i:=n >>
         else <<
          power:=power*prime;
      point:=elliptic!-multiply(point,prime,equation) >> >>;
      if ans
        then return ans
        else return nil
      end;

endmodule;

end;