[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)
  MAYBE_REDC4(x86_avx)
#endif

#if CPUFAM_AMD64
  MAYBE_REDC4(amd64_sse2)
  MAYBE_REDC4(amd64_avx)
#endif

static redccore__functype *pick_redccore(void)
{
#if CPUFAM_X86
  DISPATCH_PICK_COND(mpmont_reduce, maybe_redc4_x86_avx,
		     cpu_feature_p(CPUFEAT_X86_AVX));
  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_avx,
		     cpu_feature_p(CPUFEAT_X86_AVX));
  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)
  MAYBE_MUL4(x86_avx)
#endif

#if CPUFAM_AMD64
  MAYBE_MUL4(amd64_sse2)
  MAYBE_MUL4(amd64_avx)
#endif

static mulcore__functype *pick_mulcore(void)
{
#if CPUFAM_X86
  DISPATCH_PICK_COND(mpmont_mul, maybe_mul4_x86_avx,
		     cpu_feature_p(CPUFEAT_X86_AVX));
  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_avx,
		     cpu_feature_p(CPUFEAT_X86_AVX));
  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

#ifdef ENABLE_ASM_DEBUG
#  include "regdump.h"
#endif

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