pd.if.org Git - zpackage/blob - tomsfastmath/src/numtheory/fp_invmod.c

   1 /* TomsFastMath, a fast ISO C bignum library.
   2  *
   3  * This project is meant to fill in where LibTomMath
   4  * falls short.  That is speed ;-)
   5  *
   6  * This project is public domain and free for all purposes.
   7  *
   8  * Tom St Denis, tomstdenis@gmail.com
   9  */
  10 #include <tfm_private.h>
  11
  12 static int fp_invmod_slow (fp_int * a, fp_int * b, fp_int * c)
  13 {
  14   fp_int  x, y, u, v, A, B, C, D;
  15   int     res;
  16
  17   /* b cannot be negative */
  18   if (b->sign == FP_NEG || fp_iszero(b) == 1) {
  19     return FP_VAL;
  20   }
  21
  22   /* init temps */
  23   fp_init(&x);    fp_init(&y);
  24   fp_init(&u);    fp_init(&v);
  25   fp_init(&A);    fp_init(&B);
  26   fp_init(&C);    fp_init(&D);
  27
  28   /* x = a, y = b */
  29   if ((res = fp_mod(a, b, &x)) != FP_OKAY) {
  30       return res;
  31   }
  32   fp_copy(b, &y);
  33
  34   /* 2. [modified] if x,y are both even then return an error! */
  35   if (fp_iseven (&x) == 1 && fp_iseven (&y) == 1) {
  36     return FP_VAL;
  37   }
  38
  39   /* 3. u=x, v=y, A=1, B=0, C=0,D=1 */
  40   fp_copy (&x, &u);
  41   fp_copy (&y, &v);
  42   fp_set (&A, 1);
  43   fp_set (&D, 1);
  44
  45 top:
  46   /* 4.  while u is even do */
  47   while (fp_iseven (&u) == 1) {
  48     /* 4.1 u = u/2 */
  49     fp_div_2 (&u, &u);
  50
  51     /* 4.2 if A or B is odd then */
  52     if (fp_isodd (&A) == 1 || fp_isodd (&B) == 1) {
  53       /* A = (A+y)/2, B = (B-x)/2 */
  54       fp_add (&A, &y, &A);
  55       fp_sub (&B, &x, &B);
  56     }
  57     /* A = A/2, B = B/2 */
  58     fp_div_2 (&A, &A);
  59     fp_div_2 (&B, &B);
  60   }
  61
  62   /* 5.  while v is even do */
  63   while (fp_iseven (&v) == 1) {
  64     /* 5.1 v = v/2 */
  65     fp_div_2 (&v, &v);
  66
  67     /* 5.2 if C or D is odd then */
  68     if (fp_isodd (&C) == 1 || fp_isodd (&D) == 1) {
  69       /* C = (C+y)/2, D = (D-x)/2 */
  70       fp_add (&C, &y, &C);
  71       fp_sub (&D, &x, &D);
  72     }
  73     /* C = C/2, D = D/2 */
  74     fp_div_2 (&C, &C);
  75     fp_div_2 (&D, &D);
  76   }
  77
  78   /* 6.  if u >= v then */
  79   if (fp_cmp (&u, &v) != FP_LT) {
  80     /* u = u - v, A = A - C, B = B - D */
  81     fp_sub (&u, &v, &u);
  82     fp_sub (&A, &C, &A);
  83     fp_sub (&B, &D, &B);
  84   } else {
  85     /* v - v - u, C = C - A, D = D - B */
  86     fp_sub (&v, &u, &v);
  87     fp_sub (&C, &A, &C);
  88     fp_sub (&D, &B, &D);
  89   }
  90
  91   /* if not zero goto step 4 */
  92   if (fp_iszero (&u) == 0)
  93     goto top;
  94
  95   /* now a = C, b = D, gcd == g*v */
  96
  97   /* if v != 1 then there is no inverse */
  98   if (fp_cmp_d (&v, 1) != FP_EQ) {
  99     return FP_VAL;
 100   }
 101
 102   /* if its too low */
 103   while (fp_cmp_d(&C, 0) == FP_LT) {
 104       fp_add(&C, b, &C);
 105   }
 106
 107   /* too big */
 108   while (fp_cmp_mag(&C, b) != FP_LT) {
 109       fp_sub(&C, b, &C);
 110   }
 111
 112   /* C is now the inverse */
 113   fp_copy(&C, c);
 114   return FP_OKAY;
 115 }
 116
 117 /* c = 1/a (mod b) for odd b only */
 118 int fp_invmod(fp_int *a, fp_int *b, fp_int *c)
 119 {
 120   fp_int  x, y, u, v, B, D;
 121   int     neg;
 122
 123   /* 2. [modified] b must be odd   */
 124   if (fp_iseven (b) == FP_YES) {
 125     return fp_invmod_slow(a,b,c);
 126   }
 127
 128   /* init all our temps */
 129   fp_init(&x);  fp_init(&y);
 130   fp_init(&u);  fp_init(&v);
 131   fp_init(&B);  fp_init(&D);
 132
 133   /* x == modulus, y == value to invert */
 134   fp_copy(b, &x);
 135
 136   /* we need y = |a| */
 137   fp_abs(a, &y);
 138
 139   /* 3. u=x, v=y, A=1, B=0, C=0,D=1 */
 140   fp_copy(&x, &u);
 141   fp_copy(&y, &v);
 142   fp_set (&D, 1);
 143
 144 top:
 145   /* 4.  while u is even do */
 146   while (fp_iseven (&u) == FP_YES) {
 147     /* 4.1 u = u/2 */
 148     fp_div_2 (&u, &u);
 149
 150     /* 4.2 if B is odd then */
 151     if (fp_isodd (&B) == FP_YES) {
 152       fp_sub (&B, &x, &B);
 153     }
 154     /* B = B/2 */
 155     fp_div_2 (&B, &B);
 156   }
 157
 158   /* 5.  while v is even do */
 159   while (fp_iseven (&v) == FP_YES) {
 160     /* 5.1 v = v/2 */
 161     fp_div_2 (&v, &v);
 162
 163     /* 5.2 if D is odd then */
 164     if (fp_isodd (&D) == FP_YES) {
 165       /* D = (D-x)/2 */
 166       fp_sub (&D, &x, &D);
 167     }
 168     /* D = D/2 */
 169     fp_div_2 (&D, &D);
 170   }
 171
 172   /* 6.  if u >= v then */
 173   if (fp_cmp (&u, &v) != FP_LT) {
 174     /* u = u - v, B = B - D */
 175     fp_sub (&u, &v, &u);
 176     fp_sub (&B, &D, &B);
 177   } else {
 178     /* v - v - u, D = D - B */
 179     fp_sub (&v, &u, &v);
 180     fp_sub (&D, &B, &D);
 181   }
 182
 183   /* if not zero goto step 4 */
 184   if (fp_iszero (&u) == FP_NO) {
 185     goto top;
 186   }
 187
 188   /* now a = C, b = D, gcd == g*v */
 189
 190   /* if v != 1 then there is no inverse */
 191   if (fp_cmp_d (&v, 1) != FP_EQ) {
 192     return FP_VAL;
 193   }
 194
 195   /* b is now the inverse */
 196   neg = a->sign;
 197   while (D.sign == FP_NEG) {
 198     fp_add (&D, b, &D);
 199   }
 200   fp_copy (&D, c);
 201   c->sign = neg;
 202   return FP_OKAY;
 203 }
 204
 205 /* $Source$ */
 206 /* $Revision$ */
 207 /* $Date$ */