chiark - git - mdw - catacomb/blob - mpbarrett.c

   1 /* -*-c-*-
   2  *
   3  * $Id: mpbarrett.c,v 1.3 1999/12/12 15:08:52 mdw Exp $
   4  *
   5  * Barrett modular reduction
   6  *
   7  * (c) 1999 Straylight/Edgeware
   8  */
   9
  10 /*----- Licensing notice --------------------------------------------------*
  11  *
  12  * This file is part of Catacomb.
  13  *
  14  * Catacomb is free software; you can redistribute it and/or modify
  15  * it under the terms of the GNU Library General Public License as
  16  * published by the Free Software Foundation; either version 2 of the
  17  * License, or (at your option) any later version.
  18  *
  19  * Catacomb is distributed in the hope that it will be useful,
  20  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  21  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  22  * GNU Library General Public License for more details.
  23  *
  24  * You should have received a copy of the GNU Library General Public
  25  * License along with Catacomb; if not, write to the Free
  26  * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
  27  * MA 02111-1307, USA.
  28  */
  29
  30 /*----- Revision history --------------------------------------------------*
  31  *
  32  * $Log: mpbarrett.c,v $
  33  * Revision 1.3  1999/12/12 15:08:52  mdw
  34  * Don't bother shifting %$q$% in @mpbarrett_reduce@, just skip the least
  35  * significant digits.
  36  *
  37  * Revision 1.2  1999/12/11 01:50:56  mdw
  38  * Improve initialization slightly.
  39  *
  40  * Revision 1.1  1999/12/10 23:21:59  mdw
  41  * Barrett reduction support: works with even moduli.
  42  *
  43  */
  44
  45 /*----- Header files ------------------------------------------------------*/
  46
  47 #include "mp.h"
  48 #include "mpbarrett.h"
  49
  50 /*----- Main code ---------------------------------------------------------*/
  51
  52 /* --- @mpbarrett_create@ --- *
  53  *
  54  * Arguments:   @mpbarrett *mb@ = pointer to Barrett reduction context
  55  *              @mp *m@ = modulus to work to
  56  *
  57  *
  58  * Returns:     ---
  59  *
  60  * Use:         Initializes a Barrett reduction context ready for use.
  61  */
  62
  63 void mpbarrett_create(mpbarrett *mb, mp *m)
  64 {
  65   mp *b;
  66
  67   /* --- Validate the arguments --- */
  68
  69   assert(((void)"Barrett modulus must be positive", (m->f & MP_NEG) == 0));
  70
  71   /* --- Compute %$\mu$% --- */
  72
  73   mp_shrink(m);
  74   mb->k = MP_LEN(m);
  75   mb->m = MP_COPY(m);
  76   b = mp_create(2 * mb->k + 1);
  77   MPX_ZERO(b->v, b->vl - 1);
  78   b->vl[-1] = 1;
  79   mp_div(&b, 0, b, m);
  80   mb->mu = b;
  81 }
  82
  83 /* --- @mpbarrett_destroy@ --- *
  84  *
  85  * Arguments:   @mpbarrett *mb@ = pointer to Barrett reduction context
  86  *
  87  * Returns:     ---
  88  *
  89  * Use:         Destroys a Barrett reduction context releasing any resources
  90  *              claimed.
  91  */
  92
  93 void mpbarrett_destroy(mpbarrett *mb)
  94 {
  95   mp_drop(mb->m);
  96   mp_drop(mb->mu);
  97 }
  98
  99 /* --- @mpbarrett_reduce@ --- *
 100  *
 101  * Arguments:   @mpbarrett *mb@ = pointer to Barrett reduction context
 102  *              @mp *d@ = destination for result
 103  *              @mp *m@ = number to reduce
 104  *
 105  * Returns:     The residue of @m@ modulo the number in the reduction
 106  *              context.
 107  *
 108  * Use:         Performs an efficient modular reduction.  The argument is
 109  *              assumed to be positive.
 110  */
 111
 112 mp *mpbarrett_reduce(mpbarrett *mb, mp *d, mp *m)
 113 {
 114   mp *q;
 115   size_t k = mb->k;
 116
 117   /* --- Special case if @m@ is too small --- */
 118
 119   if (MP_LEN(m) < k) {
 120     m = MP_COPY(m);
 121     MP_DROP(d);
 122     return (m);
 123   }
 124
 125   /* --- First stage --- */
 126
 127   {
 128     mp qq;
 129     mp_build(&qq, m->v + (k - 1), m->vl);
 130     q = mp_mul(MP_NEW, &qq, mb->mu);
 131     if (MP_LEN(q) <= k) {
 132       m = MP_COPY(m);
 133       MP_DROP(d);
 134       return (m);
 135     }
 136   }
 137
 138   /* --- Second stage --- */
 139
 140   {
 141     mp *r;
 142     mpw *mvl;
 143
 144     MP_COPY(m);
 145     if (MP_LEN(m) <= k + 1)
 146       mvl = m->vl;
 147     else
 148       mvl = m->v + k + 1;
 149     r = mp_create(k + 1);
 150     mpx_umul(r->v, r->vl, q->v + k + 1, q->vl, mb->m->v, mb->m->vl);
 151     r->f = (q->f | mb->m->f) & MP_BURN;
 152     MP_MODIFY(d, k + 1);
 153     mpx_usub(d->v, d->vl, m->v, mvl, r->v, r->vl);
 154     d->f = (m->f | r->f) & MP_BURN;
 155     MP_DROP(r);
 156     MP_DROP(q);
 157     MP_DROP(m);
 158   }
 159
 160   /* --- Final stage --- */
 161
 162   MP_SHRINK(d);
 163   while (MPX_UCMP(d->v, d->vl, >=, mb->m->v, mb->m->vl))
 164     mpx_usub(d->v, d->vl, d->v, d->vl, mb->m->v, mb->m->vl);
 165
 166   MP_SHRINK(d);
 167   return (d);
 168 }
 169
 170 /* --- @mpbarrett_exp@ --- *
 171  *
 172  * Arguments:   @mpbarrett *mb@ = pointer to Barrett reduction context
 173  *              @mp *d@ = fake destination
 174  *              @mp *a@ = base
 175  *              @mp *e@ = exponent
 176  *
 177  * Returns:     Result, %$a^e \bmod m$%.
 178  */
 179
 180 mp *mpbarrett_exp(mpbarrett *mb, mp *d, mp *a, mp *e)
 181 {
 182   mpscan sc;
 183   mp *x = MP_ONE;
 184   mp *spare = MP_NEW;
 185
 186   a = MP_COPY(a);
 187   mp_scan(&sc, e);
 188   if (MP_STEP(&sc)) {
 189     size_t sq = 0;
 190     for (;;) {
 191       mp *dd;
 192       if (MP_BIT(&sc)) {
 193         while (sq) {
 194           dd = mp_sqr(spare, a);
 195           dd = mpbarrett_reduce(mb, dd, dd);
 196           spare = a; a = dd;
 197           sq--;
 198         }
 199         dd = mp_mul(spare, x, a);
 200         dd = mpbarrett_reduce(mb, dd, dd);
 201         spare = x; x = dd;
 202       }
 203       sq++;
 204       if (!MP_STEP(&sc))
 205         break;
 206     }
 207   }
 208
 209   MP_DROP(a);
 210   if (spare != MP_NEW)
 211     MP_DROP(spare);
 212   if (d != MP_NEW)
 213     MP_DROP(d);
 214   return (x);
 215 }
 216
 217 /*----- Test rig ----------------------------------------------------------*/
 218
 219 #ifdef TEST_RIG
 220
 221 static int vmod(dstr *v)
 222 {
 223   mp *x = *(mp **)v[0].buf;
 224   mp *n = *(mp **)v[1].buf;
 225   mp *r = *(mp **)v[2].buf;
 226   mp *s;
 227   mpbarrett mb;
 228   int ok = 1;
 229
 230   mpbarrett_create(&mb, n);
 231   s = mpbarrett_reduce(&mb, MP_NEW, x);
 232
 233   if (MP_CMP(s, !=, r)) {
 234     fputs("\n*** barrett reduction failure\n", stderr);
 235     fputs("x = ", stderr); mp_writefile(x, stderr, 10); fputc('\n', stderr);
 236     fputs("n = ", stderr); mp_writefile(n, stderr, 10); fputc('\n', stderr);
 237     fputs("r = ", stderr); mp_writefile(r, stderr, 10); fputc('\n', stderr);
 238     fputs("s = ", stderr); mp_writefile(s, stderr, 10); fputc('\n', stderr);
 239     ok = 0;
 240   }
 241
 242   mpbarrett_destroy(&mb);
 243   mp_drop(x);
 244   mp_drop(n);
 245   mp_drop(r);
 246   mp_drop(s);
 247   assert(mparena_count(MPARENA_GLOBAL) == 0);
 248   return (ok);
 249 }
 250
 251 static int vexp(dstr *v)
 252 {
 253   mp *m = *(mp **)v[0].buf;
 254   mp *a = *(mp **)v[1].buf;
 255   mp *b = *(mp **)v[2].buf;
 256   mp *r = *(mp **)v[3].buf;
 257   mp *mr;
 258   int ok = 1;
 259
 260   mpbarrett mb;
 261   mpbarrett_create(&mb, m);
 262
 263   mr = mpbarrett_exp(&mb, MP_NEW, a, b);
 264
 265   if (MP_CMP(mr, !=, r)) {
 266     fputs("\n*** barrett modexp failed", stderr);
 267     fputs("\n m = ", stderr); mp_writefile(m, stderr, 10);
 268     fputs("\n a = ", stderr); mp_writefile(a, stderr, 10);
 269     fputs("\n e = ", stderr); mp_writefile(b, stderr, 10);
 270     fputs("\n r = ", stderr); mp_writefile(r, stderr, 10);
 271     fputs("\nmr = ", stderr); mp_writefile(mr, stderr, 10);
 272     fputc('\n', stderr);
 273     ok = 0;
 274   }
 275
 276   mp_drop(m);
 277   mp_drop(a);
 278   mp_drop(b);
 279   mp_drop(r);
 280   mp_drop(mr);
 281   mpbarrett_destroy(&mb);
 282   assert(mparena_count(MPARENA_GLOBAL) == 0);
 283   return ok;
 284 }
 285
 286 static test_chunk tests[] = {
 287   { "mpbarrett-reduce", vmod, { &type_mp, &type_mp, &type_mp, 0 } },
 288   { "mpbarrett-exp", vexp, { &type_mp, &type_mp, &type_mp, &type_mp, 0 } },
 289   { 0, 0, { 0 } }
 290 };
 291
 292 int main(int argc, char *argv[])
 293 {
 294   sub_init();
 295   test_run(argc, argv, tests, SRCDIR "/tests/mpbarrett");
 296   return (0);
 297 }
 298
 299 #endif
 300
 301 /*----- That's all, folks -------------------------------------------------*/