s_mp_invmod_odd.c

   1 #include "tommath_private.h"
   2 #ifdef S_MP_INVMOD_ODD_C
   3 /* LibTomMath, multiple-precision integer library -- Tom St Denis */
   4 /* SPDX-License-Identifier: Unlicense */
   5
   6 /* computes the modular inverse via binary extended euclidean algorithm,
   7  * that is c = 1/a mod b
   8  *
   9  * Based on slow invmod except this is optimized for the case where b is
  10  * odd as per HAC Note 14.64 on pp. 610
  11  */
  12 mp_err s_mp_invmod_odd(const mp_int *a, const mp_int *b, mp_int *c)
  13 {
  14    mp_int  x, y, u, v, B, D;
  15    mp_err  err;
  16
  17    /* 2. [modified] b must be odd   */
  18    if (mp_iseven(b)) {
  19       return MP_VAL;
  20    }
  21
  22    /* init all our temps */
  23    if ((err = mp_init_multi(&x, &y, &u, &v, &B, &D, NULL)) != MP_OKAY) {
  24       return err;
  25    }
  26
  27    /* x == modulus, y == value to invert */
  28    if ((err = mp_copy(b, &x)) != MP_OKAY)                         goto LBL_ERR;
  29
  30    /* y needs to be positive but the remainder d of mp_div(a,b,c,d) might be negative */
  31    if ((err = mp_mod(a, b, &y)) != MP_OKAY)                       goto LBL_ERR;
  32
  33    /* if one of x,y is zero return an error! */
  34    if (mp_iszero(&x) || mp_iszero(&y)) {
  35       err = MP_VAL;
  36       goto LBL_ERR;
  37    }
  38
  39    /* 3. u=x, v=y, A=1, B=0, C=0,D=1 */
  40    if ((err = mp_copy(&x, &u)) != MP_OKAY)                        goto LBL_ERR;
  41    if ((err = mp_copy(&y, &v)) != MP_OKAY)                        goto LBL_ERR;
  42    mp_set(&D, 1uL);
  43
  44    do {
  45       /* 4.  while u is even do */
  46       while (mp_iseven(&u)) {
  47          /* 4.1 u = u/2 */
  48          if ((err = mp_div_2(&u, &u)) != MP_OKAY)                    goto LBL_ERR;
  49
  50          /* 4.2 if B is odd then */
  51          if (mp_isodd(&B)) {
  52             if ((err = mp_sub(&B, &x, &B)) != MP_OKAY)               goto LBL_ERR;
  53          }
  54          /* B = B/2 */
  55          if ((err = mp_div_2(&B, &B)) != MP_OKAY)                    goto LBL_ERR;
  56       }
  57
  58       /* 5.  while v is even do */
  59       while (mp_iseven(&v)) {
  60          /* 5.1 v = v/2 */
  61          if ((err = mp_div_2(&v, &v)) != MP_OKAY)                    goto LBL_ERR;
  62
  63          /* 5.2 if D is odd then */
  64          if (mp_isodd(&D)) {
  65             /* D = (D-x)/2 */
  66             if ((err = mp_sub(&D, &x, &D)) != MP_OKAY)               goto LBL_ERR;
  67          }
  68          /* D = D/2 */
  69          if ((err = mp_div_2(&D, &D)) != MP_OKAY)                    goto LBL_ERR;
  70       }
  71
  72       /* 6.  if u >= v then */
  73       if (mp_cmp(&u, &v) != MP_LT) {
  74          /* u = u - v, B = B - D */
  75          if ((err = mp_sub(&u, &v, &u)) != MP_OKAY)                  goto LBL_ERR;
  76
  77          if ((err = mp_sub(&B, &D, &B)) != MP_OKAY)                  goto LBL_ERR;
  78       } else {
  79          /* v - v - u, D = D - B */
  80          if ((err = mp_sub(&v, &u, &v)) != MP_OKAY)                  goto LBL_ERR;
  81
  82          if ((err = mp_sub(&D, &B, &D)) != MP_OKAY)                  goto LBL_ERR;
  83       }
  84
  85       /* if not zero goto step 4 */
  86    } while (!mp_iszero(&u));
  87
  88    /* now a = C, b = D, gcd == g*v */
  89
  90    /* if v != 1 then there is no inverse */
  91    if (mp_cmp_d(&v, 1uL) != MP_EQ) {
  92       err = MP_VAL;
  93       goto LBL_ERR;
  94    }
  95
  96    /* b is now the inverse */
  97    while (mp_isneg(&D)) {
  98       if ((err = mp_add(&D, b, &D)) != MP_OKAY)                   goto LBL_ERR;
  99    }
 100
 101    /* too big */
 102    while (mp_cmp_mag(&D, b) != MP_LT) {
 103       if ((err = mp_sub(&D, b, &D)) != MP_OKAY)                   goto LBL_ERR;
 104    }
 105
 106    mp_exch(&D, c);
 107    err = MP_OKAY;
 108
 109 LBL_ERR:
 110    mp_clear_multi(&x, &y, &u, &v, &B, &D, NULL);
 111    return err;
 112 }
 113 #endif