clear;

Fq:=GF(503);
a:=0; b:=1;
E:=EllipticCurve([Fq|a,b]);

P:=E![433,98];
Q:=E![127,258];
R:=E![339,199];
S:=E![124,42];
T:=E![147,443];

DP:=Divisor(P);
DQ:=Divisor(Q);

F<x,y>:=FunctionField(E);

lambda:=(3*P[1]^2+a)/(2*P[2]);
nu:=P[2]-lambda*P[1];
f2P:=y-(lambda*x+nu);

lambda:=(3*Q[1]^2+a)/(2*Q[2]);
nu:=Q[2]-lambda*Q[1];
f2Q:=y-(lambda*x+nu);

f:=f2P/f2Q;

g<x>:=Interpolation([R[1],S[1],T[1]],[R[2],S[2],T[2]]);
g:=y-g;

Proj<X,Y,Z>:=ProjectiveSpace(Fq,2);
INT:=Curve(P,[Y*Z-(252*X^2+374*X*Z+302*Z^2)]);
EC:=Curve(P,Y^2*Z-(X^3+a*X*Z^2+b*Z^3));

IntersectionPoints(INT,EC);
U:=E![413,369]; //The new point of intersection 

Support(Divisor(f));
Support(Divisor(g));

Evaluate(f,R)*Evaluate(f,S)*Evaluate(f,T)*Evaluate(f,U)/(Evaluate(f,Identity(E)))^4;

Evaluate(g,P)^2*Evaluate(g,(-2*P))/(Evaluate(g,Q)^2*Evaluate(g,(-2*Q)));

