Math 640:348 Prof. Kontorovich
Spring 2015, 3/31 lecture

We will use the index calculus method to solve a Discrete Log Problem

First lets get a “safe” prime

p=Prime[11113]

FactorInteger[(p-1)/2]

q=(p-1)/2

Good, so q and p are a “Sophie Germain prime pair” - in the literature q is often called a Sophie Germain prime. This means the Discrete Log Problem (DLP) mod p is not susceptible to Pohlig-Hellman attacks. 

Next we need a way of testing whether a number is smooth, say 5-smooth (meaning all its prime factors are 2, 3, or 5). Here is a number called “smooth” with the property that: any n<117000 which is 5-smooth has gcd(n,smooth)=n. That is, all prime power factors (r^e for prime r) of any 5-smooth number less than p are also factors of “smooth”:

smooth=2^Floor[Log[2,p]]3^Floor[Log[3,p]]5^Floor[Log[5,p]]

Then “is5smooth” returns “True” if n is 5-smooth, and “False” otherwise. (There are surely better ways to implement this...)

is5smooth[n_]:=(GCD[n,smooth]==n);

Now let’s get a random primitive root, say,

g=7;
MultiplicativeOrder[g,p]==p-1

In class, we chose a random value of a, say

a=1919

Ok, now we’re ready to solve the DLP. We need to find x (remember that x is only determined mod p-1, which is the same as mod 2q) so that g^x = a (mod p).

Step 1: solve the DLP for small primes.

Let’s randomly test some values of j, hoping to find a smooth value for g^j(mod p):

For[j=5000,j<6000,j++,
If[
is5smooth[PowerMod[g,j,p]]
,
Print["g^"<>ToString[j]<>" mod p is smooth"]
];
]

We found three values,

j1=5189;
j2=5664;
j3=5838;

How do the values of g^j(mod p) factor?

FactorInteger[PowerMod[g,j1,p]]

FactorInteger[PowerMod[g,j2,p]]

FactorInteger[PowerMod[g,j3,p]]

Writing ell2 for the exponent of g which gives 2, and similarly ell3, ell5, we have the system of equations (mod 2q):

3 ell2 + ell3 + 2 ell5 = j1
2 ell2 + 4 ell3 + ell5 = j2
8 ell2 + 2 ell3 + 2 ell5 = j3,

or in “augmented” matrix form: 
(in Mathematica, use “MatrixForm” to make it look like a matrix, instead of a sequence)

mat={
{3,1,2,j1},
{2,4,1,j2},
{8,2,2,j3}
};
MatrixForm[mat]

Now we just need to Gaussian Eliminate this matrix to solve for ell2, ell3, and ell5. But there’s a catch! The coefficients will in general not be invertible mod 2q! 

So you should first solve the equations mod 2, then mod q, and then “Chinese Remainder Theorem” them back together.

First let’s look mod 2:

MatrixForm[Mod[mat,2]]

So our equations are: 

ell2+ell3 = 1 (mod 2)
ell5 = 0 (mod 2)

We do not have a full-rank matrix, so can’t solve exactly; that’s ok, in the end we’ll have to guess whether ell2=0 or 1(mod 2), from which everything else will be determined...

Now let’s look mod q, and Gaussian Eliminate:
Again, the matrix is:

MatrixForm[mat]

So if we want to use the first row to eliminate coefficients below the “3”, we first need to invert the top row. First we need to know the inverse of 3 mod q:

inverse3modq=PowerMod[3,-1,q]

And now we multiply the top row by 3^(-1), that is, multiply the matrix “mat” on the left by a 3x3 diagonal matrix with diagonal entries 3^(-1), 1, and 1. And of course reduce everything mod q:

mat1=Mod[
DiagonalMatrix[{inverse3modq,1,1}].mat
,q];
MatrixForm[mat1]

Next subtract off 2x(top row) from the second row, and 8x(top row) from the third row, as always, reducing mod q:

mat2=Mod[{{1,0,0},{-2,1,0},{-8,0,1}}.mat1,q];
MatrixForm[mat2]

We continue Gaussian Elimination; now we need to turn that “19633” in the middle into a “1”, so we need its inverse mod q

inverse19633=PowerMod[19633,-1,q]

Multiply the second row by this number, and reduce mod q

mat3=Mod[DiagonalMatrix[{1,inverse19633,1}].mat2,q];
MatrixForm[mat3]

Now subtract 19630x(second row) from the first row, and 19629x(second row) from the third row:

mat4=Mod[{{1,-19630,0},{0,1,0},{0,-19629,1}}.mat3,q];
MatrixForm[mat4]

Again invert the last diagonal “35330”:

inverse35330=PowerMod[35330,-1,q]

Multiply through

mat5=Mod[DiagonalMatrix[{1,1,inverse35330}].mat4,q];
MatrixForm[mat5]

And subtract

mat6=Mod[{{1,0,-41223},{0,1,-53000},{0,0,1}}.mat5,q];
MatrixForm[mat6]

Yay! Now we know that

ell2q=55378;
ell3q=18204;
ell5q=57648;

I called these ell2“q”, etc., with q’s at the end are because these are the values mod q, not mod 2q=p-1. Here comes the Chinese Remainder Theorem step:

We already know that ell5=0(mod 2), so also knowing its value mod q determines its value mod 2q:

ell5 = ell5q

(I trust you can figure out why I did that.) 

Let’s check to make sure we didn’t make any mistakes - does raising g to this power mod p give us 5?

PowerMod[g,ell5,p]

Great!

For ell2 and ell3, we do not know their values mod 2, but we do know that 
ell2+ell3=1(mod 2), which means they’re different (one odd, one even). 
Let’s guess that ell2 is even. If that really was the case, then ell2 would be the same as ell2q (which is already even). So we test:

PowerMod[g,ell2q,p]

Aha! This is *not* 2, so ell2q is not ell2! That means that ell2 is in fact odd. So:

ell2=ell2q+q

(Again I trust you can figure out why I did that!...)

Let’s test if we got it right:

PowerMod[g,ell2,p]

Two down, one to go! Since ell2 is odd, we know that ell3 is even, so

ell3=ell3q

As always, check your work:

PowerMod[g,ell3,p]

Now Step 1 is done - we have solved the DLP for the small primes 2, 3, and 5.

Extra credit question: what would happen if we just row reduced the original matrix, and interpreted fractions as inverses mod 2q? Would we find the values of ell2, ell3, and ell5 directly? Why or why not?...

Step 2: Find smooth values of a g^-j

The next step is to find some random j so that a g^(-j) (mod p) is also smooth; then using the DLP solution above, we should be able to solve DLP for “a”. Again, we loop over possible j values:

For[j=5000,j<6000,j++,
If[
is5smooth[Mod[ a PowerMod[g,-j,p],p]]
,
Print["a g^-"<>ToString[j]<>" mod p is smooth"]
];
]

Great! We found three, but should only need one j. Let’s play with the value

j0=5057

Then a g^(-j0) (mod p) is

Mod[ a PowerMod[g,-j0,p],p]

which factors as

FactorInteger[59049]

Solving for “a” in the equation

a g^(-j0) = 3^10 = (g^ell3)^10 = g^(10 ell3)

gives

a = g^(j0 + 10 ell3),

so

x = Mod[ j0 + 10 ell3, p-1]

should be our desired exponent! (Note that we reduced mod p-1.) Did it work?

PowerMod[g,x,p]

a

The index calculus wins again...