epoints := proc(ec, x, ub, p)
    local ecurve, z, pct, k, i;
    pct := 0;
    for k from 0 to p-1 while pct <= ub do
        z := subs(x=k, ec) mod p;
        if z = 0 then
           pct := pct + 1;
           ecurve[pct] := [k, z];
        fi:
        if z &^ ((p-1)/2) mod p = 1 then
           z := z &^ ((p+1)/4) mod p;
           ecurve[pct+1] := [k, z];
           ecurve[pct+2] := [k, -z mod p];
           pct := pct + 2;
        fi:
    od:
    if pct > ub then
       pct := ub:
    fi:
    seq(ecurve[i], i = 1..pct):
end: