Moved @mpmont_factor@ to <mp.h>. Documented interface change to

[u/mdw/catacomb] / mpmont-mexp.c

/* -*-c-*-
 *
 * $Id: mpmont-mexp.c,v 1.5 2000/10/08 12:11:22 mdw Exp $
 *
 * Multiple simultaneous exponentiations
 *
 * (c) 1999 Straylight/Edgeware
 */

/*----- Licensing notice --------------------------------------------------* 
 *
 * This file is part of Catacomb.
 *
 * Catacomb is free software; you can redistribute it and/or modify
 * it under the terms of the GNU Library General Public License as
 * published by the Free Software Foundation; either version 2 of the
 * License, or (at your option) any later version.
 * 
 * Catacomb is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Library General Public License for more details.
 * 
 * You should have received a copy of the GNU Library General Public
 * License along with Catacomb; if not, write to the Free
 * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 * MA 02111-1307, USA.
 */

/*----- Revision history --------------------------------------------------* 
 *
 * $Log: mpmont-mexp.c,v $
 * Revision 1.5  2000/10/08 12:11:22  mdw
 * Use @MP_EQ@ instead of @MP_CMP@.
 *
 * Revision 1.4  2000/06/17 11:45:09  mdw
 * Major memory management overhaul.  Added arena support.  Use the secure
 * arena for secret integers.  Replace and improve the MP management macros
 * (e.g., replace MP_MODIFY by MP_DEST).
 *
 * Revision 1.3  1999/12/10 23:18:39  mdw
 * Change interface for suggested destinations.
 *
 * Revision 1.2  1999/11/21 11:35:10  mdw
 * Performance improvement: use @mp_sqr@ and @mpmont_reduce@ instead of
 * @mpmont_mul@ for squaring in exponentiation.
 *
 * Revision 1.1  1999/11/19 13:19:29  mdw
 * Simultaneous exponentiation support.
 *
 */

/*----- Header files ------------------------------------------------------*/

#include "mp.h"
#include "mpmont.h"

/*----- Main code ---------------------------------------------------------*/

/* --- @mpmont_mexpr@ --- *
 *
 * Arguments:	@mpmont *mm@ = pointer to Montgomery reduction context
 *		@mp *d@ = fake destination
 *		@mpmont_factor *f@ = pointer to array of factors
 *		@size_t n@ = number of factors supplied
 *
 * Returns:	If the bases are %$g_0, g_1, \ldots, g_{n-1}$% and the
 *		exponents are %$e_0, e_1, \ldots, e_{n-1}$% then the result
 *		is:
 *
 *		%$g_0^{e_0} g_1^{e_1} \ldots g_{n-1}^{e_{n-1}} R \bmod m$%
 */

typedef struct scan {
  size_t len;
  mpw w;
} scan;

mp *mpmont_mexpr(mpmont *mm, mp *d, mpmont_factor *f, size_t n)
{
  size_t vn = 1 << n;
  mp **v = xmalloc(vn * sizeof(mp *));
  scan *s;
  size_t o;
  unsigned b;
  mp *a = MP_COPY(mm->r);
  mp *spare = MP_NEW;

  /* --- Perform the precomputation --- */

  {
    size_t i, j;
    size_t mask;

    /* --- Fill in the rest of the array --- *
     *
     * Zero never gets used.
     */

    j = 0;
    mask = 0;
    for (i = 1; i < vn; i++) {

      /* --- Check for a new bit entering --- *
       *
       * If a bit gets set that wasn't set before, then all the lower bits
       * are zeroes and I've got to introduce a new base into the array.
       */

      if ((i & mask) == 0) {
	v[i] = mpmont_mul(mm, MP_NEW, f[j++].base, mm->r2);
	mask = i;
      }

      /* --- Otherwise I can get away with a single multiplication --- *
       *
       * In particular, if %$i$% has more than one bit set, then I only need
       * to calculate %$v_i = v_{\mathit{mask}} v_{i - \mathit{mask}}$%.
       * Since both are less than %$i$%, they must have already been
       * computed.
       */

      else
	v[i] = mpmont_mul(mm, MP_NEW, v[mask], v[i & ~mask]);
    }
  }

  /* --- Set up the bitscanners --- *
   *
   * I must scan the exponents from left to right, which is a shame.  It
   * means that I can't use the standard @mpscan@ stuff, in particular.
   *
   * If any of the exponents are considered secret then make the accumulator
   * automatically set the secret bit.
   */

  {
    size_t i;

    s = xmalloc(n * sizeof(scan));
    o = 0;
    for (i = 0; i < n; i++) {
      s[i].len = MP_LEN(f[i].exp);
      if (s[i].len > o)
	o = s[i].len;
      if (f[i].exp->f & MP_BURN)
	spare = MP_NEWSEC;
    }
    b = 0;
  }

  /* --- Now do the actual calculation --- */

  b = 0;
  for (;;) {
    size_t i;
    size_t j;
    mp *dd;

    /* --- If no more bits, get some more --- */

    if (!b) {
      if (!o)
	break;
      o--;
      b = MPW_BITS;
    }

    /* --- Work out the next index --- */

    j = 0;
    b--;
    for (i = 0; i < n; i++) {
      if (o < s[i].len)
	j |= (((f[i].exp->v[o] >> b) & 1) << i);
    }

    /* --- Accumulate the result --- */

    if (spare) {
      dd = mp_sqr(spare, a);
      dd = mpmont_reduce(mm, dd, dd);
      spare = a;
      a = dd;
    }
    
    if (j) {
      dd = mpmont_mul(mm, spare, a, v[j]);
      spare = a;
      a = dd;
    }
  }

  /* --- Tidy up afterwards --- */

  {
    size_t i;
    for (i = 1; i < vn; i++)
      MP_DROP(v[i]);
    if (spare)
      MP_DROP(spare);
    free(v);
    free(s);
  }

  if (d != MP_NEW)
    MP_DROP(d);

  return (a);
}

/* --- @mpmont_mexp@ --- *
 *
 * Arguments:	@mpmont *mm@ = pointer to Montgomery reduction context
 *		@mp *d@ = fake destination
 *		@mpmont_factor *f@ = pointer to array of factors
 *		@size_t n@ = number of factors supplied
 *
 * Returns:	Product of bases raised to exponents, all mod @m@.
 *
 * Use:		Convenient interface over @mpmont_mexpr@.
 */

mp *mpmont_mexp(mpmont *mm, mp *d, mpmont_factor *f, size_t n)
{
  d = mpmont_mexpr(mm, d, f, n);
  d = mpmont_reduce(mm, d, d);
  return (d);
}

/*----- Test rig ----------------------------------------------------------*/

#ifdef TEST_RIG

#include <mLib/testrig.h>

static int verify(size_t n, dstr *v)
{
  mp *m = *(mp **)v[0].buf;
  mpmont_factor *f = xmalloc(n * sizeof(*f));
  mp *r, *rr;
  size_t i, j;
  mpmont mm;
  int ok = 1;

  j = 1;
  for (i = 0; i < n; i++) {
    f[i].base = *(mp **)v[j++].buf;
    f[i].exp = *(mp **)v[j++].buf;
  }

  rr = *(mp **)v[j].buf;
  mpmont_create(&mm, m);
  r = mpmont_mexp(&mm, MP_NEW, f, n);
  if (!MP_EQ(r, rr)) {
    fputs("\n*** mexp failed\n", stderr);
    fputs("m = ", stderr); mp_writefile(m, stderr, 10);
    for (i = 0; i < n; i++) {
      fprintf(stderr, "\ng_%u = ", i);
      mp_writefile(f[i].base, stderr, 10);
      fprintf(stderr, "\ne_%u = ", i);
      mp_writefile(f[i].exp, stderr, 10);
    }
    fputs("\nr = ", stderr); mp_writefile(r, stderr, 10);
    fputs("\nR = ", stderr); mp_writefile(rr, stderr, 10);
    fputc('\n', stderr);
    ok = 0;
  }

  for (i = 0; i < n; i++) {
    MP_DROP(f[i].base);
    MP_DROP(f[i].exp);
  }
  MP_DROP(m);
  MP_DROP(r);
  MP_DROP(rr);
  mpmont_destroy(&mm);
  assert(mparena_count(MPARENA_GLOBAL) == 0);
  return (ok);
}

static int t1(dstr *v) { return verify(1, v); }
static int t2(dstr *v) { return verify(2, v); }
static int t3(dstr *v) { return verify(3, v); }
static int t4(dstr *v) { return verify(4, v); }
static int t5(dstr *v) { return verify(5, v); }

static test_chunk tests[] = {
  { "mexp-1", t1, { &type_mp,
		    &type_mp, &type_mp,
		    &type_mp, 0 } },
  { "mexp-2", t2, { &type_mp,
		    &type_mp, &type_mp,
		    &type_mp, &type_mp,
		    &type_mp, 0 } },
  { "mexp-3", t3, { &type_mp,
		    &type_mp, &type_mp,
		    &type_mp, &type_mp,
		    &type_mp, &type_mp,
		    &type_mp, 0 } },
  { "mexp-4", t4, { &type_mp,
		    &type_mp, &type_mp,
		    &type_mp, &type_mp,
		    &type_mp, &type_mp,
		    &type_mp, &type_mp,
		    &type_mp, 0 } },
  { "mexp-5", t5, { &type_mp,
		    &type_mp, &type_mp,
		    &type_mp, &type_mp,
		    &type_mp, &type_mp,
		    &type_mp, &type_mp,
		    &type_mp, &type_mp,
		    &type_mp, 0 } },
  { 0, 0, { 0 } }
};

int main(int argc, char *argv[])
{
  sub_init();
  test_run(argc, argv, tests, SRCDIR "/tests/mpmont");
  return (0);
}   

#endif

/*----- That's all, folks -------------------------------------------------*/