math/scaf.c: Fix trivial typo.

[catacomb] / math / mp-jacobi.c

/* -*-c-*-
 *
 * 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
 *		@mp *n@ = another integer
 *
 * Returns:	@-1@, @0@ or @1@ -- the Jacobi symbol %$J(a, n)$%.
 *
 * Use:		Computes the Kronecker symbol %$\jacobi{a}{n}$%.  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.
 *
 *		If @n@ is composite, then this computes the Kronecker symbol
 *
 *		  %$\jacobi{a}{n}=\jacobi{a}{u}\prod_i\jacobi{a}{p_i}^{e_i}$%
 *
 *		where %$n = u p_0^{e_0} \ldots p_{n-1}^{e_{n-1}}$% is the
 *		prime factorization of %$n$%.  The missing bits are:
 *
 *		  * %$\jacobi{a}{1} = 1$%;
 *		  * %$\jacobi{a}{-1} = 1$% if @a@ is negative, or 1 if
 *		    positive;
 *		  * %$\jacobi{a}{0} = 0$%;
 *		  * %$\jacobi{a}{2}$ is 0 if @a@ is even, 1 if @a@ is
 *		    congruent to 1 or 7 (mod 8), or %$-1$% otherwise.
 *
 *		If %$n$% is positive and odd, then this is the Jacobi
 *		symbol.  (The Kronecker symbol is a consistent domain
 *		extension; the Jacobi symbol was implemented first, and the
 *		name stuck.)
 */

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

  /* --- Handle zero specially --- *
   *
   * I can't find any specific statement for what to do when %$n = 0$%; PARI
   * opts to set %$\jacobi{\pm1}{0} = \pm 1$% and %$\jacobi{a}{0} = 0$% for
   * other %$a$%.
   */

  if (MP_ZEROP(n)) {
    if (MP_EQ(a, MP_ONE)) return (+1);
    else if (MP_EQ(a, MP_MONE)) return (-1);
    else return (0);
  }

  /* --- Deal with powers of two --- *
   *
   * This implicitly takes a copy of %$n$%.  Copy %$a$% at the same time to
   * make cleanup easier.
   */

  MP_COPY(a);
  n = mp_odd(MP_NEW, n, &p2);
  if (p2) {
    if (MP_EVENP(a)) {
      s = 0;
      goto done;
    } else if ((p2 & 1) && ((a->v[0] & 7) == 3 || (a->v[0] & 7) == 5))
      s = -s;
  }

  /* --- Deal with negative %$n$% --- */

  if (MP_NEGP(n)) {
    n = mp_neg(n, n);
    if (MP_NEGP(a))
      s = -s;
  }

  /* --- Check for unit %$n$% --- */

  if (MP_EQ(n, MP_ONE))
    goto done;

  /* --- Reduce %$a$% modulo %$n$% --- */

  if (MP_NEGP(a) || MP_CMP(a, >=, n))
    mp_div(0, &a, a, n);

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

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

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

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

    /* --- Strip out powers of two from %$a$% --- */

    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;

    /* --- Reduce and swap, applying quadratic residuosity --- */

    if ((nn & 3) == 3 && (a->v[0] & 3) == 3)
      s = -s;
    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 "/t/mp");
  return (0);
}

#endif

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