Commit | Line | Data |
---|---|---|
5b00a0ea | 1 | /* -*-c-*- |
2 | * | |
a69a3efd | 3 | * $Id$ |
5b00a0ea | 4 | * |
5 | * Compute Jacobi symbol | |
6 | * | |
7 | * (c) 1999 Straylight/Edgeware | |
8 | */ | |
9 | ||
45c0fd36 | 10 | /*----- Licensing notice --------------------------------------------------* |
5b00a0ea | 11 | * |
12 | * This file is part of Catacomb. | |
13 | * | |
14 | * Catacomb is free software; you can redistribute it and/or modify | |
15 | * it under the terms of the GNU Library General Public License as | |
16 | * published by the Free Software Foundation; either version 2 of the | |
17 | * License, or (at your option) any later version. | |
45c0fd36 | 18 | * |
5b00a0ea | 19 | * Catacomb is distributed in the hope that it will be useful, |
20 | * but WITHOUT ANY WARRANTY; without even the implied warranty of | |
21 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | |
22 | * GNU Library General Public License for more details. | |
45c0fd36 | 23 | * |
5b00a0ea | 24 | * You should have received a copy of the GNU Library General Public |
25 | * License along with Catacomb; if not, write to the Free | |
26 | * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, | |
27 | * MA 02111-1307, USA. | |
28 | */ | |
29 | ||
5b00a0ea | 30 | /*----- Header files ------------------------------------------------------*/ |
31 | ||
32 | #include "mp.h" | |
33 | ||
34 | /*----- Main code ---------------------------------------------------------*/ | |
35 | ||
36 | /* --- @mp_jacobi@ --- * | |
37 | * | |
6791ed17 MW |
38 | * Arguments: @mp *a@ = an integer |
39 | * @mp *n@ = another integer | |
5b00a0ea | 40 | * |
41 | * Returns: @-1@, @0@ or @1@ -- the Jacobi symbol %$J(a, n)$%. | |
42 | * | |
6791ed17 MW |
43 | * Use: Computes the Kronecker symbol %$\jacobi{a}{n}$%. If @n@ is |
44 | * prime, this is the Legendre symbol and is equal to 1 if and | |
45 | * only if @a@ is a quadratic residue mod @n@. The result is | |
46 | * zero if and only if @a@ and @n@ have a common factor greater | |
47 | * than one. | |
48 | * | |
49 | * If @n@ is composite, then this computes the Kronecker symbol | |
50 | * | |
51 | * %$\jacobi{a}{n}=\jacobi{a}{u}\prod_i\jacobi{a}{p_i}^{e_i}$% | |
52 | * | |
53 | * where %$n = u p_0^{e_0} \ldots p_{n-1}^{e_{n-1}}$% is the | |
54 | * prime factorization of %$n$%. The missing bits are: | |
55 | * | |
56 | * * %$\jacobi{a}{1} = 1$%; | |
57 | * * %$\jacobi{a}{-1} = 1$% if @a@ is negative, or 1 if | |
58 | * positive; | |
59 | * * %$\jacobi{a}{0} = 0$%; | |
60 | * * %$\jacobi{a}{2}$ is 0 if @a@ is even, 1 if @a@ is | |
61 | * congruent to 1 or 7 (mod 8), or %$-1$% otherwise. | |
62 | * | |
63 | * If %$n$% is positive and odd, then this is the Jacobi | |
64 | * symbol. (The Kronecker symbol is a consistant domain | |
65 | * extension; the Jacobi symbol was implemented first, and the | |
66 | * name stuck.) | |
5b00a0ea | 67 | */ |
68 | ||
69 | int mp_jacobi(mp *a, mp *n) | |
70 | { | |
71 | int s = 1; | |
6791ed17 MW |
72 | size_t p2; |
73 | ||
74 | /* --- Handle zero specially --- * | |
75 | * | |
76 | * I can't find any specific statement for what to do when %$n = 0$%; PARI | |
77 | * opts to set %$\jacobi{\pm1}{0} = \pm 1$% and %$\jacobi{a}{0} = 0$% for | |
78 | * other %$a$%. | |
79 | */ | |
80 | ||
81 | if (MP_ZEROP(n)) { | |
82 | if (MP_EQ(a, MP_ONE)) return (+1); | |
83 | else if (MP_EQ(a, MP_MONE)) return (-1); | |
84 | else return (0); | |
85 | } | |
86 | ||
87 | /* --- Deal with powers of two --- * | |
88 | * | |
89 | * This implicitly takes a copy of %$n$%. Copy %$a$% at the same time to | |
90 | * make cleanup easier. | |
91 | */ | |
92 | ||
93 | MP_COPY(a); | |
94 | n = mp_odd(MP_NEW, n, &p2); | |
95 | if (p2) { | |
96 | if (MP_EVENP(a)) { | |
97 | s = 0; | |
98 | goto done; | |
99 | } else if ((p2 & 1) && ((a->v[0] & 7) == 3 || (a->v[0] & 7) == 5)) | |
100 | s = -s; | |
101 | } | |
102 | ||
103 | /* --- Deal with negative %$n$% --- */ | |
104 | ||
105 | if (MP_NEGP(n)) { | |
106 | n = mp_neg(n, n); | |
107 | if (MP_NEGP(a)) | |
108 | s = -s; | |
109 | } | |
110 | ||
111 | /* --- Check for unit %$n$% --- */ | |
5b00a0ea | 112 | |
6791ed17 MW |
113 | if (MP_EQ(n, MP_ONE)) |
114 | goto done; | |
c76161cc | 115 | |
6791ed17 | 116 | /* --- Reduce %$a$% modulo %$n$% --- */ |
5b00a0ea | 117 | |
6791ed17 MW |
118 | if (MP_NEGP(a) || MP_CMP(a, >=, n)) |
119 | mp_div(0, &a, a, n); | |
5b00a0ea | 120 | |
121 | /* --- Main recursive mess, flattened out into something nice --- */ | |
122 | ||
123 | for (;;) { | |
774d32a5 | 124 | mpw nn; |
125 | size_t e; | |
5b00a0ea | 126 | |
127 | /* --- Some simple special cases --- */ | |
128 | ||
129 | MP_SHRINK(a); | |
a69a3efd | 130 | if (MP_ZEROP(a)) { |
5b00a0ea | 131 | s = 0; |
132 | goto done; | |
133 | } | |
134 | ||
774d32a5 | 135 | /* --- Main case with powers of two --- */ |
5b00a0ea | 136 | |
774d32a5 | 137 | a = mp_odd(a, a, &e); |
138 | nn = n->v[0] & 7; | |
139 | if ((e & 1) && (nn == 3 || nn == 5)) | |
140 | s = -s; | |
141 | if (MP_LEN(a) == 1 && a->v[0] == 1) | |
142 | goto done; | |
143 | if ((nn & 3) == 3 && (a->v[0] & 3) == 3) | |
144 | s = -s; | |
5b00a0ea | 145 | |
146 | /* --- Reduce and swap --- */ | |
147 | ||
148 | mp_div(0, &n, n, a); | |
149 | { mp *t = n; n = a; a = t; } | |
150 | } | |
151 | ||
152 | /* --- Wrap everything up --- */ | |
153 | ||
154 | done: | |
155 | MP_DROP(a); | |
156 | MP_DROP(n); | |
157 | return (s); | |
158 | } | |
159 | ||
160 | /*----- Test rig ----------------------------------------------------------*/ | |
161 | ||
162 | #ifdef TEST_RIG | |
163 | ||
164 | #include <mLib/testrig.h> | |
165 | ||
166 | static int verify(dstr *v) | |
167 | { | |
168 | mp *a = *(mp **)v[0].buf; | |
169 | mp *n = *(mp **)v[1].buf; | |
170 | int s = *(int *)v[2].buf; | |
171 | int j = mp_jacobi(a, n); | |
172 | int ok = 1; | |
173 | ||
174 | if (s != j) { | |
175 | fputs("\n*** fail", stderr); | |
176 | fputs("a = ", stderr); mp_writefile(a, stderr, 10); fputc('\n', stderr); | |
177 | fputs("n = ", stderr); mp_writefile(n, stderr, 10); fputc('\n', stderr); | |
178 | fprintf(stderr, "s = %i\n", s); | |
179 | fprintf(stderr, "j = %i\n", j); | |
180 | ok = 0; | |
181 | } | |
182 | ||
183 | mp_drop(a); | |
184 | mp_drop(n); | |
0e895689 | 185 | assert(mparena_count(MPARENA_GLOBAL) == 0); |
5b00a0ea | 186 | return (ok); |
187 | } | |
188 | ||
189 | static test_chunk tests[] = { | |
190 | { "jacobi", verify, { &type_mp, &type_mp, &type_int, 0 } }, | |
191 | { 0, 0, { 0 } } | |
192 | }; | |
193 | ||
194 | int main(int argc, char *argv[]) | |
195 | { | |
196 | sub_init(); | |
197 | test_run(argc, argv, tests, SRCDIR "/tests/mp"); | |
198 | return (0); | |
199 | } | |
200 | ||
201 | #endif | |
202 | ||
203 | /*----- That's all, folks -------------------------------------------------*/ |