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

   1 /* -*-c-*-
   2  *
   3  * $Id: mpbarrett.c,v 1.2 1999/12/11 01:50:56 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.2  1999/12/11 01:50:56  mdw
  34  * Improve initialization slightly.
  35  *
  36  * Revision 1.1  1999/12/10 23:21:59  mdw
  37  * Barrett reduction support: works with even moduli.
  38  *
  39  */
  40
  41 /*----- Header files ------------------------------------------------------*/
  42
  43 #include "mp.h"
  44 #include "mpbarrett.h"
  45
  46 /*----- Main code ---------------------------------------------------------*/
  47
  48 /* --- @mpbarrett_create@ --- *
  49  *
  50  * Arguments:   @mpbarrett *mb@ = pointer to Barrett reduction context
  51  *              @mp *m@ = modulus to work to
  52  *
  53  *
  54  * Returns:     ---
  55  *
  56  * Use:         Initializes a Barrett reduction context ready for use.
  57  */
  58
  59 void mpbarrett_create(mpbarrett *mb, mp *m)
  60 {
  61   mp *b;
  62
  63   /* --- Validate the arguments --- */
  64
  65   assert(((void)"Barrett modulus must be positive", (m->f & MP_NEG) == 0));
  66
  67   /* --- Compute %$\mu$% --- */
  68
  69   mp_shrink(m);
  70   mb->k = MP_LEN(m);
  71   mb->m = MP_COPY(m);
  72   b = mp_create(2 * mb->k + 1);
  73   MPX_ZERO(b->v, b->vl - 1);
  74   b->vl[-1] = 1;
  75   mp_div(&b, 0, b, m);
  76   mb->mu = b;
  77 }
  78
  79 /* --- @mpbarrett_destroy@ --- *
  80  *
  81  * Arguments:   @mpbarrett *mb@ = pointer to Barrett reduction context
  82  *
  83  * Returns:     ---
  84  *
  85  * Use:         Destroys a Barrett reduction context releasing any resources
  86  *              claimed.
  87  */
  88
  89 void mpbarrett_destroy(mpbarrett *mb)
  90 {
  91   mp_drop(mb->m);
  92   mp_drop(mb->mu);
  93 }
  94
  95 /* --- @mpbarrett_reduce@ --- *
  96  *
  97  * Arguments:   @mpbarrett *mb@ = pointer to Barrett reduction context
  98  *              @mp *d@ = destination for result
  99  *              @mp *m@ = number to reduce
 100  *
 101  * Returns:     The residue of @m@ modulo the number in the reduction
 102  *              context.
 103  *
 104  * Use:         Performs an efficient modular reduction.  The argument is
 105  *              assumed to be positive.
 106  */
 107
 108 mp *mpbarrett_reduce(mpbarrett *mb, mp *d, mp *m)
 109 {
 110   mp *q;
 111   size_t k = mb->k;
 112
 113   /* --- Special case if @m@ is too small --- */
 114
 115   if (MP_LEN(m) < k) {
 116     m = MP_COPY(m);
 117     MP_DROP(d);
 118     return (m);
 119   }
 120
 121   /* --- First stage --- */
 122
 123   {
 124     mp qq;
 125     mp_build(&qq, m->v + (k - 1), m->vl);
 126     q = mp_mul(MP_NEW, &qq, mb->mu);
 127     q = mp_lsr(q, q, MPW_BITS * (k + 1));
 128   }
 129
 130   /* --- Second stage --- */
 131
 132   {
 133     mp *r;
 134     mpw *mvl;
 135
 136     MP_COPY(m);
 137     if (MP_LEN(m) <= k + 1)
 138       mvl = m->vl;
 139     else
 140       mvl = m->v + k + 1;
 141     r = mp_create(k + 1);
 142     mpx_umul(r->v, r->vl, q->v, q->vl, mb->m->v, mb->m->vl);
 143     r->f = (q->f | mb->m->f) & MP_BURN;
 144     MP_MODIFY(d, k + 1);
 145     mpx_usub(d->v, d->vl, m->v, mvl, r->v, r->vl);
 146     d->f = (m->f | r->f) & MP_BURN;
 147     MP_DROP(r);
 148     MP_DROP(q);
 149     MP_DROP(m);
 150   }
 151
 152   /* --- Final stage --- */
 153
 154   MP_SHRINK(d);
 155   while (MPX_UCMP(d->v, d->vl, >=, mb->m->v, mb->m->vl))
 156     mpx_usub(d->v, d->vl, d->v, d->vl, mb->m->v, mb->m->vl);
 157
 158   MP_SHRINK(d);
 159   return (d);
 160 }
 161
 162 /* --- @mpbarrett_exp@ --- *
 163  *
 164  * Arguments:   @mpbarrett *mb@ = pointer to Barrett reduction context
 165  *              @mp *d@ = fake destination
 166  *              @mp *a@ = base
 167  *              @mp *e@ = exponent
 168  *
 169  * Returns:     Result, %$a^e \bmod m$%.
 170  */
 171
 172 mp *mpbarrett_exp(mpbarrett *mb, mp *d, mp *a, mp *e)
 173 {
 174   mpscan sc;
 175   mp *x = MP_ONE;
 176   mp *spare = MP_NEW;
 177
 178   a = MP_COPY(a);
 179   mp_scan(&sc, e);
 180   if (MP_STEP(&sc)) {
 181     size_t sq = 0;
 182     for (;;) {
 183       mp *dd;
 184       if (MP_BIT(&sc)) {
 185         while (sq) {
 186           dd = mp_sqr(spare, a);
 187           dd = mpbarrett_reduce(mb, dd, dd);
 188           spare = a; a = dd;
 189           sq--;
 190         }
 191         dd = mp_mul(spare, x, a);
 192         dd = mpbarrett_reduce(mb, dd, dd);
 193         spare = x; x = dd;
 194       }
 195       sq++;
 196       if (!MP_STEP(&sc))
 197         break;
 198     }
 199   }
 200
 201   MP_DROP(a);
 202   if (spare != MP_NEW)
 203     MP_DROP(spare);
 204   if (d != MP_NEW)
 205     MP_DROP(d);
 206   return (x);
 207 }
 208
 209 /*----- Test rig ----------------------------------------------------------*/
 210
 211 #ifdef TEST_RIG
 212
 213 static int vmod(dstr *v)
 214 {
 215   mp *x = *(mp **)v[0].buf;
 216   mp *n = *(mp **)v[1].buf;
 217   mp *r = *(mp **)v[2].buf;
 218   mp *s;
 219   mpbarrett mb;
 220   int ok = 1;
 221
 222   mpbarrett_create(&mb, n);
 223   s = mpbarrett_reduce(&mb, MP_NEW, x);
 224
 225   if (MP_CMP(s, !=, r)) {
 226     fputs("\n*** barrett reduction failure\n", stderr);
 227     fputs("x = ", stderr); mp_writefile(x, stderr, 10); fputc('\n', stderr);
 228     fputs("n = ", stderr); mp_writefile(n, stderr, 10); fputc('\n', stderr);
 229     fputs("r = ", stderr); mp_writefile(r, stderr, 10); fputc('\n', stderr);
 230     fputs("s = ", stderr); mp_writefile(s, stderr, 10); fputc('\n', stderr);
 231     ok = 0;
 232   }
 233
 234   mpbarrett_destroy(&mb);
 235   mp_drop(x);
 236   mp_drop(n);
 237   mp_drop(r);
 238   mp_drop(s);
 239   assert(mparena_count(MPARENA_GLOBAL) == 0);
 240   return (ok);
 241 }
 242
 243 static int vexp(dstr *v)
 244 {
 245   mp *m = *(mp **)v[0].buf;
 246   mp *a = *(mp **)v[1].buf;
 247   mp *b = *(mp **)v[2].buf;
 248   mp *r = *(mp **)v[3].buf;
 249   mp *mr;
 250   int ok = 1;
 251
 252   mpbarrett mb;
 253   mpbarrett_create(&mb, m);
 254
 255   mr = mpbarrett_exp(&mb, MP_NEW, a, b);
 256
 257   if (MP_CMP(mr, !=, r)) {
 258     fputs("\n*** barrett modexp failed", stderr);
 259     fputs("\n m = ", stderr); mp_writefile(m, stderr, 10);
 260     fputs("\n a = ", stderr); mp_writefile(a, stderr, 10);
 261     fputs("\n e = ", stderr); mp_writefile(b, stderr, 10);
 262     fputs("\n r = ", stderr); mp_writefile(r, stderr, 10);
 263     fputs("\nmr = ", stderr); mp_writefile(mr, stderr, 10);
 264     fputc('\n', stderr);
 265     ok = 0;
 266   }
 267
 268   mp_drop(m);
 269   mp_drop(a);
 270   mp_drop(b);
 271   mp_drop(r);
 272   mp_drop(mr);
 273   mpbarrett_destroy(&mb);
 274   assert(mparena_count(MPARENA_GLOBAL) == 0);
 275   return ok;
 276 }
 277
 278 static test_chunk tests[] = {
 279   { "mpbarrett-reduce", vmod, { &type_mp, &type_mp, &type_mp, 0 } },
 280   { "mpbarrett-exp", vexp, { &type_mp, &type_mp, &type_mp, &type_mp, 0 } },
 281   { 0, 0, { 0 } }
 282 };
 283
 284 int main(int argc, char *argv[])
 285 {
 286   sub_init();
 287   test_run(argc, argv, tests, SRCDIR "/tests/mpbarrett");
 288   return (0);
 289 }
 290
 291 #endif
 292
 293 /*----- That's all, folks -------------------------------------------------*/