#OK to post homework
#Ramesh Balaji,4/7/2024,hw21

with(NumberTheory);
RP:=proc(a,b) local x: x:=rand(0..200)()/100: [x,sqrt(x^3+a*x+b)]:end:

# Part 3: RandPt
RandPt := proc(a,b,p):
    y:=0.5;
    while not type(y, integer) do:
        x:=rand(0..p-1)();
        y:=sqrt(x^3 + a*x + b) mod p;
    od:

    return [x,y];
end:

RandPt(1,1,10);

# Part 4: Alice and Bob

AliceToBobEC := proc(a,b,p,g):
    A:=rand(0..p-1)();
    return A,(A*g) mod p;
end:

BobToAliceEC := proc(a,b,p,g):
    B:=rand(0..p-1)();
    return B,(B*g) mod p;
end:


for i from 0 to 10 do:
    a:=rand(0..10)();
    b:=rand(0..10)();
    p:=ithprime(rand(1..10)());
    g:=RandPt(a,b,p);
    A,Ag:=AliceToBobEC(a,b,p,g);
    B,Bg:=BobToAliceEC(a,b,p,g);

    # all OK printed.

    if B*Ag mod p <> A*Bg mod p then:
        print("FAIL");
    else
        print("OK");
    fi:
end:

# Part 5:

Hasse := proc(a,b,K):
    npts := table();
    i:=3;
    while i <= K do:
        # number of points on the elliptic curves
        npts_i := 0;

        for j from 0 to i-1 do:
            if type(sqrt(j^3 + a*j + b) mod p, integer) then:
                if type(npts[i], integer) then:
                    npts[i] += 1;
                else:
                    npts[i] := 0;
                fi:
            fi:
        od:

        if not type(npts[i], integer) then:
            npts[i] := 0;
        fi:

        i:=nextprime(i);
    od:

    maxp := 3;
    minp := 3;
    for i, j in eval(npts) do:
        # i is the prime

        newj:=(j - (i+1)) / (sqrt(i));

        if evalf(newj) < evalf(npts[minp]) then:
            minp:=i;
        fi:

        if evalf(newj) > evalf(npts[maxp]) then:
            maxp:=i;
        fi:

    od:

    return [minp, maxp, (npts[minp] - (minp+1)) / (sqrt(minp)), (npts[maxp] - (maxp+1)) / (sqrt(maxp))];
end:


for a from 1 to 10 do:
    for b from 1 to 10 do:
        h:=Hasse(a,b,500);
        print(cat("a ", a, " b ", b, " h ", h))
    od:
od:

# They all return [227 (min),3 (max)]. I'm not sure if this is right...
# it is very weird. Maybe I did something wrong in calculating the min/max? Is
# evalf rounding, because Maple did not let me compare the rational values?