[u/mdw/catacomb] / mpmont.c

/* -*-c-*-
 *
 * $Id: mpmont.c,v 1.8 1999/12/22 15:55:00 mdw Exp $
 *
 * 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.
 */

/*----- Revision history --------------------------------------------------* 
 *
 * $Log: mpmont.c,v $
 * Revision 1.8  1999/12/22 15:55:00  mdw
 * Adjust Karatsuba parameters.
 *
 * Revision 1.7  1999/12/11 01:51:14  mdw
 * Use a Karatsuba-based reduction for large moduli.
 *
 * Revision 1.6  1999/12/10 23:18:39  mdw
 * Change interface for suggested destinations.
 *
 * Revision 1.5  1999/11/22 13:58:40  mdw
 * Add an option to disable Montgomery reduction, so that performance
 * comparisons can be done.
 *
 * Revision 1.4  1999/11/21 12:27:06  mdw
 * Remove a division from the Montgomery setup by calculating
 * %$R^2 \bmod m$% first and then %$R \bmod m$% by Montgomery reduction of
 * %$R^2$%.
 *
 * Revision 1.3  1999/11/21 11:35:10  mdw
 * Performance improvement: use @mp_sqr@ and @mpmont_reduce@ instead of
 * @mpmont_mul@ for squaring in exponentiation.
 *
 * Revision 1.2  1999/11/19 13:17:26  mdw
 * Add extra interface to exponentiation which returns a Montgomerized
 * result.
 *
 * Revision 1.1  1999/11/17 18:02:16  mdw
 * New multiprecision integer arithmetic suite.
 *
 */

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

#include "mp.h"
#include "mpmont.h"

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

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

/* #define MPMONT_DISABLE */

/*----- Main code ---------------------------------------------------------*/

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

#ifdef MPMONT_DISABLE

void 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;
}

#else

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

  /* --- Validate the arguments --- */

  assert(((void)"Montgomery modulus must be positive",
	  (m->f & MP_NEG) == 0));
  assert(((void)"Montgomery modulus must be odd", m->v[0] & 1));

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

  mp_shrink(m);
  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_NEW;
  mp_gcd(0, 0, &mm->mi, &r, m);
  mm->mi = mp_sub(mm->mi, &r, mm->mi);

  /* --- 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);
}

#endif

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

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

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

#ifdef MPMONT_DISABLE

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

#else

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

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

  if (n > KARATSUBA_CUTOFF * 3) {
    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_drop(u);
  }

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

  else {
    mpw *dv, *dvl;
    mpw *mv, *mvl;
    mpw mi;
    size_t k = n;

    /* --- Initial conditioning of the arguments --- */

    if (d == a)
      MP_MODIFY(d, 2 * n + 1);
    else {
      MP_MODIFY(d, 2 * n + 1);
      MPX_COPY(d->v, d->vl, a->v, a->vl);
    }
    
    dv = d->v; dvl = d->vl;
    mv = mm->m->v; mvl = mm->m->vl;

    /* --- Let's go to work --- */

    mi = mm->mi->v[0];
    while (k--) {
      mpw u = MPW(*dv * mi);
      MPX_UMLAN(dv, dvl, mv, mvl, u);
      dv++;
    }
  }

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

  d->f = a->f & MP_BURN;
  memmove(d->v, d->v + n, MPWS(MP_LEN(d) - n));
  d->vl -= n;
  if (MP_CMP(d, >=, mm->m))
    d = mp_sub(d, d, mm->m);
  MP_SHRINK(d);
  return (d);
}

#endif

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

#ifdef MPMONT_DISABLE

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

#else

mp *mpmont_mul(mpmont *mm, mp *d, mp *a, mp *b)
{
  if (mm->n > KARATSUBA_CUTOFF * 3) {
    d = mp_mul(d, a, b);
    d = mpmont_reduce(mm, d, d);
  } else {
    mpw *dv, *dvl;
    mpw *av, *avl;
    mpw *bv, *bvl;
    mpw *mv, *mvl;
    mpw y;
    size_t n, i;
    mpw mi;

    /* --- Initial conditioning of the arguments --- */

    if (MP_LEN(a) > MP_LEN(b)) {
      mp *t = a; a = b; b = t;
    }
    n = MP_LEN(mm->m);

    a = MP_COPY(a);
    b = MP_COPY(b);
    MP_MODIFY(d, 2 * n + 1);
    dv = d->v; dvl = d->vl;
    MPX_ZERO(dv, dvl);
    av = a->v; avl = a->vl;
    bv = b->v; bvl = b->vl;
    mv = mm->m->v; mvl = mm->m->vl;
    y = *bv;

    /* --- Montgomery multiplication phase --- */

    i = 0;
    mi = mm->mi->v[0];
    while (i < n && av < avl) {
      mpw x = *av++;
      mpw u = MPW((*dv + x * y) * mi);
      MPX_UMLAN(dv, dvl, bv, bvl, x);
      MPX_UMLAN(dv, dvl, mv, mvl, u);
      dv++;
      i++;
    }

    /* --- Simpler Montgomery reduction phase --- */

    while (i < n) {
      mpw u = MPW(*dv * mi);
      MPX_UMLAN(dv, dvl, mv, mvl, u);
      dv++;
      i++;
    }

    /* --- Done --- */

    memmove(d->v, dv, MPWS(dvl - dv));
    d->vl -= dv - d->v;
    MP_SHRINK(d);
    d->f = (a->f | b->f) & MP_BURN;
    if (MP_CMP(d, >=, mm->m))
      d = mp_sub(d, d, mm->m);
    MP_DROP(a);
    MP_DROP(b);
  }

  return (d);
}

#endif

/* --- @mpmont_expr@ --- *
 *
 * Arguments:	@mpmont *mm@ = pointer to Montgomery reduction context
 *		@mp *d@ = fake destination
 *		@mp *a@ = base
 *		@mp *e@ = exponent
 *
 * Returns:	Result, %$a^e R \bmod m$%.
 */

mp *mpmont_expr(mpmont *mm, mp *d, mp *a, mp *e)
{
  mpscan sc;
  mp *ar = mpmont_mul(mm, MP_NEW, a, mm->r2);
  mp *x = MP_COPY(mm->r);
  mp *spare = MP_NEW;

  mp_scan(&sc, e);

  if (MP_STEP(&sc)) {
    size_t sq = 0;
    for (;;) {
      mp *dd;
      if (MP_BIT(&sc)) {
	while (sq) {
	  dd = mp_sqr(spare, ar);
	  dd = mpmont_reduce(mm, dd, dd);
	  spare = ar; ar = dd;
	  sq--;
	}
	dd = mpmont_mul(mm, spare, x, ar);
	spare = x; x = dd;
      }
      sq++;
      if (!MP_STEP(&sc))
	break;
    }
  }
  MP_DROP(ar);
  if (spare != MP_NEW)
    MP_DROP(spare);
  if (d != MP_NEW)
    MP_DROP(d);
  return (x);
}

/* --- @mpmont_exp@ --- *
 *
 * Arguments:	@mpmont *mm@ = pointer to Montgomery reduction context
 *		@mp *d@ = fake destination
 *		@mp *a@ = base
 *		@mp *e@ = exponent
 *
 * Returns:	Result, %$a^e \bmod m$%.
 */

mp *mpmont_exp(mpmont *mm, mp *d, mp *a, mp *e)
{
  d = mpmont_expr(mm, d, a, e);
  d = mpmont_reduce(mm, d, d);
  return (d);
}

/*----- 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_CMP(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_CMP(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_CMP(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_CMP(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 int texp(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;
  mp *mr;
  int ok = 1;

  mpmont mm;
  mpmont_create(&mm, m);

  mr = mpmont_exp(&mm, MP_NEW, a, b);

  if (MP_CMP(mr, !=, r)) {
    fputs("\n*** montgomery modexp failed", stderr);
    fputs("\n m = ", stderr); mp_writefile(m, stderr, 10);
    fputs("\n a = ", stderr); mp_writefile(a, stderr, 10);
    fputs("\n e = ", 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(m);
  MP_DROP(a);
  MP_DROP(b);
  MP_DROP(r);
  MP_DROP(mr);
  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 } },
  { "exp", texp, { &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 "/tests/mpmont");
  return (0);
}

#endif

/*----- That's all, folks -------------------------------------------------*/
Commit	Line	Data
d3409d5e	1	/* --c--
d3409d5e	2	*
01f6ed1a	3	* $Id: mpmont.c,v 1.8 1999/12/22 15:55:00 mdw Exp $
d3409d5e	4	*
	5	* Montgomery reduction
	6	*
	7	* (c) 1999 Straylight/Edgeware
	8	*/
	9
	10	/----- Licensing notice --------------------------------------------------
	11	*
	12	* This file is part of Catacomb.
	13	*
	14	* Catacomb is free software; you can redistribute it and/or modify
	15	* it under the terms of the GNU Library General Public License as
	16	* published by the Free Software Foundation; either version 2 of the
	17	* License, or (at your option) any later version.
	18	*
	19	* Catacomb is distributed in the hope that it will be useful,
	20	* but WITHOUT ANY WARRANTY; without even the implied warranty of
	21	* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
	22	* GNU Library General Public License for more details.
	23	*
	24	* You should have received a copy of the GNU Library General Public
	25	* License along with Catacomb; if not, write to the Free
	26	* Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
	27	* MA 02111-1307, USA.
	28	*/
	29
	30	/----- Revision history --------------------------------------------------
	31	*
	32	* $Log: mpmont.c,v $
01f6ed1a	33	* Revision 1.8 1999/12/22 15:55:00 mdw
	34	* Adjust Karatsuba parameters.
	35	*
f5f35081	36	* Revision 1.7 1999/12/11 01:51:14 mdw
	37	* Use a Karatsuba-based reduction for large moduli.
	38	*
ef5f4810	39	* Revision 1.6 1999/12/10 23:18:39 mdw
	40	* Change interface for suggested destinations.
	41	*
52e4b041	42	* Revision 1.5 1999/11/22 13:58:40 mdw
	43	* Add an option to disable Montgomery reduction, so that performance
	44	* comparisons can be done.
	45	*
93feaa6e	46	* Revision 1.4 1999/11/21 12:27:06 mdw
	47	* Remove a division from the Montgomery setup by calculating
	48	* %$R^2 \bmod m$% first and then %$R \bmod m$% by Montgomery reduction of
	49	* %$R^2$%.
	50	*
79a34029	51	* Revision 1.3 1999/11/21 11:35:10 mdw
	52	* Performance improvement: use @mp_sqr@ and @mpmont_reduce@ instead of
	53	* @mpmont_mul@ for squaring in exponentiation.
	54	*
17ad212e	55	* Revision 1.2 1999/11/19 13:17:26 mdw
	56	* Add extra interface to exponentiation which returns a Montgomerized
	57	* result.
	58	*
d3409d5e	59	* Revision 1.1 1999/11/17 18:02:16 mdw
	60	* New multiprecision integer arithmetic suite.
	61	*
	62	*/
	63
	64	/----- Header files ------------------------------------------------------/
	65
	66	#include "mp.h"
	67	#include "mpmont.h"
	68
52e4b041	69	/----- Tweakables --------------------------------------------------------/
	70
	71	/* --- @MPMONT_DISABLE@ --- *
	72	*
	73	* Replace all the clever Montgomery reduction with good old-fashioned long
	74	* division.
	75	*/
	76
	77	/* #define MPMONT_DISABLE */
	78
d3409d5e	79	/----- Main code ---------------------------------------------------------/
	80
	81	/* --- @mpmont_create@ --- *
	82	*
	83	* Arguments: @mpmont *mm@ = pointer to Montgomery reduction context
	84	* @mp *m@ = modulus to use
	85	*
	86	* Returns: ---
	87	*
	88	* Use: Initializes a Montgomery reduction context ready for use.
ef5f4810	89	* The argument @m@ must be a positive odd integer.
d3409d5e	90	*/
d3409d5e	91
52e4b041	92	#ifdef MPMONT_DISABLE
	93
	94	void mpmont_create(mpmont mm, mp m)
	95	{
	96	mp_shrink(m);
	97	mm->m = MP_COPY(m);
	98	mm->r = MP_ONE;
	99	mm->r2 = MP_ONE;
f5f35081	100	mm->mi = MP_ONE;
52e4b041	101	}
	102
	103	#else
	104
d3409d5e	105	void mpmont_create(mpmont mm, mp m)
d3409d5e	106	{
f5f35081	107	size_t n = MP_LEN(m);
	108	mp r2 = mp_create(2 n + 1);
	109	mp r;
	110
ef5f4810	111	/* --- Validate the arguments --- */
	112
	113	assert(((void)"Montgomery modulus must be positive",
	114	(m->f & MP_NEG) == 0));
	115	assert(((void)"Montgomery modulus must be odd", m->v[0] & 1));
	116
d3409d5e	117	/* --- Take a copy of the modulus --- */
	118
	119	mp_shrink(m);
	120	mm->m = MP_COPY(m);
	121
f5f35081	122	/* --- Determine %$R^2$% --- */
d3409d5e	123
f5f35081	124	mm->n = n;
	125	MPX_ZERO(r2->v, r2->vl - 1);
	126	r2->vl[-1] = 1;
d3409d5e	127
f5f35081	128	/* --- Find the magic value @mi@ --- */
	129
	130	mp_build(&r, r2->v + n, r2->vl);
	131	mm->mi = MP_NEW;
	132	mp_gcd(0, 0, &mm->mi, &r, m);
	133	mm->mi = mp_sub(mm->mi, &r, mm->mi);
d3409d5e	134
	135	/* --- Discover the values %$R \bmod m$% and %$R^2 \bmod m$% --- */
	136
f5f35081	137	mm->r2 = MP_NEW;
	138	mp_div(0, &mm->r2, r2, m);
	139	mm->r = mpmont_reduce(mm, MP_NEW, mm->r2);
	140	MP_DROP(r2);
d3409d5e	141	}
d3409d5e	142
52e4b041	143	#endif
52e4b041	144
d3409d5e	145	/* --- @mpmont_destroy@ --- *
	146	*
	147	* Arguments: @mpmont *mm@ = pointer to a Montgomery reduction context
	148	*
	149	* Returns: ---
	150	*
	151	* Use: Disposes of a context when it's no longer of any use to
	152	* anyone.
	153	*/
	154
	155	void mpmont_destroy(mpmont *mm)
	156	{
	157	MP_DROP(mm->m);
	158	MP_DROP(mm->r);
	159	MP_DROP(mm->r2);
f5f35081	160	MP_DROP(mm->mi);
d3409d5e	161	}
	162
	163	/* --- @mpmont_reduce@ --- *
	164	*
	165	* Arguments: @mpmont *mm@ = pointer to Montgomery reduction context
	166	* @mp *d@ = destination
ef5f4810	167	* @mp *a@ = source, assumed positive
d3409d5e	168	*
	169	* Returns: Result, %$a R^{-1} \bmod m$%.
	170	*/
	171
52e4b041	172	#ifdef MPMONT_DISABLE
52e4b041	173
ef5f4810	174	mp mpmont_reduce(mpmont mm, mp d, mp a)
52e4b041	175	{
	176	mp_div(0, &d, a, mm->m);
	177	return (d);
	178	}
	179
	180	#else
	181
ef5f4810	182	mp mpmont_reduce(mpmont mm, mp d, mp a)
d3409d5e	183	{
f5f35081	184	size_t n = mm->n;
	185
	186	/* --- Check for serious Karatsuba reduction --- */
	187
01f6ed1a	188	if (n > KARATSUBA_CUTOFF * 3) {
f5f35081	189	mp al;
	190	mpw *vl;
	191	mp *u;
	192
	193	if (MP_LEN(a) >= n)
	194	vl = a->v + n;
	195	else
	196	vl = a->vl;
	197	mp_build(&al, a->v, vl);
	198	u = mp_mul(MP_NEW, &al, mm->mi);
	199	if (MP_LEN(u) > n)
	200	u->vl = u->v + n;
	201	u = mp_mul(u, u, mm->m);
	202	d = mp_add(d, a, u);
	203	mp_drop(u);
	204	}
d3409d5e	205
f5f35081	206	/* --- Otherwise do it the hard way --- */
d3409d5e	207
d3409d5e	208	else {
f5f35081	209	mpw dv, dvl;
	210	mpw mv, mvl;
	211	mpw mi;
	212	size_t k = n;
	213
	214	/* --- Initial conditioning of the arguments --- */
	215
	216	if (d == a)
	217	MP_MODIFY(d, 2 * n + 1);
	218	else {
	219	MP_MODIFY(d, 2 * n + 1);
	220	MPX_COPY(d->v, d->vl, a->v, a->vl);
	221	}
d3409d5e	222
f5f35081	223	dv = d->v; dvl = d->vl;
f5f35081	224	mv = mm->m->v; mvl = mm->m->vl;
d3409d5e	225
f5f35081	226	/* --- Let's go to work --- */
d3409d5e	227
f5f35081	228	mi = mm->mi->v[0];
	229	while (k--) {
	230	mpw u = MPW(dv mi);
	231	MPX_UMLAN(dv, dvl, mv, mvl, u);
	232	dv++;
	233	}
d3409d5e	234	}
d3409d5e	235
f5f35081	236	/* --- Wrap everything up --- */
d3409d5e	237
d3409d5e	238	d->f = a->f & MP_BURN;
f5f35081	239	memmove(d->v, d->v + n, MPWS(MP_LEN(d) - n));
f5f35081	240	d->vl -= n;
17ad212e	241	if (MP_CMP(d, >=, mm->m))
17ad212e	242	d = mp_sub(d, d, mm->m);
f5f35081	243	MP_SHRINK(d);
d3409d5e	244	return (d);
	245	}
	246
52e4b041	247	#endif
52e4b041	248
d3409d5e	249	/* --- @mpmont_mul@ --- *
	250	*
	251	* Arguments: @mpmont *mm@ = pointer to Montgomery reduction context
	252	* @mp *d@ = destination
ef5f4810	253	* @mp a, b@ = sources, assumed positive
d3409d5e	254	*
	255	* Returns: Result, %$a b R^{-1} \bmod m$%.
	256	*/
	257
52e4b041	258	#ifdef MPMONT_DISABLE
52e4b041	259
ef5f4810	260	mp mpmont_mul(mpmont mm, mp d, mp a, mp *b)
52e4b041	261	{
	262	d = mp_mul(d, a, b);
	263	mp_div(0, &d, d, mm->m);
	264	return (d);
	265	}
	266
	267	#else
	268
ef5f4810	269	mp mpmont_mul(mpmont mm, mp d, mp a, mp *b)
d3409d5e	270	{
01f6ed1a	271	if (mm->n > KARATSUBA_CUTOFF * 3) {
ef5f4810	272	d = mp_mul(d, a, b);
	273	d = mpmont_reduce(mm, d, d);
	274	} else {
	275	mpw dv, dvl;
	276	mpw av, avl;
	277	mpw bv, bvl;
	278	mpw mv, mvl;
	279	mpw y;
	280	size_t n, i;
f5f35081	281	mpw mi;
ef5f4810	282
	283	/* --- Initial conditioning of the arguments --- */
	284
	285	if (MP_LEN(a) > MP_LEN(b)) {
	286	mp *t = a; a = b; b = t;
	287	}
	288	n = MP_LEN(mm->m);
d3409d5e	289
ef5f4810	290	a = MP_COPY(a);
	291	b = MP_COPY(b);
	292	MP_MODIFY(d, 2 * n + 1);
	293	dv = d->v; dvl = d->vl;
	294	MPX_ZERO(dv, dvl);
	295	av = a->v; avl = a->vl;
	296	bv = b->v; bvl = b->vl;
	297	mv = mm->m->v; mvl = mm->m->vl;
	298	y = *bv;
	299
	300	/* --- Montgomery multiplication phase --- */
	301
	302	i = 0;
f5f35081	303	mi = mm->mi->v[0];
ef5f4810	304	while (i < n && av < avl) {
ef5f4810	305	mpw x = *av++;
f5f35081	306	mpw u = MPW((dv + x y) * mi);
ef5f4810	307	MPX_UMLAN(dv, dvl, bv, bvl, x);
	308	MPX_UMLAN(dv, dvl, mv, mvl, u);
	309	dv++;
	310	i++;
	311	}
d3409d5e	312
ef5f4810	313	/* --- Simpler Montgomery reduction phase --- */
d3409d5e	314
ef5f4810	315	while (i < n) {
f5f35081	316	mpw u = MPW(dv mi);
ef5f4810	317	MPX_UMLAN(dv, dvl, mv, mvl, u);
	318	dv++;
	319	i++;
	320	}
d3409d5e	321
ef5f4810	322	/* --- Done --- */
d3409d5e	323
ef5f4810	324	memmove(d->v, dv, MPWS(dvl - dv));
	325	d->vl -= dv - d->v;
	326	MP_SHRINK(d);
	327	d->f = (a->f \| b->f) & MP_BURN;
	328	if (MP_CMP(d, >=, mm->m))
	329	d = mp_sub(d, d, mm->m);
	330	MP_DROP(a);
	331	MP_DROP(b);
d3409d5e	332	}
d3409d5e	333
d3409d5e	334	return (d);
	335	}
	336
52e4b041	337	#endif
52e4b041	338
17ad212e	339	/* --- @mpmont_expr@ --- *
d3409d5e	340	*
d3409d5e	341	* Arguments: @mpmont *mm@ = pointer to Montgomery reduction context
ef5f4810	342	* @mp *d@ = fake destination
	343	* @mp *a@ = base
	344	* @mp *e@ = exponent
d3409d5e	345	*
17ad212e	346	* Returns: Result, %$a^e R \bmod m$%.
d3409d5e	347	*/
d3409d5e	348
ef5f4810	349	mp mpmont_expr(mpmont mm, mp d, mp a, mp *e)
d3409d5e	350	{
	351	mpscan sc;
	352	mp *ar = mpmont_mul(mm, MP_NEW, a, mm->r2);
ef5f4810	353	mp *x = MP_COPY(mm->r);
17ad212e	354	mp *spare = MP_NEW;
d3409d5e	355
	356	mp_scan(&sc, e);
	357
	358	if (MP_STEP(&sc)) {
17ad212e	359	size_t sq = 0;
d3409d5e	360	for (;;) {
	361	mp *dd;
	362	if (MP_BIT(&sc)) {
17ad212e	363	while (sq) {
79a34029	364	dd = mp_sqr(spare, ar);
79a34029	365	dd = mpmont_reduce(mm, dd, dd);
17ad212e	366	spare = ar; ar = dd;
	367	sq--;
	368	}
ef5f4810	369	dd = mpmont_mul(mm, spare, x, ar);
ef5f4810	370	spare = x; x = dd;
d3409d5e	371	}
17ad212e	372	sq++;
d3409d5e	373	if (!MP_STEP(&sc))
d3409d5e	374	break;
d3409d5e	375	}
	376	}
	377	MP_DROP(ar);
17ad212e	378	if (spare != MP_NEW)
17ad212e	379	MP_DROP(spare);
ef5f4810	380	if (d != MP_NEW)
	381	MP_DROP(d);
	382	return (x);
17ad212e	383	}
	384
	385	/* --- @mpmont_exp@ --- *
	386	*
	387	* Arguments: @mpmont *mm@ = pointer to Montgomery reduction context
ef5f4810	388	* @mp *d@ = fake destination
	389	* @mp *a@ = base
	390	* @mp *e@ = exponent
17ad212e	391	*
	392	* Returns: Result, %$a^e \bmod m$%.
	393	*/
	394
ef5f4810	395	mp mpmont_exp(mpmont mm, mp d, mp a, mp *e)
17ad212e	396	{
ef5f4810	397	d = mpmont_expr(mm, d, a, e);
17ad212e	398	d = mpmont_reduce(mm, d, d);
17ad212e	399	return (d);
d3409d5e	400	}
	401
	402	/----- Test rig ----------------------------------------------------------/
	403
	404	#ifdef TEST_RIG
	405
	406	static int tcreate(dstr *v)
	407	{
	408	mp m = (mp **)v[0].buf;
	409	mp mi = (mp **)v[1].buf;
	410	mp r = (mp **)v[2].buf;
	411	mp r2 = (mp **)v[3].buf;
	412
	413	mpmont mm;
	414	int ok = 1;
	415
	416	mpmont_create(&mm, m);
	417
f5f35081	418	if (mm.mi->v[0] != mi->v[0]) {
d3409d5e	419	fprintf(stderr, "\n*** bad mi: found %lu, expected %lu",
f5f35081	420	(unsigned long)mm.mi->v[0], (unsigned long)mi->v[0]);
d3409d5e	421	fputs("\nm = ", stderr); mp_writefile(m, stderr, 10);
	422	fputc('\n', stderr);
	423	ok = 0;
	424	}
	425
	426	if (MP_CMP(mm.r, !=, r)) {
	427	fputs("\n*** bad r", stderr);
	428	fputs("\nm = ", stderr); mp_writefile(m, stderr, 10);
	429	fputs("\nexpected ", stderr); mp_writefile(r, stderr, 10);
17ad212e	430	fputs("\n found ", stderr); mp_writefile(mm.r, stderr, 10);
d3409d5e	431	fputc('\n', stderr);
	432	ok = 0;
	433	}
	434
	435	if (MP_CMP(mm.r2, !=, r2)) {
	436	fputs("\n*** bad r2", stderr);
	437	fputs("\nm = ", stderr); mp_writefile(m, stderr, 10);
	438	fputs("\nexpected ", stderr); mp_writefile(r2, stderr, 10);
17ad212e	439	fputs("\n found ", stderr); mp_writefile(mm.r2, stderr, 10);
d3409d5e	440	fputc('\n', stderr);
	441	ok = 0;
	442	}
	443
	444	MP_DROP(m);
	445	MP_DROP(mi);
	446	MP_DROP(r);
	447	MP_DROP(r2);
	448	mpmont_destroy(&mm);
ef5f4810	449	assert(mparena_count(MPARENA_GLOBAL) == 0);
d3409d5e	450	return (ok);
	451	}
	452
	453	static int tmul(dstr *v)
	454	{
	455	mp m = (mp **)v[0].buf;
	456	mp a = (mp **)v[1].buf;
	457	mp b = (mp **)v[2].buf;
	458	mp r = (mp **)v[3].buf;
d3409d5e	459	int ok = 1;
	460
	461	mpmont mm;
	462	mpmont_create(&mm, m);
	463
	464	{
79a34029	465	mp *qr = mp_mul(MP_NEW, a, b);
	466	mp_div(0, &qr, qr, m);
	467
	468	if (MP_CMP(qr, !=, r)) {
	469	fputs("\n*** classical modmul failed", stderr);
	470	fputs("\n m = ", stderr); mp_writefile(m, stderr, 10);
	471	fputs("\n a = ", stderr); mp_writefile(a, stderr, 10);
	472	fputs("\n b = ", stderr); mp_writefile(b, stderr, 10);
	473	fputs("\n r = ", stderr); mp_writefile(r, stderr, 10);
	474	fputs("\nqr = ", stderr); mp_writefile(qr, stderr, 10);
	475	fputc('\n', stderr);
	476	ok = 0;
	477	}
	478
	479	mp_drop(qr);
	480	}
	481
	482	{
d3409d5e	483	mp *ar = mpmont_mul(&mm, MP_NEW, a, mm.r2);
d3409d5e	484	mp *br = mpmont_mul(&mm, MP_NEW, b, mm.r2);
79a34029	485	mp *mr = mpmont_mul(&mm, MP_NEW, ar, br);
d3409d5e	486	mr = mpmont_reduce(&mm, mr, mr);
79a34029	487	if (MP_CMP(mr, !=, r)) {
	488	fputs("\n*** montgomery modmul failed", stderr);
	489	fputs("\n m = ", stderr); mp_writefile(m, stderr, 10);
	490	fputs("\n a = ", stderr); mp_writefile(a, stderr, 10);
	491	fputs("\n b = ", stderr); mp_writefile(b, stderr, 10);
	492	fputs("\n r = ", stderr); mp_writefile(r, stderr, 10);
	493	fputs("\nmr = ", stderr); mp_writefile(mr, stderr, 10);
	494	fputc('\n', stderr);
	495	ok = 0;
	496	}
d3409d5e	497	MP_DROP(ar); MP_DROP(br);
79a34029	498	mp_drop(mr);
d3409d5e	499	}
d3409d5e	500
d3409d5e	501
	502	MP_DROP(m);
	503	MP_DROP(a);
	504	MP_DROP(b);
	505	MP_DROP(r);
d3409d5e	506	mpmont_destroy(&mm);
ef5f4810	507	assert(mparena_count(MPARENA_GLOBAL) == 0);
d3409d5e	508	return ok;
	509	}
	510
	511	static int texp(dstr *v)
	512	{
	513	mp m = (mp **)v[0].buf;
	514	mp a = (mp **)v[1].buf;
	515	mp b = (mp **)v[2].buf;
	516	mp r = (mp **)v[3].buf;
	517	mp *mr;
	518	int ok = 1;
	519
	520	mpmont mm;
	521	mpmont_create(&mm, m);
	522
ef5f4810	523	mr = mpmont_exp(&mm, MP_NEW, a, b);
d3409d5e	524
	525	if (MP_CMP(mr, !=, r)) {
	526	fputs("\n*** montgomery modexp failed", stderr);
	527	fputs("\n m = ", stderr); mp_writefile(m, stderr, 10);
	528	fputs("\n a = ", stderr); mp_writefile(a, stderr, 10);
	529	fputs("\n e = ", stderr); mp_writefile(b, stderr, 10);
	530	fputs("\n r = ", stderr); mp_writefile(r, stderr, 10);
	531	fputs("\nmr = ", stderr); mp_writefile(mr, stderr, 10);
	532	fputc('\n', stderr);
	533	ok = 0;
	534	}
	535
	536	MP_DROP(m);
	537	MP_DROP(a);
	538	MP_DROP(b);
	539	MP_DROP(r);
	540	MP_DROP(mr);
	541	mpmont_destroy(&mm);
ef5f4810	542	assert(mparena_count(MPARENA_GLOBAL) == 0);
d3409d5e	543	return ok;
	544	}
	545
	546
	547	static test_chunk tests[] = {
ef5f4810	548	{ "create", tcreate, { &type_mp, &type_mp, &type_mp, &type_mp, 0 } },
	549	{ "mul", tmul, { &type_mp, &type_mp, &type_mp, &type_mp, 0 } },
	550	{ "exp", texp, { &type_mp, &type_mp, &type_mp, &type_mp, 0 } },
d3409d5e	551	{ 0, 0, { 0 } },
	552	};
	553
	554	int main(int argc, char *argv[])
	555	{
	556	sub_init();
	557	test_run(argc, argv, tests, SRCDIR "/tests/mpmont");
	558	return (0);
	559	}
	560
	561	#endif
	562
	563	/----- That's all, folks -------------------------------------------------/