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

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