[u/mdw/catacomb] / math / mpbarrett-mexp.c

/* -*-c-*-
 *
 * 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 "mpbarrett.h"

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

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

/* --- @mpbarrett_mexp@ --- *
 *
 * Arguments:	@mpbarrett *mb@ = pointer to Barrett 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$%
 */

mp *mpbarrett_mexp(mpbarrett *mb, mp *d, const mp_expfactor *f, size_t n)
{
  mp_expfactor *ff = xmalloc(n * sizeof(mp_expfactor));
  mp *a = MP_ONE;
  mp *spare;
  mp *g = MP_NEW;
  size_t i;

  spare = MP_NEW;
  for (i = 0; i < n; i++) {
    if (f[i].exp->f & MP_BURN)
      spare = MP_NEWSEC;
    if (MP_NEGP(f[i].exp))
      ff[i].base = mp_modinv(MP_NEW, f[i].base, mb->m);
    else
      ff[i].base = MP_COPY(f[i].base);
    ff[i].exp = f[i].exp;
  }
  mp_drop(g);
  EXP_SIMUL(a, ff, n);
  mp_drop(d);
  mp_drop(spare);
  for (i = 0; i < n; i++)
    mp_drop(ff[i].base);
  xfree(ff);
  return (a);
}

/*----- 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;
  mpbarrett mb;
  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;
  mpbarrett_create(&mb, m);
  r = mpbarrett_mexp(&mb, 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_%lu = ", (unsigned long)i);
      mp_writefile(f[i].base, stderr, 10);
      fprintf(stderr, "\ne_%lu = ", (unsigned long)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);
  mpbarrett_destroy(&mb);
  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 "/t/mpbarrett");
  return (0);
}

#endif

/*----- That's all, folks -------------------------------------------------*/
Commit	Line	Data
ab0350a6	1	/* --c--
ab0350a6	2	*
ab0350a6	3	* Multiple simultaneous exponentiations
	4	*
	5	* (c) 1999 Straylight/Edgeware
	6	*/
	7
45c0fd36	8	/----- Licensing notice --------------------------------------------------
ab0350a6	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.
45c0fd36	16	*
ab0350a6	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.
45c0fd36	21	*
ab0350a6	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
ab0350a6	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: @mpbarrett *mb@ = pointer to Barrett reduction context
	41	* @mp *d@ = fake destination
34e4f738	42	* @const mp_expfactor *f@ = pointer to array of factors
ab0350a6	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
34e4f738	52	mp mpbarrett_mexp(mpbarrett mb, mp d, const mp_expfactor f, size_t n)
ab0350a6	53	{
34e4f738	54	mp_expfactor ff = xmalloc(n sizeof(mp_expfactor));
ab0350a6	55	mp *a = MP_ONE;
ab0350a6	56	mp *spare;
34e4f738	57	mp *g = MP_NEW;
ab0350a6	58	size_t i;
	59
	60	spare = MP_NEW;
	61	for (i = 0; i < n; i++) {
34e4f738	62	if (f[i].exp->f & MP_BURN)
ab0350a6	63	spare = MP_NEWSEC;
a69a3efd	64	if (MP_NEGP(f[i].exp))
b817bfc6	65	ff[i].base = mp_modinv(MP_NEW, f[i].base, mb->m);
a69a3efd	66	else
a69a3efd	67	ff[i].base = MP_COPY(f[i].base);
34e4f738	68	ff[i].exp = f[i].exp;
ab0350a6	69	}
34e4f738	70	mp_drop(g);
34e4f738	71	EXP_SIMUL(a, ff, n);
ab0350a6	72	mp_drop(d);
ab0350a6	73	mp_drop(spare);
34e4f738	74	for (i = 0; i < n; i++)
	75	mp_drop(ff[i].base);
	76	xfree(ff);
ab0350a6	77	return (a);
	78	}
	79
	80	/----- Test rig ----------------------------------------------------------/
	81
	82	#ifdef TEST_RIG
	83
	84	#include <mLib/testrig.h>
	85
	86	static int verify(size_t n, dstr *v)
	87	{
	88	mp m = (mp **)v[0].buf;
	89	mp_expfactor f = xmalloc(n sizeof(*f));
	90	mp r, rr;
	91	size_t i, j;
	92	mpbarrett mb;
	93	int ok = 1;
	94
	95	j = 1;
	96	for (i = 0; i < n; i++) {
	97	f[i].base = (mp *)v[j++].buf;
	98	f[i].exp = (mp *)v[j++].buf;
	99	}
	100
	101	rr = (mp *)v[j].buf;
	102	mpbarrett_create(&mb, m);
	103	r = mpbarrett_mexp(&mb, MP_NEW, f, n);
	104	if (!MP_EQ(r, rr)) {
	105	fputs("\n*** mexp failed\n", stderr);
	106	fputs("m = ", stderr); mp_writefile(m, stderr, 10);
	107	for (i = 0; i < n; i++) {
bb77b1d1	108	fprintf(stderr, "\ng_%lu = ", (unsigned long)i);
ab0350a6	109	mp_writefile(f[i].base, stderr, 10);
bb77b1d1	110	fprintf(stderr, "\ne_%lu = ", (unsigned long)i);
ab0350a6	111	mp_writefile(f[i].exp, stderr, 10);
	112	}
	113	fputs("\nr = ", stderr); mp_writefile(r, stderr, 10);
	114	fputs("\nR = ", stderr); mp_writefile(rr, stderr, 10);
	115	fputc('\n', stderr);
	116	ok = 0;
	117	}
	118
	119	for (i = 0; i < n; i++) {
	120	MP_DROP(f[i].base);
	121	MP_DROP(f[i].exp);
	122	}
	123	MP_DROP(m);
	124	MP_DROP(r);
	125	MP_DROP(rr);
	126	mpbarrett_destroy(&mb);
	127	assert(mparena_count(MPARENA_GLOBAL) == 0);
	128	return (ok);
	129	}
	130
	131	static int t1(dstr *v) { return verify(1, v); }
	132	static int t2(dstr *v) { return verify(2, v); }
	133	static int t3(dstr *v) { return verify(3, v); }
	134	static int t4(dstr *v) { return verify(4, v); }
	135	static int t5(dstr *v) { return verify(5, v); }
	136
	137	static test_chunk tests[] = {
	138	{ "mexp-1", t1, { &type_mp,
	139	&type_mp, &type_mp,
	140	&type_mp, 0 } },
	141	{ "mexp-2", t2, { &type_mp,
	142	&type_mp, &type_mp,
	143	&type_mp, &type_mp,
	144	&type_mp, 0 } },
	145	{ "mexp-3", t3, { &type_mp,
	146	&type_mp, &type_mp,
	147	&type_mp, &type_mp,
	148	&type_mp, &type_mp,
	149	&type_mp, 0 } },
	150	{ "mexp-4", t4, { &type_mp,
	151	&type_mp, &type_mp,
	152	&type_mp, &type_mp,
	153	&type_mp, &type_mp,
	154	&type_mp, &type_mp,
	155	&type_mp, 0 } },
	156	{ "mexp-5", t5, { &type_mp,
	157	&type_mp, &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, 0 } },
	163	{ 0, 0, { 0 } }
	164	};
	165
	166	int main(int argc, char *argv[])
	167	{
	168	sub_init();
0f00dc4c	169	test_run(argc, argv, tests, SRCDIR "/t/mpbarrett");
ab0350a6	170	return (0);
45c0fd36	171	}
ab0350a6	172
	173	#endif
	174
	175	/----- That's all, folks -------------------------------------------------/