1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
/*
Copyright (C) 2014 William Hart
This file is part of FLINT.
FLINT is free software: you can redistribute it and/or modify it under
the terms of the GNU Lesser General Public License (LGPL) as published
by the Free Software Foundation; either version 3 of the License, or
(at your option) any later version. See <https://www.gnu.org/licenses/>.
*/
#include "fmpz.h"
int fmpz_is_strong_probabprime(const fmpz_t n, const fmpz_t base)
{
fmpz_t a, nm1, t, y;
int res = 0;
if (fmpz_cmp_ui(n, 1) <= 0)
return 0;
fmpz_init(a);
fmpz_init(t);
fmpz_init(nm1);
fmpz_sub_ui(nm1, n, 1);
if (fmpz_cmp(base, n) >= 0)
fmpz_mod(a, base, n);
else
fmpz_set(a, base);
if (fmpz_is_one(a) || fmpz_equal(a, nm1) || fmpz_is_zero(a))
res = 1;
else
{
slong s = 0;
fmpz_init(y);
s = fmpz_val2(nm1);
fmpz_tdiv_q_2exp(t, nm1, s);
fmpz_powm(y, a, t, n);
if (fmpz_is_one(y))
res = 1;
else
{
res = fmpz_equal(y, nm1);
for (s--; s > 0 && !res; s--)
{
fmpz_mul(t, y, y);
fmpz_mod(y, t, n);
res = fmpz_equal(y, nm1);
}
}
fmpz_clear(y);
}
fmpz_clear(nm1);
fmpz_clear(a);
fmpz_clear(t);
return res;
}