[u/mdw/catacomb] / exp.h

/* -*-c-*-
 *
 * $Id$
 *
 * Generalized exponentiation
 *
 * (c) 2001 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.
 */

#ifdef CATACOMB_EXP_H
#  error "Multiple inclusion of <catacomb/exp.h>"
#endif

#define CATACOMB_EXP_H

#ifdef __cplusplus
  extern "C" {
#endif

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

#include <stddef.h>

#include <mLib/alloc.h>

#ifndef CATACOMB_MP_H
#  include "mp.h"
#endif

/*----- Data structures ---------------------------------------------------*/

typedef struct exp_simulscan {
  mpw w;
  size_t len;
  const mpw *v;
} exp_simulscan;

typedef struct exp_simul {
  unsigned b;
  size_t o, n;
  exp_simulscan *s;
} exp_simul;

/*----- Macros provided ---------------------------------------------------*/

/* --- Parameters --- */

#ifndef EXP_WINSZ			/* Sliding window size */
#  define EXP_WINSZ 4			/* Predefine if you need to */
#endif

/* --- These are determined from the window size --- *
 *
 * Given a %$k$%-bit exponent, I expect to do %$k/2$% multiplies if I use the
 * simple way.  If I use an n-bit sliding window, then I do %$2^n$%
 * multiplies up front, but I only do %$(2^n - 1)/2^n k/n$% multiplies for
 * the exponentiation.  This is a win when
 *
 *   %$k \ge \frac{n 2^{n+1}}{n - 2}$%
 */

#define EXP_TABSZ (1 << EXP_WINSZ)
#define EXP_THRESH							\
    ((EXP_WINSZ * (2 << EXP_WINSZ))/((EXP_WINSZ - 2) * MPW_BITS))

/* --- Required operations --- *
 *
 * The macros here are independent of the underlying group elements.  You
 * must provide the necessary group operations and other definitions.  The
 * group operation is assumed to be written multiplicatively.
 *
 * @EXP_TYPE@			The type of a group element, e.g., @mp *@.
 *
 * @EXP_COPY(d, x)@		Makes @d@ be a copy of @x@.
 *
 * @EXP_DROP(x)@		Discards the element @x@, reclaiming any
 *				memory it used.
 *
 * @EXP_MUL(a, x)@		Multiplies @a@ by @x@ (writing the result
 *				back to @a@).
 *
 * @EXP_FIX(x)@			Makes @x@ be a canonical representation of
 *				its value.  All multiplications have the
 *				right argument canonical.
 *
 * @EXP_SQR(a)@			Multiplies @a@ by itself.
 *
 * @EXP_SETMUL(d, x, y)@	Sets @d@ to be the product of @x@ and @y@.
 *				The value @d@ has not been initialized.
 *
 * @EXP_SETSQR(d, x)@		Sets @d@ to be the square of @x@.
 *
 * Only @EXP_TYPE@, @EXP_MUL@ and @EXP_SQR@ are required for simple
 * exponentation.  Sliding window and simultaneous exponentation require all
 * of the operations.
 */

#ifndef EXP_TYPE
#  error "EXP_TYPE not defined for <catacomb/exp.h>"
#endif

/* --- @EXP_SIMPLE@ --- *
 *
 * Arguments:	@a@ = the result object, initially a multiplicative identity
 *		@g@ = the object to exponentiate
 *		@x@ = the exponent, as a multiprecision integer
 *
 * Use:		Performs a simple left-to-right exponentiation.  At the end
 *		of the code, the answer is left in @a@; @g@ and @x@ are
 *		unchanged.
 */

#define EXP_SIMPLE(a, g, x) do {					\
  mpscan sc;								\
  unsigned sq = 0;							\
									\
  /* --- Begin scanning --- */						\
									\
  mp_rscan(&sc, x);							\
  if (!MP_RSTEP(&sc))							\
    goto exp_simple_exit;						\
  while (!MP_RBIT(&sc))							\
    MP_RSTEP(&sc);							\
									\
  /* --- Do the main body of the work --- */				\
									\
  EXP_FIX(g);								\
  for (;;) {								\
    EXP_MUL(a, g);							\
    sq = 0;								\
    for (;;) {								\
      if (!MP_RSTEP(&sc))						\
	goto exp_simple_done;						\
      sq++;								\
      if (MP_RBIT(&sc))							\
	break;								\
    }									\
    while (sq--) EXP_SQR(a);						\
  }									\
									\
  /* --- Do a final round of squaring --- */				\
									\
exp_simple_done:							\
  while (sq--) EXP_SQR(a);						\
exp_simple_exit:;							\
} while (0)

/* --- @EXP_WINDOW@ --- *
 *
 * Arguments:	@a@ = the result object, initially a multiplicative identity
 *		@g@ = the object to exponentiate
 *		@x@ = the exponent, as a multiprecision integer
 *
 * Use:		Performs a sliding-window exponentiation.  At the end of the
 *		code, the answer is left in @a@; @g@ and @x@ are unchanged.
 */

#define EXP_WINDOW(a, g, x) do {					\
  EXP_TYPE *v;								\
  EXP_TYPE g2;								\
  unsigned i, sq = 0;							\
  mpscan sc;								\
									\
  /* --- Get going --- */						\
									\
  mp_rscan(&sc, x);							\
  if (!MP_RSTEP(&sc))							\
    goto exp_window_exit;						\
									\
  /* --- Do the precomputation --- */					\
									\
  EXP_FIX(g);								\
  EXP_SETSQR(g2, g);							\
  EXP_FIX(g2);								\
  v = xmalloc(EXP_TABSZ * sizeof(EXP_TYPE));				\
  EXP_COPY(v[0], g);							\
  for (i = 1; i < EXP_TABSZ; i++) {					\
    EXP_SETMUL(v[i], v[i - 1], g2);					\
    EXP_FIX(v[i]);							\
  }									\
  EXP_DROP(g2);								\
									\
  /* --- Skip top-end zero bits --- *					\
   *									\
   * If the initial step worked, there must be a set bit somewhere, so	\
   * keep stepping until I find it.					\
   */									\
									\
  while (!MP_RBIT(&sc))							\
    MP_RSTEP(&sc);							\
									\
  /* --- Now for the main work --- */					\
									\
  for (;;) {								\
    unsigned l = 1;							\
    unsigned z = 0;							\
									\
    /* --- The next bit is set, so read a window index --- *		\
     *									\
     * Reset @i@ to zero and increment @sq@.  Then, until either I read	\
     * @WINSZ@ bits or I run out of bits, scan in a bit: if it's clear,	\
     * bump the @z@ counter; if it's set, push a set bit into @i@,	\
     * shift it over by @z@ bits, bump @sq@ by @z + 1@ and clear @z@.	\
     * By the end of this palaver, @i@ is an index to the precomputed	\
     * value in @v@.							\
     */									\
									\
    i = 0;								\
    sq++;								\
    while (l < EXP_WINSZ && MP_RSTEP(&sc)) {				\
      l++;								\
      if (!MP_RBIT(&sc))						\
	z++;								\
      else {								\
	i = ((i << 1) | 1) << z;					\
	sq += z + 1;							\
	z = 0;								\
      }									\
    }									\
									\
    /* --- Do the squaring --- *					\
     *									\
     * Remember that @sq@ carries over from the zero-skipping stuff	\
     * below.								\
     */									\
									\
    while (sq--) EXP_SQR(a);						\
									\
    /* --- Do the multiply --- */					\
									\
    EXP_MUL(a, v[i]);							\
									\
    /* --- Now grind along through the rest of the bits --- */		\
									\
    sq = z;								\
    for (;;) {								\
      if (!MP_RSTEP(&sc))						\
	goto exp_window_done;						\
      if (MP_RBIT(&sc))							\
	break;								\
      sq++;								\
    }									\
  }									\
									\
  /* --- Do a final round of squaring --- */				\
									\
exp_window_done:							\
  while (sq--) EXP_SQR(a);						\
  for (i = 0; i < EXP_TABSZ; i++)					\
    EXP_DROP(v[i]);							\
  xfree(v);								\
exp_window_exit:;							\
} while (0)

/* --- @EXP_SIMUL@ --- *
 *
 * Arguments:	@a@ = the result object, initially a multiplicative identity
 *		@f@ = pointer to a vector of base/exp pairs
 *		@n@ = the number of base/exp pairs
 *
 * Use:		Performs a simultaneous sliding-window exponentiation.  The
 *		@f@ table is an array of structures containing members @base@
 *		of type @EXP_TYPE@, and @exp@ of type @mp *@.
 */

#define EXP_SIMUL(a, f, n) do {						\
  size_t i, j, jj, k;							\
  size_t vn = 1 << (EXP_WINSZ * n), m = (1 << n) - 1;			\
  EXP_TYPE *v = xmalloc(vn * sizeof(EXP_TYPE));				\
  exp_simul e;								\
  unsigned sq = 0;							\
									\
  /* --- Fill in the precomputed table --- */				\
									\
  j = 1;								\
  for (i = 0; i < n; i++) {						\
    EXP_COPY(v[j], f[n - 1 - i].base);					\
    EXP_FIX(v[j]);							\
    j <<= 1;								\
  }									\
  k = n * EXP_WINSZ;							\
  jj = 1;								\
  for (; i < k; i++) {							\
    EXP_SETSQR(v[j], v[jj]);						\
    EXP_FIX(v[j]);							\
    j <<= 1; jj <<= 1;							\
  }									\
  for (i = 1; i < vn; i <<= 1) {					\
    for (j = 1; j < i; j++) {						\
      EXP_SETMUL(v[j + i], v[j], v[i]);					\
      EXP_FIX(v[j + i]);						\
    }									\
  }									\
									\
  /* --- Set up the bitscanners --- *					\
   *									\
   * Got to use custom scanners, to keep them all in sync.		\
   */									\
									\
  e.n = n;								\
  e.b = 0;								\
  e.s = xmalloc(n * sizeof(*e.s));					\
  e.o = 0;								\
  for (i = 0; i < n; i++) {						\
    MP_SHRINK(f[i].exp);						\
    e.s[i].len = MP_LEN(f[i].exp);					\
    e.s[i].v = f[i].exp->v;						\
    if (e.s[i].len > e.o)						\
      e.o = e.s[i].len;							\
  }									\
									\
  /* --- Skip as far as a nonzero column in the exponent matrix --- */	\
									\
  do {									\
    if (!e.o && !e.b)							\
      goto exp_simul_done;						\
    i = exp_simulnext(&e, 0);						\
  } while (!(i & m));							\
									\
  /* --- Now for the main work --- */					\
									\
  for (;;) {								\
    unsigned l = 1;							\
    unsigned z = 0;							\
									\
    /* --- Just read a nonzero column, so read a window index --- *	\
     *									\
     * Clear high bits of @i@ and increment @sq@.  Then, until either I	\
     * read @WINSZ@ columns or I run out, scan in a column and append	\
     * it to @i@.  If it's zero, bump the @z@ counter; if it's nonzero,	\
     * bump @sq@ by @z + 1@ and clear @z@.  By the end of this palaver,	\
     * @i@ is an index to the precomputed value in @v@, followed by	\
     * @n * z@ zero bits.						\
     */									\
									\
    sq++;								\
    while (l < EXP_WINSZ && (e.o || e.b)) {				\
      l++;								\
      i = exp_simulnext(&e, i);						\
      if (!(i & m))							\
	z++;								\
      else {								\
	sq += z + 1;							\
	z = 0;								\
      }									\
    }									\
									\
    /* --- Do the squaring --- *					\
     *									\
     * Remember that @sq@ carries over from the zero-skipping stuff	\
     * below.								\
     */									\
									\
    while (sq--) EXP_SQR(a);						\
									\
    /* --- Do the multiply --- */					\
									\
    i >>= (z * n);							\
    EXP_MUL(a, v[i]);							\
									\
    /* --- Now grind along through the rest of the bits --- */		\
									\
    sq = z;								\
    for (;;) {								\
      if (!e.o && !e.b)							\
	goto exp_simul_done;						\
      if ((i = exp_simulnext(&e, 0)) != 0)				\
	break;								\
      sq++;								\
    }									\
  }									\
									\
  /* --- Do a final round of squaring --- */				\
									\
exp_simul_done:								\
  while (sq--) EXP_SQR(a);						\
  for (i = 1; i < vn; i++)						\
    EXP_DROP(v[i]);							\
  xfree(v);								\
  xfree(e.s);								\
} while (0)

/*----- Functions provided ------------------------------------------------*/

/* --- @exp_simulnext@ --- *
 *
 * Arguments:	@exp_simul *e@ = pointer to state structure
 *		@size_t x@ = a current accumulator
 *
 * Returns:	The next column of bits.
 *
 * Use:		Scans the next column of bits for a simultaneous
 *		exponentiation.
 */

extern size_t exp_simulnext(exp_simul */*e*/, size_t /*x*/);

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

#ifdef __cplusplus
  }
#endif
Commit	Line	Data
7765e926	1	/* --c--
7765e926	2	*
ef13e9a4	3	* $Id$
7765e926	4	*
	5	* Generalized exponentiation
	6	*
	7	* (c) 2001 Straylight/Edgeware
	8	*/
	9
45c0fd36	10	/----- Licensing notice --------------------------------------------------
7765e926	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	*
7765e926	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	*
7765e926	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
7765e926	30	#ifdef CATACOMB_EXP_H
	31	# error "Multiple inclusion of <catacomb/exp.h>"
	32	#endif
	33
	34	#define CATACOMB_EXP_H
	35
	36	#ifdef __cplusplus
	37	extern "C" {
	38	#endif
	39
	40	/----- Header files ------------------------------------------------------/
	41
	42	#include <stddef.h>
	43
	44	#include <mLib/alloc.h>
	45
	46	#ifndef CATACOMB_MP_H
	47	# include "mp.h"
	48	#endif
	49
	50	/----- Data structures ---------------------------------------------------/
	51
	52	typedef struct exp_simulscan {
	53	mpw w;
	54	size_t len;
	55	const mpw *v;
45c0fd36	56	} exp_simulscan;
7765e926	57
	58	typedef struct exp_simul {
	59	unsigned b;
	60	size_t o, n;
	61	exp_simulscan *s;
	62	} exp_simul;
	63
	64	/----- Macros provided ---------------------------------------------------/
	65
	66	/* --- Parameters --- */
	67
	68	#ifndef EXP_WINSZ /* Sliding window size */
	69	# define EXP_WINSZ 4 /* Predefine if you need to */
	70	#endif
	71
391faf42	72	/* --- These are determined from the window size --- *
	73	*
	74	* Given a %$k$%-bit exponent, I expect to do %$k/2$% multiplies if I use the
	75	* simple way. If I use an n-bit sliding window, then I do %$2^n$%
	76	* multiplies up front, but I only do %$(2^n - 1)/2^n k/n$% multiplies for
	77	* the exponentiation. This is a win when
	78	*
	79	* %$k \ge \frac{n 2^{n+1}}{n - 2}$%
	80	*/
7765e926	81
7765e926	82	#define EXP_TABSZ (1 << EXP_WINSZ)
391faf42	83	#define EXP_THRESH \
391faf42	84	((EXP_WINSZ * (2 << EXP_WINSZ))/((EXP_WINSZ - 2) * MPW_BITS))
7765e926	85
	86	/* --- Required operations --- *
	87	*
	88	* The macros here are independent of the underlying group elements. You
	89	* must provide the necessary group operations and other definitions. The
	90	* group operation is assumed to be written multiplicatively.
	91	*
	92	* @EXP_TYPE@ The type of a group element, e.g., @mp *@.
	93	*
	94	* @EXP_COPY(d, x)@ Makes @d@ be a copy of @x@.
	95	*
	96	* @EXP_DROP(x)@ Discards the element @x@, reclaiming any
	97	* memory it used.
	98	*
	99	* @EXP_MUL(a, x)@ Multiplies @a@ by @x@ (writing the result
	100	* back to @a@).
	101	*
8823192f	102	* @EXP_FIX(x)@ Makes @x@ be a canonical representation of
	103	* its value. All multiplications have the
	104	* right argument canonical.
	105	*
7765e926	106	* @EXP_SQR(a)@ Multiplies @a@ by itself.
	107	*
	108	* @EXP_SETMUL(d, x, y)@ Sets @d@ to be the product of @x@ and @y@.
	109	* The value @d@ has not been initialized.
	110	*
	111	* @EXP_SETSQR(d, x)@ Sets @d@ to be the square of @x@.
	112	*
	113	* Only @EXP_TYPE@, @EXP_MUL@ and @EXP_SQR@ are required for simple
	114	* exponentation. Sliding window and simultaneous exponentation require all
	115	* of the operations.
	116	*/
	117
	118	#ifndef EXP_TYPE
	119	# error "EXP_TYPE not defined for <catacomb/exp.h>"
	120	#endif
	121
	122	/* --- @EXP_SIMPLE@ --- *
	123	*
	124	* Arguments: @a@ = the result object, initially a multiplicative identity
	125	* @g@ = the object to exponentiate
	126	* @x@ = the exponent, as a multiprecision integer
	127	*
	128	* Use: Performs a simple left-to-right exponentiation. At the end
	129	* of the code, the answer is left in @a@; @g@ and @x@ are
	130	* unchanged.
	131	*/
	132
	133	#define EXP_SIMPLE(a, g, x) do { \
	134	mpscan sc; \
	135	unsigned sq = 0; \
	136	\
	137	/* --- Begin scanning --- */ \
	138	\
	139	mp_rscan(&sc, x); \
	140	if (!MP_RSTEP(&sc)) \
	141	goto exp_simple_exit; \
	142	while (!MP_RBIT(&sc)) \
	143	MP_RSTEP(&sc); \
	144	\
	145	/* --- Do the main body of the work --- */ \
	146	\
8823192f	147	EXP_FIX(g); \
7765e926	148	for (;;) { \
	149	EXP_MUL(a, g); \
	150	sq = 0; \
	151	for (;;) { \
	152	if (!MP_RSTEP(&sc)) \
	153	goto exp_simple_done; \
	154	sq++; \
	155	if (MP_RBIT(&sc)) \
	156	break; \
	157	} \
	158	while (sq--) EXP_SQR(a); \
	159	} \
	160	\
	161	/* --- Do a final round of squaring --- */ \
	162	\
	163	exp_simple_done: \
	164	while (sq--) EXP_SQR(a); \
	165	exp_simple_exit:; \
	166	} while (0)
	167
	168	/* --- @EXP_WINDOW@ --- *
	169	*
	170	* Arguments: @a@ = the result object, initially a multiplicative identity
	171	* @g@ = the object to exponentiate
	172	* @x@ = the exponent, as a multiprecision integer
	173	*
	174	* Use: Performs a sliding-window exponentiation. At the end of the
	175	* code, the answer is left in @a@; @g@ and @x@ are unchanged.
	176	*/
	177
	178	#define EXP_WINDOW(a, g, x) do { \
	179	EXP_TYPE *v; \
	180	EXP_TYPE g2; \
	181	unsigned i, sq = 0; \
	182	mpscan sc; \
	183	\
	184	/* --- Get going --- */ \
	185	\
	186	mp_rscan(&sc, x); \
	187	if (!MP_RSTEP(&sc)) \
	188	goto exp_window_exit; \
	189	\
	190	/* --- Do the precomputation --- */ \
	191	\
8823192f	192	EXP_FIX(g); \
7765e926	193	EXP_SETSQR(g2, g); \
8823192f	194	EXP_FIX(g2); \
7765e926	195	v = xmalloc(EXP_TABSZ * sizeof(EXP_TYPE)); \
7765e926	196	EXP_COPY(v[0], g); \
8823192f	197	for (i = 1; i < EXP_TABSZ; i++) { \
7765e926	198	EXP_SETMUL(v[i], v[i - 1], g2); \
8823192f	199	EXP_FIX(v[i]); \
8823192f	200	} \
7765e926	201	EXP_DROP(g2); \
	202	\
	203	/* --- Skip top-end zero bits --- * \
	204	* \
	205	* If the initial step worked, there must be a set bit somewhere, so \
	206	* keep stepping until I find it. \
	207	*/ \
	208	\
	209	while (!MP_RBIT(&sc)) \
	210	MP_RSTEP(&sc); \
	211	\
	212	/* --- Now for the main work --- */ \
	213	\
	214	for (;;) { \
	215	unsigned l = 1; \
	216	unsigned z = 0; \
	217	\
	218	/* --- The next bit is set, so read a window index --- * \
	219	* \
	220	* Reset @i@ to zero and increment @sq@. Then, until either I read \
	221	* @WINSZ@ bits or I run out of bits, scan in a bit: if it's clear, \
	222	* bump the @z@ counter; if it's set, push a set bit into @i@, \
	223	* shift it over by @z@ bits, bump @sq@ by @z + 1@ and clear @z@. \
	224	* By the end of this palaver, @i@ is an index to the precomputed \
	225	* value in @v@. \
	226	*/ \
	227	\
	228	i = 0; \
	229	sq++; \
	230	while (l < EXP_WINSZ && MP_RSTEP(&sc)) { \
	231	l++; \
	232	if (!MP_RBIT(&sc)) \
	233	z++; \
	234	else { \
	235	i = ((i << 1) \| 1) << z; \
	236	sq += z + 1; \
	237	z = 0; \
	238	} \
	239	} \
	240	\
	241	/* --- Do the squaring --- * \
	242	* \
	243	* Remember that @sq@ carries over from the zero-skipping stuff \
	244	* below. \
	245	*/ \
	246	\
	247	while (sq--) EXP_SQR(a); \
	248	\
	249	/* --- Do the multiply --- */ \
	250	\
	251	EXP_MUL(a, v[i]); \
	252	\
	253	/* --- Now grind along through the rest of the bits --- */ \
	254	\
	255	sq = z; \
	256	for (;;) { \
	257	if (!MP_RSTEP(&sc)) \
	258	goto exp_window_done; \
	259	if (MP_RBIT(&sc)) \
	260	break; \
	261	sq++; \
	262	} \
	263	} \
	264	\
265	/* --- Do a final round of squaring --- */ \
266	\
267	exp_window_done: \
268	while (sq--) EXP_SQR(a); \
269	for (i = 0; i < EXP_TABSZ; i++) \
270	EXP_DROP(v[i]); \
271	xfree(v); \
272	exp_window_exit:; \
273	} while (0)
274
275	/* --- @EXP_SIMUL@ --- *
276	*
277	* Arguments: @a@ = the result object, initially a multiplicative identity
278	* @f@ = pointer to a vector of base/exp pairs
279	* @n@ = the number of base/exp pairs
280	*
281	* Use: Performs a simultaneous sliding-window exponentiation. The
282	* @f@ table is an array of structures containing members @base@
283	* of type @EXP_TYPE@, and @exp@ of type @mp *@.
284	*/
285
286	#define EXP_SIMUL(a, f, n) do { \
287	size_t i, j, jj, k; \
288	size_t vn = 1 << (EXP_WINSZ * n), m = (1 << n) - 1; \
289	EXP_TYPE v = xmalloc(vn sizeof(EXP_TYPE)); \
290	exp_simul e; \
291	unsigned sq = 0; \
292	\
293	/* --- Fill in the precomputed table --- */ \
294	\
295	j = 1; \
296	for (i = 0; i < n; i++) { \
297	EXP_COPY(v[j], f[n - 1 - i].base); \
8823192f	298	EXP_FIX(v[j]); \
7765e926	299	j <<= 1; \
	300	} \
	301	k = n * EXP_WINSZ; \
	302	jj = 1; \
	303	for (; i < k; i++) { \
	304	EXP_SETSQR(v[j], v[jj]); \
8823192f	305	EXP_FIX(v[j]); \
7765e926	306	j <<= 1; jj <<= 1; \
	307	} \
	308	for (i = 1; i < vn; i <<= 1) { \
8823192f	309	for (j = 1; j < i; j++) { \
7765e926	310	EXP_SETMUL(v[j + i], v[j], v[i]); \
8823192f	311	EXP_FIX(v[j + i]); \
8823192f	312	} \
7765e926	313	} \
	314	\
	315	/* --- Set up the bitscanners --- * \
	316	* \
	317	* Got to use custom scanners, to keep them all in sync. \
	318	*/ \
	319	\
	320	e.n = n; \
	321	e.b = 0; \
	322	e.s = xmalloc(n * sizeof(*e.s)); \
	323	e.o = 0; \
	324	for (i = 0; i < n; i++) { \
	325	MP_SHRINK(f[i].exp); \
	326	e.s[i].len = MP_LEN(f[i].exp); \
	327	e.s[i].v = f[i].exp->v; \
	328	if (e.s[i].len > e.o) \
	329	e.o = e.s[i].len; \
	330	} \
	331	\
	332	/* --- Skip as far as a nonzero column in the exponent matrix --- */ \
	333	\
	334	do { \
	335	if (!e.o && !e.b) \
	336	goto exp_simul_done; \
	337	i = exp_simulnext(&e, 0); \
	338	} while (!(i & m)); \
	339	\
	340	/* --- Now for the main work --- */ \
	341	\
	342	for (;;) { \
	343	unsigned l = 1; \
	344	unsigned z = 0; \
	345	\
	346	/* --- Just read a nonzero column, so read a window index --- * \
	347	* \
	348	* Clear high bits of @i@ and increment @sq@. Then, until either I \
	349	* read @WINSZ@ columns or I run out, scan in a column and append \
	350	* it to @i@. If it's zero, bump the @z@ counter; if it's nonzero, \
	351	* bump @sq@ by @z + 1@ and clear @z@. By the end of this palaver, \
	352	* @i@ is an index to the precomputed value in @v@, followed by \
	353	* @n * z@ zero bits. \
	354	*/ \
	355	\
	356	sq++; \
	357	while (l < EXP_WINSZ && (e.o \|\| e.b)) { \
	358	l++; \
	359	i = exp_simulnext(&e, i); \
	360	if (!(i & m)) \
	361	z++; \
	362	else { \
	363	sq += z + 1; \
	364	z = 0; \
	365	} \
	366	} \
	367	\
	368	/* --- Do the squaring --- * \
	369	* \
	370	* Remember that @sq@ carries over from the zero-skipping stuff \
	371	* below. \
	372	*/ \
	373	\
	374	while (sq--) EXP_SQR(a); \
	375	\
	376	/* --- Do the multiply --- */ \
377	\
378	i >>= (z * n); \
379	EXP_MUL(a, v[i]); \
380	\
381	/* --- Now grind along through the rest of the bits --- */ \
382	\
383	sq = z; \
384	for (;;) { \
385	if (!e.o && !e.b) \
386	goto exp_simul_done; \
387	if ((i = exp_simulnext(&e, 0)) != 0) \
388	break; \
389	sq++; \
390	} \
391	} \
392	\
393	/* --- Do a final round of squaring --- */ \
394	\
395	exp_simul_done: \
396	while (sq--) EXP_SQR(a); \
8823192f	397	for (i = 1; i < vn; i++) \
7765e926	398	EXP_DROP(v[i]); \
7765e926	399	xfree(v); \
ef13e9a4	400	xfree(e.s); \
7765e926	401	} while (0)
	402
	403	/----- Functions provided ------------------------------------------------/
	404
	405	/* --- @exp_simulnext@ --- *
	406	*
	407	* Arguments: @exp_simul *e@ = pointer to state structure
	408	* @size_t x@ = a current accumulator
	409	*
	410	* Returns: The next column of bits.
	411	*
	412	* Use: Scans the next column of bits for a simultaneous
	413	* exponentiation.
	414	*/
	415
	416	extern size_t exp_simulnext(exp_simul /e/, size_t /x*/);
	417
	418	/----- That's all, folks -------------------------------------------------/
	419
	420	#ifdef __cplusplus
	421	}
	422	#endif