/////////////////////////////////////////////////////////////////////////////

/*
  
  General test implementation of the Weil pairing.
   No fancy speedups included.

  written by Florian Hess, version 20050829
  
*/

/////////////////////////////////////////////////////////////////////////////
  
  
// The coordinate functions

x := function(P)
   assert not IsZero(P);
   return Eltseq(P)[1];
end function;
   
y := function(P)   
   assert not IsZero(P);
   return Eltseq(P)[2];
end function;   
   
  
f := function(A, B, C)  
   
   // Q may be zero then problem
   if IsZero(C) then
      error "C is zero in f";
   end if;
   
   if IsZero(A) and IsZero(B) then
      return BaseRing(Curve(A))!1;
   elif IsZero(A) then
      return x(C) - x(B);
   elif IsZero(B) then
      return x(C) - x(A);
   elif x(A) ne x(B) then
      return (x(B) - x(A)) * y(C) - (y(B) - y(A)) * x(C) 
             - x(B) * y(A) + x(A) * y(B);
   elif y(A) ne y(B) then   
      return x(C) - x(A);
   else
      a := Eltseq(Curve(A));
      return ( a[1] * y(A) - 3 * x(A)^2 - 2 * a[2] * x(A) - a[4] ) *    
             ( x(C) - x(A) ) + 
             ( 2 * y(A) + a[1] * x(A) + a[3] ) * ( y(C) - y(A) );
   end if;
   
end function;
   
  
g := function(D, A, B, C)
   
   // muss sicherstellen dass C nicht O
   // wird in w getan  
     
   num := f(A, B, C);
   
   if IsZero(num) then
      //assert C in { A, B, A+B, -A-B };
      //assert Log(D, C) ne -1;
      return false, _;
   end if;
   
   S := A + B;
   den := f(S, -S, C);
   
   if IsZero(den) then
      //assert C in { A, B, A+B, -A-B };
      //assert Log(D, C) ne -1;
      return false, _;
   end if;
   
   return true, num/den;
   
end function;
  
  
  
h := function(D, C1, C2, l)

   n := [];
   ll := l;
   while not IsZero(ll) do
      Append(~n, IsOdd(ll) select 1 else 0);
      ll := ll div 2; 
   end while;
   t := #n - 1;
   assert &+ [ n[i+1] * 2^i : i in [0..t] ] eq l;
   
   A := D;
   fr1 := BaseRing(Curve(D))!1;
   fr2 := BaseRing(Curve(D))!1;
   // (*)
   r := 1;  
   //assert A eq r*D;
   
   for i in Reverse([0..t-1]) do
      
      //print "loop in h ", i;
      
      ok1, g1 := g(D, A, A, C1);
      if not ok1 then
         //assert Log(D, C1) ne -1;
         return false, _, _;
      end if;
      ok2, g2 := g(D, A, A, C2);
      if not ok2 then
         //assert Log(D, C2) ne -1;
         return false, _, _;
      end if;
         
      fr1 := fr1 * fr1 * g1;
      fr2 := fr2 * fr2 * g2;   // could in fact be combined into the
                               // single variable fr1/fr2
      A := A + A;
      
      // (*)   
      r := 2 * r;
      //assert A eq r*D;
      
      if not IsZero(n[i+1]) then 
         
         ok1, g1 := g(D, A, D, C1);
         if not ok1 then
            //assert Log(D, C1) ne -1;
            return false, _, _;
         end if;
         ok2, g2 := g(D, A, D, C2);
         if not ok1 then
            //assert Log(D, C1) ne -1;
            return false, _, _;
         end if;
         
         fr1 := fr1 * g1;
         fr2 := fr2 * g2;
         A := A + D;
         
         // (*)
         r := r + 1;
         //assert A eq r*D;
         
      end if;
      
   end for;
   
   assert r eq l;
   if not IsZero(A) then
      error "A not in E[l]";
   end if;
  
   return true, fr1, fr2;
   
end function;


w := function(P, Q, l)
   
   assert IsOdd(l);
   
   E := Curve(P);
   
   if IsZero(P) or IsZero(Q) or x(P) eq x(Q) then
      return BaseRing(Curve(P))!1;
   end if;
   
   T := P - Q;
   
   PT := P + T;
   if IsZero(PT) then
      return BaseRing(Curve(P))!1;
   end if;
   
   ok1, num1, den1 := h(Q, PT, T, l);  
   //print "ok1 is ", ok1;
   if not ok1 then
      //assert Log(Q, PT) ne -1 or Log(Q, T) ne -1;
      return BaseRing(Curve(P))!1;
   end if;
   
   QT := Q - T;
   if IsZero(QT) then
      return BaseRing(Curve(P))!1;
   end if;
   
   ok2, num2, den2 := h(P, QT, -T, l);
   //print "ok2 is ", ok2;
   if not ok2 then
      //assert Log(P, QT) ne -1 or Log(P, -T) ne -1;
      return BaseRing(Curve(P))!1;
   end if;

   return num1 * den2 / (num2 * den1);
   
   // If ok1 or ok2 is false then
   // we have lam Q = 2P or lam Q = P or lam P = 2Q or lam P = Q
   // Since l is odd it follows that the pairing is one.  
     
   // If l is not odd: If lam is odd the pairing = 1,   
   // otherwise cen reduce to a l=2 case ... see manuscript
         
end function;
   
   
   
///////////////////////////////////////////////////////////////////////////
// Another implementation. This requires product reps in an essential way     
//  for large examples.

  
w_general := function(P, Q, l)  
   
   E := Curve(P);
   
   // if we don't do this for Q then the next while loop 
   // necessarily is infinite
   if IsZero(P) or IsZero(Q) then
      return BaseRing(Curve(P))!1;
   end if;
   
   // if the until-condition is not satisfied then
   // there will be zero numerator or denominator later on
   repeat
      T := Random(E);
      PT := P + T;   
   until #{ Identity(E), Q, T, PT} eq 4;
   
   // Everything must work from here  
     
   ok, ff := IsPrincipal( l * (Place(Q) - Place(Identity(E))) );  
   assert ok;
   
   ok, gg := IsPrincipal( l * (Place(PT) - Place(T)) );  
   assert ok;
   
   num := Evaluate( ff, Place(PT) ) / Evaluate( ff, Place(T) );
   den := Evaluate( gg, Place(Q) ) / Evaluate( gg, Place(Identity(E)) );
   
   return num/den;
   
end function;

  
   
///////////////////////////////////////////////////////////////////////////
// Test


F := GF( 3^97 );
R<x> := PolynomialRing( F );
F2 := ext< F |  x^3 - x - 1 >;
rho := F2.1;
F6 := ext< F2 | x^2 + 1 >;
sigma := F6.1;
E := EllipticCurve( [ 0, 0,0, F ! -1, F ! 1 ] );
P := Random(E );
N := #E;
time N*P;
E := EllipticCurve( [ 0, 0,0, F6 ! -1, F6 ! 1 ] );
P1 := E ! P;
// Image under distortion map
P2 := E ! [ rho - P1[1], sigma* P1[2] ];
time N*P1;
time N*P2;

print "Computing WeilPairing ...";
time a := WeilPairing( P1, P2, N );
print "Computing w Pairing ...";
time b := w(P1, P2, N);

assert a eq 1/b;