pb_div.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
  15 /* Divides a by b (uses the char of b as the char of the result) only for
  16  * polynomials over GF(p^k)[x]
  17  */
  18 int pb_div(pb_poly *a, pb_poly *b, pb_poly *c, pb_poly *d)
  19 {
  20    int     err, x;
  21    pb_poly p, q, r;
  22    mp_int  tmp, tmp2;
  23
  24    /* if ||b|| > ||a|| it's a simple copy over */
  25    if (b->used > a->used) {
  26       if (d != NULL) {
  27          pb_copy(a, d);
  28       }
  29       if (c != NULL) {
  30          pb_zero(c);
  31       }
  32    }
  33
  34    /* zero divisor */
  35    if (b->used == 0) {
  36       return MP_VAL;
  37    }
  38
  39    /* compare chars */
  40    if (mp_cmp(&(a->characteristic), &(b->characteristic)) != MP_EQ) {
  41       return MP_VAL;
  42    }
  43
  44    /* we don't support char zero! */
  45    if (mp_iszero(&(b->characteristic)) == MP_YES) {
  46       return MP_VAL;
  47    }
  48
  49    /* get a copy of "a" to work with */
  50    if ((err = pb_init_copy(&p, a)) != MP_OKAY) {
  51       return err;
  52    }
  53
  54    /* init a temp quotient */
  55    if ((err = pb_init_size(&q, &(b->characteristic), a->used - b->used + 1)) != MP_OKAY) {
  56       goto __P;
  57    }
  58
  59    /* init a temp polynomial we can work with */
  60    if ((err = pb_init_size(&r, &(b->characteristic), a->used)) != MP_OKAY) {
  61       goto __Q;
  62    }
  63
  64    /* now init an mp_int we can work with */
  65    if ((err = mp_init(&tmp)) != MP_OKAY) {
  66       goto __R;
  67    }
  68
  69    /* and the inverse of the leading digit of b(x) */
  70    if ((err = mp_init(&tmp2)) != MP_OKAY) {
  71       goto __TMP;
  72    }
  73    if ((err = mp_invmod(&(b->terms[b->used - 1]), &(b->characteristic), &tmp2)) != MP_OKAY) {
  74       goto __TMP2;
  75    }
  76
  77    /* now let's reduce us some digits */
  78    for (x = a->used - 1; x >= b->used-1; x--) {
  79       /* skip zero digits */
  80       if (mp_iszero(&(p.terms[x])) == MP_YES) {
  81           continue;
  82       }
  83
  84       /* now we have a leading digit of p(x), call it A
  85          and a leading digit of b(x), call it B
  86
  87          Now we want to reduce it so we need A + CB = 0 which turns into
  88          A + CB = 0
  89          CB = -A
  90          C = -A/B
  91
  92          So we multiply b(x) by C * x^k [where k = ||p(x)|| - ||b(x)||], add that to p(x) and
  93          we must reduced one digit
  94       */
  95
  96       /* multiply 1/B [tmp2] by A */
  97       if ((err = mp_mulmod(&tmp2, &(p.terms[x]), &(b->characteristic), &tmp)) != MP_OKAY)    { goto __TMP2; }
  98
  99       /* tmp is now a term of the quotient */
 100       if ((err = mp_copy(&tmp, &(q.terms[x - b->used + 1]))) != MP_OKAY)                     { goto __TMP2; }
 101
 102       /* create r(x) = C */
 103       pb_zero(&r);
 104       if ((err = mp_copy(&tmp, &(r.terms[0]))) != MP_OKAY)                                   { goto __TMP2; }
 105       r.used = 1;
 106
 107       /* now multiply r(x) by b(x)*x^k and subtract from p(x) */
 108       if ((err = pb_mul(b, &r, &r)) != MP_OKAY)                                              { goto __TMP2; }
 109       if ((err = pb_lshd(&r, x - b->used + 1)) != MP_OKAY)                                   { goto __TMP2; }
 110       if ((err = pb_sub(&p, &r, &p)) != MP_OKAY)                                             { goto __TMP2; }
 111    }
 112
 113    /* setup q properly (so q.used is valid) */
 114    q.used = a->used;
 115    pb_clamp(&q);
 116
 117    /* ok so now p(x) is the remainder and q(x) is the quotient */
 118    if (c != NULL) {
 119       pb_exch(&q, c);
 120    }
 121    if (d != NULL) {
 122       pb_exch(&p, d);
 123    }
 124
 125    err = MP_OKAY;
 126
 127 __TMP2: mp_clear(&tmp2);
 128 __TMP : mp_clear(&tmp);
 129 __R   : pb_clear(&r);
 130 __Q   : pb_clear(&q);
 131 __P   : pb_clear(&p);
 132    return err;
 133 }
 134
 135