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

   1 /* -*-c-*-
   2  *
   3  * $Id: mpbarrett.c,v 1.9 2004/04/01 12:50:09 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.9  2004/04/01 12:50:09  mdw
  34  * Add cyclic group abstraction, with test code.  Separate off exponentation
  35  * functions for better static linking.  Fix a buttload of bugs on the way.
  36  * Generally ensure that negative exponents do inversion correctly.  Add
  37  * table of standard prime-field subgroups.  (Binary field subgroups are
  38  * currently unimplemented but easy to add if anyone ever finds a good one.)
  39  *
  40  * Revision 1.8  2001/06/16 13:00:20  mdw
  41  * Use the generic exponentiation functions.
  42  *
  43  * Revision 1.7  2001/04/19 18:25:26  mdw
  44  * Use sliding-window exponentiation.
  45  *
  46  * Revision 1.6  2000/10/08 12:03:44  mdw
  47  * (mpbarrett_reduce): Cope with negative numbers.
  48  *
  49  * Revision 1.5  2000/07/29 17:04:33  mdw
  50  * Change to use left-to-right bitwise exponentiation.  This will improve
  51  * performance when the base is small.
  52  *
  53  * Revision 1.4  2000/06/17 11:45:09  mdw
  54  * Major memory management overhaul.  Added arena support.  Use the secure
  55  * arena for secret integers.  Replace and improve the MP management macros
  56  * (e.g., replace MP_MODIFY by MP_DEST).
  57  *
  58  * Revision 1.3  1999/12/12 15:08:52  mdw
  59  * Don't bother shifting %$q$% in @mpbarrett_reduce@, just skip the least
  60  * significant digits.
  61  *
  62  * Revision 1.2  1999/12/11 01:50:56  mdw
  63  * Improve initialization slightly.
  64  *
  65  * Revision 1.1  1999/12/10 23:21:59  mdw
  66  * Barrett reduction support: works with even moduli.
  67  *
  68  */
  69
  70 /*----- Header files ------------------------------------------------------*/
  71
  72 #include "mp.h"
  73 #include "mpbarrett.h"
  74
  75 /*----- Main code ---------------------------------------------------------*/
  76
  77 /* --- @mpbarrett_create@ --- *
  78  *
  79  * Arguments:   @mpbarrett *mb@ = pointer to Barrett reduction context
  80  *              @mp *m@ = modulus to work to
  81  *
  82  *
  83  * Returns:     ---
  84  *
  85  * Use:         Initializes a Barrett reduction context ready for use.
  86  */
  87
  88 void mpbarrett_create(mpbarrett *mb, mp *m)
  89 {
  90   mp *b;
  91
  92   /* --- Validate the arguments --- */
  93
  94   assert(((void)"Barrett modulus must be positive", (m->f & MP_NEG) == 0));
  95
  96   /* --- Compute %$\mu$% --- */
  97
  98   mp_shrink(m);
  99   mb->k = MP_LEN(m);
 100   mb->m = MP_COPY(m);
 101   b = mp_new(2 * mb->k + 1, 0);
 102   MPX_ZERO(b->v, b->vl - 1);
 103   b->vl[-1] = 1;
 104   mp_div(&b, 0, b, m);
 105   mb->mu = b;
 106 }
 107
 108 /* --- @mpbarrett_destroy@ --- *
 109  *
 110  * Arguments:   @mpbarrett *mb@ = pointer to Barrett reduction context
 111  *
 112  * Returns:     ---
 113  *
 114  * Use:         Destroys a Barrett reduction context releasing any resources
 115  *              claimed.
 116  */
 117
 118 void mpbarrett_destroy(mpbarrett *mb)
 119 {
 120   mp_drop(mb->m);
 121   mp_drop(mb->mu);
 122 }
 123
 124 /* --- @mpbarrett_reduce@ --- *
 125  *
 126  * Arguments:   @mpbarrett *mb@ = pointer to Barrett reduction context
 127  *              @mp *d@ = destination for result
 128  *              @mp *m@ = number to reduce
 129  *
 130  * Returns:     The residue of @m@ modulo the number in the reduction
 131  *              context.
 132  *
 133  * Use:         Performs an efficient modular reduction.
 134  */
 135
 136 mp *mpbarrett_reduce(mpbarrett *mb, mp *d, mp *m)
 137 {
 138   mp *q;
 139   size_t k = mb->k;
 140
 141   /* --- Special case if @m@ is too small --- */
 142
 143   if (MP_LEN(m) < k) {
 144     m = MP_COPY(m);
 145     if (d)
 146       MP_DROP(d);
 147     return (m);
 148   }
 149
 150   /* --- First stage --- */
 151
 152   {
 153     mp qq;
 154     mp_build(&qq, m->v + (k - 1), m->vl);
 155     q = mp_mul(MP_NEW, &qq, mb->mu);
 156     if (MP_LEN(q) <= k) {
 157       m = MP_COPY(m);
 158       if (d)
 159         MP_DROP(d);
 160       return (m);
 161     }
 162   }
 163
 164   /* --- Second stage --- */
 165
 166   {
 167     mp *r;
 168     mpw *mvl;
 169
 170     MP_COPY(m);
 171     if (MP_LEN(m) <= k + 1)
 172       mvl = m->vl;
 173     else
 174       mvl = m->v + k + 1;
 175     r = mp_new(k + 1, (q->f | mb->m->f) & MP_BURN);
 176     mpx_umul(r->v, r->vl, q->v + k + 1, q->vl, mb->m->v, mb->m->vl);
 177     MP_DEST(d, k + 1, r->f);
 178     mpx_usub(d->v, d->vl, m->v, mvl, r->v, r->vl);
 179     d->f = (m->f | r->f) & (MP_BURN | MP_NEG);
 180     MP_DROP(r);
 181     MP_DROP(q);
 182     MP_DROP(m);
 183   }
 184
 185   /* --- Final stage --- */
 186
 187   MP_SHRINK(d);
 188   while (MPX_UCMP(d->v, d->vl, >=, mb->m->v, mb->m->vl))
 189     mpx_usub(d->v, d->vl, d->v, d->vl, mb->m->v, mb->m->vl);
 190
 191   /* --- Fix up the sign --- */
 192
 193   if (d->f & MP_NEG) {
 194     mpx_usub(d->v, d->vl, mb->m->v, mb->m->vl, d->v, d->vl);
 195     d->f &= ~MP_NEG;
 196   }
 197
 198   MP_SHRINK(d);
 199   return (d);
 200 }
 201
 202 /*----- Test rig ----------------------------------------------------------*/
 203
 204 #ifdef TEST_RIG
 205
 206 static int vmod(dstr *v)
 207 {
 208   mp *x = *(mp **)v[0].buf;
 209   mp *n = *(mp **)v[1].buf;
 210   mp *r = *(mp **)v[2].buf;
 211   mp *s;
 212   mpbarrett mb;
 213   int ok = 1;
 214
 215   mpbarrett_create(&mb, n);
 216   s = mpbarrett_reduce(&mb, MP_NEW, x);
 217
 218   if (!MP_EQ(s, r)) {
 219     fputs("\n*** barrett reduction failure\n", stderr);
 220     fputs("x = ", stderr); mp_writefile(x, stderr, 10); fputc('\n', stderr);
 221     fputs("n = ", stderr); mp_writefile(n, stderr, 10); fputc('\n', stderr);
 222     fputs("r = ", stderr); mp_writefile(r, stderr, 10); fputc('\n', stderr);
 223     fputs("s = ", stderr); mp_writefile(s, stderr, 10); fputc('\n', stderr);
 224     ok = 0;
 225   }
 226
 227   mpbarrett_destroy(&mb);
 228   mp_drop(x);
 229   mp_drop(n);
 230   mp_drop(r);
 231   mp_drop(s);
 232   assert(mparena_count(MPARENA_GLOBAL) == 0);
 233   return (ok);
 234 }
 235
 236 static test_chunk tests[] = {
 237   { "mpbarrett-reduce", vmod, { &type_mp, &type_mp, &type_mp, 0 } },
 238   { 0, 0, { 0 } }
 239 };
 240
 241 int main(int argc, char *argv[])
 242 {
 243   sub_init();
 244   test_run(argc, argv, tests, SRCDIR "/tests/mpbarrett");
 245   return (0);
 246 }
 247
 248 #endif
 249
 250 /*----- That's all, folks -------------------------------------------------*/