pb_exteuclid.c

   1 /* LibTomPoly, Polynomial Basis Math -- Tom St Denis
   2  *
   3  * LibTomPoly is a public domain library that provides
   4  * polynomial basis arithmetic support.  It relies on
   5  * LibTomMath for large integer support.
   6  *
   7  * This library is free for all purposes without any
   8  * express guarantee that it works.
   9  *
  10  * Tom St Denis, tomstdenis@iahu.ca, http://poly.libtomcrypt.org
  11  */
  12 #include <tompoly.h>
  13
  14 /* Extended euclidean algorithm of (a, b) produces
  15    a*u1 + b*u2 = u3
  16  */
  17 int pb_exteuclid(pb_poly *a, pb_poly *b, pb_poly *U1, pb_poly *U2, pb_poly *U3)
  18 {
  19    pb_poly u1,u2,u3,v1,v2,v3,t1,t2,t3,q,tmp;
  20    int err;
  21
  22    if ((err = pb_init_multi(&(b->characteristic),
  23         &u1, &u2, &u3, &v1, &v2, &v3, &t1, &t2, &t3, &q, &tmp, NULL)) != MP_OKAY) {
  24       return err;
  25    }
  26
  27    /* initialize, (u1,u2,u3) = (1,0,a) */
  28    mp_set(&(u1.terms[0]), 1); u1.used = 1;
  29    if ((err = pb_copy(a, &u3)) != MP_OKAY)                                        { goto _ERR; }
  30
  31    /* initialize, (v1,v2,v3) = (0,1,b) */
  32    mp_set(&(v2.terms[0]), 1); v2.used = 1;
  33    if ((err = pb_copy(b, &v3)) != MP_OKAY)                                        { goto _ERR; }
  34
  35    /* loop while v3 != 0 */
  36    while (v3.used > 0) {
  37        /* q = u3/v3 */
  38        if ((err = pb_div(&u3, &v3, &q, NULL)) != MP_OKAY)                         { goto _ERR; }
  39
  40        /* (t1,t2,t3) = (u1,u2,u3) - (v1,v2,v3)q */
  41        if ((err = pb_mul(&v1, &q, &tmp)) != MP_OKAY)                              { goto _ERR; }
  42        if ((err = pb_sub(&u1, &tmp, &t1)) != MP_OKAY)                             { goto _ERR; }
  43        if ((err = pb_mul(&v2, &q, &tmp)) != MP_OKAY)                              { goto _ERR; }
  44        if ((err = pb_sub(&u2, &tmp, &t2)) != MP_OKAY)                             { goto _ERR; }
  45        if ((err = pb_mul(&v3, &q, &tmp)) != MP_OKAY)                              { goto _ERR; }
  46        if ((err = pb_sub(&u3, &tmp, &t3)) != MP_OKAY)                             { goto _ERR; }
  47
  48        /* (u1,u2,u3) = (v1,v2,v3) */
  49        if ((err = pb_copy(&v1, &u1)) != MP_OKAY)                                  { goto _ERR; }
  50        if ((err = pb_copy(&v2, &u2)) != MP_OKAY)                                  { goto _ERR; }
  51        if ((err = pb_copy(&v3, &u3)) != MP_OKAY)                                  { goto _ERR; }
  52
  53        /* (v1,v2,v3) = (t1,t2,t3) */
  54        if ((err = pb_copy(&t1, &v1)) != MP_OKAY)                                  { goto _ERR; }
  55        if ((err = pb_copy(&t2, &v2)) != MP_OKAY)                                  { goto _ERR; }
  56        if ((err = pb_copy(&t3, &v3)) != MP_OKAY)                                  { goto _ERR; }
  57    }
  58
  59    /* reduce U3 to monic but we have to mangle all three to remain consistent */
  60    pb_zero(&v1);
  61    if ((err = mp_copy(&(u3.terms[u3.used-1]), &(v1.terms[0]))) != MP_OKAY)        { goto _ERR; }
  62    v1.used = 1;
  63
  64    /* copy result out */
  65    if (U1 != NULL) { if ((err = pb_div(&u1, &v1, U1, NULL)) != MP_OKAY)           { goto _ERR; }}
  66    if (U2 != NULL) { if ((err = pb_div(&u2, &v1, U2, NULL)) != MP_OKAY)           { goto _ERR; }}
  67    if (U3 != NULL) { if ((err = pb_div(&u3, &v1, U3, NULL)) != MP_OKAY)           { goto _ERR; }}
  68
  69    err = MP_OKAY;
  70 _ERR: pb_clear_multi(&u1, &u2, &u3, &v1, &v2, &v3, &t1, &t2, &t3, &q, &tmp, NULL);
  71    return err;
  72 }
  73
  74