solving third degree modular equations in matlab?

hey

i am trying to write a program to find x and y for equations of the form :

y^2 mod p = (x^3+ax+b)mod p

y--> unknown. range from 0 to p-1. x-->unknown. range from 0 to p-1. p--> a large prime number.

i am currently using plain brute force to solve this problem. that is, for every value of x, y is varied from 0 to p-1 till lhs equals rhs.

while this works fine for smaller numbers (uptil 5 digit primes), for larger primes the program takes too long (over 3 hours) to give me the answer. is there any other, faster way to solve the above equation?

[Muthu Annamalai]

You can vectorize it and gain some speedup by calculating,

[xx,yy]=meshgrid(x,y);

z=( yy.^2 mod p ) + ( xx.^3 + ax + b mod p)

idx=find(z == 0);

% your solutions are,
yy(idx),xx(idx)

[Adam Danz]

Muthu's solution will be MUCH faster and two nested loops. The x and y inputs to meshgrid will be vectors of every possible value to be tested.

Answer 1

[John D'Errico]

Far too late for an answer. Sorry, as I just saw this question. But it is also perhaps quite a bit easier than you may think. A virtue of this solution is I show it is quite easy to solve, as well as pretty fast to compute. A solution will exist for roughly 50% or all x for any given p, a, and b (and when that solution does esist, there will be two such solutions. Thus if (x,y) is a solution pair, then (x,p-y) will also be a solution pair.

This last claim falls out of the fact that

mod((p-y)^2,p) == mod(p^2 - 2*p*y + y^2,p) == mod(y^2,p)

I'll use the symbolic toolbox, as the computations willl overflow a double when I start squaring things for sufficiently large p. I could also be using my own VPI tools, though they are slower then the symbolic TB.

Let me pick some large prime p, as well, values for a and b.

p = nextprime(sym(1e12))
p =
1000000000039

a = 539;
b = 100000001;

I hope that is sufficiently large for p.

Now, just choose any integer value for x. I'm not being very creative here. Sorry about that.

x = 1e5;

rhs = mod(sym(x)^3 + a*x + b,p);
rhs =
153861001

What good does that do us? A it turns out, lots. The rest is now easy, as all we need to do is compute the modular square root of 153861001, with respect to modulus p.

That is, we need to find a value oy y, such that y^2 is congruent to rhs, modulo p.

Hmm. My vpi toolbox has a modular square root functionality. And mupad has numlib::sqrtmodp. But the docs tell me that will be removed. Sigh. Even the Java.Math.BigInteger class offers a modular square root, but Java is also going away from MATLAB before long.

This next call should work as long as the mupad functionality remains, although it seems to take really long for large p. Its not at all bad for p on the order of 1e12 though.

y = feval(symengine, 'numlib::sqrtmodp', rhs, p)
y =
299284269818

Is that a solution to the problem? Yes, it must be so.

mod(y^2 - (x^3 +a*x + b),p) 
ans =
0

The point of all this is the solution is easy, as long as you can compute a modular square root. Pick any value of x. It will generate a result for the right hand side. And with p a prime, roughly 50% of the time, a modular square root will exist. That is true because for some values of x, the right hand side number will not have a square root, modulo p.

Consider a different value for x:

x
x =
      100003
rhs = mod(sym(x)^3 + a*x + b,p)
rhs =
90156562645

y = feval(symengine, 'numlib::sqrtmodp', rhs, p)
y =
660417945355

mod(y^2 - (x^3 +a*x + b),p)
ans =
819686874749

The failure here arises because 90156562645 does not have a square root mod p. And that will happen roughly 50% of the time. The call to sqrtmodp failed, but it did not return a warning. We could have predicted that failure using this call:

feval(symengine, 'numlib::isquadres', rhs, p)
ans =
FALSE

Oh well. Don't blame me. I should probably provide a code based on Shanks-Tonelli, that will run for doubles or in the symbolic toolbox. Its not that difficult as I recall. When p is a prime of the form 4*k+3, the answer is trivial. It is only when p is of the form 4*k+1 there is any problem at all.

Finally, can we express all possible solutions to the problem for a given a, b, and prime modulus p?

As a simple example, consider p=17. (My personally favorite prime.) And still not feeling very creative, I'll pick a and b randomly. What really matters is that b is not a multiple of p.

p = 17;
a = randi(p-1)
a =
     2
     
b = randi(p-1)
b =
     5

So as long as p is prime, and b was less than p, We are ok in that respect.

We wish to investigate all solutions of the problem

mod(y^2 - (x^3 + a*x + b),p) == 0

Its sort of a modified Pell equation. So, first, consider what happens when x exceeds p.

x = reshape(0:3*p-1,p,3)
x =
     0    17    34
     1    18    35
     2    19    36
     3    20    37
     4    21    38
     5    22    39
     6    23    40
     7    24    41
     8    25    42
     9    26    43
    10    27    44
    11    28    45
    12    29    46
    13    30    47
    14    31    48
    15    32    49
    16    33    50
    
rhs = mod(x.^3 + a*x + b,p)
rhs =
     5     5     5
     8     8     8
     0     0     0
     4     4     4
     9     9     9
     4     4     4
    12    12    12
     5     5     5
     6     6     6
     4     4     4
     5     5     5
    15    15    15
     6     6     6
     1     1     1
     6     6     6
    10    10    10
     2     2     2

You should recognize what happens. Consider the value x0<p. What happens to the modular form of the cubic polynomial (x^3 + a*x + b), for x=x0+p? As we might expect, the values of this polynomial will periodically repeat.

It is important to recognize that when the polynomial in x returns a zero solution, that the only solution for y is then y==0. However, numlib::sqrtmodp seems to fail for that case. (I really do need to write that code. Sigh.)

So as long as we are concerned only with unique solutions to the problam at hand, we need consider only x in the set [0:p-1].

x= 0:p-1;
rhs = mod(x.^3 + a*x + b,p);

Now, unfortunately, the mupad tools do not seem to be vectorized, but a loop will suffice here.

y = nan(size(x));
for i = 1:numel(x)
  if rhs(i) == 0
    y(i) = 0;
  elseif feval(symengine, 'numlib::isquadres', rhs(i), p)
    y(i) = feval(symengine, 'numlib::sqrtmodp', rhs(i), p);
  end
end

[x,rhs,y,p-y]
ans =
     0     5   NaN   NaN
     1     8    12     5
     2     0     0    17
     3     4     2    15
     4     9    14     3
     5     4     2    15
     6    12   NaN   NaN
     7     5   NaN   NaN
     8     6   NaN   NaN
     9     4     2    15
    10     5   NaN   NaN
    11    15     7    10
    12     6   NaN   NaN
    13     1     1    16
    14     6   NaN   NaN
    15    10   NaN   NaN
    16     2     6    11

Don't forget that if (x,y) is a solution, then (x,p-y) is also a solution, with the one special case when the modular polynomial in x results in zero.

So when p==17 with a==2 and b==5, we have found all the solutions, shown as columns [1,3,4] of the displayed array. As I claimed in the beginning of my answer, we will expect this to happen for essentially 50% of the possible values of x<p.

Sadly, I'm now done kicking this problem virtually to death. Why am I sad? Because I only noticed the question almost 9 years after it had been asked, so my answer is not terribly a timely one.