sys/src/libmp/port/mpmodop.c

   1 #include "os.h"
   2 #include <mp.h>
   3
   4 /* operands need to have m->top+1 digits of space and satisfy 0 ≤ a ≤ m-1 */
   5 static mpint*
   6 modarg(mpint *a, mpint *m)
   7 {
   8         if(a->size <= m->top || a->sign < 0 || mpmagcmp(a, m) >= 0){
   9                 a = mpcopy(a);
  10                 mpmod(a, m, a);
  11                 mpbits(a, Dbits*(m->top+1));
  12                 a->top = m->top;
  13         } else if(a->top < m->top){
  14                 memset(&a->p[a->top], 0, (m->top - a->top)*Dbytes);
  15         }
  16         return a;
  17 }
  18
  19 void
  20 mpmodadd(mpint *b1, mpint *b2, mpint *m, mpint *sum)
  21 {
  22         mpint *a, *b;
  23         mpdigit d;
  24         int i, j;
  25
  26         a = modarg(b1, m);
  27         b = modarg(b2, m);
  28
  29         sum->flags |= (a->flags | b->flags) & MPtimesafe;
  30         mpbits(sum, Dbits*2*(m->top+1));
  31
  32         mpvecadd(a->p, m->top, b->p, m->top, sum->p);
  33         mpvecsub(sum->p, m->top+1, m->p, m->top, sum->p+m->top+1);
  34
  35         d = sum->p[2*m->top+1];
  36         for(i = 0, j = m->top+1; i < m->top; i++, j++)
  37                 sum->p[i] = (sum->p[i] & d) | (sum->p[j] & ~d);
  38
  39         sum->top = m->top;
  40         sum->sign = 1;
  41         mpnorm(sum);
  42
  43         if(a != b1)
  44                 mpfree(a);
  45         if(b != b2)
  46                 mpfree(b);
  47 }
  48
  49 void
  50 mpmodsub(mpint *b1, mpint *b2, mpint *m, mpint *diff)
  51 {
  52         mpint *a, *b;
  53         mpdigit d;
  54         int i, j;
  55
  56         a = modarg(b1, m);
  57         b = modarg(b2, m);
  58
  59         diff->flags |= (a->flags | b->flags) & MPtimesafe;
  60         mpbits(diff, Dbits*2*(m->top+1));
  61
  62         a->p[m->top] = 0;
  63         mpvecsub(a->p, m->top+1, b->p, m->top, diff->p);
  64         mpvecadd(diff->p, m->top, m->p, m->top, diff->p+m->top+1);
  65
  66         d = ~diff->p[m->top];
  67         for(i = 0, j = m->top+1; i < m->top; i++, j++)
  68                 diff->p[i] = (diff->p[i] & d) | (diff->p[j] & ~d);
  69
  70         diff->top = m->top;
  71         diff->sign = 1;
  72         mpnorm(diff);
  73
  74         if(a != b1)
  75                 mpfree(a);
  76         if(b != b2)
  77                 mpfree(b);
  78 }
  79
  80 void
  81 mpmodmul(mpint *b1, mpint *b2, mpint *m, mpint *prod)
  82 {
  83         mpint *a, *b;
  84
  85         a = modarg(b1, m);
  86         b = modarg(b2, m);
  87
  88         mpmul(a, b, prod);
  89         mpmod(prod, m, prod);
  90
  91         if(a != b1)
  92                 mpfree(a);
  93         if(b != b2)
  94                 mpfree(b);
  95 }