[catacomb] / math / mpmont.c

/* -*-c-*-
 *
 * Montgomery reduction
 *
 * (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 "config.h"
#include "dispatch.h"
#include "mp.h"
#include "mpmont.h"

/*----- Tweakables --------------------------------------------------------*/

/* --- @MPMONT_DISABLE@ --- *
 *
 * Replace all the clever Montgomery reduction with good old-fashioned long
 * division.
 */

/* #define MPMONT_DISABLE */

#define MPMONT_KTHRESH (16*MPK_THRESH)

/*----- Low-level implementation ------------------------------------------*/

#ifndef MPMONT_DISABLE

/* --- @redccore@ --- *
 *
 * Arguments:	@mpw *dv, *dvl@ = base and limit of source/destination
 *		@const mpw *mv@ = base of modulus %$m$%
 *		@size_t n@ = length of modulus
 *		@const mpw *mi@ = base of REDC coefficient %$m'$%
 *
 * Returns:	---
 *
 * Use:		Let %$a$% be the input operand.  Store in %$d$% the value
 *		%$a + (m' a \bmod R) m$%.  The destination has space for at
 *		least %$2 n + 1$% words of result.
 */

CPU_DISPATCH(static, (void), void, redccore,
	     (mpw *dv, mpw *dvl, const mpw *mv, size_t n, const mpw *mi),
	     (dv, dvl, mv, n, mi), pick_redccore, simple_redccore);

static void simple_redccore(mpw *dv, mpw *dvl, const mpw *mv,
			    size_t n, const mpw *mi)
{
  mpw mi0 = *mi;
  size_t i;

  for (i = 0; i < n; i++) {
    MPX_UMLAN(dv, dvl, mv, mv + n, MPW(*dv*mi0));
    dv++;
  }
}

#define MAYBE_REDC4(impl)						\
  extern void mpxmont_redc4_##impl(mpw *dv, mpw *dvl, const mpw *mv,	\
				   size_t n, const mpw *mi);		\
  static void maybe_redc4_##impl(mpw *dv, mpw *dvl, const mpw *mv,	\
				 size_t n, const mpw *mi)		\
  {									\
    if (n%4) simple_redccore(dv, dvl, mv, n, mi);			\
    else mpxmont_redc4_##impl(dv, dvl, mv, n, mi);			\
  }

#if CPUFAM_X86
  MAYBE_REDC4(x86_sse2)
#endif

#if CPUFAM_AMD64
  MAYBE_REDC4(amd64_sse2)
#endif

static redccore__functype *pick_redccore(void)
{
#if CPUFAM_X86
  DISPATCH_PICK_COND(mpmont_reduce, maybe_redc4_x86_sse2,
		     cpu_feature_p(CPUFEAT_X86_SSE2));
#endif
#if CPUFAM_AMD64
  DISPATCH_PICK_COND(mpmont_reduce, maybe_redc4_amd64_sse2,
		     cpu_feature_p(CPUFEAT_X86_SSE2));
#endif
  DISPATCH_PICK_FALLBACK(mpmont_reduce, simple_redccore);
}

/* --- @redccore@ --- *
 *
 * Arguments:	@mpw *dv, *dvl@ = base and limit of source/destination
 *		@const mpw *av, *avl@ = base and limit of first multiplicand
 *		@const mpw *bv, *bvl@ = base and limit of second multiplicand
 *		@const mpw *mv@ = base of modulus %$m$%
 *		@size_t n@ = length of modulus
 *		@const mpw *mi@ = base of REDC coefficient %$m'$%
 *
 * Returns:	---
 *
 * Use:		Let %$a$% and %$b$% be the multiplicands.  Let %$w = a b$%.
 *		Store in %$d$% the value %$a b + (m' a b \bmod R) m$%.
 */

CPU_DISPATCH(static, (void), void, mulcore,
	     (mpw *dv, mpw *dvl, const mpw *av, const mpw *avl,
	      const mpw *bv, const mpw *bvl, const mpw *mv,
	      size_t n, const mpw *mi),
	     (dv, dvl, av, avl, bv, bvl, mv, n, mi),
	     pick_mulcore, simple_mulcore);

static void simple_mulcore(mpw *dv, mpw *dvl,
			   const mpw *av, const mpw *avl,
			   const mpw *bv, const mpw *bvl,
			   const mpw *mv, size_t n, const mpw *mi)
{
  mpw ai, b0, y, mi0 = *mi;
  const mpw *tv, *tvl;
  const mpw *mvl = mv + n;
  size_t i = 0;

  /* --- Initial setup --- */

  MPX_ZERO(dv, dvl);
  if (avl - av > bvl - bv) {
    tv = av; av = bv; bv = tv;
    tvl = avl; avl = bvl; bvl = tvl;
  }
  b0 = *bv;

  /* --- Multiply, until we run out of multiplicand --- */

  while (i < n && av < avl) {
    ai = *av++;
    y = MPW((*dv + ai*b0)*mi0);
    MPX_UMLAN(dv, dvl, bv, bvl, ai);
    MPX_UMLAN(dv, dvl, mv, mvl, y);
    dv++; i++;
  }

  /* --- Continue reducing until we run out of modulus --- */

  while (i < n) {
    y = MPW(*dv*mi0);
    MPX_UMLAN(dv, dvl, mv, mvl, y);
    dv++; i++;
  }
}

#define MAYBE_MUL4(impl)						\
  extern void mpxmont_mul4_##impl(mpw *dv,				\
				  const mpw *av, const mpw *bv,		\
				  const mpw *mv,			\
				  size_t n, const mpw *mi);		\
  static void maybe_mul4_##impl(mpw *dv, mpw *dvl,			\
			   const mpw *av, const mpw *avl,		\
			   const mpw *bv, const mpw *bvl,		\
			   const mpw *mv, size_t n, const mpw *mi)	\
  {									\
    size_t an = avl - av, bn = bvl - bv;				\
    if (n%4 || an != n || bn != n)					\
      simple_mulcore(dv, dvl, av, avl, bv, bvl, mv, n, mi);		\
    else {								\
      mpxmont_mul4_##impl(dv, av, bv, mv, n, mi);			\
      MPX_ZERO(dv + 2*n + 1, dvl);					\
    }									\
  }

#if CPUFAM_X86
  MAYBE_MUL4(x86_sse2)
#endif

#if CPUFAM_AMD64
  MAYBE_MUL4(amd64_sse2)
#endif

static mulcore__functype *pick_mulcore(void)
{
#if CPUFAM_X86
  DISPATCH_PICK_COND(mpmont_mul, maybe_mul4_x86_sse2,
		     cpu_feature_p(CPUFEAT_X86_SSE2));
#endif
#if CPUFAM_AMD64
  DISPATCH_PICK_COND(mpmont_mul, maybe_mul4_amd64_sse2,
		     cpu_feature_p(CPUFEAT_X86_SSE2));
#endif
  DISPATCH_PICK_FALLBACK(mpmont_mul, simple_mulcore);
}

/* --- @finish@ --- *
 *
 * Arguments:	@mpmont *mm@ = pointer to a Montgomery reduction context
 *		*mp *d@ = pointer to mostly-reduced operand
 *
 * Returns:	---
 *
 * Use:		Applies the finishing touches to Montgomery reduction.  The
 *		operand @d@ is a multiple of %$R%$ at this point, so it needs
 *		to be shifted down; the result might need a further
 *		subtraction to get it into the right interval; and we may
 *		need to do an additional subtraction if %$d$% is negative.
 */

static void finish(mpmont *mm, mp *d)
{
  mpw *dv = d->v, *dvl = d->vl;
  size_t n = mm->n;

  memmove(dv, dv + n, MPWS(dvl - (dv + n)));
  dvl -= n;

  if (MPX_UCMP(dv, dvl, >=, mm->m->v, mm->m->vl))
    mpx_usub(dv, dvl, dv, dvl, mm->m->v, mm->m->vl);

  if (d->f & MP_NEG) {
    mpx_usub(dv, dvl, mm->m->v, mm->m->vl, dv, dvl);
    d->f &= ~MP_NEG;
  }

  d->vl = dvl;
  MP_SHRINK(d);
}

#endif

/*----- Reduction and multiplication --------------------------------------*/

/* --- @mpmont_create@ --- *
 *
 * Arguments:	@mpmont *mm@ = pointer to Montgomery reduction context
 *		@mp *m@ = modulus to use
 *
 * Returns:	Zero on success, nonzero on error.
 *
 * Use:		Initializes a Montgomery reduction context ready for use.
 *		The argument @m@ must be a positive odd integer.
 */

#ifdef MPMONT_DISABLE

int mpmont_create(mpmont *mm, mp *m)
{
  mp_shrink(m);
  mm->m = MP_COPY(m);
  mm->r = MP_ONE;
  mm->r2 = MP_ONE;
  mm->mi = MP_ONE;
  return (0);
}

#else

int mpmont_create(mpmont *mm, mp *m)
{
  size_t n = MP_LEN(m);
  mp *r2 = mp_new(2 * n + 1, 0);
  mp r;

  /* --- Take a copy of the modulus --- */

 if (!MP_POSP(m) || !MP_ODDP(m))
   return (-1);
  mm->m = MP_COPY(m);

  /* --- Determine %$R^2$% --- */

  mm->n = n;
  MPX_ZERO(r2->v, r2->vl - 1);
  r2->vl[-1] = 1;

  /* --- Find the magic value @mi@ --- */

  mp_build(&r, r2->v + n, r2->vl);
  mm->mi = mp_modinv(MP_NEW, m, &r);
  mm->mi = mp_sub(mm->mi, &r, mm->mi);
  MP_ENSURE(mm->mi, n);

  /* --- Discover the values %$R \bmod m$% and %$R^2 \bmod m$% --- */

  mm->r2 = MP_NEW;
  mp_div(0, &mm->r2, r2, m);
  mm->r = mpmont_reduce(mm, MP_NEW, mm->r2);
  MP_DROP(r2);
  return (0);
}

#endif

/* --- @mpmont_destroy@ --- *
 *
 * Arguments:	@mpmont *mm@ = pointer to a Montgomery reduction context
 *
 * Returns:	---
 *
 * Use:		Disposes of a context when it's no longer of any use to
 *		anyone.
 */

void mpmont_destroy(mpmont *mm)
{
  MP_DROP(mm->m);
  MP_DROP(mm->r);
  MP_DROP(mm->r2);
  MP_DROP(mm->mi);
}

/* --- @mpmont_reduce@ --- *
 *
 * Arguments:	@mpmont *mm@ = pointer to Montgomery reduction context
 *		@mp *d@ = destination
 *		@mp *a@ = source, assumed positive
 *
 * Returns:	Result, %$a R^{-1} \bmod m$%.
 */

#ifdef MPMONT_DISABLE

mp *mpmont_reduce(mpmont *mm, mp *d, mp *a)
{
  mp_div(0, &d, a, mm->m);
  return (d);
}

#else

mp *mpmont_reduce(mpmont *mm, mp *d, mp *a)
{
  size_t n = mm->n;

  /* --- Check for serious Karatsuba reduction --- */

  if (n > MPMONT_KTHRESH) {
    mp al;
    mpw *vl;
    mp *u;

    if (MP_LEN(a) >= n) vl = a->v + n;
    else vl = a->vl;
    mp_build(&al, a->v, vl);
    u = mp_mul(MP_NEW, &al, mm->mi);
    if (MP_LEN(u) > n) u->vl = u->v + n;
    u = mp_mul(u, u, mm->m);
    d = mp_add(d, a, u);
    MP_ENSURE(d, n);
    mp_drop(u);
  }

  /* --- Otherwise do it the hard way --- */

  else {
    a = MP_COPY(a);
    if (d) MP_DROP(d);
    d = a;
    MP_DEST(d, 2*mm->n + 1, a->f);
    redccore(d->v, d->vl, mm->m->v, mm->n, mm->mi->v);
  }

  /* --- Wrap everything up --- */

  finish(mm, d);
  return (d);
}

#endif

/* --- @mpmont_mul@ --- *
 *
 * Arguments:	@mpmont *mm@ = pointer to Montgomery reduction context
 *		@mp *d@ = destination
 *		@mp *a, *b@ = sources, assumed positive
 *
 * Returns:	Result, %$a b R^{-1} \bmod m$%.
 */

#ifdef MPMONT_DISABLE

mp *mpmont_mul(mpmont *mm, mp *d, mp *a, mp *b)
{
  d = mp_mul(d, a, b);
  mp_div(0, &d, d, mm->m);
  return (d);
}

#else

mp *mpmont_mul(mpmont *mm, mp *d, mp *a, mp *b)
{
  size_t n = mm->n;

  if (n > MPMONT_KTHRESH) {
    d = mp_mul(d, a, b);
    d = mpmont_reduce(mm, d, d);
  } else {
    a = MP_COPY(a); b = MP_COPY(b);
    MP_DEST(d, 2*n + 1, a->f | b->f | MP_UNDEF);
    mulcore(d->v, d->vl, a->v, a->vl, b->v, b->vl,
	    mm->m->v, mm->n, mm->mi->v);
    d->f = ((a->f | b->f) & MP_BURN) | ((a->f ^ b->f) & MP_NEG);
    finish(mm, d);
    MP_DROP(a); MP_DROP(b);
  }

  return (d);
}

#endif

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

#ifdef TEST_RIG

static int tcreate(dstr *v)
{
  mp *m = *(mp **)v[0].buf;
  mp *mi = *(mp **)v[1].buf;
  mp *r = *(mp **)v[2].buf;
  mp *r2 = *(mp **)v[3].buf;

  mpmont mm;
  int ok = 1;

  mpmont_create(&mm, m);

  if (mm.mi->v[0] != mi->v[0]) {
    fprintf(stderr, "\n*** bad mi: found %lu, expected %lu",
	    (unsigned long)mm.mi->v[0], (unsigned long)mi->v[0]);
    fputs("\nm = ", stderr); mp_writefile(m, stderr, 10);
    fputc('\n', stderr);
    ok = 0;
  }

  if (!MP_EQ(mm.r, r)) {
    fputs("\n*** bad r", stderr);
    fputs("\nm = ", stderr); mp_writefile(m, stderr, 10);
    fputs("\nexpected ", stderr); mp_writefile(r, stderr, 10);
    fputs("\n	found ", stderr); mp_writefile(mm.r, stderr, 10);
    fputc('\n', stderr);
    ok = 0;
  }

  if (!MP_EQ(mm.r2, r2)) {
    fputs("\n*** bad r2", stderr);
    fputs("\nm = ", stderr); mp_writefile(m, stderr, 10);
    fputs("\nexpected ", stderr); mp_writefile(r2, stderr, 10);
    fputs("\n	found ", stderr); mp_writefile(mm.r2, stderr, 10);
    fputc('\n', stderr);
    ok = 0;
  }

  MP_DROP(m);
  MP_DROP(mi);
  MP_DROP(r);
  MP_DROP(r2);
  mpmont_destroy(&mm);
  assert(mparena_count(MPARENA_GLOBAL) == 0);
  return (ok);
}

static int tmul(dstr *v)
{
  mp *m = *(mp **)v[0].buf;
  mp *a = *(mp **)v[1].buf;
  mp *b = *(mp **)v[2].buf;
  mp *r = *(mp **)v[3].buf;
  int ok = 1;

  mpmont mm;
  mpmont_create(&mm, m);

  {
    mp *qr = mp_mul(MP_NEW, a, b);
    mp_div(0, &qr, qr, m);

    if (!MP_EQ(qr, r)) {
      fputs("\n*** classical modmul failed", stderr);
      fputs("\n m = ", stderr); mp_writefile(m, stderr, 10);
      fputs("\n a = ", stderr); mp_writefile(a, stderr, 10);
      fputs("\n b = ", stderr); mp_writefile(b, stderr, 10);
      fputs("\n r = ", stderr); mp_writefile(r, stderr, 10);
      fputs("\nqr = ", stderr); mp_writefile(qr, stderr, 10);
      fputc('\n', stderr);
      ok = 0;
    }

    mp_drop(qr);
  }

  {
    mp *ar = mpmont_mul(&mm, MP_NEW, a, mm.r2);
    mp *br = mpmont_mul(&mm, MP_NEW, b, mm.r2);
    mp *mr = mpmont_mul(&mm, MP_NEW, ar, br);
    mr = mpmont_reduce(&mm, mr, mr);
    if (!MP_EQ(mr, r)) {
      fputs("\n*** montgomery modmul failed", stderr);
      fputs("\n m = ", stderr); mp_writefile(m, stderr, 10);
      fputs("\n a = ", stderr); mp_writefile(a, stderr, 10);
      fputs("\n b = ", stderr); mp_writefile(b, stderr, 10);
      fputs("\n r = ", stderr); mp_writefile(r, stderr, 10);
      fputs("\nmr = ", stderr); mp_writefile(mr, stderr, 10);
      fputc('\n', stderr);
      ok = 0;
    }
    MP_DROP(ar); MP_DROP(br);
    mp_drop(mr);
  }


  MP_DROP(m);
  MP_DROP(a);
  MP_DROP(b);
  MP_DROP(r);
  mpmont_destroy(&mm);
  assert(mparena_count(MPARENA_GLOBAL) == 0);
  return ok;
}

static test_chunk tests[] = {
  { "create", tcreate, { &type_mp, &type_mp, &type_mp, &type_mp, 0 } },
  { "mul", tmul, { &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/mpmont");
  return (0);
}

#endif

/*----- That's all, folks -------------------------------------------------*/
Commit	Line	Data
d3409d5e	1	/* --c--
d3409d5e	2	*
d3409d5e	3	* Montgomery reduction
	4	*
	5	* (c) 1999 Straylight/Edgeware
	6	*/
	7
45c0fd36	8	/----- Licensing notice --------------------------------------------------
d3409d5e	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	*
d3409d5e	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	*
d3409d5e	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
d3409d5e	28	/----- Header files ------------------------------------------------------/
d3409d5e	29
444083ae MW	30	#include "config.h"
444083ae MW	31	#include "dispatch.h"
d3409d5e	32	#include "mp.h"
	33	#include "mpmont.h"
	34
52e4b041	35	/----- Tweakables --------------------------------------------------------/
	36
	37	/* --- @MPMONT_DISABLE@ --- *
	38	*
	39	* Replace all the clever Montgomery reduction with good old-fashioned long
	40	* division.
	41	*/
	42
	43	/* #define MPMONT_DISABLE */
	44
d6b9dc04 MW	45	#define MPMONT_KTHRESH (16*MPK_THRESH)
d6b9dc04 MW	46
0e70bd46 MW	47	/----- Low-level implementation ------------------------------------------/
	48
	49	#ifndef MPMONT_DISABLE
	50
	51	/* --- @redccore@ --- *
	52	*
	53	* Arguments: @mpw dv, dvl@ = base and limit of source/destination
	54	* @const mpw *mv@ = base of modulus %$m$%
	55	* @size_t n@ = length of modulus
	56	* @const mpw *mi@ = base of REDC coefficient %$m'$%
	57	*
	58	* Returns: ---
	59	*
	60	* Use: Let %$a$% be the input operand. Store in %$d$% the value
	61	* %$a + (m' a \bmod R) m$%. The destination has space for at
	62	* least %$2 n + 1$% words of result.
	63	*/
	64
444083ae MW	65	CPU_DISPATCH(static, (void), void, redccore,
	66	(mpw dv, mpw dvl, const mpw mv, size_t n, const mpw mi),
	67	(dv, dvl, mv, n, mi), pick_redccore, simple_redccore);
	68
	69	static void simple_redccore(mpw dv, mpw dvl, const mpw *mv,
	70	size_t n, const mpw *mi)
0e70bd46 MW	71	{
	72	mpw mi0 = *mi;
	73	size_t i;
	74
	75	for (i = 0; i < n; i++) {
	76	MPX_UMLAN(dv, dvl, mv, mv + n, MPW(dvmi0));
	77	dv++;
	78	}
	79	}
	80
444083ae MW	81	#define MAYBE_REDC4(impl) \
	82	extern void mpxmont_redc4_##impl(mpw dv, mpw dvl, const mpw *mv, \
	83	size_t n, const mpw *mi); \
	84	static void maybe_redc4_##impl(mpw dv, mpw dvl, const mpw *mv, \
	85	size_t n, const mpw *mi) \
	86	{ \
	87	if (n%4) simple_redccore(dv, dvl, mv, n, mi); \
	88	else mpxmont_redc4_##impl(dv, dvl, mv, n, mi); \
	89	}
	90
	91	#if CPUFAM_X86
	92	MAYBE_REDC4(x86_sse2)
	93	#endif
	94
3119b3ae MW	95	#if CPUFAM_AMD64
	96	MAYBE_REDC4(amd64_sse2)
	97	#endif
	98
444083ae MW	99	static redccore__functype *pick_redccore(void)
	100	{
	101	#if CPUFAM_X86
	102	DISPATCH_PICK_COND(mpmont_reduce, maybe_redc4_x86_sse2,
	103	cpu_feature_p(CPUFEAT_X86_SSE2));
	104	#endif
3119b3ae MW	105	#if CPUFAM_AMD64
	106	DISPATCH_PICK_COND(mpmont_reduce, maybe_redc4_amd64_sse2,
	107	cpu_feature_p(CPUFEAT_X86_SSE2));
	108	#endif
444083ae MW	109	DISPATCH_PICK_FALLBACK(mpmont_reduce, simple_redccore);
	110	}
	111
0e70bd46 MW	112	/* --- @redccore@ --- *
	113	*
	114	* Arguments: @mpw dv, dvl@ = base and limit of source/destination
	115	* @const mpw av, avl@ = base and limit of first multiplicand
	116	* @const mpw bv, bvl@ = base and limit of second multiplicand
	117	* @const mpw *mv@ = base of modulus %$m$%
	118	* @size_t n@ = length of modulus
	119	* @const mpw *mi@ = base of REDC coefficient %$m'$%
	120	*
	121	* Returns: ---
	122	*
	123	* Use: Let %$a$% and %$b$% be the multiplicands. Let %$w = a b$%.
	124	* Store in %$d$% the value %$a b + (m' a b \bmod R) m$%.
	125	*/
	126
444083ae MW	127	CPU_DISPATCH(static, (void), void, mulcore,
	128	(mpw dv, mpw dvl, const mpw av, const mpw avl,
	129	const mpw bv, const mpw bvl, const mpw *mv,
	130	size_t n, const mpw *mi),
	131	(dv, dvl, av, avl, bv, bvl, mv, n, mi),
	132	pick_mulcore, simple_mulcore);
	133
	134	static void simple_mulcore(mpw dv, mpw dvl,
	135	const mpw av, const mpw avl,
	136	const mpw bv, const mpw bvl,
	137	const mpw mv, size_t n, const mpw mi)
0e70bd46 MW	138	{
	139	mpw ai, b0, y, mi0 = *mi;
	140	const mpw tv, tvl;
	141	const mpw *mvl = mv + n;
	142	size_t i = 0;
	143
	144	/* --- Initial setup --- */
	145
	146	MPX_ZERO(dv, dvl);
	147	if (avl - av > bvl - bv) {
	148	tv = av; av = bv; bv = tv;
	149	tvl = avl; avl = bvl; bvl = tvl;
	150	}
	151	b0 = *bv;
	152
	153	/* --- Multiply, until we run out of multiplicand --- */
	154
	155	while (i < n && av < avl) {
	156	ai = *av++;
	157	y = MPW((dv + aib0)*mi0);
	158	MPX_UMLAN(dv, dvl, bv, bvl, ai);
	159	MPX_UMLAN(dv, dvl, mv, mvl, y);
	160	dv++; i++;
	161	}
	162
	163	/* --- Continue reducing until we run out of modulus --- */
	164
	165	while (i < n) {
	166	y = MPW(dvmi0);
	167	MPX_UMLAN(dv, dvl, mv, mvl, y);
	168	dv++; i++;
	169	}
	170	}
	171
444083ae MW	172	#define MAYBE_MUL4(impl) \
	173	extern void mpxmont_mul4_##impl(mpw *dv, \
	174	const mpw av, const mpw bv, \
	175	const mpw *mv, \
	176	size_t n, const mpw *mi); \
	177	static void maybe_mul4_##impl(mpw dv, mpw dvl, \
	178	const mpw av, const mpw avl, \
	179	const mpw bv, const mpw bvl, \
	180	const mpw mv, size_t n, const mpw mi) \
	181	{ \
	182	size_t an = avl - av, bn = bvl - bv; \
	183	if (n%4 \|\| an != n \|\| bn != n) \
	184	simple_mulcore(dv, dvl, av, avl, bv, bvl, mv, n, mi); \
	185	else { \
	186	mpxmont_mul4_##impl(dv, av, bv, mv, n, mi); \
	187	MPX_ZERO(dv + 2*n + 1, dvl); \
	188	} \
	189	}
	190
	191	#if CPUFAM_X86
	192	MAYBE_MUL4(x86_sse2)
	193	#endif
	194
3119b3ae MW	195	#if CPUFAM_AMD64
	196	MAYBE_MUL4(amd64_sse2)
	197	#endif
	198
444083ae MW	199	static mulcore__functype *pick_mulcore(void)
	200	{
	201	#if CPUFAM_X86
	202	DISPATCH_PICK_COND(mpmont_mul, maybe_mul4_x86_sse2,
	203	cpu_feature_p(CPUFEAT_X86_SSE2));
	204	#endif
3119b3ae MW	205	#if CPUFAM_AMD64
	206	DISPATCH_PICK_COND(mpmont_mul, maybe_mul4_amd64_sse2,
	207	cpu_feature_p(CPUFEAT_X86_SSE2));
	208	#endif
444083ae MW	209	DISPATCH_PICK_FALLBACK(mpmont_mul, simple_mulcore);
	210	}
	211
0e70bd46 MW	212	/* --- @finish@ --- *
	213	*
	214	* Arguments: @mpmont *mm@ = pointer to a Montgomery reduction context
	215	* mp d@ = pointer to mostly-reduced operand
	216	*
	217	* Returns: ---
	218	*
	219	* Use: Applies the finishing touches to Montgomery reduction. The
	220	* operand @d@ is a multiple of %$R%$ at this point, so it needs
	221	* to be shifted down; the result might need a further
	222	* subtraction to get it into the right interval; and we may
	223	* need to do an additional subtraction if %$d$% is negative.
	224	*/
	225
	226	static void finish(mpmont mm, mp d)
	227	{
	228	mpw dv = d->v, dvl = d->vl;
	229	size_t n = mm->n;
	230
	231	memmove(dv, dv + n, MPWS(dvl - (dv + n)));
	232	dvl -= n;
	233
	234	if (MPX_UCMP(dv, dvl, >=, mm->m->v, mm->m->vl))
	235	mpx_usub(dv, dvl, dv, dvl, mm->m->v, mm->m->vl);
	236
	237	if (d->f & MP_NEG) {
	238	mpx_usub(dv, dvl, mm->m->v, mm->m->vl, dv, dvl);
	239	d->f &= ~MP_NEG;
	240	}
	241
	242	d->vl = dvl;
	243	MP_SHRINK(d);
	244	}
	245
	246	#endif
	247
4640a0dd	248	/----- Reduction and multiplication --------------------------------------/
d3409d5e	249
	250	/* --- @mpmont_create@ --- *
	251	*
	252	* Arguments: @mpmont *mm@ = pointer to Montgomery reduction context
	253	* @mp *m@ = modulus to use
	254	*
f4535c64	255	* Returns: Zero on success, nonzero on error.
d3409d5e	256	*
d3409d5e	257	* Use: Initializes a Montgomery reduction context ready for use.
ef5f4810	258	* The argument @m@ must be a positive odd integer.
d3409d5e	259	*/
d3409d5e	260
52e4b041	261	#ifdef MPMONT_DISABLE
52e4b041	262
f4535c64	263	int mpmont_create(mpmont mm, mp m)
52e4b041	264	{
	265	mp_shrink(m);
	266	mm->m = MP_COPY(m);
	267	mm->r = MP_ONE;
	268	mm->r2 = MP_ONE;
f5f35081	269	mm->mi = MP_ONE;
f4535c64	270	return (0);
52e4b041	271	}
	272
	273	#else
	274
f4535c64	275	int mpmont_create(mpmont mm, mp m)
d3409d5e	276	{
f5f35081	277	size_t n = MP_LEN(m);
d34decd2	278	mp r2 = mp_new(2 n + 1, 0);
f5f35081	279	mp r;
f5f35081	280
d3409d5e	281	/* --- Take a copy of the modulus --- */
d3409d5e	282
f4535c64	283	if (!MP_POSP(m) \|\| !MP_ODDP(m))
f4535c64	284	return (-1);
d3409d5e	285	mm->m = MP_COPY(m);
d3409d5e	286
f5f35081	287	/* --- Determine %$R^2$% --- */
d3409d5e	288
f5f35081	289	mm->n = n;
	290	MPX_ZERO(r2->v, r2->vl - 1);
	291	r2->vl[-1] = 1;
d3409d5e	292
f5f35081	293	/* --- Find the magic value @mi@ --- */
	294
	295	mp_build(&r, r2->v + n, r2->vl);
b817bfc6	296	mm->mi = mp_modinv(MP_NEW, m, &r);
f5f35081	297	mm->mi = mp_sub(mm->mi, &r, mm->mi);
362c3d18	298	MP_ENSURE(mm->mi, n);
d3409d5e	299
	300	/* --- Discover the values %$R \bmod m$% and %$R^2 \bmod m$% --- */
	301
f5f35081	302	mm->r2 = MP_NEW;
	303	mp_div(0, &mm->r2, r2, m);
	304	mm->r = mpmont_reduce(mm, MP_NEW, mm->r2);
	305	MP_DROP(r2);
f4535c64	306	return (0);
d3409d5e	307	}
d3409d5e	308
52e4b041	309	#endif
52e4b041	310
d3409d5e	311	/* --- @mpmont_destroy@ --- *
	312	*
	313	* Arguments: @mpmont *mm@ = pointer to a Montgomery reduction context
	314	*
	315	* Returns: ---
	316	*
	317	* Use: Disposes of a context when it's no longer of any use to
	318	* anyone.
	319	*/
	320
	321	void mpmont_destroy(mpmont *mm)
	322	{
	323	MP_DROP(mm->m);
	324	MP_DROP(mm->r);
	325	MP_DROP(mm->r2);
f5f35081	326	MP_DROP(mm->mi);
d3409d5e	327	}
	328
	329	/* --- @mpmont_reduce@ --- *
	330	*
	331	* Arguments: @mpmont *mm@ = pointer to Montgomery reduction context
	332	* @mp *d@ = destination
ef5f4810	333	* @mp *a@ = source, assumed positive
d3409d5e	334	*
	335	* Returns: Result, %$a R^{-1} \bmod m$%.
	336	*/
	337
52e4b041	338	#ifdef MPMONT_DISABLE
52e4b041	339
ef5f4810	340	mp mpmont_reduce(mpmont mm, mp d, mp a)
52e4b041	341	{
	342	mp_div(0, &d, a, mm->m);
	343	return (d);
	344	}
	345
	346	#else
	347
ef5f4810	348	mp mpmont_reduce(mpmont mm, mp d, mp a)
d3409d5e	349	{
f5f35081	350	size_t n = mm->n;
	351
	352	/* --- Check for serious Karatsuba reduction --- */
	353
d6b9dc04	354	if (n > MPMONT_KTHRESH) {
f5f35081	355	mp al;
	356	mpw *vl;
	357	mp *u;
	358
0e70bd46 MW	359	if (MP_LEN(a) >= n) vl = a->v + n;
0e70bd46 MW	360	else vl = a->vl;
f5f35081	361	mp_build(&al, a->v, vl);
f5f35081	362	u = mp_mul(MP_NEW, &al, mm->mi);
0e70bd46	363	if (MP_LEN(u) > n) u->vl = u->v + n;
f5f35081	364	u = mp_mul(u, u, mm->m);
f5f35081	365	d = mp_add(d, a, u);
366e9701	366	MP_ENSURE(d, n);
f5f35081	367	mp_drop(u);
f5f35081	368	}
d3409d5e	369
f5f35081	370	/* --- Otherwise do it the hard way --- */
d3409d5e	371
d3409d5e	372	else {
d34decd2	373	a = MP_COPY(a);
0e70bd46	374	if (d) MP_DROP(d);
d34decd2	375	d = a;
0e70bd46 MW	376	MP_DEST(d, 2*mm->n + 1, a->f);
0e70bd46 MW	377	redccore(d->v, d->vl, mm->m->v, mm->n, mm->mi->v);
d3409d5e	378	}
d3409d5e	379
f5f35081	380	/* --- Wrap everything up --- */
d3409d5e	381
0e70bd46	382	finish(mm, d);
d3409d5e	383	return (d);
	384	}
	385
52e4b041	386	#endif
52e4b041	387
d3409d5e	388	/* --- @mpmont_mul@ --- *
	389	*
	390	* Arguments: @mpmont *mm@ = pointer to Montgomery reduction context
	391	* @mp *d@ = destination
ef5f4810	392	* @mp a, b@ = sources, assumed positive
d3409d5e	393	*
	394	* Returns: Result, %$a b R^{-1} \bmod m$%.
	395	*/
	396
52e4b041	397	#ifdef MPMONT_DISABLE
52e4b041	398
ef5f4810	399	mp mpmont_mul(mpmont mm, mp d, mp a, mp *b)
52e4b041	400	{
	401	d = mp_mul(d, a, b);
	402	mp_div(0, &d, d, mm->m);
	403	return (d);
	404	}
	405
	406	#else
	407
ef5f4810	408	mp mpmont_mul(mpmont mm, mp d, mp a, mp *b)
d3409d5e	409	{
444083ae MW	410	size_t n = mm->n;
444083ae MW	411
d6b9dc04	412	if (n > MPMONT_KTHRESH) {
ef5f4810	413	d = mp_mul(d, a, b);
	414	d = mpmont_reduce(mm, d, d);
	415	} else {
444083ae MW	416	a = MP_COPY(a); b = MP_COPY(b);
444083ae MW	417	MP_DEST(d, 2*n + 1, a->f \| b->f \| MP_UNDEF);
0e70bd46 MW	418	mulcore(d->v, d->vl, a->v, a->vl, b->v, b->vl,
	419	mm->m->v, mm->n, mm->mi->v);
	420	d->f = ((a->f \| b->f) & MP_BURN) \| ((a->f ^ b->f) & MP_NEG);
	421	finish(mm, d);
	422	MP_DROP(a); MP_DROP(b);
d3409d5e	423	}
d3409d5e	424
d3409d5e	425	return (d);
	426	}
	427
52e4b041	428	#endif
52e4b041	429
d3409d5e	430	/----- Test rig ----------------------------------------------------------/
	431
	432	#ifdef TEST_RIG
	433
	434	static int tcreate(dstr *v)
	435	{
	436	mp m = (mp **)v[0].buf;
	437	mp mi = (mp **)v[1].buf;
	438	mp r = (mp **)v[2].buf;
	439	mp r2 = (mp **)v[3].buf;
	440
	441	mpmont mm;
	442	int ok = 1;
	443
	444	mpmont_create(&mm, m);
	445
f5f35081	446	if (mm.mi->v[0] != mi->v[0]) {
d3409d5e	447	fprintf(stderr, "\n*** bad mi: found %lu, expected %lu",
f5f35081	448	(unsigned long)mm.mi->v[0], (unsigned long)mi->v[0]);
d3409d5e	449	fputs("\nm = ", stderr); mp_writefile(m, stderr, 10);
	450	fputc('\n', stderr);
	451	ok = 0;
	452	}
	453
032099d1	454	if (!MP_EQ(mm.r, r)) {
d3409d5e	455	fputs("\n*** bad r", stderr);
	456	fputs("\nm = ", stderr); mp_writefile(m, stderr, 10);
	457	fputs("\nexpected ", stderr); mp_writefile(r, stderr, 10);
45c0fd36	458	fputs("\n found ", stderr); mp_writefile(mm.r, stderr, 10);
d3409d5e	459	fputc('\n', stderr);
	460	ok = 0;
	461	}
	462
032099d1	463	if (!MP_EQ(mm.r2, r2)) {
d3409d5e	464	fputs("\n*** bad r2", stderr);
	465	fputs("\nm = ", stderr); mp_writefile(m, stderr, 10);
	466	fputs("\nexpected ", stderr); mp_writefile(r2, stderr, 10);
45c0fd36	467	fputs("\n found ", stderr); mp_writefile(mm.r2, stderr, 10);
d3409d5e	468	fputc('\n', stderr);
	469	ok = 0;
	470	}
	471
	472	MP_DROP(m);
	473	MP_DROP(mi);
	474	MP_DROP(r);
	475	MP_DROP(r2);
	476	mpmont_destroy(&mm);
ef5f4810	477	assert(mparena_count(MPARENA_GLOBAL) == 0);
d3409d5e	478	return (ok);
	479	}
	480
	481	static int tmul(dstr *v)
	482	{
	483	mp m = (mp **)v[0].buf;
	484	mp a = (mp **)v[1].buf;
	485	mp b = (mp **)v[2].buf;
	486	mp r = (mp **)v[3].buf;
d3409d5e	487	int ok = 1;
	488
	489	mpmont mm;
	490	mpmont_create(&mm, m);
	491
	492	{
79a34029	493	mp *qr = mp_mul(MP_NEW, a, b);
	494	mp_div(0, &qr, qr, m);
	495
032099d1	496	if (!MP_EQ(qr, r)) {
79a34029	497	fputs("\n*** classical modmul failed", stderr);
	498	fputs("\n m = ", stderr); mp_writefile(m, stderr, 10);
	499	fputs("\n a = ", stderr); mp_writefile(a, stderr, 10);
	500	fputs("\n b = ", stderr); mp_writefile(b, stderr, 10);
	501	fputs("\n r = ", stderr); mp_writefile(r, stderr, 10);
	502	fputs("\nqr = ", stderr); mp_writefile(qr, stderr, 10);
	503	fputc('\n', stderr);
	504	ok = 0;
	505	}
	506
	507	mp_drop(qr);
	508	}
	509
	510	{
d3409d5e	511	mp *ar = mpmont_mul(&mm, MP_NEW, a, mm.r2);
d3409d5e	512	mp *br = mpmont_mul(&mm, MP_NEW, b, mm.r2);
79a34029	513	mp *mr = mpmont_mul(&mm, MP_NEW, ar, br);
d3409d5e	514	mr = mpmont_reduce(&mm, mr, mr);
032099d1	515	if (!MP_EQ(mr, r)) {
79a34029	516	fputs("\n*** montgomery modmul failed", stderr);
	517	fputs("\n m = ", stderr); mp_writefile(m, stderr, 10);
	518	fputs("\n a = ", stderr); mp_writefile(a, stderr, 10);
	519	fputs("\n b = ", stderr); mp_writefile(b, stderr, 10);
	520	fputs("\n r = ", stderr); mp_writefile(r, stderr, 10);
	521	fputs("\nmr = ", stderr); mp_writefile(mr, stderr, 10);
	522	fputc('\n', stderr);
	523	ok = 0;
	524	}
d3409d5e	525	MP_DROP(ar); MP_DROP(br);
79a34029	526	mp_drop(mr);
d3409d5e	527	}
d3409d5e	528
d3409d5e	529
	530	MP_DROP(m);
	531	MP_DROP(a);
	532	MP_DROP(b);
	533	MP_DROP(r);
d3409d5e	534	mpmont_destroy(&mm);
ef5f4810	535	assert(mparena_count(MPARENA_GLOBAL) == 0);
d3409d5e	536	return ok;
	537	}
	538
d3409d5e	539	static test_chunk tests[] = {
ef5f4810	540	{ "create", tcreate, { &type_mp, &type_mp, &type_mp, &type_mp, 0 } },
ef5f4810	541	{ "mul", tmul, { &type_mp, &type_mp, &type_mp, &type_mp, 0 } },
d3409d5e	542	{ 0, 0, { 0 } },
	543	};
	544
	545	int main(int argc, char *argv[])
	546	{
	547	sub_init();
0f00dc4c	548	test_run(argc, argv, tests, SRCDIR "/t/mpmont");
d3409d5e	549	return (0);
	550	}
	551
	552	#endif
	553
	554	/----- That's all, folks -------------------------------------------------/