[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

#if CPUFAM_ARMEL
  MAYBE_REDC4(arm_neon)
#endif

#if CPUFAM_ARM64
  MAYBE_REDC4(arm64_simd)
#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
#if CPUFAM_ARMEL
  DISPATCH_PICK_COND(mpmont_reduce, maybe_redc4_arm_neon,
		     cpu_feature_p(CPUFEAT_ARM_NEON));
#endif
#if CPUFAM_ARM64
  DISPATCH_PICK_COND(mpmont_reduce, maybe_redc4_arm64_simd,
		     cpu_feature_p(CPUFEAT_ARM_NEON));
#endif
  DISPATCH_PICK_FALLBACK(mpmont_reduce, simple_redccore);
}

/* --- @mulcore@ --- *
 *
 * 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

#if CPUFAM_ARMEL
  MAYBE_MUL4(arm_neon)
#endif

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