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       return (m);
 123     }
 124   }
 125
 126   /* --- Second stage --- */
 127
 128   {
 129     mp *r;
 130     mpw *mvl;
 131
 132     MP_COPY(m);
 133     if (MP_LEN(m) <= k + 1)
 134       mvl = m->vl;
 135     else
 136       mvl = m->v + k + 1;
 137     r = mp_new(k + 1, (q->f | mb->m->f) & MP_BURN);
 138     mpx_umul(r->v, r->vl, q->v + k + 1, q->vl, mb->m->v, mb->m->vl);
 139     MP_DEST(d, k + 1, r->f);
 140     mpx_usub(d->v, d->vl, m->v, mvl, r->v, r->vl);
 141     d->f = (m->f | r->f) & (MP_BURN | MP_NEG);
 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   /* --- Fix up the sign --- */
 154
 155   if (d->f & MP_NEG) {
 156     mpx_usub(d->v, d->vl, mb->m->v, mb->m->vl, d->v, d->vl);
 157     d->f &= ~MP_NEG;
 158   }
 159
 160   MP_SHRINK(d);
 161   return (d);
 162 }
 163
 164 /*----- Test rig ----------------------------------------------------------*/
 165
 166 #ifdef TEST_RIG
 167
 168 static int vmod(dstr *v)
 169 {
 170   mp *x = *(mp **)v[0].buf;
 171   mp *n = *(mp **)v[1].buf;
 172   mp *r = *(mp **)v[2].buf;
 173   mp *s;
 174   mpbarrett mb;
 175   int ok = 1;
 176
 177   mpbarrett_create(&mb, n);
 178   s = mpbarrett_reduce(&mb, MP_NEW, x);
 179
 180   if (!MP_EQ(s, r)) {
 181     fputs("\n*** barrett reduction failure\n", stderr);
 182     fputs("x = ", stderr); mp_writefile(x, stderr, 10); fputc('\n', stderr);
 183     fputs("n = ", stderr); mp_writefile(n, stderr, 10); fputc('\n', stderr);
 184     fputs("r = ", stderr); mp_writefile(r, stderr, 10); fputc('\n', stderr);
 185     fputs("s = ", stderr); mp_writefile(s, stderr, 10); fputc('\n', stderr);
 186     ok = 0;
 187   }
 188
 189   mpbarrett_destroy(&mb);
 190   mp_drop(x);
 191   mp_drop(n);
 192   mp_drop(r);
 193   mp_drop(s);
 194   assert(mparena_count(MPARENA_GLOBAL) == 0);
 195   return (ok);
 196 }
 197
 198 static test_chunk tests[] = {
 199   { "mpbarrett-reduce", vmod, { &type_mp, &type_mp, &type_mp, 0 } },
 200   { 0, 0, { 0 } }
 201 };
 202
 203 int main(int argc, char *argv[])
 204 {
 205   sub_init();
 206   test_run(argc, argv, tests, SRCDIR "/tests/mpbarrett");
 207   return (0);
 208 }
 209
 210 #endif
 211
 212 /*----- That's all, folks -------------------------------------------------*/