Fast Exponentiation to find all points on Elliptic Curve over a Finite Field showing incorrect results for some primes.

Function to find the points:

function EC_Points = ECPoints(a, b, p)
% Input:  Two integers a, b, and a prime p > 3.  
% 
%  It is assumed that 4a^3 + 27b^2 is non-zero, 
% so that the associated elliptic curve
%         y^2 = x^3 + ax + b
%is nonsingular.
% Output:  Points, a 2 column vector listing all points belonging to this
% elliptic curve mod p, including the point at infinity.  Each row of
% points will be a point on this elliptic curve.
% The method will simply run through all possible values for x (mod p),
% compute the right side r = x^3 + ax + b, and check to see whether
% r^(p+1)/4 (mod p) is a square root of r.
PointIndex = 1;
for x = 0:p-1
    Y = mod(x^3 + a*x + b, p);
    rpower = FastExp(Y,(p+1)/4,p);
    if mod(rpower^2,p)==Y %+/- rpower are sqare roots:
        if Y == 0  %only one square root
            EC_Points(PointIndex, :) = [x 0];  PointIndex = PointIndex + 1; 
        else
            EC_Points(PointIndex, :) = [x rpower];  PointIndex = PointIndex + 1; 
            EC_Points(PointIndex, :) = [x p-rpower];  PointIndex = PointIndex + 1; 
        end
    end
end
EC_Points(PointIndex, :) = [Inf Inf];

Function for the FastExp

function c = FastExp(a,b,p)
% Input:  a, b, < p positive integers, p = the modulus
% Output: c = a^b (mod p) 
% Program uses the fast exponentiation method to calculate c = a^b (mod p).  
% Beginning with factor = a^1: 
% on each pass, it is determined whether this is a desired factor.  If so, then 
% it is multiplied into a variable containing a cumulative product of desired 
% factors.  Then the factor is squared, and the process repeated.

expon = b;
c = ones(size(a));
factor = a;

while(expon > 0)
    if (mod(expon,2) == 1)
        c = mod(c .* factor,p);
    end
    expon = floor(expon/2);
    factor = mod(factor.*factor,p);
end

RESULTS

Running the program for a=-1, b=188, and p=11369

EC_Points=ECPoints(-1, 188, 11369);

EC_Points =

209 1

209 11368

3676 1

3676 11368

5926 0

7484 1

7484 11368

Inf Inf

This is not correct as the curve has a lot of other points; code works with most primes p but not some such as p=11369 as shown above.

Bruteforce method to find all points on the curve takes a lot of time and thus, I tried implementing with the fast exponentiation method. If anyone has any better way to implement the fast exponentiation method or can find flaws in my code that leads to such incorrect result for some primes, please do share/let me know.

Thank you in advance!

Answer 1

[Umang Pandey]

Hi Gazi,

The issue seems to be in the code you are using to check whether the number is square in the finite field. Before computing the fourth root of the "Y", you have to verify if "Y" is a quadratic residue, i.e., it has a square and a square root.

You can modify the for loop as follows to check if "Y" is a quadratic residue modulo "p" before computing the square root:

for x = 0:p-1 
    Y = mod (x^3 + a*x + b, p); 
    if FastExp(Y, (p-1)/2, p) == 1 % Check if Y is a square 
        rpower = FastExp(Y, (p+1)/4, p); % Compute the square root of Y 
        if Y == 0 %only one square root 
            EC_Points(PointIndex, :) = [x 0]; PointIndex = PointIndex + 1;  
        else 
            EC_Points(PointIndex, :) = [x rpower]; PointIndex = PointIndex + 1;  
            EC_Points(PointIndex, :) = [x p-rpower]; PointIndex = PointIndex + 1;  
        end 
    end 
end 

Best,

Umang