fn factorize_smooth(
x: u64,
row: &mut [u64],
br_primes: &[BarrettReduction<u64>],
) -> boolExamples found in repository?
crates/competitive/src/math/discrete_logarithm.rs (line 210)
184fn index_calculus_for_primitive_root(
185 p: u64,
186 ord: u64,
187 br_primes: &[BarrettReduction<u64>],
188 prec: &QdrtPowPrec,
189) -> Vec<u64> {
190 let br_ord = BarrettReduction::<u128>::new(ord as u128);
191 let mul = |x: u64, y: u64| br_ord.rem(x as u128 * y as u128) as u64;
192 let sub = |x: u64, y: u64| if x < y { x + ord - y } else { x - y };
193
194 let pc = br_primes.len();
195 let mut mat: Vec<Vec<u64>> = vec![];
196 let mut rows: Vec<Vec<u64>> = vec![];
197
198 let mut rng = Xorshift::default();
199 let br = BarrettReduction::<u128>::new(p as u128);
200
201 for i in 0..pc {
202 for ri in 0usize.. {
203 let mut row = vec![0u64; pc + 1];
204 let mut kk = rng.rand(ord - 1) + 1;
205 let mut gkk = prec.pow(kk, &br);
206 let mut k = kk;
207 let mut gk = gkk;
208 while ri >= rows.len() {
209 row[pc] = k;
210 if factorize_smooth(gk, &mut row, br_primes) {
211 rows.push(row);
212 break;
213 }
214 if k + kk < ord {
215 k += kk;
216 gk = br.rem(gk as u128 * gkk as u128) as u64;
217 } else {
218 kk = rng.rand(ord - 1) + 1;
219 gkk = prec.pow(kk, &br);
220 k = kk;
221 gk = gkk;
222 }
223 }
224 let row = &mut rows[ri];
225 for j in 0..i {
226 if row[j] != 0 {
227 let b = mul(modinv(mat[j][j], ord), row[j]);
228 for (r, a) in row[j..].iter_mut().zip(&mat[j][j..]) {
229 *r = sub(*r, mul(*a, b));
230 }
231 }
232 assert_eq!(row[j], 0);
233 }
234 if gcd(row[i], ord) == 1 {
235 let last = rows.len() - 1;
236 rows.swap(ri, last);
237 mat.push(rows.pop().unwrap());
238 break;
239 }
240 }
241 }
242 for i in (0..pc).rev() {
243 for j in i + 1..pc {
244 mat[i][pc] = sub(mat[i][pc], mul(mat[i][j], mat[j][pc]));
245 }
246 mat[i][pc] = mul(mat[i][pc], modinv(mat[i][i], ord));
247 }
248 (0..pc).map(|i| mat[i][pc]).collect()
249}
250
251#[derive(Debug)]
252struct IndexCalculusWithPrimitiveRoot {
253 p: u64,
254 ord: u64,
255 prec: QdrtPowPrec,
256 coeff: Vec<u64>,
257}
258
259impl IndexCalculusWithPrimitiveRoot {
260 fn new(p: u64, br_primes: &[BarrettReduction<u64>]) -> Self {
261 let ord = p - 1;
262 let g = primitive_root(p);
263 let br = BarrettReduction::<u128>::new(p as u128);
264 let prec = QdrtPowPrec::new(g, ord, &br);
265 let coeff = index_calculus_for_primitive_root(p, ord, br_primes, &prec);
266 Self {
267 p,
268 ord,
269 prec,
270 coeff,
271 }
272 }
273 fn index_calculus(&self, a: u64, br_primes: &[BarrettReduction<u64>]) -> Option<u64> {
274 let p = self.p;
275 let ord = self.ord;
276 let br = BarrettReduction::<u128>::new(p as u128);
277 let a = br.rem(a as _) as u64;
278 if a == 1 {
279 return Some(0);
280 }
281 if p == 2 {
282 return None;
283 }
284
285 let mut rng = Xorshift::new();
286 let mut row = vec![0u64; br_primes.len()];
287 let mut kk = rng.rand(ord - 1) + 1;
288 let mut gkk = self.prec.pow(kk, &br);
289 let mut k = kk;
290 let mut gk = br.rem(gkk as u128 * a as u128) as u64;
291 loop {
292 if factorize_smooth(gk, &mut row, br_primes) {
293 let mut res = ord - k;
294 for (&c, &r) in self.coeff.iter().zip(&row) {
295 for _ in 0..r {
296 res += c;
297 if res >= ord {
298 res -= ord;
299 }
300 }
301 }
302 return Some(res);
303 }
304 if k + kk < ord {
305 k += kk;
306 gk = br.rem(gk as u128 * gkk as u128) as u64;
307 } else {
308 kk = rng.rand(ord - 1) + 1;
309 gkk = self.prec.pow(kk, &br);
310 k = kk;
311 gk = br.rem(gkk as u128 * a as u128) as u64;
312 }
313 }
314 }