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

/* -*-c-*-
 *
 * $Id$
 *
 * 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.
 */

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

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

#define EXP_WINSZ 3
#include "mpmont-exp.h"

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

/* --- @mpmont_mexpr@ --- *
 *
 * Arguments:	@mpmont *mm@ = pointer to Montgomery reduction context
 *		@mp *d@ = fake destination
 *		@const mp_expfactor *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}} \bmod m$%
 *
 *		except that the %$g_i$% and result are in Montgomery form.
 */

static mp *mexpr(mpmont *mm, mp *d, mp_expfactor *f, size_t n)
{
  mp *a = MP_COPY(mm->r);
  mp *spare = MP_NEW;
  size_t i;

  for (i = 0; i < n; i++) {
    mp *t;
    if (f[i].exp->f & MP_BURN)
      spare = MP_NEWSEC;
    if (MP_NEGP(f[i].exp)) {
      t = mpmont_reduce(mm, f[i].base, f[i].base);
      t = mp_modinv(t, t, mm->m);
      f[i].base = mpmont_mul(mm, t, t, mm->r2);
    }
  }
  EXP_SIMUL(a, f, n);
  mp_drop(d);
  mp_drop(spare);
  for (i = 0; i < n; i++)
    MP_DROP(f[i].base);
  xfree(f);
  return (a);
}

mp *mpmont_mexpr(mpmont *mm, mp *d, const mp_expfactor *f, size_t n)
{
  mp_expfactor *ff = xmalloc(n * sizeof(mp_expfactor));
  size_t i;

  for (i = 0; i < n; i++) {
    ff[i].base = MP_COPY(f[i].base);
    ff[i].exp = f[i].exp;
  }
  return (mexpr(mm, d, ff, n));
}

/* --- @mpmont_mexp@ --- *
 *
 * Arguments:	@mpmont *mm@ = pointer to Montgomery reduction context
 *		@mp *d@ = fake destination
 *		@const mp_expfactor *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, const mp_expfactor *f, size_t n)
{
  mp_expfactor *ff = xmalloc(n * sizeof(mp_expfactor));
  size_t i;

  for (i = 0; i < n; i++) {
    ff[i].base = mpmont_mul(mm, MP_NEW, f[i].base, mm->r2);
    ff[i].exp = f[i].exp;
  }
  d = mexpr(mm, d, ff, n);
  return (mpmont_reduce(mm, d, d));
}

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

#ifdef TEST_RIG

#include <mLib/testrig.h>

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