X-Git-Url: https://git.distorted.org.uk/u/mdw/catacomb/blobdiff_plain/ba6e6b64033b1f9de49feccb5c9cd438354481f7..0f00dc4c8eb47e67bc0f148c2dd109f73a451e0a:/mp-jacobi.c?ds=sidebyside diff --git a/mp-jacobi.c b/mp-jacobi.c deleted file mode 100644 index 3674f22..0000000 --- a/mp-jacobi.c +++ /dev/null @@ -1,203 +0,0 @@ -/* -*-c-*- - * - * $Id$ - * - * 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 consistant 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; - } - - /* --- 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 - -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 -------------------------------------------------*/