Makefile: Link tests against stuff like -lm.

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

/* -*-c-*-
 *
 * $Id$
 *
 * 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 "/tests/bbs");
  return (0);
}

#endif

/*----- That's all, folks -------------------------------------------------*/