General utilities cleanup. Add signature support to catcrypt. Throw in

[u/mdw/catacomb] / mp-jacobi.c

/* -*-c-*-
 *
 * $Id: mp-jacobi.c,v 1.5 2004/04/08 01:36:15 mdw Exp $
 *
 * Compute Jacobi symbol
 *
 * (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 "mp.h"

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

/* --- @mp_jacobi@ --- *
 *
 * Arguments:	@mp *a@ = an integer less than @n@
 *		@mp *n@ = an odd integer
 *
 * Returns:	@-1@, @0@ or @1@ -- the Jacobi symbol %$J(a, n)$%.
 *
 * Use:		Computes the Jacobi symbol.  If @n@ is prime, this is the
 *		Legendre symbol and is equal to 1 if and only if @a@ is a
 *		quadratic residue mod @n@.  The result is zero if and only if
 *		@a@ and @n@ have a common factor greater than one.
 */

int mp_jacobi(mp *a, mp *n)
{
  int s = 1;

  assert(((void)"n must be odd in mp_jacobi", MP_LEN(n) && (n->v[0] & 1)));

  /* --- Take copies of the arguments --- */

  a = MP_COPY(a);
  n = MP_COPY(n);

  /* --- Main recursive mess, flattened out into something nice --- */

  for (;;) {
    mpw nn;
    size_t e;

    /* --- Some simple special cases --- */

    MP_SHRINK(a);
    if (MP_LEN(a) == 0) {
      s = 0;
      goto done;
    }

    /* --- Main case with powers of two --- */

    a = mp_odd(a, a, &e);
    nn = n->v[0] & 7;
    if ((e & 1) && (nn == 3 || nn == 5))
      s = -s;
    if (MP_LEN(a) == 1 && a->v[0] == 1)
      goto done;
    if ((nn & 3) == 3 && (a->v[0] & 3) == 3)
      s = -s;

    /* --- Reduce and swap --- */

    mp_div(0, &n, n, a);
    { mp *t = n; n = a; a = t; }
  }

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

done:
  MP_DROP(a);
  MP_DROP(n);
  return (s);
}

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

#ifdef TEST_RIG

#include <mLib/testrig.h>

static int verify(dstr *v)
{
  mp *a = *(mp **)v[0].buf;
  mp *n = *(mp **)v[1].buf;
  int s = *(int *)v[2].buf;
  int j = mp_jacobi(a, n);
  int ok = 1;

  if (s != j) {
    fputs("\n*** fail", stderr);
    fputs("a = ", stderr); mp_writefile(a, stderr, 10); fputc('\n', stderr);
    fputs("n = ", stderr); mp_writefile(n, stderr, 10); fputc('\n', stderr);
    fprintf(stderr, "s = %i\n", s);
    fprintf(stderr, "j = %i\n", j);
    ok = 0;
  }

  mp_drop(a);
  mp_drop(n);
  assert(mparena_count(MPARENA_GLOBAL) == 0);
  return (ok);
}

static test_chunk tests[] = {
  { "jacobi", verify, { &type_mp, &type_mp, &type_int, 0 } },
  { 0, 0, { 0 } }
};

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

#endif

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