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

   1 /* -*-c-*-
   2  *
   3  * $Id$
   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 /*----- Header files ------------------------------------------------------*/
  31
  32 #include "mp.h"
  33 #include "mpbarrett.h"
  34
  35 /*----- Main code ---------------------------------------------------------*/
  36
  37 /* --- @mpbarrett_create@ --- *
  38  *
  39  * Arguments:   @mpbarrett *mb@ = pointer to Barrett reduction context
  40  *              @mp *m@ = modulus to work to
  41  *
  42  *
  43  * Returns:     Zero on success, nonzero on error.
  44  *
  45  * Use:         Initializes a Barrett reduction context ready for use.
  46  */
  47
  48 int mpbarrett_create(mpbarrett *mb, mp *m)
  49 {
  50   mp *b;
  51
  52   /* --- Validate the arguments --- */
  53
  54   if (!MP_POSP(m))
  55     return (-1);
  56
  57   /* --- Compute %$\mu$% --- */
  58
  59   mp_shrink(m);
  60   mb->k = MP_LEN(m);
  61   mb->m = MP_COPY(m);
  62   b = mp_new(2 * mb->k + 1, 0);
  63   MPX_ZERO(b->v, b->vl - 1);
  64   b->vl[-1] = 1;
  65   mp_div(&b, 0, b, m);
  66   mb->mu = b;
  67   return (0);
  68 }
  69
  70 /* --- @mpbarrett_destroy@ --- *
  71  *
  72  * Arguments:   @mpbarrett *mb@ = pointer to Barrett reduction context
  73  *
  74  * Returns:     ---
  75  *
  76  * Use:         Destroys a Barrett reduction context releasing any resources
  77  *              claimed.
  78  */
  79
  80 void mpbarrett_destroy(mpbarrett *mb)
  81 {
  82   mp_drop(mb->m);
  83   mp_drop(mb->mu);
  84 }
  85
  86 /* --- @mpbarrett_reduce@ --- *
  87  *
  88  * Arguments:   @mpbarrett *mb@ = pointer to Barrett reduction context
  89  *              @mp *d@ = destination for result
  90  *              @mp *m@ = number to reduce
  91  *
  92  * Returns:     The residue of @m@ modulo the number in the reduction
  93  *              context.
  94  *
  95  * Use:         Performs an efficient modular reduction.
  96  */
  97
  98 mp *mpbarrett_reduce(mpbarrett *mb, mp *d, mp *m)
  99 {
 100   mp *q;
 101   size_t k = mb->k;
 102
 103   /* --- Special case if @m@ is too small --- */
 104
 105   if (MP_LEN(m) < k) {
 106     m = MP_COPY(m);
 107     if (d)
 108       MP_DROP(d);
 109     return (m);
 110   }
 111
 112   /* --- First stage --- */
 113
 114   {
 115     mp qq;
 116     mp_build(&qq, m->v + (k - 1), m->vl);
 117     q = mp_mul(MP_NEW, &qq, mb->mu);
 118     if (MP_LEN(q) <= k) {
 119       m = MP_COPY(m);
 120       if (d)
 121         MP_DROP(d);
 122       MP_DROP(q);
 123       return (m);
 124     }
 125   }
 126
 127   /* --- Second stage --- */
 128
 129   {
 130     mp *r;
 131     mpw *mvl;
 132
 133     MP_COPY(m);
 134     if (MP_LEN(m) <= k + 1)
 135       mvl = m->vl;
 136     else
 137       mvl = m->v + k + 1;
 138     r = mp_new(k + 1, (q->f | mb->m->f) & MP_BURN);
 139     mpx_umul(r->v, r->vl, q->v + k + 1, q->vl, mb->m->v, mb->m->vl);
 140     MP_DEST(d, k + 1, r->f | MP_UNDEF);
 141     mpx_usub(d->v, d->vl, m->v, mvl, r->v, r->vl);
 142     d->f = (m->f | r->f) & (MP_BURN | MP_NEG);
 143     MP_DROP(r);
 144     MP_DROP(q);
 145     MP_DROP(m);
 146   }
 147
 148   /* --- Final stage --- */
 149
 150   MP_SHRINK(d);
 151   while (MPX_UCMP(d->v, d->vl, >=, mb->m->v, mb->m->vl))
 152     mpx_usub(d->v, d->vl, d->v, d->vl, mb->m->v, mb->m->vl);
 153
 154   /* --- Fix up the sign --- */
 155
 156   if (d->f & MP_NEG) {
 157     mpx_usub(d->v, d->vl, mb->m->v, mb->m->vl, d->v, d->vl);
 158     d->f &= ~MP_NEG;
 159   }
 160
 161   MP_SHRINK(d);
 162   return (d);
 163 }
 164
 165 /*----- Test rig ----------------------------------------------------------*/
 166
 167 #ifdef TEST_RIG
 168
 169 static int vmod(dstr *v)
 170 {
 171   mp *x = *(mp **)v[0].buf;
 172   mp *n = *(mp **)v[1].buf;
 173   mp *r = *(mp **)v[2].buf;
 174   mp *s;
 175   mpbarrett mb;
 176   int ok = 1;
 177
 178   mpbarrett_create(&mb, n);
 179   s = mpbarrett_reduce(&mb, MP_NEW, x);
 180
 181   if (!MP_EQ(s, r)) {
 182     fputs("\n*** barrett reduction failure\n", stderr);
 183     fputs("x = ", stderr); mp_writefile(x, stderr, 10); fputc('\n', stderr);
 184     fputs("n = ", stderr); mp_writefile(n, stderr, 10); fputc('\n', stderr);
 185     fputs("r = ", stderr); mp_writefile(r, stderr, 10); fputc('\n', stderr);
 186     fputs("s = ", stderr); mp_writefile(s, stderr, 10); fputc('\n', stderr);
 187     ok = 0;
 188   }
 189
 190   mpbarrett_destroy(&mb);
 191   mp_drop(x);
 192   mp_drop(n);
 193   mp_drop(r);
 194   mp_drop(s);
 195   assert(mparena_count(MPARENA_GLOBAL) == 0);
 196   return (ok);
 197 }
 198
 199 static test_chunk tests[] = {
 200   { "mpbarrett-reduce", vmod, { &type_mp, &type_mp, &type_mp, 0 } },
 201   { 0, 0, { 0 } }
 202 };
 203
 204 int main(int argc, char *argv[])
 205 {
 206   sub_init();
 207   test_run(argc, argv, tests, SRCDIR "/tests/mpbarrett");
 208   return (0);
 209 }
 210
 211 #endif
 212
 213 /*----- That's all, folks -------------------------------------------------*/