1 #include "cado.h" // IWYU pragma: keep
2 #include "mod64.hpp"
3 #include "modredc64.hpp" // IWYU pragma: keep
4
5 typedef Modulus64 Modulus;
6 #include "mod64_common.cpp"
7 #include "macros.h"
8
9 /* Put 1/s (mod t) in r and return 1 if s is invertible,
10 or set r to 0 and return 0 if s is not invertible */
11
12 bool
inv(Residue & r,const Residue & sp) const13 Modulus64::inv (Residue &r, const Residue &sp) const
14 {
15 int64_t u1, v1;
16 uint64_t u2, v2, s, t;
17 #ifndef NDEBUG
18 Residue tmp(*this);
19 #endif
20
21 #ifndef NDEBUG
22 /* Remember input in case r overwrites it */
23 set (tmp, sp);
24 #endif
25
26 s = get_u64 (sp);
27 t = getmod_u64 ();
28
29 ASSERT (t > 0);
30 ASSERT (s < t);
31
32 if (s == 0)
33 {
34 set0(r); /* Not invertible */
35 return 0;
36 }
37
38 if (s == 1)
39 {
40 set1(r);
41 return 1;
42 }
43
44 u1 = 1;
45 u2 = s;
46 v1 = - (int64_t) (t / s); /* No overflow, since s >= 2 */
47 v2 = t % s;
48
49 if (v2 == 1)
50 {
51 u1 = v1 + t;
52 }
53 else
54 {
55 while (v2 != 0)
56 {
57 uint64_t q;
58 /* unroll twice and swap u/v */
59 q = u2 / v2;
60 ASSERT_EXPENSIVE (q <= (uint64_t) INT64_MAX);
61 u1 = u1 - (int64_t) q * v1;
62 u2 = u2 - q * v2;
63
64 if (u2 == 0)
65 {
66 u1 = v1;
67 u2 = v2;
68 break;
69 }
70
71 q = v2 / u2;
72 ASSERT_EXPENSIVE (q <= (uint64_t) INT64_MAX);
73 v1 = v1 - (int64_t) q * u1;
74 v2 = v2 - q * u2;
75 }
76
77 if (u2 != 1)
78 {
79 /* printf ("s=%lu t=%lu found %lu\n", s[0], t[0], u2); */
80 set0(r); /* non-trivial gcd */
81 return 0;
82 }
83
84 if (u1 < 0)
85 u1 = u1 + t;
86 }
87
88 ASSERT ((uint64_t) u1 < t);
89
90 set (r, (uint64_t) u1);
91
92 #ifndef NDEBUG
93 mul (tmp, tmp, r);
94 ASSERT(is1 (tmp));
95 #endif
96
97 return 1;
98 }
99
100 /* even_inv_lookup_table[i] is 1/(2*i+1) mod 128 */
101 static uint64_t even_inv_lookup_table[64] = {
102 1, 43, 77, 55, 57, 35, 69, 111, 113, 27, 61, 39, 41, 19, 53, 95, 97, 11, 45,
103 23, 25, 3, 37, 79, 81, 123, 29, 7, 9, 115, 21, 63, 65, 107, 13, 119, 121, 99,
104 5, 47, 49, 91, 125, 103, 105, 83, 117, 31, 33, 75, 109, 87, 89, 67, 101, 15,
105 17, 59, 93, 71, 73, 51, 85, 127 } ;
106
107
108 /* Faster modul_inv for the case where m = 2^k */
109 bool
inv_powerof2(Residue & r,const Residue & A) const110 Modulus64::inv_powerof2 (Residue &r, const Residue &A) const
111 {
112 uint64_t x = getmod_u64(), y = get_u64(A);
113
114 ASSERT (!(x & (x-1))); /* assert that x is a power of 2 */
115 ASSERT (y < x);
116 if (!(y & 1))
117 return 0;
118 else
119 {
120 if (!(x >> 4)) /* x = 2, 4 or 8 */
121 set(r, y);
122 else if (!(x >> 8)) /* x = 16, 32, 64, or 128 */
123 set(r, even_inv_lookup_table[(y-1) >> 1] & (x-1));
124 else
125 {
126 uint64_t h = x >> (u64arith_ctz(x) >> 1);
127 Modulus64 m2(h);
128 Residue B(*this);
129 m2.set_reduced (B, (y & (h-1)));
130
131 m2.inv_powerof2 (r, B);
132 uint64_t t = get_u64(r);
133 set(r, (2 * t - t*t*y) & (x-1));
134 }
135 return 1;
136 }
137 }
138
139 /* Faster modul_inv for the case where m is odd */
140 bool
inv_odd(Residue & r,const Residue & A) const141 Modulus64::inv_odd (Residue &r, const Residue &A) const
142 {
143 /* TODO: do REDC trick like in ModulusREDC64 */
144 return inv(r, A);
145 }
146