NOTE: This site is obsolete. See book draft (in PDF):

The Laws of Cryptography with Java Code.

---

The Extended Euclidean Algorithm.

The previous section introduced the field known as the integers mod p, denoted Z_p or GF(p). Most of the field operations are straightforward, since they are just the ordinary arithmetic operations followed by remainder on division by p. However the multiplicative inverse is not intuitive and requires some theory to compute. If a is a non-zero element of the field, then a^-1 can be computed eficiently using what is known as the extended Euclidean algorithm, which gives the greatest common divisor (gcd) along with other integers that give the inverse. It is the topic of the remainder of this section.

Law GCD1: The cryptographer's first and oldest favorite algorithm is the extended Euclidean algorithm, which computes the greatest common divisor of two positive integers a and b and also supplies integers x and y such that

x*a + y*b = gcd(a, b).

The Basic Euclidean Algorithm to give the gcd: Consider the calculation of the greatest common divisior (gcd) of 819 and 462. One could factor the numbers as: 819 = 3*3*7*13 and 462 = 2*3*7*11, to see immediately that the gcd is 21 = 3*7. The problem with this method is that there is no efficient algorithm to factor integers. In fact, the security of the RSA cryptosystem relies on the difficulty of factoring, and we need an extended gcd algorithm to implement RSA. It turns out that there is another better algorithm for the gcd -- developed 2500 years ago by Euclid (or mathematicians before him), called (surprise) the Euclidean algorithm.

The algorithm is simple: just repeatedly divide the larger one by the smaller, and write an equation of the form larger = smaller * quotient + remainder. Then repeat using the two numbers smaller and remainder. When you get a 0 remainder, then you have the gcd of the original two numbers. Here is the sequence of calculations for the same example as before:

819   =  462 * 1  +  357   (Step  0)
462   =  357 * 1  +  105   (Step  1)
357   =  105 * 3  +   42   (Step  2)
105   =  42  * 2  +   21   (Step  3, giving the GCD: 21)
42    =  21  * 2  +    0   (Step  4)

The proof that this works is simple: a common divisor of the first two numbers must also be a common divisor of all three numbers all the way down. (Any number is a divisor of 0 (sort of on an honorary basis).) One also has to argue that the algorithm will terminate and not go on forever, but this is clear since the remainders must be smaller at each stage.

Here is Java code for two versions of the GCD algorithm: one iterative and one recursive. (There is also a more complicated binary version that is efficient and does not require division.)

public static long gcd1(long x, long y) {
   if (y == 0) return x;
   return gcd1(y, x%y);
}

public static long gcd2(long x, long y) {
   while (y != 0) {
      long r = x % y;
      x = y; y = r;
   }
   return x;
}

Here is a complete Java program using the above two functions: gcd calculation.

The Extended GCD Algorithm: Given the two positive integers 819 and 462, the extended Euclidean algorithm finds unique integers a and b so that a*819 + b*462 = gcd(819, 462) = 21. In this case, (-9)*819 + 16*462 = 21.

For this example, to calculate the magic a and b, just work backwards through the original equations, from step 3 back to step 0 (see above). In red and green below are equations, where each shows two numbers x and y from a step of the original algorithm, and corresponding integers a and b (in green), such that a*x + b*y = gcd(x,y). Between each pair of equations is an equation in black that leads to the next equation.

1*105 + (-2)* 42 = 21                            (from Step 3 above)
(-2)*357 + (-2)(-3)*105 = (-2)*42 = (-1)*105 + 21  (Step 2 times -2)
(-2)*357 + 7*105 = 21       (Combine and simplify previous equation)
7 *462 + (7)(-1)* 357 = 7*105 = 2*357 + 21          (Step 1 times 7)
7*462 + (-9)*357 = 21       (Combine and simplify previous equation)
(-9)*819 + (-9)(-1)*462 = (-9)*357 = (-7)*462 + 21   (Step 0 * (-9))
(-9)*819 + 16*462 =  21                (Simplfy -- the final answer)

It's possible to code the extended gcd algorithm following the model above, first using a loop to calculate the gcd, while saving the quotients at each stage, and then using a second loop as above to work back through the equations, solving for the gcd in terms of the original two numbers. However, there is a much shorter, tricky version of the extended gcd algorithm, adapted from D. Knuth.

public static long[] GCD(long x, long y) {
   long[] u = {1, 0, x}, v = {0, 1, y}, t = new long[3];
   while (v[2] != 0) {
      long q = u[2]/v[2];
      for (int i = 0; i < 3; i++) {
         t[i] = u[i] - v[i]*q; u[i] = v[i]; v[i] = t[i];
      }
   }
   return u;
}

A complete Java program using the above function is here.

The above code rather hides what is going on, so with additional comments and checks, the code is rewritten below. Note that at every stage of the algorithm below, if w stands for any of the three vectors u, v or t, then one has x*w[0] + y*w[1] = w[2]. The function check verifies that this condition is met, checking in each case the vector that has just been changed. Since at the end, u[2] is the gcd, u[0] and u[1] must be the desired integers.

public static long[] GCD(long x, long y) {
   long[] u = new long[3];
   long[] v = new long[3];
   long[] t = new long[3];
   // at all stages, if w is any of the 3 vectors u, v or t, then
   //   x*w[0] + y*w[1] = w[2]  (this is verified by "check" below)
   // vector initializations: u = {1, 0, u}; v = {0, 1, v};
   u[0] = 1; u[1] = 0; u[2] = x; v[0] = 0; v[1] = 1; v[2] = y;
   System.out.println("q\tu[0]\tu[1]\tu[2]\tv[0]\tv[1]\tv[2]");

   while (v[2] != 0) {
      long q = u[2]/v[2];
      // vector equation:  t = u - v*q
      t[0] = u[0] - v[0]*q; t[1] = u[1] - v[1]*q; t[2] = u[2] - v[2]*q;
      check(x, y, t);
      // vector equation:  u = v;
      u[0] = v[0]; u[1] = v[1]; u[2] = v[2]; check(x, y, u);
      // vector equation:  v = t;
      v[0] = t[0]; v[1] = t[1]; v[2] = t[2]; check(x, y, v);
      System.out.println(q + "\t"+ u[0] + "\t" + u[1] + "\t" + u[2] +
                             "\t"+ v[0] + "\t" + v[1] + "\t" + v[2]);
   }
   return u;
}

public static void check(long x, long y, long[] w) {
   if (x*w[0] + y*w[1] != w[2]) {
      System.out.println("*** Check fails: " + x + " " + y);
      System.exit(1);
   }
}

Here is the result of a run with the data shown above:

q       u[0]    u[1]    u[2]    v[0]    v[1]    v[2]
 
1       0       1       462     1       -1      357
1       1       -1      357     -1      2       105
3       -1      2       105     4       -7      42
2       4       -7      42      -9      16      21
2       -9      16      21      22      -39     0

gcd(819, 462) = 21
(-9)*819 + (16)*462 = 21

Here is a run starting with 40902 and 24140:

q       u[0]    u[1]    u[2]    v[0]    v[1]    v[2]
 
1       0       1       24140   1       -1      16762
1       1       -1      16762   -1      2       7378
2       -1      2       7378    3       -5      2006
3       3       -5      2006    -10     17      1360
1       -10     17      1360    13      -22     646
2       13      -22     646     -36     61      68
9       -36     61      68      337     -571    34
2       337     -571    34      -710    1203    0

gcd(40902, 24140) = 34
(337)*40902 + (-571)*24140 = 34

Here is a complete Java program with the above functions, along with other example runs: Complete Extended GCD Program.

Exercise: Prove that the long (debug) version of the Extended GCD Algorithm works. First show that u[2] is the gcd by throwing out all references to array indexes 0 and 1, leaving just u[2], v[2], and t[2]. Show that this still terminates and just calculates the simple gcd, without reference to the other array indexes. (This shows that at the end of the complicated algorithm, u[2] actually is the gcd.)

Next show mathematically that the three special equations are true at the start of the algorithm, and that each stage of the algorithm leaves them true. (One says that they are left invariant.)

Finally deduce that algorithm is correct.

---

Fast Integer Exponentiation (Raise to a Power).

A number of cryptosystems require arithmetic on large integers. For example, the RSA public key cryptosystem uses integers that are at least 1024 bits long. An essential part of many of the algorithms involved is to raise an integer to another integer power, modulo an integer (taking the remainder on division).

Law EXP1: Many cryptosystems in modern cryptography depend on a fast algorithm to perform integer exponentiation.

It comes as a surprise to some people that in a reasonable amount of time one can raise a 1024 bit integer to a similar-sized power modulo an integer of the same size. (This calculation can be done on a modern workstation in a fraction of a second.) In fact, if one wants to calculate x¹⁰²⁴ (a 10-bit exponent), there is no need to multiply x by itself 1024 times, but one only needs to square x and keep squaring the result 10 times. Similarly, 20 squarings yields x^1048576 (a 20-bit exponent), and an exponent with 1024 bits requires only that many squarings if it is an exact power of 2. Intermediate powers come from saving intermediate results and multiplying them in. RSA would be useless if there were no efficient exponentiation algorithm.

There are two distinct fast algorithms for raising a number to an integer power. Here is pseudo-code for raising an integer x to power an integer Y:

Inputs: integers x, Y = Yk Yk-1 ... Y1 Y0  (in binary)
Output: integer z
Algorithm:
   int exp(int x, int Y[], int k) {
      int y = 0, z = 1;
      for (int i = k; i >= 0; i--) {
         y = 2*y;
         z = z*z;
         if (Y[i] == 1) {
            y++;
            z = z*x;
         }
      }
      return z;
   }

The variable y is only present to give a loop invariant, since at the beginning and end of each loop, as well as just before the if statement, the invariant x^y = z holds, and after the loop terminates y = Y is also true, so at the end, z = x^Y. To find x^y mod n one should add a remainder on division by n to the two lines that calculate z.

Here is the other exponentiation algorithm. It is very similar to the previous algorithm, but differs in processing the binary bits of the exponent in the opposite order. It also creates those bits as it goes, while the other assumes they are given. (In the HTML below, it is hard to tell the capital X and Y from the lower case, so I have put the capitals in red.)

Inputs: integers X, Y
Output: integer z
Algorithm:
   int exp(int X, int Y) {
      int x = X, y = Y, z = 1;
      while (y > 0) {
         while (y%2 == 0) {
            x = x*x;
            y = y/2;
         }
         z = z*x;
         y = y - 1;
      }
      return z;
   }

The loop invariant at each stage and after the each iteration of the inner while statement is:

z*xy = XY.

Here is a mathematical proof that the second algorithm actually calculates X^Y. Just before the while loop starts, since x = X, y = Y, and z = 1, it is obvious that the loop invariant is true. (In these equations, the = is a mathematical equals, not an assignment.)

Now suppose that at the start of one of the iterations of the while loop, the invariant holds. Use x', y', and z' for the new values of x, y, and z after executing the statements inside one iteration of the inner while statement. Notice that this assumes that y is even. Then the following are true:

x' = x*x
y' = y/2  (an exact integer without truncation because y is even)
z' = z
z'*x'y' = z*(x*x)y/2  = z*xy = XY.

This means that the loop invariant holds at the end of each iteration of the inner while statement for the new values of x, y, and z. Similarly, use the same prime notation for the variables at the end of the while loop.

x' = x
y' = y - 1
z' = z*x
z'*x'y' = z*x*(x)y - 1  = z*xy = XY.

So once again the loop invariant holds. After the loop terminates, the variable y must be 0, so that the loop invariant equation says:

XY = z*xy = z*x0 = z.

For a complete proof, one must also carefully argue that the loop will always terminate.

Here is a test of the two exponentiation functions implemented in Java: Exponentiation functions in Java.

---

Checking that an Integer is Probably a Prime Number.

For 2500 years mathematicians studied prime numbers just because they were interesting, without any idea they would have practical applications. Where do prime numbers come up in the real world? Well, there's always the 7-Up soft drink, and there are sometimes a prime number of ball bearings arranged in a circle, to cut down on periodic wear. Now finally, in cryptography, prime numbers have come into their own.

Law PRIME1: A source of large random prime integers is an essential part of many current cryptosystems.

Usually large random primes are created (or found) by starting with a random integer n, and checking each successive integer after that point to see if it is prime. The present situation is interesting: there are reasonabe algorithms to check that a large integer is prime, but these algorithms are not very efficient. On the other hand, it is very quick to check that an integer is ``probably'' prime. It is not very satisfactory to know that an integer is only probably prime, but if the chances of making a mistake about the number being a prime are reduced to a quantity close enough to zero, the users can just discount the chances of such a mistake.

Tests to check if a number is probably prime are called pseudo-prime tests. A number of such tests are available, but most use mathematical overkill. Anyway, one starts with a property of a prime number, such as Fermat's Theorem, mentioned in the previous chapter: if p is a prime and a is any non-zero number less than p, then a^p-1 mod p = 1. If one can find a number a for which Fermat's Theorem does not hold, then the number p in the theorem is definitely not a prime. If the theorem holds, then p is called a pseudo-prime with repect to a, and it might actually be a prime.

So the simplest possible pseudo-prime test would just take a small value of a, say 2 or 3, and check if Fermat's Theorem is true.

Simple Pseudo-prime Test: If a very large random integer p (100 decimal digits or more) is not divisible by a small prime, and if 3^p-1 mod p = 1, then the number is prime except for a vanishingly small probability, which one can ignore.

One could just repeat the test for other integers besides 3 as the base, but unfortunately there are non-primes (called Carmichael numbers) that satisfy Fermat's theorem for all values of a even though they are not prime. The smallest such number is ??? = ??*??, but these numbers become very rare in the larger range, such as 100 digit numbers. Corman et al. estimate that the chances of a mistake with just the above simple test are less than ???, although in practice commercial cryptosystems use better tests for which there is a proof of the low probability. Such better tests are not really needed, since even if the almost inconceivable happened and a mistake were made, the cryptosystem wouldn't work, and one could just choose another pseudo-prime.

Law PRIME2: Just one simple pseudo-prime test is enough to test that a very large random integer is probably prime.