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

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