load "curve.magma";
  // The finite field is k = F_p, where p = 10^8 + 7 is the
  // first prime greater than 10^8.
  // FG, FG2 in the source file represent two sample elements
  // of the Jacobian, corresponding to divisors D and D2
  // (more accurately, to the degree zero divisor classes
  // [D-3P_infty] and [D2-3P_infty] ).
  //

FG;
FG2;
  //The values of a,...,f and the inverse ainv=1/a for each
  //divisor.
  //

FGaddition := makeAddition(FG,FG2);
FGdoubling := makeDoubling(FG);
FGaddition;
FGdoubling;
  //Just to show the form of the output.
  //
  //

  //Now start the verification in terms of polynomials,
  //using the routines for ideals in polynomial rings
  //that are built in to Magma.  Along the way, we verify
  //that various intermediate results are correct. 
F:=getF(FG);
G:=getG(FG);
F2:=getF(FG2);
G2:=getG(FG2);
Faddition:=getF(FGaddition);
Gaddition:=getG(FGaddition);
Fdoubling:=getF(FGdoubling);
Gdoubling:=getG(FGdoubling);
F,G;
F2,G2;
Faddition,Gaddition;
Fdoubling,Gdoubling;



I:=makeIdeal(FG);
I2:=makeIdeal(FG2);
  //These are the ideals in the polynomial ring RR = k[x,y]
  //generated by {F,G} and {F2,G2} respectively.
  //Their images in the coordinate algebra RR/<ff> of the
  //affine open set C-{P_infty} are the elements vanishing
  //on the divisors D and D2, respectively.
  //

IC := ideal<RR|ff>;
  //Thus the coordinate algebra is RR/IC, and all the ideals
  //that we consider should contain IC.
  //

Iaddition:=makeIdeal(FGaddition);
Idoubling:=makeIdeal(FGdoubling);
Iinvaddition:=makeIdeal(makeInverse(FGaddition));
Iinvdoubling:=makeIdeal(makeInverse(FGdoubling));
  //Make the ideals representing the sum and the double
  //in the Jacobian.  Also make the ideals of their 
  //inverses (this will be used for checking the result).
  //For convenience later, let us write Daddition and Ddoubling
  //for the degree three divisors representing the sum and
  //the doubling in the Jacobian, respectively.  Thus
  // [Daddition-3P_infty] = [D-3P_infty] + [D2-3P_infty],
  // [Ddoubling-3P_infty] = 2[D-3P_infty].
  //Similarly, we refer to the divisors Dinvaddition and
  //Dinvdoubling, for which
  // [Dinvaddition-3P_infty] = -[Daddition-3P_infty],
  // [Dinvdoubling-3P_infty] = -[Ddoubling-3P_infty].
  // 

ff in I;
ff in I2;
Dimension(RR/I);
Dimension(RR/I2);
  //check that the divisors D and D2 actually lie on the curve
  //and are of degree 3.

ff in I*I2;
ff in I*I;
Dimension(RR/(I*I2));
Dimension(RR/(I*I));
  //Here I*I2 represents D+D2, which lies on the curve even
  //when viewed as a dimension zero subscheme of the plane.
  //This is because D and D2 are disjoint, as expected for
  //typical divisors.
  //However, I*I represents 2D as a dimension zero 
  //subscheme of the plane, and this is NOT on the curve.
  //The reason is that the "double points" of D, when viewed
  //on the plane, are "thicker" than the double points
  //when viewed on C itself, which only extend along the
  //tangent directions to C.  (Double points in the plane
  //extend along all tangent directions, i.e. with two 
  //"extra" dimensions of vanishing per double point.)
  //

ItimesI2 := I*I2;
ItimesI := I*I + IC;
Dimension(RR/ItimesI);
  //The correct answer is that one should consider the
  //product (I/IC)*(I/IC) of ideals in the Dedekind
  //domain RR/IC.  The pullback to RR is then
  //I*I + IC.  On the other hand, I*I2 already contains IC.
  //

MprimeFGplusFG2:=makeMprimeForAddition(FG,FG2);
MprimeFGplusFG:=makeMprimeForDoubling(FG);
MprimeFGplusFG2;
v1v2addition:=kernelOfMprime(MprimeFGplusFG2);
v1v2doubling:=kernelOfMprime(MprimeFGplusFG);
v1v2addition;
staddition:=makeCoefficientsOfst(FG,v1v2addition);
stdoubling:=makeCoefficientsOfst(FG,v1v2doubling);
staddition;
  //Make the pairs {s,t} in each case of addition
  //and of doubling.  We have shown the intermediate
  //results for the case of addition.
  //

saddition:=gets(staddition);
taddition:=gett(staddition);
saddition, taddition;
saddition in ItimesI2;
taddition in ItimesI2;
ideal<RR|saddition,taddition,ff> eq ItimesI2;
  //Check that the polynomials s,t actually belong
  //to ItimesI2, and in fact that they generate ItimesI2
  //when we work in the quotient RR/IC.
  //

sdoubling:=gets(stdoubling);
tdoubling:=gett(stdoubling);
sdoubling, tdoubling;
sdoubling in ItimesI;
tdoubling in ItimesI;
ideal<RR|sdoubling,tdoubling,ff> eq ItimesI;
  //Same for sdoubling and tdoubling with ItimesI.
  //
  //

  //Now we verify the correctness of our answers.
(I*I2*Iinvaddition + IC) eq ideal<RR|saddition,ff>;
  //Check that the principal ideal in RR/IC generated
  //by saddition corresponds to the sum of divisors
  // D + D2 + Dinvaddition.
  //Thus the principal divisor of saddition on the
  //curve C is
  //  (saddition) = D+D2+Dinvaddition - 9P_infty.
  //

(I*I*Iinvdoubling + IC) eq ideal<RR|sdoubling,ff>;
  //Do the same for sdoubling.
  //

  //Also verify that Daddition and Dinvaddition
  //represent inverses of each other in the Jacobian.
  //Similarly for doubling.
(Iaddition*Iinvaddition + IC) eq ideal<RR|Faddition,ff>;
(Idoubling*Iinvdoubling + IC) eq ideal<RR|Fdoubling,ff>;
  //This is because the principal ideal in RR/IC generated
  //by Faddition corresponds to Daddition + Dinvaddition.
  //In terms of principal divisors on C,
  // (Faddition) = Daddition + Dinvaddition - 6P_infty.