fp_invmod.c

00001 /* TomsFastMath, a fast ISO C bignum library.
00002  * 
00003  * This project is meant to fill in where LibTomMath
00004  * falls short.  That is speed ;-)
00005  *
00006  * This project is public domain and free for all purposes.
00007  * 
00008  * Tom St Denis, tomstdenis@gmail.com
00009  */
00010 #include <tfm.h>
00011 
00012 static int fp_invmod_slow (fp_int * a, fp_int * b, fp_int * c)
00013 {
00014   fp_int  x, y, u, v, A, B, C, D;
00015   int     res;
00016 
00017   /* b cannot be negative */
00018   if (b->sign == FP_NEG || fp_iszero(b) == 1) {
00019     return FP_VAL;
00020   }
00021 
00022   /* init temps */
00023   fp_init(&x);    fp_init(&y);
00024   fp_init(&u);    fp_init(&v);
00025   fp_init(&A);    fp_init(&B);
00026   fp_init(&C);    fp_init(&D);
00027 
00028   /* x = a, y = b */
00029   if ((res = fp_mod(a, b, &x)) != FP_OKAY) {
00030       return res;
00031   }
00032   fp_copy(b, &y);
00033 
00034   /* 2. [modified] if x,y are both even then return an error! */
00035   if (fp_iseven (&x) == 1 && fp_iseven (&y) == 1) {
00036     return FP_VAL;
00037   }
00038 
00039   /* 3. u=x, v=y, A=1, B=0, C=0,D=1 */
00040   fp_copy (&x, &u);
00041   fp_copy (&y, &v);
00042   fp_set (&A, 1);
00043   fp_set (&D, 1);
00044 
00045 top:
00046   /* 4.  while u is even do */
00047   while (fp_iseven (&u) == 1) {
00048     /* 4.1 u = u/2 */
00049     fp_div_2 (&u, &u);
00050 
00051     /* 4.2 if A or B is odd then */
00052     if (fp_isodd (&A) == 1 || fp_isodd (&B) == 1) {
00053       /* A = (A+y)/2, B = (B-x)/2 */
00054       fp_add (&A, &y, &A);
00055       fp_sub (&B, &x, &B);
00056     }
00057     /* A = A/2, B = B/2 */
00058     fp_div_2 (&A, &A);
00059     fp_div_2 (&B, &B);
00060   }
00061 
00062   /* 5.  while v is even do */
00063   while (fp_iseven (&v) == 1) {
00064     /* 5.1 v = v/2 */
00065     fp_div_2 (&v, &v);
00066 
00067     /* 5.2 if C or D is odd then */
00068     if (fp_isodd (&C) == 1 || fp_isodd (&D) == 1) {
00069       /* C = (C+y)/2, D = (D-x)/2 */
00070       fp_add (&C, &y, &C);
00071       fp_sub (&D, &x, &D);
00072     }
00073     /* C = C/2, D = D/2 */
00074     fp_div_2 (&C, &C);
00075     fp_div_2 (&D, &D);
00076   }
00077 
00078   /* 6.  if u >= v then */
00079   if (fp_cmp (&u, &v) != FP_LT) {
00080     /* u = u - v, A = A - C, B = B - D */
00081     fp_sub (&u, &v, &u);
00082     fp_sub (&A, &C, &A);
00083     fp_sub (&B, &D, &B);
00084   } else {
00085     /* v - v - u, C = C - A, D = D - B */
00086     fp_sub (&v, &u, &v);
00087     fp_sub (&C, &A, &C);
00088     fp_sub (&D, &B, &D);
00089   }
00090 
00091   /* if not zero goto step 4 */
00092   if (fp_iszero (&u) == 0)
00093     goto top;
00094 
00095   /* now a = C, b = D, gcd == g*v */
00096 
00097   /* if v != 1 then there is no inverse */
00098   if (fp_cmp_d (&v, 1) != FP_EQ) {
00099     return FP_VAL;
00100   }
00101 
00102   /* if its too low */
00103   while (fp_cmp_d(&C, 0) == FP_LT) {
00104       fp_add(&C, b, &C);
00105   }
00106   
00107   /* too big */
00108   while (fp_cmp_mag(&C, b) != FP_LT) {
00109       fp_sub(&C, b, &C);
00110   }
00111   
00112   /* C is now the inverse */
00113   fp_copy(&C, c);
00114   return FP_OKAY;
00115 }
00116 
00117 /* c = 1/a (mod b) for odd b only */
00118 int fp_invmod(fp_int *a, fp_int *b, fp_int *c)
00119 {
00120   fp_int  x, y, u, v, B, D;
00121   int     neg;
00122 
00123   /* 2. [modified] b must be odd   */
00124   if (fp_iseven (b) == FP_YES) {
00125     return fp_invmod_slow(a,b,c);
00126   }
00127 
00128   /* init all our temps */
00129   fp_init(&x);  fp_init(&y);
00130   fp_init(&u);  fp_init(&v);
00131   fp_init(&B);  fp_init(&D);
00132 
00133   /* x == modulus, y == value to invert */
00134   fp_copy(b, &x);
00135 
00136   /* we need y = |a| */
00137   fp_abs(a, &y);
00138 
00139   /* 3. u=x, v=y, A=1, B=0, C=0,D=1 */
00140   fp_copy(&x, &u);
00141   fp_copy(&y, &v);
00142   fp_set (&D, 1);
00143 
00144 top:
00145   /* 4.  while u is even do */
00146   while (fp_iseven (&u) == FP_YES) {
00147     /* 4.1 u = u/2 */
00148     fp_div_2 (&u, &u);
00149 
00150     /* 4.2 if B is odd then */
00151     if (fp_isodd (&B) == FP_YES) {
00152       fp_sub (&B, &x, &B);
00153     }
00154     /* B = B/2 */
00155     fp_div_2 (&B, &B);
00156   }
00157 
00158   /* 5.  while v is even do */
00159   while (fp_iseven (&v) == FP_YES) {
00160     /* 5.1 v = v/2 */
00161     fp_div_2 (&v, &v);
00162 
00163     /* 5.2 if D is odd then */
00164     if (fp_isodd (&D) == FP_YES) {
00165       /* D = (D-x)/2 */
00166       fp_sub (&D, &x, &D);
00167     }
00168     /* D = D/2 */
00169     fp_div_2 (&D, &D);
00170   }
00171 
00172   /* 6.  if u >= v then */
00173   if (fp_cmp (&u, &v) != FP_LT) {
00174     /* u = u - v, B = B - D */
00175     fp_sub (&u, &v, &u);
00176     fp_sub (&B, &D, &B);
00177   } else {
00178     /* v - v - u, D = D - B */
00179     fp_sub (&v, &u, &v);
00180     fp_sub (&D, &B, &D);
00181   }
00182 
00183   /* if not zero goto step 4 */
00184   if (fp_iszero (&u) == FP_NO) {
00185     goto top;
00186   }
00187 
00188   /* now a = C, b = D, gcd == g*v */
00189 
00190   /* if v != 1 then there is no inverse */
00191   if (fp_cmp_d (&v, 1) != FP_EQ) {
00192     return FP_VAL;
00193   }
00194 
00195   /* b is now the inverse */
00196   neg = a->sign;
00197   while (D.sign == FP_NEG) {
00198     fp_add (&D, b, &D);
00199   }
00200   fp_copy (&D, c);
00201   c->sign = neg;
00202   return FP_OKAY;
00203 }
00204 
00205 /* $Source: /cvs/libtom/tomsfastmath/src/numtheory/fp_invmod.c,v $ */
00206 /* $Revision: 1.1 $ */
00207 /* $Date: 2007/01/24 21:25:19 $ */