Abstract:
A method for calculating greatest common divisors uses an approximate division in its reduction step. The result of this approximate division is then compared to determine if it is valid. If not, then the method applies a correction to the first approximate division to determine corrected values that have a reduced number of bits. If, during this correction step, the result is again not valid, then another method is applied to reduce the number of bits in the values. The approximate division is applied only when the number of significant bits in the two values differ by at least a predetermined number. When the number of bits in the two values differ by less than this number, an alternative GCD algorithm is applied but only to reduce the number of bits in the intermediate values.

Description:
BACKGROUND OF THE INVENTION 
     The present invention is directed to greatest common divisor algorithms and in particular to such an algorithm that uses multi-precision arithmetic and iteratively reduces the number of bits in each of the values it is testing either by performing an approximate division reduction step or by performing a novel inter-reduction step. 
     Modular arithmetic is used in many applications such as encryption algorithms and the generation of error correcting codes. When the modulus is a prime number, the operations of addition, subtraction, multiplication and division as well as the associative, commutative and distributive properties are defined over the set of numbers from zero to p. This set of numbers defines a Finite field modulo p, Fp. 
     GCD algorithms are commonly used to screen candidate large prime numbers. If a number passes the GCD test a more rigorous test is applied to determine if the number is prime. As set forth above, large prime number are important because they define Finite fields which are of interest for encryption purposes. One type of encryption algorithm encrypts data using exponentiation, relying on the inherent difficulty of the reverse of exponentiation, the discrete logarithm problem, to hold the data secure. Because encryption performed on a large Finite field is more secure than encryption performed on a small field, these algorithms work best when the exponentiation operations are carried out over a large Finite field. Thus, there is a need to identify large prime numbers. One problem with using large fields, however, is the size of the numbers being processed. Typical numbers used in data encryption have hundreds of bits. These numbers are too large to be easily handled by commonly available microprocessors that are limited to 32 or 64-bit arithmetic. As described below, calculations used in the encryption process are typically handled using multiprecision arithmetic. 
     FIG. 1 is a flow chart diagram of an algorithm that screens candidate prime numbers. The GCD algorithm is relatively quick and can find weed-out many candidate prime numbers having relatively small factors. The algorithm shown in FIG. 1 tests a number P to determine if it is prime relative to the first N thousand prime numbers where N is an integer. Prior to performing the algorithm shown in FIG. 1, an array THOUSAND_PRIMES [I] is prepared such that THOUSAND_PRIMES [ 1 ] contains the product of the first thousand prime numbers, THOUSAND_PRIMES [ 2 ] contains the product of the second thousand prime numbers and so on. The process begins at step  102  by obtaining the number P which is to be tested and the number N defining the number of thousand primes that P is to be tested against. At step  104  an index variable I is set to zero. At step  106 , I is incremented by one and a value G is calculated where G is the greatest common divisor of P and THOUSAND_ PRIMES [I]. If, at step  108 , G is greater than one then P is divisible by one of the factors of THOUSAND_PRIMES [I] and is not prime. Accordingly, control transfers to step  110  where the value FALSE is returned. 
     If, however, at step  108  G is equal to one then P is not divisible by any of the factors of THOUSAND_PRIMES [I] and control transfers to step  112 . At step  112  the variable I is compared against N. If I is less than N control transfers to step  106  where I is incremented and the greatest common divisor of P and THOUSAND_PRIMES [I] is calculated for the incremented value of I. When, at step  112 , I is equal to N, the number P has been tested against the N thousand prime numbers and has been found to be relatively prime to each of these products. Thus, at step  114 , the algorithm shown in FIG. 1 returns a value TRUE indicating that P is a good candidate to be a prime number. After performing the test shown in FIG. 1, the candidate prime number may be processed using a probabilistic primality testing routine, such as the Miller-Rabin algorithm to determine if it is prime. 
     A key element of the program shown in FIG. 1 is the greatest common divisor (GCD) algorithm. FIG. 2 is a flow chart diagram which illustrates a greatest common divisor algorithm known as Euclid&#39;s algorithm. The algorithm shown in FIG. 2 calculates the greatest common divisor of U and V where U is greater than V. The algorithm shown in FIG. 2 also calculates the inverse of V modulo U. This algorithm relies on the property that if U and V have a common divisor D so does U-V, U- 2 V and so on. Thus, using only subtraction, one can calculate the GCD of U and V. GCD algorithms run faster when they can combine multiple subtractions in a single step. When calculating the inverse, The GCD algorithm operates by maintaining the linear equations (1), (2) and (3) 
     
       
           U   1   P+U   2   X=U   3   (1) 
       
     
     
       
           V   1   P+V   2   X=V   3   (2) 
       
     
     
       
           T   1   P+T   2   X=T   3   (3) 
       
     
     Where U≧V and U 3  and V 3  are initially assigned the values of U and V, respectively. Upon termination, U 3 =GCD(P, X)=1 and U 2 =X −1  MOD P. In general terms, GCD algorithms operate by repetitively reducing the number of bits in the larger value, U 3 , and switching the two values U 3  and V 3 . Thus, the algorithm successively reduces the values of U 3  and V 3  while maintaining the linear equations. Because it maintains the values U 2  and V 2 , the algorithm shown in FIG. 2 not only calculates the greatest common divisor of U and V but also calculates the inverse of V modulo U. That is to say, V −1  where V*V −1 =1 modulo U. As described below, when U is a prime number, this inverse may be used to simplify division calculations performed within the Finite fields Fu. Furthermore, the variables U 1 , V 1  and T 1  do not need to be maintained because they can be determined from the other variables, for example, U 1  can be determined from U 2  and U 3  by the identity U 1 =(U 3 −U 2 V)/U. 
     The Euclid GCD algorithm shown in FIG. 2 begins at step  210  by obtaining the values U and V. Next, at step  212  a temporary variable U 3  is set equal to U and a temporary variable V 3  is set equal to V. Also in step  212  variable U 2  is set to zero and V 2  is set one. Next, at step  214  the variable V 3  is tested to determine if it is equal to zero. If not, step  218  is executed. This step calculates a value Q which is equal to the greatest integer less than U 3 /V 3 . Next, it determines a new value for V 3  by storing the quantity, U 3 −Q*V 3  into a temporary variable T 3 . Also at step  218 , the process stores the value U 2 −Q*V 2  into a temporary variable T 2 . This value is the new value for the variable V 2 . The value in T 3  represents a reduced value of U 3 . Because V 3  has not been reduced, however, T 3  is less than V 3 . Consequently, the value in V 3  is assigned to U 3  and the value in T 3  is assigned to V 3 . A similar operation is performed on the variables U 2 , V 2  and T 2 . After step  218 , control returns to step  214  where V 3  is once again compared to determine if it is equal to zero. When V 3  is equal to zero the process returns the values U 3  and U 2 . The value U 3  is the greatest common divisor of U and V while the value U 2  is the inverse of V modulo U (i.e. V −1  MOD U). 
     As described above, a modular arithmetic operation that is used in many encryption algorithms is raising a large number M, representing a message, to a power E modulo P. In a typical encryption example the values M, E and P are each hundreds of bits in length. Operations on numbers of this size—and especially division operations—are difficult because common computer processing hardware supports at most sixty-four bit mathematical operations. Values greater than 64 bits in length are typically calculated using multiprecision arithmetic. A number, a, raised to an exponent e can always be calculated by multiplying that number by itself the number of time represented by the exponent, or in mathematical terms: 
     
       
           a   e   =a*a*a . . . e  number of times. 
       
     
     Another method, which is significantly faster, is the multiply chain algorithm. In this case, let e=e n−1 e n−2  . . . e 1 e 0  be an n-bit exponent e i ε{0, 1}, 0≦i≦n−1 and e n−1 =1. The algorithm starts with p 1 =a, then 
     
       
           P   i+1   =P   i   2  if  e   n−1-i =0 or  a*p   i   2  if  e   n−1-i =1, where 1 ≦i≦n− 2. 
       
     
     Several methods are known in the art to reduce either the number of multiplications used to produce efficient exponentiation of the base value. 
     One method known to reduce the number of multiplication is the signed digit algorithm. Using this method, the exponent is represented as a string of bits comprising the values 0 and 1. Within the bit string, sequences (or “runs”) of 1&#39;s are replaced by 0&#39;s, with a 1 being placed in the next higher bit position to the most significant bit (MSB) position of the run, and “−1” being inserted in the least significant bit (LSB) position of the run. By thus efficiently recoding the exponent bit string, the expected number of multiplication operations is reduced from n/2 to n/3. 
     Although the signed digit algorithm reduces the number of multiplication operations needed to calculated the exponential value M e  modulo P, it replaces some of the multiplication operations by division operations. As described above, division may be implemented as the multiplication by the inverse. 
     FIG. 3 is a flow chart diagram which illustrates an exemplary method for calculating M E  modulo P. The first step in this algorithm, step  310  gets the values P, M and E. P is a large prime number, M is the message to be encoded and E is an exponent value to be used in order to encode the message. The next step in the process is to determine the inverse of M modulo P (i.e. M −1  MOD P). This is accomplished by calculating the greatest common divisor of P and M. As P is a prime number by definition, the value U 2  returned by the GCD algorithm is equal to M −1  MOD P. Next, at step  314 , the signed digit algorithm is applied to the exponent E. Finally, at step  316 , the multiply-square algorithm is applied to calculate M E  modulo P. As described above, the multiply-square algorithm both multiplies and divides by M. The value M −1  MOD P is used to convert the division operation into a multiplication operation. This follows from the following equations (a/b*b −1 /b −1 ) mod p=(a*b −1 )/(b*b −1 ) mod p=(a*b −1 )/1 mod p. 
     The Euclid algorithm is not efficient for large numbers that require multiprecision arithmetic. Thus, there is a need for an efficient GCD algorithm that calculates inverses, and operates efficiently on numbers that require multiprecision arithmetic. 
     SUMMARY OF THE INVENTION 
     The present invention is embodied in a method for calculating greatest common divisors that uses an approximate division in its reduction step. The result of this approximate division is then compared to determine if the result is negative. If it is, then an alternate method is used to reduce the number of bits in the values while maintaining the linear formulas. 
     According to another aspect of the invention, if the first approximate division produces a negative result, the method applies a correction to the first approximate division to determine corrected values that have a reduced number of bits. If, during this correction step, the result is still negative, then yet another method is applied to reduce the number of bits in the values. 
     According to another aspect of the invention, the approximate division is applied only when the number of bits in the two values differ by at least a predetermined value. When the number of bits in the two values differ by less than this number, an alternative inter-reduction step is applied. 
    
    
     BRIEF DESCRIPTION OF THE DRAWINGS 
     FIG. 1 (Prior Art) is a flow chart diagram that is useful for describing an algorithm that screens a number to determine if it is not prime in order to identify candidate numbers for more rigorous prime testing. 
     FIG. 2 (Prior Art) is a flow chart diagram of the Euclid GCD algorithm. 
     FIG. 3 (Prior Art) is a flow chart diagram which is useful for describing a method for calculating a value M e  modulo P. 
     FIG. 4 is a flow chart diagram which illustrates an exemplary GCD algorithm according to the present invention. 
     FIG. 4A is a flow-chart diagram of a TOP 2 M routine suitable for use in the algorithm shown in FIG.  4 . 
     FIG. 5 is a flow-chart diagram of a CALCULATE-XY process that is used in the CALCULATE XA+YB&lt;2 M  step of FIG.  4 . 
     FIG. 6 is a flow-chart diagram of a modified Lehmer GCD algorithm which may be used by the CALCULATE-XY process shown in FIG.  5 . 
     FIG. 7 is a flow-chart diagram of a modified Euclid GCD algorithm that may be used by the CALCULATE-XY process shown in FIG.  5 . 
     FIG. 8 is a flow chart diagram which illustrates the GCD algorithm shown in FIG. 4 including an inverse calculation. 
    
    
     DETAILED DESCRIPTION 
     The present invention is embodied in a greatest common divisor (GCD) method that is fast even when multi-precision arithmetic is used and is also efficient at computing inverses. 
     Some GCD algorithms, including the algorithm disclosed by T. Jebelean, in an article entitled “A Generalization of the Binary GCD Algorithm,”  ACM Symposium on Symbolic and Algebraic Computation  ( ISSAC ) July, 1993, pp 111-116 gain efficiency by calculating factors A and B that zero the N low-order bits of U 3  and V 3  when U 3  and V 3  are multiplied by these factors. The bit-reduced value for U 3  is then calculated as (AU 3 +BV 3 )/2 N , where division by 2 N  is implemented as a right-shift by N bits. When U 3  and V 3  are multiplied by A and B, respectively, the number of bits in the new value of U 3  increases but by less than N bits. Thus, the bit shift produces a net reduction in the number of bits. One problem with algorithms of this type is that they calculate GCD (XU 3 , V 3 ) instead of GCD (U 3 , V 3 ). Consequently, when V 3  equals zero, the returned value of the GCD may be greater than 1 even though U and V are relatively prime. This value represents the value 1 multiplied by spurious factors. 
     A second problem with algorithms of this type concerns the calculation of the inverse. The Jebelean algorithm uses a multi-precision addition/subtraction and shift operation to maintain the value of V 2  for each iteration in which V 2  is odd, greatly increasing the execution time of the algorithm. This problem requires a relatively expensive fix: to maintain equation (3) one would initially calculate T 1 U+T 2 V=T 3 /K X , where it is expensive to correct T 2  when accounting for the K X  term. In the Jebelean algorithm, K equals 2. To maintain equation (2) one must transform U′ 1 U+U′ 2 V=GCD′ into U 1 U+U 2 V=GCD where GCD|GCD′ but A=(GCD′/GCD) need not divide U′ 1  or U′ 2 . Thus, the calculation of the inverse is not a straightforward process using the Jebelean GCD algorithm. 
     The algorithm according to the present invention shares the first problem with the Jebelean algorithm but avoids the second problem by working with the top 2M bits of U and V. The reduction step employs approximate division if the difference in the number of bits in the values U and V is greater than a threshold and uses a GCD algorithm to reduce the number of bits in the top 2M bits of U and V by calculating factors A and B such that A(TOP 2 M(U))+B(TOP 2 M(V))&lt;2 M , if the difference in the number of bits in the values U and V is less than the threshold. The inventor has determined the threshold value experimentally. This process, which uses division instead of shifting, removes the K X  term and, thus, eliminates the second problem. 
     FIG. 4 is a flow-chart diagram of a GCD algorithm that includes an embodiment of present invention. The algorithm begins at step  410  when it receives the values U, V and M. The next step,  412  sets the temporary variables U 3  and V 3  equal to U and V, respectively. At step  414 , the method tests the value of M. In the exemplary embodiment of the invention, M may have a value of 16 (indicating that method uses 32-bit arithmetic) or 32 (indicating that the method uses 64-bit arithmetic). Responsive to M having a value of 32 or 16, steps  416  and  418  set a temporary variable D to 8 or 5, respectively. D is the threshold value that determines which reduction step is used. The inventor has determined that, when using 32 (64) bit arithmetic and the number of bits in the values U 3  and V 3  differs by 5 (8) or more, it is more efficient to use approximate division in the reduction step. When this difference is less than 5 (8), however, it is more efficient to use a modified GCD inter-reduction step in the algorithm. 
     At step  419 , the process sets temporary variables A and B to be the top 2M bits of U 3  and V 3 , respectively by invoking the routine TOP 2 M. This routine is described below with reference to FIG.  4 A. This operation is done before branching between the reduction step and the inter-reduction step because the 2M most significant bits of U 3  and V 3  are used in both branches. 
     Referring to FIG. 4A, the first step in the routine TOP 2 M, step  450 , gets the values of U 3 , V 3  and M. Next, at step  452 , the process sets a temporary variable PB 1  to the bit position of the most significant non-zero bit in the value U 3 . Next, at step  454 , the process sets a temporary variable BP 2  to BP 1  minus 2M plus one. Finally, at step  456 , the routine returns U 3  [BP 1 :BP 2 ], the top 2M bits of U 3 , and V 3  [BP 1 :BP 2 ], the corresponding bits in the value V 3 . This routine assumes that U 3  is greater than V 3 . 
     Referring again to FIG. 4, after step  419 , the process determines, at step  420 , if the bit position of the most significant non-zero bit in V 3  differs from the position of the most significant non-zero bit in U 3  by D or less. If it does not, then Step  434  performs an inter-reduction step by using the COMPUTE-XY routine to calculate values X and Y such that XA+YB&lt;2 M . As described below, with reference to FIGS. 5,  6  and  7 , this computation employs one of two modified GCD algorithms. Each of the GCD algorithms is modified to terminate when the GCD value is less than 2 M . After the values X and Y have been obtained, the algorithm, at step  436 , the algorithm calculates the new value for V 3  by setting a temporary variable T 3  to (XU 3 +YV 3 ), setting U 3  equal to V 3  and then setting V 3  equal to the absolute value of T 3 . 
     If, at step  420 , the number of bits in U 3  and V 3  differ by more than D, the GCD algorithm according to the subject invention, at step  422 , “approximately” divides U 3  and V 3  using only the top 2M bits of U 3  (A) and V 3  (B) to generate an estimated value Q 1  that is approximately equal to the greatest integer less than U 3  divided by V 3 . Also at step  422 , the variable T 3  is set to (U 3 −Q 1 *V 3 ). The inventor has determined that, very often, the value Q 1  is the same as would be obtained by fully calculating └U 3 /V 3 ┘ using all of the bits of these values. After step  422 , the next step in the algorithm, step  424 , tests T 3  to determine if it is less than zero. If it is, then a correction step is implemented because T 3  needs to be positive (i.e. Q 1  does not equal └U 3 /V 3 ┘). 
     In this instance, the algorithm performs a correction step by first determining, at step  425 , if the absolute value of T 3  is greater than V 3 . If it is, then at step  426 , the process calls TOP 2 M to store the top 2M bits of the absolute value of T 3  and the corresponding 2M bits of V 3  into the respective temporary variables A and B. Also at step  426 , the process calculates a new value Q 2  which is equal to Q 1 −(└A/B)┘+1). 
     If, at step  425 , the absolute value of T 3  is not greater than V 3 , then, at step  427 , the value Q 2  is calculated by simply subtracting one from Q 1  at step  427 . This correction results from the correction performed in step  426  because, if the absolute value of T 3  is less than V 3 , then └|T 3 |/V 3 ┘ is zero. 
     After step  426  or  427 , step  428  computes T 3  as U 3 −Q 2 *V 3  and, at step  429 , determines if the newly calculated value of T 3  is less than zero. Because the difference in the number of bits in U 3  and the number of bits in V 3  is almost always less than 16 when M=16, and less than 32 when M=32. The inventor has determined empirically that when Q 1  does not equal └U 3 /V 3 ┘, T 3  is greater than zero and, thus, Q 2  equals └U 3 /V 3 ┘ over 99% of the time. 
     For the rare instances, at step  429 , when T 3  is still negative, the algorithm, at step  430  sets the variable Q equal to └U 3 /V 3 ┘, calculated using all of the bits in U 3  and V 3 , and then sets T 3  equal to U 3  minus the quantity Q times V 3 . In the last step,  432 , of this branch of the reduction step, the algorithm sets U 3  equal to V 3  and V 3  equal to T 3 . 
     After step  436  or  432 , the algorithm tests V 3  to determine if it is equal to zero. If it is not, control returns to step  419  to extract the most significant non-zero 2M bits of the new values of U 3  and V 3  and perform the next reduction step or inter-reduction step, as described above. 
     When, at step  438 , V 3  is equal zero, step  440  is executed to determine if U 3  is greater than 1. If U 3  is greater than 1, then the result may include spurious factors. The algorithm removes any factors at step  444  by setting a temporary variable Y equal to GCD (U 3 , U mod U 3 ) and then setting Y equal to the GCD (Y, V mod Y). The GCD algorithms used in this step may be any GCD algorithm because the inventors have determined that, most of the time, the number of bits in U 3  is less than 2M, thus, these GCD calculations only rarely need to use multi-precision arithmetic. One exemplary GCD algorithm that may be used is the Euclid algorithm described in the above with reference to FIG.  2 . 
     If, at step  440 , U 3  is found to be equal to 1, then, at step  442 , the temporary variable Y is set to 1. At step  446 , after step  442  or step  444 , the algorithm returns the value Y as the result of the GCD calculation. 
     In step  434 , the calculation of X and Y such that XA+YB&lt;2 M  is performed using the routine COMPUTE-XY, shown in FIG.  5 . As described in more detail below, the routine COMPUTE-XY invokes one of two GCD algorithms. If M equals 16, a modified version of the Euclid algorithm, shown in FIG. 7, is used. If M equals 32, a modified Lehmer algorithm, shown in FIG. 6, is used. 
     The inventor has determined that using the Lehmer algorithm when M equals 32 has a significant effect on the speed at which the GCD and inverse (if the inventive GCD routine computes inverses) is calculated. The speedup, relative to using the Euclid algorithm, varies from 20-39% when U and V are 160 bit numbers and roughly 16-26% when U and V are  256  bit numbers. 
     Referring to FIG. 5, the routine COMPUTE-XY begins at step  510  by obtaining the values A, B and M. As set forth above, this routine determines values for variables X and Y such that XA+YB&lt;2 M . In the exemplary embodiment of the invention, M may be either 16 or 32. The next step in the process, step  512 , determines if M is equal to 32. If it is then, at step  514  the process calls the modified Lehmer GCD routine described below with reference to FIG. 6 to obtain values for X and Y. If M is equal to 16 at step  512 , then the process, at step  516 , calls the modified Euclid GCD routine to obtain the values of X and Y. This routine is described below with reference to FIG.  7 . After step  514  or step  516 , the process returns the values of X and Y at step  518 . 
     The modified Lehmer GCD algorithm is derived from the basic Lehmer-Euclid algorithm, which is described in an article by J. Sorenson entitled “An Analysis of Lehmer&#39;s Euclidean GCD Algorithm,”  ACM International Symposium on Symbolic and Algebraic Computation  ( ISSAC ), July 1995, pp 254-258. The algorithm described in the article calculates the GCD of U and V. It is modified for the subject invention to terminate as soon as U 3  is less than 2 M . FIG. 6 is a flow-chart diagram of this modified GCD algorithm. 
     At step  610  of FIG. 6, the process obtains the values U and V. Because the process shown in FIG. 6 is invoked from the CALCULATE X-Y routine (shown in FIG. 5) which is, itself invoked from the inventive GCD algorithm (shown in FIG.  4 ), the values U and V are the top 2M bits of the values being processed by the algorithm shown in FIG.  4 . The modified Lehmer algorithm is used only when M equals 32 and, so, 2M equals 64. For the sake of consistency, the variable names used in the modified Lehmer algorithm are the same as those used in the Euclid algorithm (described above with reference to FIG. 2) and the same as those used in the inventive algorithm (described above with reference to FIG.  4 ). These variables, however, are local to the Lehmer algorithm; they are not the same as the variables of the same name that are used in the algorithm shown in FIG.  4 . 
     After step  610 . the next step in the modified Lehmer GCD process, step  612 , assigns the values U and V to temporary variables U 3  and V 3 , respectively. Also at step  612 , the temporary variables U 2  and V 2  are initialized to zero and one, respectively. At step  614 , the routine determines if the most significant non-zero bit in the value U 3  differs from the most significant non-zero bit in the value V 3  by more than 32. If this condition is true, control transfers to step  630 , described below. Otherwise, control transfers to step  616 . At step  616 , the process initializes an inter-reduction step used by the Lehmer algorithm. This step sets temporary variables SU and SV to the top  32  bits of the current values of U 3  and V 3 , respectively; sets a temporary variable D to the bit position of the most significant non-zero bit in U 3  minus  32  (i.e. log 2  (U 3 )−32); sets an index variable J to zero, a condition variable DONE to false; and initializes the first two entries in the arrays X and Y, setting X[ 0 ] and Y[ 1 ] to one and X[ 1 ] and Y[ 0 ] to zero. 
     At step  618 , the process calculates a value Q as the greatest integer less than SU divided by SV. Next, step  618  assigns the value X[J]−X[J+1] to the array element X[J+2] and assigns the value Y[J]−Y[J−1] to the array element Y[J+1]. Also at step  618 , the routine reduces SU by the inter-reduction number of bits (assigned to a temporary variable T) and switches the values of SU and SV. Finally at step  618 , the routine increments the index variable J by one. 
     At step  620 , the routine determines if J is even or odd. If it is even then the routine, at step  622 , calculates a value for the Boolean variable DONE as (SU&lt;−X[J+1]) or (SU−SV&lt;Y[J−1]−Y[J]) or (SU&lt;2 D ). This last condition, (SU&lt;2 D ), is not a part of the Lehmer Euclid algorithm but is added to terminate the inter-reduction step when U 3  is less than 2 M . If J is even at step  620  then the process, at step  624 , calculates a value for the Boolean variable DONE as (SV&lt;−Y[J+1]) or (SU−SV&lt;X[J−1]−X[J]) or (SU&lt;2 D ). After either step  622  or  624 , the process tests the Boolean variable DONE. If it is true then control transfers to step  628 , described below. If it is false control transfers to step  618  to perform another inter-reduction step. 
     Step  628  calculates updated values for U 3  and U 2 , according to equations (7) and (8) assigning the values to temporary variables T 3  and T 2 . 
     
       
           T   3 = X[J −1 ]*U   3 + Y[J −1 ]*V   3   (7) 
       
     
     
       
           T   2 = X[J −1 ]*U   2 + Y[J −1 ]*V   2   (8) 
       
     
     The routine also calculates new values for V 3  and V 2  according to equations (9) and (10). 
     
       
           V   3 = X[J]*U   3 + Y[J]*V   3   (9) 
       
     
     
       
           V   2 = X[J]*U   2 + Y[J]*V   2   (10) 
       
     
     Finally, the values in the temporary variables T 3  and T 2  are assigned to U 3  and U 2 , respectively. 
     After step  628  or after step  614  if the difference in the number of bits in U 3  and V 3  is greater than 32, step  630  is executed. This step performs the same reduction step as the Euclid algorithm and then switches the values of U 3  and V 3  after reducing the number of bits in U 3 . Next, step  632  determines if U 3  is less than 2 32 . If it is, the routine is done and the values U 3  and U 2  are returned as the values X and Y. If U 3  is not less than 2 32 , then control transfers to step  614  to continue the inter-reduction steps until the test at step  632  is satisfied. When U 3  is less than 2 32  at step  632 , step  634  is executed which calculates a value for U 1  as (U 3 −U 2 *V)/U and returns the values U 1  and U 2  as the values X and Y. This is a relatively economical division as U exactly divides (U 3 −U 2 *V). 
     As set forth above, the modified Lehmer algorithm is used only when M equals 32 at step  512  of the CALCULATE-XY routine described above with reference to FIG.  5 . When M equals 16, the inventive GCD algorithm uses a modified version of the Euclid GCD algorithm. This modified Euclid algorithm is shown in FIG.  7 . Steps  710 ,  712  and  718  of this algorithm are identical to the respective steps  210 ,  212 ,  216  and  218 , described above with reference to FIG.  2 . The only differences between the two algorithms are in the termination condition at step  714  and the calculation of U 1  and return of the values U 1  and U 2  at step  716 . The modified Euclid GCD algorithm, terminates when U 3  is less than 2 16  rather than when V 3  equals zero as in the algorithm shown in FIG.  2 . 
     Thus, the GCD algorithm according to the present invention performs an approximate calculation using only the most significant 2M bits of U 3  and V 3  if the number of bits in U 3  and V 3  differ by more than a threshold value. If this condition is not met then the algorithm reduces the number of bits in U 3  by invoking a modified GCD algorithm to calculate values X and Y such that XU 3 +YV 3  is less than 2 M . This step uses a simple GCD routine, such as the Euclid routine (modified to terminate when U 3 &lt;2 M ) when M equals 16 or a more complex GCD routine, such as the Lehmer method (also modified to terminate when U 3 &lt;2 M ) when M equals 32. The inventors have determined that this algorithm produces good results on values having a number of bits greater than or equal to 96. 
     One advantage that the inventive GCD algorithm has over other GCD algorithms is the speed at which it calculates inverses. The inventive algorithm described above with reference to FIGS. 4-7 does not calculate inverses. It only calculates GCD&#39;s. FIG. 8 shows a version of the inventive algorithm that also calculates inverses. 
     In FIG. 8, the algorithm begins at step  810  when it receives the values U, V and M. The next step,  812  sets the temporary variables U 3  and V 3  equal to U and V respectively, and initializes the variables U 2  to zero and V 2  to one. At step  814 , the method tests the value of M. If M is 32, step  816  sets a temporary variable D to 8. If M is 16, step  818  sets D to 5. 
     At step  819 , the process sets temporary variables A and B to be the top 2M bits of U 3  and V 3 , respectively by invoking the routine TOP 2 M. This routine is described above with reference to FIG.  4 A. This operation is done before branching between the reduction step and the inter-reduction step because the 2M most significant bits of U 3  and V 3  are used in both branches. 
     If, at step  820 , the number of bits in U 3  and V 3  differ by more than D, the GCD algorithm according to the subject invention, at step  832 , “approximately” calculates └A/B┘ to generate an estimated value Q 1  that is approximately equal to └U 3 /V 3 ┘. Also at step  832 , the process calculates a new value for T 3  as U 3  minus Q 1  times V 3 . After step  832 , the next step in the algorithm, step  834 , tests T 3  to determine if it is less than zero. If it is, then a correction step is implemented because T 3  needs to be positive (i.e. Q 1  does not equal └U 3 /V 3 ┘). 
     In this instance, the algorithm performs a correction step by first determining, at step  837 , if the absolute value of T 3  is greater than V 3 . If it is, then at step  838 , the process calls TOP 2 M to store the top 2M non-zero bits of the absolute value of T 3  and the corresponding 2M bits of V 3  into the respective temporary variables A and B. Also at step  838 , the process calculates a new value Q 2  as Q 2 =Q 1 −(└A/B)┘+1). 
     If, at step  837 , the absolute value of T 3  is not greater than V 3 , then, at step  839 , the value Q 2  is calculated by simply subtracting one from Q 1  at step  839 . This correction is consistent with the correction performed in step  838 , because, if the absolute value of T 3  is less than V 3 , then └T 3 |/V 3 ┘ is zero. 
     After step  838  or  839 , step  840  computes T 3  as U 3 −Q 2 *V 3  and, at step  841 , determines if the newly calculated value of T 3  is less than zero. Because the difference in the number of bits in U 3  and the number of bits in V 3  is almost always less than 16 when M=16, and less than 32 when M=32, The inventor has determined empirically that when Q 1  does not equal └U 3 /V 3 ┘, T 3  is greater than zero and, thus, Q 2  equals └U 3 /V 3 ┘ over 99% of the time. 
     As in the algorithm shown in FIG.  4 . when Q 2  does not equal └U 3 /V 3 ┘, the algorithm, at step  844  sets the variable Q equal to └U 3 /V 3 ┘, calculated using all of the bits in U 3  and V 3 , and then uses the value Q to calculate a new value for U 3 , storing it into the temporary variable T 3 . After step  834  if T 3  is not less than zero then, at step  836 , the algorithm assigns the value in Q 1  to Q. Similarly, after step  841  if T 3  is not less than zero, the process sets Q equal to Q 2 . At step  846 , the algorithm completes the task of switching U 3  and V 3  and, at the same time, calculates a new value for U 2 , using Q, and switches the values of U 2  and V 2 . 
     When, at step  820 , the number of bits in V 3  differs from the number of bits in U 3  by D or less, the algorithm according to the subject invention, performs the inter-reduction step by computing X and Y such that XA+YB&lt;2 M . At step  824 , the algorithm then sets T 3 =XU 3 +YV 3  and sets T 2 =XU 2 +YV 2 . Next, at step  826 , the algorithm tests T 3  to determine if it is less than zero. If it is, the sign of T 2  is switched at step  828 . After step  826  or  828 , step  830  is executed. This step completes the switching of U 3  and V 3  and the switching of U 2  and V 2  by assigning V 3  to U 3 , assigning the absolute value of T 3  to V 3 , assigning V 2  to U 2  and assigning T 2  to V 2 . 
     After step  830  or  846 , the algorithm, at step  848  tests V 3  to determine if it is equal to zero. If it is not, control returns to step  820  to perform the next reduction step or inter-reduction step. When, at step  848 , V 3  is equal zero, step  850  is executed to determine if U 3  is greater than 1. If U 3  is greater than 1, then the result may include spurious factors. The algorithm removes these factors at step  854  by setting a temporary variable Y equal to GCD (U 3 , U MOD U 3 ), setting a temporary variable X equal to └U 3 /Y┘ and then calculating the inverse of X with respect to U by using a GCD algorithm. Also at step  854 , the process corrects U 2  by multiplying U 2  by the X −1  modulo U. The process also removes the spurious factors from the GCD value, currently stored in the variable Y by setting Y equal to the GCD (Y, V mod Y). The GCD algorithms used in this step may be any GCD algorithm, including the inventive GCD algorithm or the Euclid algorithm described above with reference to FIG. 2 because the inventors have determined that, most of the time, the number of bits in U 3  is less than 2M, thus, these GCD calculations only rarely need to use multi-precision arithmetic. 
     The use of the approximate division in the inventive GCD algorithm provides a significant advantage over other methods. The inventor has determined that adding this step to the inventive algorithm increases its performance by approximately 35 percent when operating on 160-bit numbers and by approximately 32 percent when operating on 256-bit numbers. 
     If, at step  850 , U 3  is found to be equal to 1, then no spurious factors exist and, at step  852 , the temporary variable Y is set to 1. At step  856 , after step  852  or step  854 , the algorithm returns the value Y as the result of the GCD calculation and the value U 2  as the inverse of V with respect to U. Note that when calculating the inverse it is desirable to calculate U 2  mod U because U 2  may be greater than U upon termination. This occurs because V 2  grows by log(Y)+log(Q) bits per iteration while U 3  and V 3  are reduced by approximately (M−1+log(Q)) bits every other iteration, where log(Y) is less than or equal to M−1. Consequently, log(U 2 )&lt;2*log(U). 
     Although the inventive method for identifying the GCD of two numbers has been described in terms of applying a single algorithm, it is contemplated that two or more GCD routines may be combined in a single calculation. One such combination would be to invoke the Lehmer-Euclid routine as soon as the number of bits in U 3  becomes less than or equal to 64. 
     The GCD routines may be further optimized by analyzing the processors that execute the routines. If, for example, an Intel Pentium II™ microprocessor is used, the algorithm may be optimized based on the methods used to perform 32 and 64 bit operations and on the type of variable that is used to hold the values. For example, although the Microsoft C++ compiler version 6.0 supports 64-bit multiplications, the Pentium II processor does not, only the lower 64 bits of the result are returned. This “feature” of the Pentium II processor may be exploited to perform a multiplication modulo 2 64 . Furthermore, the inventor has determined that dividing 64-bit integers is more efficient than dividing 32-bit integers for the amount of reduction gained. Division of double precision values, however, is even more efficient. 
     As described above, the inventive GCD algorithm may be applied in cryptography and integer factoring. The GCD algorithm without the inverse calculations may be applied to find the factors of a number or to help determine if a number is prime and, thus, to define a Finite field suitable for use in a data encryption operation. The GCD algorithm with the inverse calculation may be used to identify mathematical inverses of values in the Finite field in order to replace a relatively expensive multiprecision division operation by a simpler multiplication operation. 
     Multiprecision division is used, for example in the elliptic curve encryption algorithm. In this algorithm, a decoding key is identified by repeatedly adding a large number over a group of numbers that are defined by an elliptic curve in a Finite field. Because the addition is along the curve, it is not a simple process but requires a division step. In general, the addition of values in elliptic curve cryptography is defined by equations (11) through (14). All of these operations are modulo P where P is the large prime number that defines the Finite field. 
     
       
         ( x   1 ,  y   1 )+( x   2 ,  y   2 )=( x   3 ,  y   3 )  (11) 
       
     
     
       
           x   3 = L   2   −x   1 − x   2   (12) 
       
     
     
       
           y   3 = L ( x   1 − x   3 )− y   1   (13) 
       
     
     
       
           L =( y   1 − y   2 )/( x   1 − x   2 )  (14). 
       
     
     Thus, to find an encryption key in an elliptic curve cryptographic system, one would perform a division step, as shown in equation (14). This may be simplified according to the present invention by calculating (x 1 −x 2 ) and then invoking the inventive GCD routine to determine the inverse of the calculated value, modulo P. 
     Although the invention has been described in terms of exemplary embodiments, it is contemplated that it may be practiced with variations within the scope of the appended claims.