[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:	@const 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(const 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:	@const 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(const mpmont *mm, mp *d, mp *a)
{
  mp_div(0, &d, a, mm->m);
  return (d);
}

#else

mp *mpmont_reduce(const 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:	@const 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(const 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(const 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@ --- *
0e70bd46 MW	213	*
295f4f90 MW	214	* Arguments: @const mpmont *mm@ = pointer to a Montgomery reduction
295f4f90 MW	215	* context
0e70bd46 MW	216	* mp d@ = pointer to mostly-reduced operand
	217	*
	218	* Returns: ---
	219	*
	220	* Use: Applies the finishing touches to Montgomery reduction. The
	221	* operand @d@ is a multiple of %$R%$ at this point, so it needs
	222	* to be shifted down; the result might need a further
	223	* subtraction to get it into the right interval; and we may
	224	* need to do an additional subtraction if %$d$% is negative.
	225	*/
	226
295f4f90	227	static void finish(const mpmont mm, mp d)
0e70bd46 MW	228	{
	229	mpw dv = d->v, dvl = d->vl;
	230	size_t n = mm->n;
	231
	232	memmove(dv, dv + n, MPWS(dvl - (dv + n)));
	233	dvl -= n;
	234
	235	if (MPX_UCMP(dv, dvl, >=, mm->m->v, mm->m->vl))
	236	mpx_usub(dv, dvl, dv, dvl, mm->m->v, mm->m->vl);
	237
	238	if (d->f & MP_NEG) {
	239	mpx_usub(dv, dvl, mm->m->v, mm->m->vl, dv, dvl);
	240	d->f &= ~MP_NEG;
	241	}
	242
	243	d->vl = dvl;
	244	MP_SHRINK(d);
	245	}
	246
	247	#endif
	248
4640a0dd	249	/----- Reduction and multiplication --------------------------------------/
d3409d5e	250
	251	/* --- @mpmont_create@ --- *
	252	*
	253	* Arguments: @mpmont *mm@ = pointer to Montgomery reduction context
	254	* @mp *m@ = modulus to use
	255	*
f4535c64	256	* Returns: Zero on success, nonzero on error.
d3409d5e	257	*
d3409d5e	258	* Use: Initializes a Montgomery reduction context ready for use.
ef5f4810	259	* The argument @m@ must be a positive odd integer.
d3409d5e	260	*/
d3409d5e	261
52e4b041	262	#ifdef MPMONT_DISABLE
52e4b041	263
f4535c64	264	int mpmont_create(mpmont mm, mp m)
52e4b041	265	{
	266	mp_shrink(m);
	267	mm->m = MP_COPY(m);
	268	mm->r = MP_ONE;
	269	mm->r2 = MP_ONE;
f5f35081	270	mm->mi = MP_ONE;
f4535c64	271	return (0);
52e4b041	272	}
	273
	274	#else
	275
f4535c64	276	int mpmont_create(mpmont mm, mp m)
d3409d5e	277	{
f5f35081	278	size_t n = MP_LEN(m);
d34decd2	279	mp r2 = mp_new(2 n + 1, 0);
f5f35081	280	mp r;
f5f35081	281
d3409d5e	282	/* --- Take a copy of the modulus --- */
d3409d5e	283
f4535c64	284	if (!MP_POSP(m) \|\| !MP_ODDP(m))
f4535c64	285	return (-1);
d3409d5e	286	mm->m = MP_COPY(m);
d3409d5e	287
f5f35081	288	/* --- Determine %$R^2$% --- */
d3409d5e	289
f5f35081	290	mm->n = n;
	291	MPX_ZERO(r2->v, r2->vl - 1);
	292	r2->vl[-1] = 1;
d3409d5e	293
f5f35081	294	/* --- Find the magic value @mi@ --- */
	295
	296	mp_build(&r, r2->v + n, r2->vl);
b817bfc6	297	mm->mi = mp_modinv(MP_NEW, m, &r);
f5f35081	298	mm->mi = mp_sub(mm->mi, &r, mm->mi);
362c3d18	299	MP_ENSURE(mm->mi, n);
d3409d5e	300
	301	/* --- Discover the values %$R \bmod m$% and %$R^2 \bmod m$% --- */
	302
f5f35081	303	mm->r2 = MP_NEW;
	304	mp_div(0, &mm->r2, r2, m);
	305	mm->r = mpmont_reduce(mm, MP_NEW, mm->r2);
	306	MP_DROP(r2);
f4535c64	307	return (0);
d3409d5e	308	}
d3409d5e	309
52e4b041	310	#endif
52e4b041	311
d3409d5e	312	/* --- @mpmont_destroy@ --- *
	313	*
	314	* Arguments: @mpmont *mm@ = pointer to a Montgomery reduction context
	315	*
	316	* Returns: ---
	317	*
	318	* Use: Disposes of a context when it's no longer of any use to
	319	* anyone.
	320	*/
	321
	322	void mpmont_destroy(mpmont *mm)
	323	{
	324	MP_DROP(mm->m);
	325	MP_DROP(mm->r);
	326	MP_DROP(mm->r2);
f5f35081	327	MP_DROP(mm->mi);
d3409d5e	328	}
	329
	330	/* --- @mpmont_reduce@ --- *
	331	*
295f4f90	332	* Arguments: @const mpmont *mm@ = pointer to Montgomery reduction context
d3409d5e	333	* @mp *d@ = destination
ef5f4810	334	* @mp *a@ = source, assumed positive
d3409d5e	335	*
	336	* Returns: Result, %$a R^{-1} \bmod m$%.
	337	*/
	338
52e4b041	339	#ifdef MPMONT_DISABLE
52e4b041	340
295f4f90	341	mp mpmont_reduce(const mpmont mm, mp d, mp a)
52e4b041	342	{
	343	mp_div(0, &d, a, mm->m);
	344	return (d);
	345	}
	346
	347	#else
	348
295f4f90	349	mp mpmont_reduce(const mpmont mm, mp d, mp a)
d3409d5e	350	{
f5f35081	351	size_t n = mm->n;
	352
	353	/* --- Check for serious Karatsuba reduction --- */
	354
d6b9dc04	355	if (n > MPMONT_KTHRESH) {
f5f35081	356	mp al;
	357	mpw *vl;
	358	mp *u;
	359
0e70bd46 MW	360	if (MP_LEN(a) >= n) vl = a->v + n;
0e70bd46 MW	361	else vl = a->vl;
f5f35081	362	mp_build(&al, a->v, vl);
f5f35081	363	u = mp_mul(MP_NEW, &al, mm->mi);
0e70bd46	364	if (MP_LEN(u) > n) u->vl = u->v + n;
f5f35081	365	u = mp_mul(u, u, mm->m);
f5f35081	366	d = mp_add(d, a, u);
366e9701	367	MP_ENSURE(d, n);
f5f35081	368	mp_drop(u);
f5f35081	369	}
d3409d5e	370
f5f35081	371	/* --- Otherwise do it the hard way --- */
d3409d5e	372
d3409d5e	373	else {
d34decd2	374	a = MP_COPY(a);
0e70bd46	375	if (d) MP_DROP(d);
d34decd2	376	d = a;
0e70bd46 MW	377	MP_DEST(d, 2*mm->n + 1, a->f);
0e70bd46 MW	378	redccore(d->v, d->vl, mm->m->v, mm->n, mm->mi->v);
d3409d5e	379	}
d3409d5e	380
f5f35081	381	/* --- Wrap everything up --- */
d3409d5e	382
0e70bd46	383	finish(mm, d);
d3409d5e	384	return (d);
	385	}
	386
52e4b041	387	#endif
52e4b041	388
d3409d5e	389	/* --- @mpmont_mul@ --- *
d3409d5e	390	*
295f4f90	391	* Arguments: @const mpmont *mm@ = pointer to Montgomery reduction context
d3409d5e	392	* @mp *d@ = destination
ef5f4810	393	* @mp a, b@ = sources, assumed positive
d3409d5e	394	*
	395	* Returns: Result, %$a b R^{-1} \bmod m$%.
	396	*/
	397
52e4b041	398	#ifdef MPMONT_DISABLE
52e4b041	399
295f4f90	400	mp mpmont_mul(const mpmont mm, mp d, mp a, mp *b)
52e4b041	401	{
	402	d = mp_mul(d, a, b);
	403	mp_div(0, &d, d, mm->m);
	404	return (d);
	405	}
	406
	407	#else
	408
295f4f90	409	mp mpmont_mul(const mpmont mm, mp d, mp a, mp *b)
d3409d5e	410	{
444083ae MW	411	size_t n = mm->n;
444083ae MW	412
d6b9dc04	413	if (n > MPMONT_KTHRESH) {
ef5f4810	414	d = mp_mul(d, a, b);
	415	d = mpmont_reduce(mm, d, d);
	416	} else {
444083ae MW	417	a = MP_COPY(a); b = MP_COPY(b);
444083ae MW	418	MP_DEST(d, 2*n + 1, a->f \| b->f \| MP_UNDEF);
0e70bd46 MW	419	mulcore(d->v, d->vl, a->v, a->vl, b->v, b->vl,
	420	mm->m->v, mm->n, mm->mi->v);
	421	d->f = ((a->f \| b->f) & MP_BURN) \| ((a->f ^ b->f) & MP_NEG);
	422	finish(mm, d);
	423	MP_DROP(a); MP_DROP(b);
d3409d5e	424	}
d3409d5e	425
d3409d5e	426	return (d);
	427	}
	428
52e4b041	429	#endif
52e4b041	430
d3409d5e	431	/----- Test rig ----------------------------------------------------------/
	432
	433	#ifdef TEST_RIG
	434
	435	static int tcreate(dstr *v)
	436	{
	437	mp m = (mp **)v[0].buf;
	438	mp mi = (mp **)v[1].buf;
	439	mp r = (mp **)v[2].buf;
	440	mp r2 = (mp **)v[3].buf;
	441
	442	mpmont mm;
	443	int ok = 1;
	444
	445	mpmont_create(&mm, m);
	446
f5f35081	447	if (mm.mi->v[0] != mi->v[0]) {
d3409d5e	448	fprintf(stderr, "\n*** bad mi: found %lu, expected %lu",
f5f35081	449	(unsigned long)mm.mi->v[0], (unsigned long)mi->v[0]);
d3409d5e	450	fputs("\nm = ", stderr); mp_writefile(m, stderr, 10);
	451	fputc('\n', stderr);
	452	ok = 0;
	453	}
	454
032099d1	455	if (!MP_EQ(mm.r, r)) {
d3409d5e	456	fputs("\n*** bad r", stderr);
	457	fputs("\nm = ", stderr); mp_writefile(m, stderr, 10);
	458	fputs("\nexpected ", stderr); mp_writefile(r, stderr, 10);
45c0fd36	459	fputs("\n found ", stderr); mp_writefile(mm.r, stderr, 10);
d3409d5e	460	fputc('\n', stderr);
	461	ok = 0;
	462	}
	463
032099d1	464	if (!MP_EQ(mm.r2, r2)) {
d3409d5e	465	fputs("\n*** bad r2", stderr);
	466	fputs("\nm = ", stderr); mp_writefile(m, stderr, 10);
	467	fputs("\nexpected ", stderr); mp_writefile(r2, stderr, 10);
45c0fd36	468	fputs("\n found ", stderr); mp_writefile(mm.r2, stderr, 10);
d3409d5e	469	fputc('\n', stderr);
	470	ok = 0;
	471	}
	472
	473	MP_DROP(m);
	474	MP_DROP(mi);
	475	MP_DROP(r);
	476	MP_DROP(r2);
	477	mpmont_destroy(&mm);
ef5f4810	478	assert(mparena_count(MPARENA_GLOBAL) == 0);
d3409d5e	479	return (ok);
	480	}
	481
	482	static int tmul(dstr *v)
	483	{
	484	mp m = (mp **)v[0].buf;
	485	mp a = (mp **)v[1].buf;
	486	mp b = (mp **)v[2].buf;
	487	mp r = (mp **)v[3].buf;
d3409d5e	488	int ok = 1;
	489
	490	mpmont mm;
	491	mpmont_create(&mm, m);
	492
	493	{
79a34029	494	mp *qr = mp_mul(MP_NEW, a, b);
	495	mp_div(0, &qr, qr, m);
	496
032099d1	497	if (!MP_EQ(qr, r)) {
79a34029	498	fputs("\n*** classical modmul failed", stderr);
	499	fputs("\n m = ", stderr); mp_writefile(m, stderr, 10);
	500	fputs("\n a = ", stderr); mp_writefile(a, stderr, 10);
	501	fputs("\n b = ", stderr); mp_writefile(b, stderr, 10);
	502	fputs("\n r = ", stderr); mp_writefile(r, stderr, 10);
	503	fputs("\nqr = ", stderr); mp_writefile(qr, stderr, 10);
	504	fputc('\n', stderr);
	505	ok = 0;
	506	}
	507
	508	mp_drop(qr);
	509	}
	510
	511	{
d3409d5e	512	mp *ar = mpmont_mul(&mm, MP_NEW, a, mm.r2);
d3409d5e	513	mp *br = mpmont_mul(&mm, MP_NEW, b, mm.r2);
79a34029	514	mp *mr = mpmont_mul(&mm, MP_NEW, ar, br);
d3409d5e	515	mr = mpmont_reduce(&mm, mr, mr);
032099d1	516	if (!MP_EQ(mr, r)) {
79a34029	517	fputs("\n*** montgomery modmul failed", stderr);
	518	fputs("\n m = ", stderr); mp_writefile(m, stderr, 10);
	519	fputs("\n a = ", stderr); mp_writefile(a, stderr, 10);
	520	fputs("\n b = ", stderr); mp_writefile(b, stderr, 10);
	521	fputs("\n r = ", stderr); mp_writefile(r, stderr, 10);
	522	fputs("\nmr = ", stderr); mp_writefile(mr, stderr, 10);
	523	fputc('\n', stderr);
	524	ok = 0;
	525	}
d3409d5e	526	MP_DROP(ar); MP_DROP(br);
79a34029	527	mp_drop(mr);
d3409d5e	528	}
d3409d5e	529
d3409d5e	530
	531	MP_DROP(m);
	532	MP_DROP(a);
	533	MP_DROP(b);
	534	MP_DROP(r);
d3409d5e	535	mpmont_destroy(&mm);
ef5f4810	536	assert(mparena_count(MPARENA_GLOBAL) == 0);
d3409d5e	537	return ok;
	538	}
	539
d3409d5e	540	static test_chunk tests[] = {
ef5f4810	541	{ "create", tcreate, { &type_mp, &type_mp, &type_mp, &type_mp, 0 } },
ef5f4810	542	{ "mul", tmul, { &type_mp, &type_mp, &type_mp, &type_mp, 0 } },
d3409d5e	543	{ 0, 0, { 0 } },
	544	};
	545
	546	int main(int argc, char *argv[])
	547	{
	548	sub_init();
0f00dc4c	549	test_run(argc, argv, tests, SRCDIR "/t/mpmont");
d3409d5e	550	return (0);
	551	}
	552
	553	#endif
	554
	555	/----- That's all, folks -------------------------------------------------/