chiark - git - mdw - catacomb/blob - math/mpmont-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 "mpmont.h"
  32
  33 #define EXP_WINSZ 3
  34 #include "mpmont-exp.h"
  35
  36 /*----- Main code ---------------------------------------------------------*/
  37
  38 /* --- @mpmont_mexpr@ --- *
  39  *
  40  * Arguments:   @const mpmont *mm@ = pointer to Montgomery 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  *              except that the %$g_i$% and result are in Montgomery form.
  52  */
  53
  54 static mp *mexpr(const mpmont *mm, mp *d, mp_expfactor *f, size_t n)
  55 {
  56   mp *a = MP_COPY(mm->r);
  57   mp *spare = MP_NEW;
  58   size_t i;
  59
  60   for (i = 0; i < n; i++) {
  61     mp *t;
  62     if (f[i].exp->f & MP_BURN)
  63       spare = MP_NEWSEC;
  64     if (MP_NEGP(f[i].exp)) {
  65       t = mpmont_reduce(mm, f[i].base, f[i].base);
  66       t = mp_modinv(t, t, mm->m);
  67       f[i].base = mpmont_mul(mm, t, t, mm->r2);
  68     }
  69   }
  70   EXP_SIMUL(a, f, n);
  71   mp_drop(d);
  72   mp_drop(spare);
  73   for (i = 0; i < n; i++)
  74     MP_DROP(f[i].base);
  75   xfree(f);
  76   return (a);
  77 }
  78
  79 mp *mpmont_mexpr(const mpmont *mm, mp *d, const mp_expfactor *f, size_t n)
  80 {
  81   mp_expfactor *ff = xmalloc(n * sizeof(mp_expfactor));
  82   size_t i;
  83
  84   for (i = 0; i < n; i++) {
  85     ff[i].base = MP_COPY(f[i].base);
  86     ff[i].exp = f[i].exp;
  87   }
  88   return (mexpr(mm, d, ff, n));
  89 }
  90
  91 /* --- @mpmont_mexp@ --- *
  92  *
  93  * Arguments:   @const mpmont *mm@ = pointer to Montgomery reduction context
  94  *              @mp *d@ = fake destination
  95  *              @const mp_expfactor *f@ = pointer to array of factors
  96  *              @size_t n@ = number of factors supplied
  97  *
  98  * Returns:     Product of bases raised to exponents, all mod @m@.
  99  *
 100  * Use:         Convenient interface over @mpmont_mexpr@.
 101  */
 102
 103 mp *mpmont_mexp(const mpmont *mm, mp *d, const mp_expfactor *f, size_t n)
 104 {
 105   mp_expfactor *ff = xmalloc(n * sizeof(mp_expfactor));
 106   size_t i;
 107
 108   for (i = 0; i < n; i++) {
 109     ff[i].base = mpmont_mul(mm, MP_NEW, f[i].base, mm->r2);
 110     ff[i].exp = f[i].exp;
 111   }
 112   d = mexpr(mm, d, ff, n);
 113   return (mpmont_reduce(mm, d, d));
 114 }
 115
 116 /*----- Test rig ----------------------------------------------------------*/
 117
 118 #ifdef TEST_RIG
 119
 120 #include <mLib/testrig.h>
 121
 122 static int verify(size_t n, dstr *v)
 123 {
 124   mp *m = *(mp **)v[0].buf;
 125   mp_expfactor *f = xmalloc(n * sizeof(*f));
 126   mp *r, *rr;
 127   size_t i, j;
 128   mpmont mm;
 129   int ok = 1;
 130
 131   j = 1;
 132   for (i = 0; i < n; i++) {
 133     f[i].base = *(mp **)v[j++].buf;
 134     f[i].exp = *(mp **)v[j++].buf;
 135   }
 136
 137   rr = *(mp **)v[j].buf;
 138   mpmont_create(&mm, m);
 139   r = mpmont_mexp(&mm, MP_NEW, f, n);
 140   if (!MP_EQ(r, rr)) {
 141     fputs("\n*** mexp failed\n", stderr);
 142     fputs("m = ", stderr); mp_writefile(m, stderr, 10);
 143     for (i = 0; i < n; i++) {
 144       fprintf(stderr, "\ng_%lu = ", (unsigned long)i);
 145       mp_writefile(f[i].base, stderr, 10);
 146       fprintf(stderr, "\ne_%lu = ", (unsigned long)i);
 147       mp_writefile(f[i].exp, stderr, 10);
 148     }
 149     fputs("\nr = ", stderr); mp_writefile(r, stderr, 10);
 150     fputs("\nR = ", stderr); mp_writefile(rr, stderr, 10);
 151     fputc('\n', stderr);
 152     ok = 0;
 153   }
 154
 155   for (i = 0; i < n; i++) {
 156     MP_DROP(f[i].base);
 157     MP_DROP(f[i].exp);
 158   }
 159   MP_DROP(m);
 160   MP_DROP(r);
 161   MP_DROP(rr);
 162   mpmont_destroy(&mm);
 163   assert(mparena_count(MPARENA_GLOBAL) == 0);
 164   return (ok);
 165 }
 166
 167 static int t1(dstr *v) { return verify(1, v); }
 168 static int t2(dstr *v) { return verify(2, v); }
 169 static int t3(dstr *v) { return verify(3, v); }
 170 static int t4(dstr *v) { return verify(4, v); }
 171 static int t5(dstr *v) { return verify(5, v); }
 172
 173 static test_chunk tests[] = {
 174   { "mexp-1", t1, { &type_mp,
 175                     &type_mp, &type_mp,
 176                     &type_mp, 0 } },
 177   { "mexp-2", t2, { &type_mp,
 178                     &type_mp, &type_mp,
 179                     &type_mp, &type_mp,
 180                     &type_mp, 0 } },
 181   { "mexp-3", t3, { &type_mp,
 182                     &type_mp, &type_mp,
 183                     &type_mp, &type_mp,
 184                     &type_mp, &type_mp,
 185                     &type_mp, 0 } },
 186   { "mexp-4", t4, { &type_mp,
 187                     &type_mp, &type_mp,
 188                     &type_mp, &type_mp,
 189                     &type_mp, &type_mp,
 190                     &type_mp, &type_mp,
 191                     &type_mp, 0 } },
 192   { "mexp-5", t5, { &type_mp,
 193                     &type_mp, &type_mp,
 194                     &type_mp, &type_mp,
 195                     &type_mp, &type_mp,
 196                     &type_mp, &type_mp,
 197                     &type_mp, &type_mp,
 198                     &type_mp, 0 } },
 199   { 0, 0, { 0 } }
 200 };
 201
 202 int main(int argc, char *argv[])
 203 {
 204   sub_init();
 205   test_run(argc, argv, tests, SRCDIR "/t/mpmont");
 206   return (0);
 207 }
 208
 209 #endif
 210
 211 /*----- That's all, folks -------------------------------------------------*/