[u/mdw/catacomb] / pub / bbs-jump.c

/* -*-c-*-
 *
 * Jumping around a BBS sequence
 *
 * (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 "bbs.h"
#include "mp.h"
#include "mpbarrett.h"
#include "mpcrt.h"
#include "mpint.h"

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

/* --- @jump@ --- *
 *
 * Arguments:	@bbs *b@ = pointer to BBS generator context
 *		@const bbs_priv *bp@ = pointer to BBS modulus factors
 *		@mp *n@ = number of steps to move
 *		@mp *px@ = exponent mod @p@ for a one-step jump
 *		@mp *qx@ = exponent mod @q@ for a one-step jump
 *
 * Returns:	---
 *
 * Use:		Jumps a BBS context a certain number of places (assuming the
 *		arguments are right).
 *
 *		Let the BBS modulus be %$n = pq$% and the current residue be
 *		%$x$%.  Then the computations performed are:
 *
 *		  * Calculate %$x_p = x \bmod p$% and %$x_q = x \bmod q$%.
 *
 *		  * Determine %$e_p = px^n \bmod (p - 1)$% and similarly
 *		    %$e_q = qx^n \bmod (p - 1)$%.
 *
 *		  * Calculate %$x_p' = x_p^{e_p} \bmod p$% and
 *		    %$x_q' = x_q^{e_q} \bmod q$%.
 *
 *		  * Combine %$x_p'$% and %$x_q'$% using the Chinese Remainder
 *		    Theorem.
 *
 *		If you want to step the generator forwards, simply set
 *		%$px = qx = 2$%.  If you want to step backwards, make
 *		%$px = (p + 1)/4$% and %$qx = (q + 1)/4$%.  Note that, if
 *		%$x$% is a quadratic residue mod $%p$%, then
 *
 *		%$(x^2) ^ {(p + 1)/4}$%
 *		  %${} = x^{(p + 1)/2}$%
 *		  %${} = x \cdot x^{(p - 1)/2}$%
 *		  %${} = x$%
 *
 *		Simple, no?  (Note that the division works because
 *		%$p \equiv 3 \pmod 4$%.)
 */

static void jump(bbs *b, const bbs_priv *bp, mp *n,
		 mp *px, mp *qx)
{
  mp *ep, *eq;
  mp *v[2] = { MP_NEW, MP_NEW };

  /* --- First work out the exponents --- */

  {
    mpbarrett mb;
    mp *m;

    m = mp_sub(MP_NEW, bp->p, MP_ONE);
    mpbarrett_create(&mb, m);
    ep = mpbarrett_exp(&mb, MP_NEW, px, n);
    mpbarrett_destroy(&mb);
    if (qx == px)
      eq = MP_COPY(ep);
    else {
      m = mp_sub(m, bp->q, MP_ONE);
      mpbarrett_create(&mb, m);
      eq = mpbarrett_exp(&mb, MP_NEW, qx, n);
      mpbarrett_destroy(&mb);
    }

    mp_drop(m);
  }

  /* --- Now calculate the residues of @x@ --- */

  mp_div(0, &v[0], b->x, bp->p);
  mp_div(0, &v[1], b->x, bp->q);

  /* --- Exponentiate --- */

  {
    mpbarrett mb;

    mpbarrett_create(&mb, bp->p);
    v[0] = mpbarrett_exp(&mb, v[0], v[0], ep);
    mpbarrett_destroy(&mb);

    mpbarrett_create(&mb, bp->q);
    v[1] = mpbarrett_exp(&mb, v[1], v[1], eq);
    mpbarrett_destroy(&mb);

    mp_drop(ep);
    mp_drop(eq);
  }

  /* --- Sort out the result using the Chinese Remainder Theorem --- */

  {
    mpcrt_mod mv[2];
    mpcrt c;
    int i;

    mv[0].m = MP_COPY(bp->p);
    mv[1].m = MP_COPY(bp->q);
    for (i = 0; i < 2; i++)
      mv[i].n = mv[i].ni = mv[i].nni = MP_NEW;
    mpcrt_create(&c, mv, 2, b->mb.m);
    b->x = mpcrt_solve(&c, b->x, v);
    mpcrt_destroy(&c);
  }

  /* --- Tidy away --- */

  mp_drop(v[0]);
  mp_drop(v[1]);
  b->r = b->x->v[0];
  b->b = b->k;
}

/* --- @bbs_{ff,rew}{,n}@ --- *
 *
 * Arguments:	@bbs *b@ = pointer to a BBS generator state
 *		@const bbs_priv *bp@ = pointer to BBS modulus factors
 *		@mp *n@, @unsigned long n@ = number of steps to make
 *
 * Returns:	---
 *
 * Use:		`Fast-forwards' or rewinds a Blum-Blum-Shub generator by @n@
 *		steps.  The @...n@ versions take an @unsigned long@ argument;
 *		the non-@...n@ versions a multiprecision integer.  If @n@ is
 *		negative then the generator is stepped in the reverse
 *		direction.
 */

static void ff(bbs *b, const bbs_priv *bp, mp *n)
  { jump(b, bp, n, MP_TWO, MP_TWO); }

static void rew(bbs *b, const bbs_priv *bp, mp *n)
{
  mp *px = mp_lsr(MP_NEW, bp->p, 2);
  mp *qx = mp_lsr(MP_NEW, bp->q, 2);
  px = mp_add(px, px, MP_ONE);
  qx = mp_add(qx, qx, MP_ONE);
  jump(b, bp, n, px, qx);
  mp_drop(px);
  mp_drop(qx);
}

void bbs_ff(bbs *b, const bbs_priv *bp, mp *n)
{
  if (!MP_NEGP(n))
    ff(b, bp, n);
  else {
    n = mp_neg(MP_NEW, n);
    rew(b, bp, n);
    mp_drop(n);
  }
}

void bbs_ffn(bbs *b, const bbs_priv *bp, unsigned long n)
  { mp *nn = mp_fromulong(MP_NEW, n); ff(b, bp, nn); mp_drop(nn); }

void bbs_rew(bbs *b, const bbs_priv *bp, mp *n)
{
  if (!MP_NEGP(n))
    rew(b, bp, n);
  else {
    n = mp_neg(MP_NEW, n);
    ff(b, bp, n);
    mp_drop(n);
  }
}

void bbs_rewn(bbs *b, const bbs_priv *bp, unsigned long n)
  { mp *nn = mp_fromulong(MP_NEW, n); bbs_rew(b, bp, nn); mp_drop(nn); }

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

#ifdef TEST_RIG

static int verify(dstr *v)
{
  bbs_priv bp;
  bbs b;
  mp *x;
  unsigned long n;
  int ok = 1;
  uint32 p, q, r;
  int i;

  bp.p = *(mp **)v[0].buf;
  bp.q = *(mp **)v[1].buf;
  bp.n = mp_mul(MP_NEW, bp.p, bp.q);
  x = *(mp **)v[2].buf;
  n = *(unsigned long *)v[3].buf;

  bbs_create(&b, bp.n, x);
  p = bbs_bits(&b, 32);

  bbs_seed(&b, x);
  for (i = 0; i < n; i++)
    bbs_step(&b);
  q = bbs_bits(&b, 32);
  bbs_wrap(&b);

  bbs_rewn(&b, &bp, n + (32 + b.k - 1) / b.k);
  r = bbs_bits(&b, 32);

  if (r != p) {
    fputs("\n*** bbs rewind failure\n", stderr);
    fputs("p = ", stderr); mp_writefile(bp.p, stderr, 10); fputc('\n', stderr);
    fputs("q = ", stderr); mp_writefile(bp.q, stderr, 10); fputc('\n', stderr);
    fputs("n = ", stderr); mp_writefile(bp.n, stderr, 10); fputc('\n', stderr);
    fputs("x = ", stderr); mp_writefile(x, stderr, 10); fputc('\n', stderr);
    fprintf(stderr, "stepped %lu back\n", n + (32 + b.k - 1) / b.k);
    fprintf(stderr, "expected output = %08lx, found %08lx\n",
	    (unsigned long)p, (unsigned long)r);
    ok = 0;
  }

  bbs_seed(&b, x);
  bbs_ffn(&b, &bp, n);
  r = bbs_bits(&b, 32);

  if (q != r) {
    fputs("\n*** bbs fastforward failure\n", stderr);
    fputs("p = ", stderr); mp_writefile(bp.p, stderr, 10); fputc('\n', stderr);
    fputs("q = ", stderr); mp_writefile(bp.q, stderr, 10); fputc('\n', stderr);
    fputs("n = ", stderr); mp_writefile(bp.n, stderr, 10); fputc('\n', stderr);
    fputs("x = ", stderr); mp_writefile(x, stderr, 10); fputc('\n', stderr);
    fprintf(stderr, "stepped %lu back\n", n + (32 + b.k - 1) / b.k);
    fprintf(stderr, "expected output = %08lx, found %08lx\n",
	    (unsigned long)q, (unsigned long)r);
    ok = 0;
  }

  bbs_destroy(&b);
  mp_drop(bp.p);
  mp_drop(bp.q);
  mp_drop(bp.n);
  mp_drop(x);

  assert(mparena_count(MPARENA_GLOBAL) == 0);
  return (ok);
}

static test_chunk tests[] = {
  { "bbs-jump", verify, { &type_mp, &type_mp, &type_mp, &type_ulong, 0 } },
  { 0, 0, { 0 } }
};

int main(int argc, char *argv[])
{
  sub_init();
  test_run(argc, argv, tests, SRCDIR "/t/bbs");
  return (0);
}

#endif

/*----- That's all, folks -------------------------------------------------*/
Commit	Line	Data
2c52abe6	1	/* --c--
2c52abe6	2	*
2c52abe6	3	* Jumping around a BBS sequence
	4	*
	5	* (c) 1999 Straylight/Edgeware
	6	*/
	7
45c0fd36	8	/----- Licensing notice --------------------------------------------------
2c52abe6	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	*
2c52abe6	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	*
2c52abe6	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
2c52abe6	28	/----- Header files ------------------------------------------------------/
	29
	30	#include "bbs.h"
	31	#include "mp.h"
	32	#include "mpbarrett.h"
	33	#include "mpcrt.h"
	34	#include "mpint.h"
	35
	36	/----- Main code ---------------------------------------------------------/
	37
	38	/* --- @jump@ --- *
	39	*
	40	* Arguments: @bbs *b@ = pointer to BBS generator context
2cc3d1d0	41	* @const bbs_priv *bp@ = pointer to BBS modulus factors
2cc3d1d0	42	* @mp *n@ = number of steps to move
2c52abe6	43	* @mp *px@ = exponent mod @p@ for a one-step jump
	44	* @mp *qx@ = exponent mod @q@ for a one-step jump
	45	*
	46	* Returns: ---
	47	*
	48	* Use: Jumps a BBS context a certain number of places (assuming the
	49	* arguments are right).
	50	*
	51	* Let the BBS modulus be %$n = pq$% and the current residue be
	52	* %$x$%. Then the computations performed are:
	53	*
	54	* * Calculate %$x_p = x \bmod p$% and %$x_q = x \bmod q$%.
	55	*
	56	* * Determine %$e_p = px^n \bmod (p - 1)$% and similarly
	57	* %$e_q = qx^n \bmod (p - 1)$%.
	58	*
	59	* * Calculate %$x_p' = x_p^{e_p} \bmod p$% and
	60	* %$x_q' = x_q^{e_q} \bmod q$%.
	61	*
	62	* * Combine %$x_p'$% and %$x_q'$% using the Chinese Remainder
	63	* Theorem.
	64	*
	65	* If you want to step the generator forwards, simply set
	66	* %$px = qx = 2$%. If you want to step backwards, make
38f31040	67	* %$px = (p + 1)/4$% and %$qx = (q + 1)/4$%. Note that, if
2c52abe6	68	* %$x$% is a quadratic residue mod $%p$%, then
	69	*
	70	* %$(x^2) ^ {(p + 1)/4}$%
	71	* %${} = x^{(p + 1)/2}$%
	72	* %${} = x \cdot x^{(p - 1)/2}$%
	73	* %${} = x$%
	74	*
	75	* Simple, no? (Note that the division works because
	76	* %$p \equiv 3 \pmod 4$%.)
	77	*/
	78
2cc3d1d0	79	static void jump(bbs b, const bbs_priv bp, mp *n,
2c52abe6	80	mp px, mp qx)
	81	{
	82	mp ep, eq;
	83	mp *v[2] = { MP_NEW, MP_NEW };
	84
	85	/* --- First work out the exponents --- */
	86
	87	{
	88	mpbarrett mb;
	89	mp *m;
2c52abe6	90
2c52abe6	91	m = mp_sub(MP_NEW, bp->p, MP_ONE);
2c52abe6	92	mpbarrett_create(&mb, m);
2cc3d1d0	93	ep = mpbarrett_exp(&mb, MP_NEW, px, n);
2c52abe6	94	mpbarrett_destroy(&mb);
	95	if (qx == px)
	96	eq = MP_COPY(ep);
	97	else {
	98	m = mp_sub(m, bp->q, MP_ONE);
	99	mpbarrett_create(&mb, m);
2cc3d1d0	100	eq = mpbarrett_exp(&mb, MP_NEW, qx, n);
2c52abe6	101	mpbarrett_destroy(&mb);
	102	}
	103
	104	mp_drop(m);
2c52abe6	105	}
	106
	107	/* --- Now calculate the residues of @x@ --- */
	108
	109	mp_div(0, &v[0], b->x, bp->p);
	110	mp_div(0, &v[1], b->x, bp->q);
	111
	112	/* --- Exponentiate --- */
	113
	114	{
	115	mpbarrett mb;
	116
	117	mpbarrett_create(&mb, bp->p);
	118	v[0] = mpbarrett_exp(&mb, v[0], v[0], ep);
	119	mpbarrett_destroy(&mb);
	120
	121	mpbarrett_create(&mb, bp->q);
	122	v[1] = mpbarrett_exp(&mb, v[1], v[1], eq);
	123	mpbarrett_destroy(&mb);
	124
	125	mp_drop(ep);
	126	mp_drop(eq);
	127	}
	128
	129	/* --- Sort out the result using the Chinese Remainder Theorem --- */
	130
	131	{
	132	mpcrt_mod mv[2];
	133	mpcrt c;
	134	int i;
	135
	136	mv[0].m = MP_COPY(bp->p);
	137	mv[1].m = MP_COPY(bp->q);
	138	for (i = 0; i < 2; i++)
	139	mv[i].n = mv[i].ni = mv[i].nni = MP_NEW;
	140	mpcrt_create(&c, mv, 2, b->mb.m);
	141	b->x = mpcrt_solve(&c, b->x, v);
	142	mpcrt_destroy(&c);
	143	}
	144
	145	/* --- Tidy away --- */
	146
	147	mp_drop(v[0]);
	148	mp_drop(v[1]);
	149	b->r = b->x->v[0];
	150	b->b = b->k;
	151	}
	152
2cc3d1d0	153	/* --- @bbs_{ff,rew}{,n}@ --- *
2c52abe6	154	*
2c52abe6	155	* Arguments: @bbs *b@ = pointer to a BBS generator state
2cc3d1d0	156	* @const bbs_priv *bp@ = pointer to BBS modulus factors
2cc3d1d0	157	* @mp *n@, @unsigned long n@ = number of steps to make
2c52abe6	158	*
	159	* Returns: ---
	160	*
2cc3d1d0	161	* Use: `Fast-forwards' or rewinds a Blum-Blum-Shub generator by @n@
	162	* steps. The @...n@ versions take an @unsigned long@ argument;
	163	* the non-@...n@ versions a multiprecision integer. If @n@ is
	164	* negative then the generator is stepped in the reverse
45c0fd36	165	* direction.
2c52abe6	166	*/
2c52abe6	167
2cc3d1d0	168	static void ff(bbs b, const bbs_priv bp, mp *n)
2cc3d1d0	169	{ jump(b, bp, n, MP_TWO, MP_TWO); }
2c52abe6	170
2cc3d1d0	171	static void rew(bbs b, const bbs_priv bp, mp *n)
2c52abe6	172	{
	173	mp *px = mp_lsr(MP_NEW, bp->p, 2);
	174	mp *qx = mp_lsr(MP_NEW, bp->q, 2);
	175	px = mp_add(px, px, MP_ONE);
	176	qx = mp_add(qx, qx, MP_ONE);
	177	jump(b, bp, n, px, qx);
	178	mp_drop(px);
	179	mp_drop(qx);
	180	}
	181
2cc3d1d0	182	void bbs_ff(bbs b, const bbs_priv bp, mp *n)
	183	{
	184	if (!MP_NEGP(n))
	185	ff(b, bp, n);
	186	else {
	187	n = mp_neg(MP_NEW, n);
	188	rew(b, bp, n);
	189	mp_drop(n);
	190	}
	191	}
	192
	193	void bbs_ffn(bbs b, const bbs_priv bp, unsigned long n)
	194	{ mp *nn = mp_fromulong(MP_NEW, n); ff(b, bp, nn); mp_drop(nn); }
	195
	196	void bbs_rew(bbs b, const bbs_priv bp, mp *n)
	197	{
	198	if (!MP_NEGP(n))
	199	rew(b, bp, n);
	200	else {
	201	n = mp_neg(MP_NEW, n);
	202	ff(b, bp, n);
	203	mp_drop(n);
	204	}
	205	}
	206
	207	void bbs_rewn(bbs b, const bbs_priv bp, unsigned long n)
	208	{ mp *nn = mp_fromulong(MP_NEW, n); bbs_rew(b, bp, nn); mp_drop(nn); }
	209
2c52abe6	210	/----- Test rig ----------------------------------------------------------/
	211
	212	#ifdef TEST_RIG
	213
	214	static int verify(dstr *v)
	215	{
a22bbdf6	216	bbs_priv bp;
2c52abe6	217	bbs b;
	218	mp *x;
	219	unsigned long n;
	220	int ok = 1;
	221	uint32 p, q, r;
	222	int i;
	223
	224	bp.p = (mp *)v[0].buf;
	225	bp.q = (mp *)v[1].buf;
	226	bp.n = mp_mul(MP_NEW, bp.p, bp.q);
	227	x = (mp *)v[2].buf;
	228	n = (unsigned long )v[3].buf;
	229
	230	bbs_create(&b, bp.n, x);
	231	p = bbs_bits(&b, 32);
	232
	233	bbs_seed(&b, x);
	234	for (i = 0; i < n; i++)
	235	bbs_step(&b);
	236	q = bbs_bits(&b, 32);
	237	bbs_wrap(&b);
	238
2cc3d1d0	239	bbs_rewn(&b, &bp, n + (32 + b.k - 1) / b.k);
2c52abe6	240	r = bbs_bits(&b, 32);
	241
	242	if (r != p) {
	243	fputs("\n*** bbs rewind failure\n", stderr);
	244	fputs("p = ", stderr); mp_writefile(bp.p, stderr, 10); fputc('\n', stderr);
	245	fputs("q = ", stderr); mp_writefile(bp.q, stderr, 10); fputc('\n', stderr);
	246	fputs("n = ", stderr); mp_writefile(bp.n, stderr, 10); fputc('\n', stderr);
	247	fputs("x = ", stderr); mp_writefile(x, stderr, 10); fputc('\n', stderr);
	248	fprintf(stderr, "stepped %lu back\n", n + (32 + b.k - 1) / b.k);
	249	fprintf(stderr, "expected output = %08lx, found %08lx\n",
	250	(unsigned long)p, (unsigned long)r);
	251	ok = 0;
	252	}
	253
	254	bbs_seed(&b, x);
2cc3d1d0	255	bbs_ffn(&b, &bp, n);
2c52abe6	256	r = bbs_bits(&b, 32);
	257
	258	if (q != r) {
	259	fputs("\n*** bbs fastforward failure\n", stderr);
	260	fputs("p = ", stderr); mp_writefile(bp.p, stderr, 10); fputc('\n', stderr);
	261	fputs("q = ", stderr); mp_writefile(bp.q, stderr, 10); fputc('\n', stderr);
	262	fputs("n = ", stderr); mp_writefile(bp.n, stderr, 10); fputc('\n', stderr);
	263	fputs("x = ", stderr); mp_writefile(x, stderr, 10); fputc('\n', stderr);
	264	fprintf(stderr, "stepped %lu back\n", n + (32 + b.k - 1) / b.k);
	265	fprintf(stderr, "expected output = %08lx, found %08lx\n",
	266	(unsigned long)q, (unsigned long)r);
	267	ok = 0;
	268	}
	269
	270	bbs_destroy(&b);
	271	mp_drop(bp.p);
	272	mp_drop(bp.q);
	273	mp_drop(bp.n);
	274	mp_drop(x);
	275
	276	assert(mparena_count(MPARENA_GLOBAL) == 0);
	277	return (ok);
	278	}
	279
	280	static test_chunk tests[] = {
	281	{ "bbs-jump", verify, { &type_mp, &type_mp, &type_mp, &type_ulong, 0 } },
	282	{ 0, 0, { 0 } }
	283	};
	284
	285	int main(int argc, char *argv[])
	286	{
	287	sub_init();
0f00dc4c	288	test_run(argc, argv, tests, SRCDIR "/t/bbs");
2c52abe6	289	return (0);
	290	}
	291
	292	#endif
	293
	294	/----- That's all, folks -------------------------------------------------/