1 #include "cado.h" // IWYU pragma: keep
2 #include <stdio.h>
3 #include "pp1.h"
4 #include "pm1.h"
5
6 /* Do we want backtracking when processing factors of 2 in E? */
7 #ifndef PM1_BACKTRACKING
8 /* Default is "yes." Set to 0 for "no." */
9 #define PM1_BACKTRACKING 1
10 #endif
11
12 /* #define PARI */
13
14 void
pm1_stage1(residue_t x,const unsigned long * E,const int E_nrwords,const modulus_t m)15 pm1_stage1 (residue_t x, const unsigned long *E, const int E_nrwords,
16 const modulus_t m)
17 {
18 mod_pow_mp (x, x, E, E_nrwords, m);
19 }
20
21
22 /* Looks for a factor of the modulus m, using the P-1 algorithm.
23 * The parameters of P-1 are given in plan.
24 * Returns 1 if backtracking was used, 0 otherwise.
25 */
26 int
pm1(modint_t f,const modulus_t m,const pm1_plan_t * plan)27 pm1 (modint_t f, const modulus_t m, const pm1_plan_t *plan)
28 {
29 residue_t x, t, X, one, two;
30 unsigned int i;
31 int bt = 0;
32
33 mod_init_noset0 (x, m);
34 mod_init_noset0 (one, m);
35 mod_init_noset0 (two, m);
36 mod_init_noset0 (t, m);
37 mod_set1 (one, m);
38 mod_add (two, one, one, m);
39
40 /* Stage 1, a simple exponentiation ... */
41 mod_2pow_mp (x, plan->E, plan->E_nrwords, m);
42 /* ... except for the backtracking part for the 2's in the exponent */
43 mod_set (t, x, m);
44 for (i = 0; i < plan->exp2; i++)
45 {
46 mod_mul (x, x, x, m);
47 #if PM1_BACKTRACKING
48 if (mod_is1 (x, m))
49 {
50 mod_set (x, t, m);
51 bt = 1;
52 break;
53 }
54 mod_set (t, x, m);
55 #endif
56 }
57
58 #ifdef PARI
59 printf ("E = B1_exponent (%u); x = Mod(2, %lu)^E; x == %lu /* PARI */\n",
60 plan->B1, mod_getmod_ul (m), mod_get_ul (x, m));
61 #endif
62
63 mod_sub (t, x, one, m);
64 mod_gcd (f, t, m);
65
66 if (mod_intcmp_ul (f, 1UL) > 0 || plan->B1 >= plan->stage2.B2)
67 {
68 mod_clear (one, m);
69 mod_clear (two, m);
70 mod_clear (t, m);
71 mod_clear (x, m);
72 return 0;
73 }
74
75 /* Compute X = x + 1/x */
76 mod_init_noset0 (X, m);
77 /* I think that we have the guarantee that m is odd, so that x is
78 * necessarily coprime to all prime factors of m, and is thus in
79 * (Z/mZ)^*. Whence 1/x cannot go bananas.
80 */
81 // coverity[check_return]
82 mod_inv (X, x, m);
83 mod_add (X, X, x, m);
84
85 #ifdef PARI
86 printf ("X = x+1/x; X == %lu /* PARI */\n", mod_get_ul (X, m));
87 #endif
88
89 bt = pp1_stage2 (t, X, &(plan->stage2), two, m);
90 mod_gcd (f, t, m);
91
92 mod_clear (one, m);
93 mod_clear (two, m);
94 mod_clear (t, m);
95 mod_clear (X, m);
96 mod_clear (x, m);
97 return bt;
98 }
99