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

   1 /* -*-c-*-
   2  *
   3  * Multiple simultaneous exponentiations
   4  *
   5  * (c) 1999 Straylight/Edgeware
   6  */
   7
   8 /*----- Licensing notice --------------------------------------------------*
   9  *
  10  * This file is part of Catacomb.
  11  *
  12  * Catacomb is free software; you can redistribute it and/or modify
  13  * it under the terms of the GNU Library General Public License as
  14  * published by the Free Software Foundation; either version 2 of the
  15  * License, or (at your option) any later version.
  16  *
  17  * Catacomb is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  20  * GNU Library General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Library General Public
  23  * License along with Catacomb; if not, write to the Free
  24  * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
  25  * MA 02111-1307, USA.
  26  */
  27
  28 /*----- Header files ------------------------------------------------------*/
  29
  30 #include "mp.h"
  31 #include "mpbarrett.h"
  32
  33 #define EXP_WINSZ 3
  34 #include "mpbarrett-exp.h"
  35
  36 /*----- Main code ---------------------------------------------------------*/
  37
  38 /* --- @mpbarrett_mexp@ --- *
  39  *
  40  * Arguments:   @const mpbarrett *mb@ = pointer to Barrett reduction context
  41  *              @mp *d@ = fake destination
  42  *              @const mp_expfactor *f@ = pointer to array of factors
  43  *              @size_t n@ = number of factors supplied
  44  *
  45  * Returns:     If the bases are %$g_0, g_1, \ldots, g_{n-1}$% and the
  46  *              exponents are %$e_0, e_1, \ldots, e_{n-1}$% then the result
  47  *              is:
  48  *
  49  *              %$g_0^{e_0} g_1^{e_1} \ldots g_{n-1}^{e_{n-1}} \bmod m$%
  50  */
  51
  52 mp *mpbarrett_mexp(const mpbarrett *mb, mp *d,
  53                    const mp_expfactor *f, size_t n)
  54 {
  55   mp_expfactor *ff = xmalloc(n * sizeof(mp_expfactor));
  56   mp *a = MP_ONE;
  57   mp *spare;
  58   mp *g = MP_NEW;
  59   size_t i;
  60
  61   spare = MP_NEW;
  62   for (i = 0; i < n; i++) {
  63     if (f[i].exp->f & MP_BURN)
  64       spare = MP_NEWSEC;
  65     if (MP_NEGP(f[i].exp))
  66       ff[i].base = mp_modinv(MP_NEW, f[i].base, mb->m);
  67     else
  68       ff[i].base = MP_COPY(f[i].base);
  69     ff[i].exp = f[i].exp;
  70   }
  71   mp_drop(g);
  72   EXP_SIMUL(a, ff, n);
  73   mp_drop(d);
  74   mp_drop(spare);
  75   for (i = 0; i < n; i++)
  76     mp_drop(ff[i].base);
  77   xfree(ff);
  78   return (a);
  79 }
  80
  81 /*----- Test rig ----------------------------------------------------------*/
  82
  83 #ifdef TEST_RIG
  84
  85 #include <mLib/testrig.h>
  86
  87 static int verify(size_t n, dstr *v)
  88 {
  89   mp *m = *(mp **)v[0].buf;
  90   mp_expfactor *f = xmalloc(n * sizeof(*f));
  91   mp *r, *rr;
  92   size_t i, j;
  93   mpbarrett mb;
  94   int ok = 1;
  95
  96   j = 1;
  97   for (i = 0; i < n; i++) {
  98     f[i].base = *(mp **)v[j++].buf;
  99     f[i].exp = *(mp **)v[j++].buf;
 100   }
 101
 102   rr = *(mp **)v[j].buf;
 103   mpbarrett_create(&mb, m);
 104   r = mpbarrett_mexp(&mb, MP_NEW, f, n);
 105   if (!MP_EQ(r, rr)) {
 106     fputs("\n*** mexp failed\n", stderr);
 107     fputs("m = ", stderr); mp_writefile(m, stderr, 10);
 108     for (i = 0; i < n; i++) {
 109       fprintf(stderr, "\ng_%lu = ", (unsigned long)i);
 110       mp_writefile(f[i].base, stderr, 10);
 111       fprintf(stderr, "\ne_%lu = ", (unsigned long)i);
 112       mp_writefile(f[i].exp, stderr, 10);
 113     }
 114     fputs("\nr = ", stderr); mp_writefile(r, stderr, 10);
 115     fputs("\nR = ", stderr); mp_writefile(rr, stderr, 10);
 116     fputc('\n', stderr);
 117     ok = 0;
 118   }
 119
 120   for (i = 0; i < n; i++) {
 121     MP_DROP(f[i].base);
 122     MP_DROP(f[i].exp);
 123   }
 124   MP_DROP(m);
 125   MP_DROP(r);
 126   MP_DROP(rr);
 127   mpbarrett_destroy(&mb);
 128   assert(mparena_count(MPARENA_GLOBAL) == 0);
 129   return (ok);
 130 }
 131
 132 static int t1(dstr *v) { return verify(1, v); }
 133 static int t2(dstr *v) { return verify(2, v); }
 134 static int t3(dstr *v) { return verify(3, v); }
 135 static int t4(dstr *v) { return verify(4, v); }
 136 static int t5(dstr *v) { return verify(5, v); }
 137
 138 static test_chunk tests[] = {
 139   { "mexp-1", t1, { &type_mp,
 140                     &type_mp, &type_mp,
 141                     &type_mp, 0 } },
 142   { "mexp-2", t2, { &type_mp,
 143                     &type_mp, &type_mp,
 144                     &type_mp, &type_mp,
 145                     &type_mp, 0 } },
 146   { "mexp-3", t3, { &type_mp,
 147                     &type_mp, &type_mp,
 148                     &type_mp, &type_mp,
 149                     &type_mp, &type_mp,
 150                     &type_mp, 0 } },
 151   { "mexp-4", t4, { &type_mp,
 152                     &type_mp, &type_mp,
 153                     &type_mp, &type_mp,
 154                     &type_mp, &type_mp,
 155                     &type_mp, &type_mp,
 156                     &type_mp, 0 } },
 157   { "mexp-5", t5, { &type_mp,
 158                     &type_mp, &type_mp,
 159                     &type_mp, &type_mp,
 160                     &type_mp, &type_mp,
 161                     &type_mp, &type_mp,
 162                     &type_mp, &type_mp,
 163                     &type_mp, 0 } },
 164   { 0, 0, { 0 } }
 165 };
 166
 167 int main(int argc, char *argv[])
 168 {
 169   sub_init();
 170   test_run(argc, argv, tests, SRCDIR "/t/mpbarrett");
 171   return (0);
 172 }
 173
 174 #endif
 175
 176 /*----- That's all, folks -------------------------------------------------*/