1 /*
2 Copyright (C) 2015 Vladimir Glazachev
3
4 This file is part of FLINT.
5
6 FLINT is free software: you can redistribute it and/or modify it under
7 the terms of the GNU Lesser General Public License (LGPL) as published
8 by the Free Software Foundation; either version 2.1 of the License, or
9 (at your option) any later version. See <http://www.gnu.org/licenses/>.
10 */
11
12 #include "aprcl.h"
13
14 /*
15 Returns 1 if \tau^{\sigma_n-n}(\chi)=-1; otherwise returns 0.
16 It is equal to check:
17 (\chi(-1) * q)^((n - 1) / 2) congruent -1 mod n
18
19 \tau is the Gauss sum;
20 \chi is the character of conductor q and order 2 (quadratic character);
21 \tau(\chi)^\sigma_n means \sigma_n(\tau(\chi)),
22 there \sigma_n is the ring automorphism
23 */
24 int
_aprcl_is_gausspower_2q_equal_first(ulong q,const fmpz_t n)25 _aprcl_is_gausspower_2q_equal_first(ulong q, const fmpz_t n)
26 {
27 int result;
28 fmpz_t npow, nval, ncmp;
29
30 fmpz_init_set(npow, n);
31 fmpz_init_set_ui(nval, q);
32 fmpz_init_set(ncmp, n);
33
34 /* ncmp = -1 mod n */
35 fmpz_sub_ui(ncmp, ncmp, 1);
36
37 /* nval = (\chi(-1) * q) = ((-1)^((q - 1) / 2) * q) */
38 if ((q - 1) % 2 == 1)
39 {
40 fmpz_neg(nval, nval);
41 fmpz_add(nval, nval, n);
42 }
43
44 /* npow = (n - 1) / 2 */
45 fmpz_sub_ui(npow, npow, 1);
46 fmpz_fdiv_q_2exp(npow, npow, 1);
47
48 /* nval = (\chi(-1) * q)^((n - 1) / 2) mod n */
49 fmpz_powm(nval, nval, npow, n);
50
51 result = 0;
52 if (fmpz_equal(nval, ncmp))
53 result = 1;
54
55 fmpz_clear(npow);
56 fmpz_clear(nval);
57 fmpz_clear(ncmp);
58
59 return result;
60 }
61
62 /*
63 Returns 1 if \tau^{\sigma_n-n}(\chi^{p / 2}) = -1; otherwise returns 0.
64 It is equal to check:
65 q^((n - 1) / 2) congruent -1 mod n
66
67 \tau is the Gauss sum;
68 \chi is the character of conductor q and order 2 (quadratic character);
69 \tau(\chi)^\sigma_n means \sigma_n(\tau(\chi)),
70 there \sigma_n is the ring automorphism
71 */
72 int
_aprcl_is_gausspower_2q_equal_second(ulong q,const fmpz_t n)73 _aprcl_is_gausspower_2q_equal_second(ulong q, const fmpz_t n)
74 {
75 int result;
76 fmpz_t npow, nval, ncmp;
77
78 fmpz_init_set(npow, n);
79 fmpz_init_set_ui(nval, q);
80 fmpz_init_set(ncmp, n);
81
82 /* ncmp = -1 mod n */
83 fmpz_sub_ui(ncmp, ncmp, 1);
84
85 /* nval = q^((n - 1) / 2) mod n */
86 fmpz_sub_ui(npow, npow, 1);
87 fmpz_fdiv_q_2exp(npow, npow, 1);
88 fmpz_powm(nval, nval, npow, n);
89
90 result = 0;
91 if (fmpz_equal(nval, ncmp))
92 result = 1;
93
94 fmpz_clear(npow);
95 fmpz_clear(nval);
96 fmpz_clear(ncmp);
97
98 return result;
99 }
100
101 /*
102 Returns non-negative value if \tau(\chi)^(\sigma_n - n) from <\zeta_p>
103 cyclic group. It is equal to check that for some i from 0 to p - 1
104 \tau(\chi^n) == \zeta_p^i * \tau^n(\chi). If such i exists returns i.
105 Otherwise returns -1.
106
107 \tau is the Gauss sum;
108 \chi is the character of conductor q and order r;
109 \tau(\chi)^\sigma_n means \sigma_n(\tau(\chi)),
110 there \sigma_n is the ring automorphism
111 */
112 slong
_aprcl_is_gausspower_from_unity_p(ulong q,ulong r,const fmpz_t n)113 _aprcl_is_gausspower_from_unity_p(ulong q, ulong r, const fmpz_t n)
114 {
115 slong result;
116 ulong i;
117 unity_zpq temp, gauss, gausspow, gausssigma;
118
119 unity_zpq_init(gauss, q, r, n);
120 unity_zpq_init(gausssigma, q, r, n);
121 unity_zpq_init(gausspow, q, r, n);
122 unity_zpq_init(temp, q, r, n);
123
124 /* gauss = \tau(\chi) */
125 unity_zpq_gauss_sum(gauss, q, r);
126 /* gausssigma = \tau(\chi^n) */
127 unity_zpq_gauss_sum_sigma_pow(gausssigma, q, r);
128 /* gausspow = \tau^n(\chi) */
129 unity_zpq_pow(gausspow, gauss, n);
130
131 result = -1;
132 for (i = 0; i < r; i++)
133 {
134 /* temp = \zeta_p^i * \tau^n(\chi) */
135 unity_zpq_mul_unity_p_pow(temp, gausspow, i);
136 if (unity_zpq_equal(gausssigma, temp))
137 {
138 result = i;
139 break;
140 }
141 }
142
143 unity_zpq_clear(gauss);
144 unity_zpq_clear(gausssigma);
145 unity_zpq_clear(gausspow);
146 unity_zpq_clear(temp);
147
148 return result;
149 }
150
151 primality_test_status
_aprcl_is_prime_gauss(const fmpz_t n,const aprcl_config config)152 _aprcl_is_prime_gauss(const fmpz_t n, const aprcl_config config)
153 {
154 int *lambdas;
155 ulong i, j, k, nmod4;
156 primality_test_status result;
157
158 /*
159 Condition (Lp) is satisfied iff:
160 For each p | R we must show that for all prime r | n and all
161 positive integers a there exists l such that:
162 r^(p-1) congruent n^(l(p-1)) mod p^a
163
164 If (Lp), then lambdas_p = 3.
165 */
166 lambdas = (int*) flint_malloc(sizeof(int) * config->rs.num);
167 for (i = 0; i < config->rs.num; i++)
168 lambdas[i] = 0;
169
170 result = PROBABPRIME;
171
172 /* nmod4 = n % 4 */
173 nmod4 = fmpz_tdiv_ui(n, 4);
174
175 /* for every prime q | s */
176 for (i = 0; i < config->qs->num; i++)
177 {
178 n_factor_t q_factors;
179 ulong q;
180 if (result == COMPOSITE) break;
181
182 q = fmpz_get_ui(config->qs->p + i);
183
184 /* n == q, q - prime => n - prime */
185 if (fmpz_equal_ui(n, q))
186 {
187 result = PRIME;
188 break;
189 }
190
191 /* find prime factors of q - 1 */
192 n_factor_init(&q_factors);
193 n_factor(&q_factors, q - 1, 1);
194
195 /* for every prime p | q - 1 */
196 for (j = 0; j < q_factors.num; j++)
197 {
198 int state, pind;
199 ulong p;
200 if (result == COMPOSITE) break;
201
202 p = q_factors.p[j];
203
204 pind = _aprcl_p_ind(config, p);
205 state = lambdas[pind];
206
207 /*
208 (Lp.a)
209 if p == 2 and n = 1 mod 4 then (Lp) is equal to:
210 for quadratic character \chi (\tau(\chi))^(\sigma_n-n) = -1
211 */
212 if (p == 2 && state == 0 && nmod4 == 1)
213 {
214 if (_aprcl_is_gausspower_2q_equal_first(q, n) == 1)
215 {
216 state = 3;
217 lambdas[pind] = state;
218 }
219 }
220
221 /*
222 (Lp.b)
223 if p == 2, r = 2^k >= 4 and n = 3 mod 4 then (Lp) is equal to:
224 1) for quadratic character \chi
225 (\tau(\chi^(r / 2)))^(\sigma_n-n) = -1
226 2) for character \chi = \chi_{r, q}
227 (\tau(\chi))^(\sigma_n-n) is a generator of cyclic
228 group <\zeta_r>
229
230 if 1) is true, then lambdas_p = 1
231 if 2) is true, then lambdas_p = 2
232 if 1) and 2) is true, then lambdas_p = 3
233 */
234 if (p == 2 && (state == 0 || state == 2) && nmod4 == 3)
235 {
236 if (_aprcl_is_gausspower_2q_equal_second(q, n) == 1)
237 {
238 if (state == 2)
239 state = 3;
240 else
241 state = 1;
242 lambdas[pind] = state;
243 }
244 }
245
246 /* for every prime power p^k | q - 1 */
247 for (k = 1; k <= q_factors.exp[j]; k++)
248 {
249 int unity_power;
250 ulong r;
251
252 /* r = p^k */
253 r = n_pow(p, k);
254
255 /* if gcd(q*r, n) != 1 */
256 if (aprcl_is_mul_coprime_ui_ui(q, r, n) == 0)
257 {
258 result = COMPOSITE;
259 break;
260 }
261
262 /*
263 if exists z such that \tau(\chi^n) = \zeta_r^z*\tau^n(\chi)
264 unity_power = z; otherwise unity_power = -1
265 */
266 unity_power = _aprcl_is_gausspower_from_unity_p(q, r, n);
267
268 /* if unity_power < 0 then n is composite */
269 if (unity_power < 0)
270 {
271 result = COMPOSITE;
272 break;
273 }
274
275 /*
276 (Lp.c)
277 if p > 2 then (Lp) is equal to:
278 (\tau(\chi))^(\sigma_n - n) is a generator of cyclic
279 group <\zeta_r>
280 */
281 if (p > 2 && state == 0 && unity_power > 0)
282 {
283 ulong upow = unity_power;
284 /*
285 if gcd(r, unity_power) = 1 then
286 (\tau(\chi))^(\sigma_n - n) is a generator
287 */
288 if (n_gcd(r, upow) == 1)
289 {
290 state = 3;
291 lambdas[pind] = state;
292 }
293 }
294
295 /*
296 (Lp.b)
297 check 2) of (Lp) if p == 2 and nmod4 == 3
298 */
299 if (p == 2 && unity_power > 0
300 && (state == 0 || state == 1) && nmod4 == 3)
301 {
302 ulong upow = unity_power;
303 if (n_gcd(r, upow) == 1)
304 {
305 if (state == 0)
306 {
307 state = 2;
308 lambdas[pind] = state;
309 }
310 if (state == 1)
311 {
312 state = 3;
313 lambdas[pind] = state;
314 }
315 }
316 }
317 }
318 }
319 }
320
321 /*
322 if for some p we have not proved (Lp)
323 then n can be as prime or composite
324 */
325 if (result == PROBABPRIME)
326 for (i = 0; i < config->rs.num; i++)
327 if (lambdas[i] != 3)
328 result = UNKNOWN;
329
330
331 /* if n can be prime we do final division */
332 if (result == UNKNOWN || result == PROBABPRIME)
333 {
334 int f_division;
335 f_division = aprcl_is_prime_final_division(n, config->s, config->R);
336 /* if (Lp) is true for all p | R and f_division == 1 then n - prime */
337 if (result == PROBABPRIME && f_division == 1)
338 result = PRIME;
339 /*
340 if we not prove (Lp) for some p and f_division == 1
341 then n still can be prime
342 */
343 if (result == UNKNOWN && f_division == 1)
344 result = PROBABPRIME;
345 /* if f_division == 0 so we find the divisor of n, so n - composite */
346 if (f_division == 0)
347 result = COMPOSITE;
348 }
349
350 flint_free(lambdas);
351
352 return result;
353 }
354
355 int
aprcl_is_prime_gauss_min_R(const fmpz_t n,ulong R)356 aprcl_is_prime_gauss_min_R(const fmpz_t n, ulong R)
357 {
358 primality_test_status result;
359 aprcl_config config;
360
361 aprcl_config_gauss_init_min_R(config, n, R);
362 result = _aprcl_is_prime_gauss(n, config);
363
364 aprcl_config_gauss_clear(config);
365
366 if (result == PRIME)
367 return 1;
368 return 0;
369 }
370
371 int
aprcl_is_prime_gauss(const fmpz_t n)372 aprcl_is_prime_gauss(const fmpz_t n)
373 {
374 ulong R;
375 primality_test_status result;
376 aprcl_config config;
377
378 if (fmpz_cmp_ui(n, 2) < 0)
379 return 0;
380
381 aprcl_config_gauss_init_min_R(config, n, 180);
382 result = _aprcl_is_prime_gauss(n, config);
383 R = config->R;
384 aprcl_config_gauss_clear(config);
385
386 /*
387 if result == PROBABPRIME it means that we have
388 not proved (Lp) for some p (most likely we fail L.c step);
389 we can try to use bigger R
390 */
391 if (result == PROBABPRIME)
392 {
393 R = R * 2;
394 aprcl_config_gauss_init_min_R(config, n, R);
395 result = _aprcl_is_prime_gauss(n, config);
396 aprcl_config_gauss_clear(config);
397 }
398
399 if (result == PROBABPRIME)
400 {
401 R = R * 3;
402 aprcl_config_gauss_init_min_R(config, n, R);
403 result = _aprcl_is_prime_gauss(n, config);
404 aprcl_config_gauss_clear(config);
405 }
406
407 if (result == PROBABPRIME)
408 {
409 R = R * 5;
410 aprcl_config_gauss_init_min_R(config, n, R);
411 result = _aprcl_is_prime_gauss(n, config);
412 aprcl_config_gauss_clear(config);
413 }
414
415 if (result == PROBABPRIME || result == UNKNOWN)
416 {
417 flint_printf("aprcl_is_prime_gauss: failed to prove n prime\n");
418 fmpz_print(n); flint_printf("\n");
419 flint_abort();
420 }
421
422 if (result == PRIME)
423 return 1;
424 return 0;
425 }
426