pub struct BarrettReduction<T> {
m: T,
im: T,
}Fields§
§m: T§im: TImplementations§
Source§impl<T> BarrettReduction<T>where
T: Barrettable,
impl<T> BarrettReduction<T>where
T: Barrettable,
Sourcepub fn new(m: T) -> Self
pub fn new(m: T) -> Self
Examples found in repository?
More examples
crates/competitive/src/num/mint/mint_basic.rs (line 135)
133 pub fn set_mod(m: u32) {
134 DYN_MODULUS_U32
135 .with(|cell| unsafe { *cell.get() = BarrettReduction::<u64>::new(m as u64) });
136 }
137 fn rem(x: u64) -> u64 {
138 DYN_MODULUS_U32.with(|cell| unsafe { (*cell.get()).rem(x) })
139 }
140}
141impl DynMIntU32 {
142 pub fn set_mod(m: u32) {
143 DynModuloU32::set_mod(m)
144 }
145}
146
147thread_local!(static DYN_MODULUS_U64: UnsafeCell<BarrettReduction<u128>> = const { UnsafeCell::new(BarrettReduction::<u128>::new_with_im(1_000_000_007, !0 / 1_000_000_007)) });
148impl DynModuloU64 {
149 pub fn set_mod(m: u64) {
150 DYN_MODULUS_U64
151 .with(|cell| unsafe { *cell.get() = BarrettReduction::<u128>::new(m as u128) })
152 }crates/competitive/src/math/primitive_root.rs (line 9)
3pub fn primitive_root(p: u64) -> u64 {
4 if p == 2 {
5 return 1;
6 }
7 let phi = p - 1;
8 let pf = prime_factors(phi);
9 let br = BarrettReduction::<u128>::new(p as _);
10 for g in 2..=3.min(p - 1) {
11 if check_primitive_root(g, phi, &br, &pf) {
12 return g;
13 }
14 }
15 let mut rng = Xorshift::default();
16 loop {
17 let g = ((rng.rand64() as u128 * (p - 2) as u128) >> 64) as u64 + 2;
18 if check_primitive_root(g, phi, &br, &pf) {
19 return g;
20 }
21 }
22}crates/competitive/src/math/floor_sum.rs (line 22)
18pub fn floor_sum(n: u64, a: u64, b: u64, m: u64) -> u64 {
19 let mut ans = Wrapping(0u64);
20 let (mut n, mut m, mut a, mut b) = (Wrapping(n), m, a, b);
21 loop {
22 let br = BarrettReduction::<u64>::new(m);
23 if a >= m {
24 let (q, r) = br.div_rem(a);
25 ans += choose2(n) * q;
26 a = r;
27 }
28 if b >= m {
29 let (q, r) = br.div_rem(b);
30 ans += n * q;
31 b = r;
32 }
33 let y_max = (n * a + b).0;
34 if y_max < m {
35 break;
36 }
37 let (q, r) = br.div_rem(y_max);
38 n = Wrapping(q);
39 b = r;
40 swap(&mut m, &mut a);
41 }
42 ans.0
43}crates/competitive/src/math/discrete_logarithm.rs (line 76)
67 fn discrete_logarithm(&mut self, a: u64, b: u64, p: u64) -> Option<(u64, u64)> {
68 let lim = ((((p as f64).log2() * (p as f64).log2().log2()).sqrt() / 2.0 + 1.).exp2() * 0.9)
69 as u32;
70 self.primes.reserve(lim);
71 let prime_count = self.primes.primes_lte(lim).count();
72 self.br_primes.extend(
73 self.primes
74 .primes_lte(lim)
75 .skip(self.br_primes.len())
76 .map(|p| BarrettReduction::<u64>::new(p.into())),
77 );
78 let br_primes = &self.br_primes[..prime_count];
79 self.ic
80 .entry(p)
81 .or_insert_with(|| IndexCalculusWithPrimitiveRoot::new(p, br_primes))
82 .discrete_logarithm(a, b, br_primes)
83 }
84}
85
86const A: [u32; 150] = [
87 62, 61, 60, 60, 59, 58, 58, 58, 57, 56, 56, 56, 56, 55, 55, 55, 54, 54, 54, 53, 53, 53, 53, 52,
88 52, 52, 52, 52, 52, 51, 50, 50, 50, 50, 49, 49, 49, 48, 48, 48, 48, 48, 47, 47, 47, 47, 47, 47,
89 47, 47, 47, 47, 47, 47, 47, 47, 45, 42, 42, 41, 41, 41, 41, 41, 41, 41, 40, 40, 40, 40, 40, 40,
90 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 38, 38, 38, 38, 38, 32, 32, 32, 32,
91 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 31, 31, 31, 31, 31,
92 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 22, 22, 22, 22,
93 22, 22, 22, 22, 22, 22,
94];
95
96fn factorize_smooth(mut x: u64, row: &mut [u64], br_primes: &[BarrettReduction<u64>]) -> bool {
97 for (j, (&br, r)) in br_primes.iter().zip(row).enumerate() {
98 *r = 0;
99 loop {
100 let (div, rem) = br.div_rem(x);
101 if rem != 0 {
102 break;
103 }
104 *r += 1;
105 x = div;
106 }
107 if j < 150 && x >= (1u64 << A[j]) {
108 break;
109 }
110 }
111 x == 1
112}
113
114#[derive(Debug)]
115struct QdrtPowPrec {
116 br_qdrt: BarrettReduction<u64>,
117 p0: Vec<u64>,
118 p1: Vec<u64>,
119 p2: Vec<u64>,
120 p3: Vec<u64>,
121}
122
123impl QdrtPowPrec {
124 fn new(a: u64, ord: u64, br: &BarrettReduction<u128>) -> Self {
125 let qdrt = (ord as f64).powf(0.25).ceil() as u64;
126 let br_qdrt = BarrettReduction::<u64>::new(qdrt);
127 let mut p0 = Vec::with_capacity(qdrt as usize);
128 let mut p1 = Vec::with_capacity(qdrt as usize);
129 let mut p2 = Vec::with_capacity(qdrt as usize);
130 let mut p3 = Vec::with_capacity(qdrt as usize);
131 let mut acc = 1u64;
132 for _ in 0..qdrt {
133 p0.push(acc);
134 acc = br.rem(acc as u128 * a as u128) as u64;
135 }
136 let a = acc;
137 acc = 1;
138 for _ in 0..qdrt {
139 p1.push(acc);
140 acc = br.rem(acc as u128 * a as u128) as u64;
141 }
142 let a = acc;
143 acc = 1;
144 for _ in 0..qdrt {
145 p2.push(acc);
146 acc = br.rem(acc as u128 * a as u128) as u64;
147 }
148 let a = acc;
149 acc = 1;
150 for _ in 0..qdrt {
151 p3.push(acc);
152 acc = br.rem(acc as u128 * a as u128) as u64;
153 }
154 Self {
155 br_qdrt,
156 p0,
157 p1,
158 p2,
159 p3,
160 }
161 }
162 fn pow(&self, mut k: u64, br: &BarrettReduction<u128>) -> u64 {
163 let (a, b) = self.br_qdrt.div_rem(k);
164 let mut x = self.p0[b as usize];
165 k = a;
166 if k > 0 {
167 let (a, b) = self.br_qdrt.div_rem(k);
168 x = br.rem(x as u128 * self.p1[b as usize] as u128) as u64;
169 k = a;
170 }
171 if k > 0 {
172 let (a, b) = self.br_qdrt.div_rem(k);
173 x = br.rem(x as u128 * self.p2[b as usize] as u128) as u64;
174 k = a;
175 }
176 if k > 0 {
177 let (_, b) = self.br_qdrt.div_rem(k);
178 x = br.rem(x as u128 * self.p3[b as usize] as u128) as u64;
179 }
180 x
181 }
182}
183
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 }
315 fn discrete_logarithm(
316 &self,
317 a: u64,
318 b: u64,
319 br_primes: &[BarrettReduction<u64>],
320 ) -> Option<(u64, u64)> {
321 let p = self.p;
322 let ord = self.ord;
323 let br = BarrettReduction::<u128>::new(p as u128);
324 let a = br.rem(a as _) as u64;
325 let b = br.rem(b as _) as u64;
326 if a == 0 {
327 return if b == 0 { Some((1, 1)) } else { None };
328 }
329 if b == 0 {
330 return None;
331 }
332
333 let x = self.index_calculus(a, br_primes)?;
334 let y = self.index_calculus(b, br_primes)?;
335 solve_linear_congruence(x, y, ord)
336 }
337}
338
339thread_local!(
340 static IC: UnsafeCell<IndexCalculus> = UnsafeCell::new(IndexCalculus::new());
341);
342
343pub fn discrete_logarithm_prime_mod(a: u64, b: u64, p: u64) -> Option<u64> {
344 IC.with(|ic| unsafe { &mut *ic.get() }.discrete_logarithm(a, b, p))
345 .map(|t| t.0)
346}
347
348/// a^x ≡ b (mod n), a has order p^e
349fn pohlig_hellman_prime_power_order(a: u64, b: u64, n: u64, p: u64, e: u32) -> Option<u64> {
350 let br = BarrettReduction::<u128>::new(n as u128);
351 let mul = |x: u64, y: u64| br.rem(x as u128 * y as u128) as u64;
352 let block_size = (p as f64).sqrt().ceil() as u64;
353 let mut baby = HashMap::<u64, u64>::new();
354 let g = pow(a, p.pow(e - 1), &br);
355 let mut xj = 1;
356 for j in 0..block_size {
357 baby.entry(xj).or_insert(j);
358 xj = mul(xj, g);
359 }
360 let xi = modinv(xj, n);
361 let mut t = 0u64;
362 for k in 0..e {
363 let mut h = pow(mul(modinv(pow(a, t, &br), n), b), p.pow(e - 1 - k), &br);
364 let mut ok = false;
365 for i in (0..block_size * block_size).step_by(block_size as usize) {
366 if let Some(j) = baby.get(&h) {
367 t += (i + j) * p.pow(k);
368 ok = true;
369 break;
370 }
371 h = mul(h, xi);
372 }
373 if !ok {
374 return None;
375 }
376 }
377 Some(t)
378}
379
380/// a^x ≡ b (mod p^e)
381fn discrete_logarithm_prime_power(a: u64, b: u64, p: u64, e: u32) -> Option<(u64, u64)> {
382 assert_ne!(p, 0);
383 assert_ne!(e, 0);
384 let n = p.pow(e);
385 assert!(a < n);
386 assert!(b < n);
387 assert_eq!(gcd(a, p), 1);
388 if p == 1 {
389 return Some((0, 1));
390 }
391 if a == 0 {
392 return if b == 0 { Some((1, 1)) } else { None };
393 }
394 if b == 0 {
395 return None;
396 }
397 if e == 1 {
398 return IC.with(|ic| unsafe { &mut *ic.get() }.discrete_logarithm(a, b, p));
399 }
400 let br = BarrettReduction::<u128>::new(n as _);
401 if p == 2 {
402 if e >= 3 {
403 if a % 4 == 1 && b % 4 != 1 {
404 return None;
405 }
406 let aa = if a % 4 == 1 { a } else { n - a };
407 let bb = if b % 4 == 1 { b } else { n - b };
408 let g = 5;
409 let ord = n / 4;
410 let x = pohlig_hellman_prime_power_order(g, aa, n, p, e - 2)?;
411 let y = pohlig_hellman_prime_power_order(g, bb, n, p, e - 2)?;
412 let t = solve_linear_congruence(x, y, ord)?;
413 match (a % 4 == 1, b % 4 == 1) {
414 (true, true) => Some(t),
415 (false, true) if t.0 % 2 == 0 => Some((t.0, lcm(t.1, 2))),
416 (false, false) if t.0 % 2 == 1 => Some((t.0, lcm(t.1, 2))),
417 (false, false) if a == b => Some((1, lcm(t.1, 2))),
418 _ => None,
419 }
420 } else if a == 1 {
421 if b == 1 { Some((0, 1)) } else { None }
422 } else {
423 assert_eq!(a, 3);
424 if b == 1 {
425 Some((0, 2))
426 } else if b == 3 {
427 Some((1, 2))
428 } else {
429 None
430 }
431 }
432 } else {
433 let ord = n - n / p;
434 let pf_ord = prime_factors(ord);
435 let g = (2..)
436 .find(|&g| check_primitive_root(g, ord, &br, &pf_ord))
437 .unwrap();
438 let mut pf_p = prime_factors(p - 1);
439 pf_p.push((p, e - 1));
440 let mut abm = vec![];
441 for (q, c) in pf_p {
442 let m = q.pow(c);
443 let d = ord / m;
444 let gg = pow(g, d, &br);
445 let aa = pow(a, d, &br);
446 let bb = pow(b, d, &br);
447 let x = pohlig_hellman_prime_power_order(gg, aa, n, q, c)?;
448 let y = pohlig_hellman_prime_power_order(gg, bb, n, q, c)?;
449 abm.push((x, y, m));
450 }
451 solve_linear_congruences(abm)
452 }
453}crates/competitive/src/math/prime_factors.rs (line 33)
32fn find_factor(n: u64) -> Option<u64> {
33 let br = BarrettReduction::<u128>::new(n as u128);
34 if miller_rabin_with_br(n, &br) {
35 return None;
36 }
37 let mr = MontgomeryReduction64::new(n);
38 let mut rng = Xorshift::default();
39 let (mut y0, mut c) = (0, n - 1);
40 loop {
41 let (mut x, mut y, mut ys, mut g, mut q, mut r, mut k) = (0, y0, 0, 1, 1, 1, 0);
42 while g == 1 && r <= 1 << 20 {
43 x = y;
44 while k < r && g == 1 {
45 ys = y;
46 for _ in 0..1024.min(r - k) {
47 y = mr.sub(mr.mul(y, y), c);
48 q = mr.mul(q, mr.sub(x, y));
49 }
50 g = gcd(q, n);
51 k += 1024;
52 }
53 k = r;
54 r <<= 1;
55 }
56 if g == n {
57 g = 1;
58 y = ys;
59 while g == 1 {
60 y = mr.sub(mr.mul(y, y), c);
61 g = gcd(mr.sub(x, y), n);
62 }
63 }
64 if g != 1 && g != n {
65 return Some(g);
66 }
67 y0 = ((rng.rand64() as u128 * (n - 2) as u128) >> 64) as u64 + 2;
68 c = ((rng.rand64() as u128 * (n - 1) as u128) >> 64) as u64 + 1;
69 }
70}Additional examples can be found in:
pub const fn new_with_im(m: T, im: T) -> Self
pub const fn get_mod(&self) -> T
Sourcepub fn div_rem(&self, a: T) -> (T, T)
pub fn div_rem(&self, a: T) -> (T, T)
Examples found in repository?
More examples
crates/competitive/src/math/discrete_logarithm.rs (line 100)
96fn factorize_smooth(mut x: u64, row: &mut [u64], br_primes: &[BarrettReduction<u64>]) -> bool {
97 for (j, (&br, r)) in br_primes.iter().zip(row).enumerate() {
98 *r = 0;
99 loop {
100 let (div, rem) = br.div_rem(x);
101 if rem != 0 {
102 break;
103 }
104 *r += 1;
105 x = div;
106 }
107 if j < 150 && x >= (1u64 << A[j]) {
108 break;
109 }
110 }
111 x == 1
112}
113
114#[derive(Debug)]
115struct QdrtPowPrec {
116 br_qdrt: BarrettReduction<u64>,
117 p0: Vec<u64>,
118 p1: Vec<u64>,
119 p2: Vec<u64>,
120 p3: Vec<u64>,
121}
122
123impl QdrtPowPrec {
124 fn new(a: u64, ord: u64, br: &BarrettReduction<u128>) -> Self {
125 let qdrt = (ord as f64).powf(0.25).ceil() as u64;
126 let br_qdrt = BarrettReduction::<u64>::new(qdrt);
127 let mut p0 = Vec::with_capacity(qdrt as usize);
128 let mut p1 = Vec::with_capacity(qdrt as usize);
129 let mut p2 = Vec::with_capacity(qdrt as usize);
130 let mut p3 = Vec::with_capacity(qdrt as usize);
131 let mut acc = 1u64;
132 for _ in 0..qdrt {
133 p0.push(acc);
134 acc = br.rem(acc as u128 * a as u128) as u64;
135 }
136 let a = acc;
137 acc = 1;
138 for _ in 0..qdrt {
139 p1.push(acc);
140 acc = br.rem(acc as u128 * a as u128) as u64;
141 }
142 let a = acc;
143 acc = 1;
144 for _ in 0..qdrt {
145 p2.push(acc);
146 acc = br.rem(acc as u128 * a as u128) as u64;
147 }
148 let a = acc;
149 acc = 1;
150 for _ in 0..qdrt {
151 p3.push(acc);
152 acc = br.rem(acc as u128 * a as u128) as u64;
153 }
154 Self {
155 br_qdrt,
156 p0,
157 p1,
158 p2,
159 p3,
160 }
161 }
162 fn pow(&self, mut k: u64, br: &BarrettReduction<u128>) -> u64 {
163 let (a, b) = self.br_qdrt.div_rem(k);
164 let mut x = self.p0[b as usize];
165 k = a;
166 if k > 0 {
167 let (a, b) = self.br_qdrt.div_rem(k);
168 x = br.rem(x as u128 * self.p1[b as usize] as u128) as u64;
169 k = a;
170 }
171 if k > 0 {
172 let (a, b) = self.br_qdrt.div_rem(k);
173 x = br.rem(x as u128 * self.p2[b as usize] as u128) as u64;
174 k = a;
175 }
176 if k > 0 {
177 let (_, b) = self.br_qdrt.div_rem(k);
178 x = br.rem(x as u128 * self.p3[b as usize] as u128) as u64;
179 }
180 x
181 }crates/competitive/src/math/floor_sum.rs (line 24)
18pub fn floor_sum(n: u64, a: u64, b: u64, m: u64) -> u64 {
19 let mut ans = Wrapping(0u64);
20 let (mut n, mut m, mut a, mut b) = (Wrapping(n), m, a, b);
21 loop {
22 let br = BarrettReduction::<u64>::new(m);
23 if a >= m {
24 let (q, r) = br.div_rem(a);
25 ans += choose2(n) * q;
26 a = r;
27 }
28 if b >= m {
29 let (q, r) = br.div_rem(b);
30 ans += n * q;
31 b = r;
32 }
33 let y_max = (n * a + b).0;
34 if y_max < m {
35 break;
36 }
37 let (q, r) = br.div_rem(y_max);
38 n = Wrapping(q);
39 b = r;
40 swap(&mut m, &mut a);
41 }
42 ans.0
43}crates/competitive/src/math/arbitrary_mod_binomial.rs (line 111)
102 fn combination(&self, mut n: u64, mut k: u64) -> u64 {
103 if k > n {
104 return 0;
105 }
106 assert!(self.size as u64 >= n.min(self.m - 1));
107 if self.m < 1 << 31 {
108 let mut res = 1u64;
109 if self.e == 1 {
110 while n > 0 {
111 let (nn, n0) = self.bp.div_rem(n);
112 let (nk, k0) = self.bp.div_rem(k);
113 if n0 < k0 {
114 return 0;
115 }
116 res = self.bm.rem(res * self.fact[n0 as usize]);
117 res = self.bm.rem(res * self.inv_fact[k0 as usize]);
118 res = self.bm.rem(res * self.inv_fact[(n0 - k0) as usize]);
119 n = nn;
120 k = nk;
121 }
122 } else {
123 let mut r = n - k;
124 let mut e0 = 0;
125 let mut eq = 0;
126 let mut i = 0;
127 while n > 0 {
128 res = self.bm.rem(res * self.fact[self.bm.rem(n) as usize]);
129 res = self.bm.rem(res * self.inv_fact[self.bm.rem(k) as usize]);
130 res = self.bm.rem(res * self.inv_fact[self.bm.rem(r) as usize]);
131 n = self.bp.div(n);
132 k = self.bp.div(k);
133 r = self.bp.div(r);
134 let eps = n - k - r;
135 e0 += eps;
136 if e0 >= self.e as u64 {
137 return 0;
138 }
139 i += 1;
140 if i >= self.e {
141 eq ^= eps & 1;
142 }
143 }
144 if eq == 1 {
145 res = self.bm.rem(res * self.delta);
146 }
147 res = self
148 .bm
149 .rem(res * pow32(self.p as _, e0 as _, &self.bm) as u64);
150 }
151 res
152 } else {
153 let mut res = 1u128;
154 if self.e == 1 {
155 while n > 0 {
156 let (nn, n0) = self.bp.div_rem(n);
157 let (nk, k0) = self.bp.div_rem(k);
158 if n0 < k0 {
159 return 0;
160 }
161 res = self.bm128.rem(res * self.fact[n0 as usize] as u128);
162 res = self.bm128.rem(res * self.inv_fact[k0 as usize] as u128);
163 res = self
164 .bm128
165 .rem(res * self.inv_fact[(n0 - k0) as usize] as u128);
166 n = nn;
167 k = nk;
168 }
169 } else {
170 let mut r = n - k;
171 let mut e0 = 0;
172 let mut eq = 0;
173 let mut i = 0;
174 while n > 0 {
175 res = self
176 .bm128
177 .rem(res * self.fact[self.bm.rem(n) as usize] as u128);
178 res = self
179 .bm128
180 .rem(res * self.inv_fact[self.bm.rem(k) as usize] as u128);
181 res = self
182 .bm128
183 .rem(res * self.inv_fact[self.bm.rem(r) as usize] as u128);
184 n = self.bp.div(n);
185 k = self.bp.div(k);
186 r = self.bp.div(r);
187 let eps = n - k - r;
188 e0 += eps;
189 if e0 >= self.e as u64 {
190 return 0;
191 }
192 i += 1;
193 if i >= self.e {
194 eq ^= eps & 1;
195 }
196 }
197 if eq == 1 {
198 res = self.bm128.rem(res * self.delta as u128);
199 }
200 res = self.bm128.rem(res * pow64(self.p, e0, &self.bm128) as u128);
201 }
202 res as u64
203 }
204 }Sourcepub fn div(&self, a: T) -> T
pub fn div(&self, a: T) -> T
Examples found in repository?
crates/competitive/src/math/arbitrary_mod_binomial.rs (line 131)
102 fn combination(&self, mut n: u64, mut k: u64) -> u64 {
103 if k > n {
104 return 0;
105 }
106 assert!(self.size as u64 >= n.min(self.m - 1));
107 if self.m < 1 << 31 {
108 let mut res = 1u64;
109 if self.e == 1 {
110 while n > 0 {
111 let (nn, n0) = self.bp.div_rem(n);
112 let (nk, k0) = self.bp.div_rem(k);
113 if n0 < k0 {
114 return 0;
115 }
116 res = self.bm.rem(res * self.fact[n0 as usize]);
117 res = self.bm.rem(res * self.inv_fact[k0 as usize]);
118 res = self.bm.rem(res * self.inv_fact[(n0 - k0) as usize]);
119 n = nn;
120 k = nk;
121 }
122 } else {
123 let mut r = n - k;
124 let mut e0 = 0;
125 let mut eq = 0;
126 let mut i = 0;
127 while n > 0 {
128 res = self.bm.rem(res * self.fact[self.bm.rem(n) as usize]);
129 res = self.bm.rem(res * self.inv_fact[self.bm.rem(k) as usize]);
130 res = self.bm.rem(res * self.inv_fact[self.bm.rem(r) as usize]);
131 n = self.bp.div(n);
132 k = self.bp.div(k);
133 r = self.bp.div(r);
134 let eps = n - k - r;
135 e0 += eps;
136 if e0 >= self.e as u64 {
137 return 0;
138 }
139 i += 1;
140 if i >= self.e {
141 eq ^= eps & 1;
142 }
143 }
144 if eq == 1 {
145 res = self.bm.rem(res * self.delta);
146 }
147 res = self
148 .bm
149 .rem(res * pow32(self.p as _, e0 as _, &self.bm) as u64);
150 }
151 res
152 } else {
153 let mut res = 1u128;
154 if self.e == 1 {
155 while n > 0 {
156 let (nn, n0) = self.bp.div_rem(n);
157 let (nk, k0) = self.bp.div_rem(k);
158 if n0 < k0 {
159 return 0;
160 }
161 res = self.bm128.rem(res * self.fact[n0 as usize] as u128);
162 res = self.bm128.rem(res * self.inv_fact[k0 as usize] as u128);
163 res = self
164 .bm128
165 .rem(res * self.inv_fact[(n0 - k0) as usize] as u128);
166 n = nn;
167 k = nk;
168 }
169 } else {
170 let mut r = n - k;
171 let mut e0 = 0;
172 let mut eq = 0;
173 let mut i = 0;
174 while n > 0 {
175 res = self
176 .bm128
177 .rem(res * self.fact[self.bm.rem(n) as usize] as u128);
178 res = self
179 .bm128
180 .rem(res * self.inv_fact[self.bm.rem(k) as usize] as u128);
181 res = self
182 .bm128
183 .rem(res * self.inv_fact[self.bm.rem(r) as usize] as u128);
184 n = self.bp.div(n);
185 k = self.bp.div(k);
186 r = self.bp.div(r);
187 let eps = n - k - r;
188 e0 += eps;
189 if e0 >= self.e as u64 {
190 return 0;
191 }
192 i += 1;
193 if i >= self.e {
194 eq ^= eps & 1;
195 }
196 }
197 if eq == 1 {
198 res = self.bm128.rem(res * self.delta as u128);
199 }
200 res = self.bm128.rem(res * pow64(self.p, e0, &self.bm128) as u128);
201 }
202 res as u64
203 }
204 }Sourcepub fn rem(&self, a: T) -> T
pub fn rem(&self, a: T) -> T
Examples found in repository?
crates/competitive/src/num/mint/mint_basic.rs (line 138)
137 fn rem(x: u64) -> u64 {
138 DYN_MODULUS_U32.with(|cell| unsafe { (*cell.get()).rem(x) })
139 }
140}
141impl DynMIntU32 {
142 pub fn set_mod(m: u32) {
143 DynModuloU32::set_mod(m)
144 }
145}
146
147thread_local!(static DYN_MODULUS_U64: UnsafeCell<BarrettReduction<u128>> = const { UnsafeCell::new(BarrettReduction::<u128>::new_with_im(1_000_000_007, !0 / 1_000_000_007)) });
148impl DynModuloU64 {
149 pub fn set_mod(m: u64) {
150 DYN_MODULUS_U64
151 .with(|cell| unsafe { *cell.get() = BarrettReduction::<u128>::new(m as u128) })
152 }
153 fn rem(x: u128) -> u128 {
154 DYN_MODULUS_U64.with(|cell| unsafe { (*cell.get()).rem(x) })
155 }More examples
crates/competitive/src/math/discrete_logarithm.rs (line 12)
7fn pow(x: u64, mut y: u64, br: &BarrettReduction<u128>) -> u64 {
8 let mut x = x as u128;
9 let mut z: u128 = 1;
10 while y > 0 {
11 if y & 1 == 1 {
12 z = br.rem(z * x);
13 }
14 x = br.rem(x * x);
15 y >>= 1;
16 }
17 z as u64
18}
19
20fn solve_linear_congruence(a: u64, b: u64, m: u64) -> Option<(u64, u64)> {
21 let g = gcd(a, m);
22 if !b.is_multiple_of(g) {
23 return None;
24 }
25 let (a, b, m) = (a / g, b / g, m / g);
26 Some(((b as u128 * modinv(a, m) as u128 % m as u128) as _, m))
27}
28
29fn solve_linear_congruences<I>(abm: I) -> Option<(u64, u64)>
30where
31 I: IntoIterator<Item = (u64, u64, u64)>,
32{
33 let mut x = 0u64;
34 let mut m0 = 1u64;
35 for (a, b, m) in abm {
36 let mut b = b + m - a * x % m;
37 if b >= m {
38 b -= m;
39 }
40 let a = a * m0;
41 let g = gcd(a, m);
42 if !b.is_multiple_of(g) {
43 return None;
44 }
45 let (a, b, m) = (a / g, b / g, m / g);
46 x += (b as u128 * modinv(a, m) as u128 % m as u128 * m0 as u128) as u64;
47 m0 *= m;
48 }
49 Some((x, m0))
50}
51
52#[derive(Debug)]
53struct IndexCalculus {
54 primes: PrimeList,
55 br_primes: Vec<BarrettReduction<u64>>,
56 ic: HashMap<u64, IndexCalculusWithPrimitiveRoot>,
57}
58
59impl IndexCalculus {
60 fn new() -> Self {
61 Self {
62 primes: PrimeList::new(2),
63 br_primes: Default::default(),
64 ic: Default::default(),
65 }
66 }
67 fn discrete_logarithm(&mut self, a: u64, b: u64, p: u64) -> Option<(u64, u64)> {
68 let lim = ((((p as f64).log2() * (p as f64).log2().log2()).sqrt() / 2.0 + 1.).exp2() * 0.9)
69 as u32;
70 self.primes.reserve(lim);
71 let prime_count = self.primes.primes_lte(lim).count();
72 self.br_primes.extend(
73 self.primes
74 .primes_lte(lim)
75 .skip(self.br_primes.len())
76 .map(|p| BarrettReduction::<u64>::new(p.into())),
77 );
78 let br_primes = &self.br_primes[..prime_count];
79 self.ic
80 .entry(p)
81 .or_insert_with(|| IndexCalculusWithPrimitiveRoot::new(p, br_primes))
82 .discrete_logarithm(a, b, br_primes)
83 }
84}
85
86const A: [u32; 150] = [
87 62, 61, 60, 60, 59, 58, 58, 58, 57, 56, 56, 56, 56, 55, 55, 55, 54, 54, 54, 53, 53, 53, 53, 52,
88 52, 52, 52, 52, 52, 51, 50, 50, 50, 50, 49, 49, 49, 48, 48, 48, 48, 48, 47, 47, 47, 47, 47, 47,
89 47, 47, 47, 47, 47, 47, 47, 47, 45, 42, 42, 41, 41, 41, 41, 41, 41, 41, 40, 40, 40, 40, 40, 40,
90 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 38, 38, 38, 38, 38, 32, 32, 32, 32,
91 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 31, 31, 31, 31, 31,
92 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 22, 22, 22, 22,
93 22, 22, 22, 22, 22, 22,
94];
95
96fn factorize_smooth(mut x: u64, row: &mut [u64], br_primes: &[BarrettReduction<u64>]) -> bool {
97 for (j, (&br, r)) in br_primes.iter().zip(row).enumerate() {
98 *r = 0;
99 loop {
100 let (div, rem) = br.div_rem(x);
101 if rem != 0 {
102 break;
103 }
104 *r += 1;
105 x = div;
106 }
107 if j < 150 && x >= (1u64 << A[j]) {
108 break;
109 }
110 }
111 x == 1
112}
113
114#[derive(Debug)]
115struct QdrtPowPrec {
116 br_qdrt: BarrettReduction<u64>,
117 p0: Vec<u64>,
118 p1: Vec<u64>,
119 p2: Vec<u64>,
120 p3: Vec<u64>,
121}
122
123impl QdrtPowPrec {
124 fn new(a: u64, ord: u64, br: &BarrettReduction<u128>) -> Self {
125 let qdrt = (ord as f64).powf(0.25).ceil() as u64;
126 let br_qdrt = BarrettReduction::<u64>::new(qdrt);
127 let mut p0 = Vec::with_capacity(qdrt as usize);
128 let mut p1 = Vec::with_capacity(qdrt as usize);
129 let mut p2 = Vec::with_capacity(qdrt as usize);
130 let mut p3 = Vec::with_capacity(qdrt as usize);
131 let mut acc = 1u64;
132 for _ in 0..qdrt {
133 p0.push(acc);
134 acc = br.rem(acc as u128 * a as u128) as u64;
135 }
136 let a = acc;
137 acc = 1;
138 for _ in 0..qdrt {
139 p1.push(acc);
140 acc = br.rem(acc as u128 * a as u128) as u64;
141 }
142 let a = acc;
143 acc = 1;
144 for _ in 0..qdrt {
145 p2.push(acc);
146 acc = br.rem(acc as u128 * a as u128) as u64;
147 }
148 let a = acc;
149 acc = 1;
150 for _ in 0..qdrt {
151 p3.push(acc);
152 acc = br.rem(acc as u128 * a as u128) as u64;
153 }
154 Self {
155 br_qdrt,
156 p0,
157 p1,
158 p2,
159 p3,
160 }
161 }
162 fn pow(&self, mut k: u64, br: &BarrettReduction<u128>) -> u64 {
163 let (a, b) = self.br_qdrt.div_rem(k);
164 let mut x = self.p0[b as usize];
165 k = a;
166 if k > 0 {
167 let (a, b) = self.br_qdrt.div_rem(k);
168 x = br.rem(x as u128 * self.p1[b as usize] as u128) as u64;
169 k = a;
170 }
171 if k > 0 {
172 let (a, b) = self.br_qdrt.div_rem(k);
173 x = br.rem(x as u128 * self.p2[b as usize] as u128) as u64;
174 k = a;
175 }
176 if k > 0 {
177 let (_, b) = self.br_qdrt.div_rem(k);
178 x = br.rem(x as u128 * self.p3[b as usize] as u128) as u64;
179 }
180 x
181 }
182}
183
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 }
315 fn discrete_logarithm(
316 &self,
317 a: u64,
318 b: u64,
319 br_primes: &[BarrettReduction<u64>],
320 ) -> Option<(u64, u64)> {
321 let p = self.p;
322 let ord = self.ord;
323 let br = BarrettReduction::<u128>::new(p as u128);
324 let a = br.rem(a as _) as u64;
325 let b = br.rem(b as _) as u64;
326 if a == 0 {
327 return if b == 0 { Some((1, 1)) } else { None };
328 }
329 if b == 0 {
330 return None;
331 }
332
333 let x = self.index_calculus(a, br_primes)?;
334 let y = self.index_calculus(b, br_primes)?;
335 solve_linear_congruence(x, y, ord)
336 }
337}
338
339thread_local!(
340 static IC: UnsafeCell<IndexCalculus> = UnsafeCell::new(IndexCalculus::new());
341);
342
343pub fn discrete_logarithm_prime_mod(a: u64, b: u64, p: u64) -> Option<u64> {
344 IC.with(|ic| unsafe { &mut *ic.get() }.discrete_logarithm(a, b, p))
345 .map(|t| t.0)
346}
347
348/// a^x ≡ b (mod n), a has order p^e
349fn pohlig_hellman_prime_power_order(a: u64, b: u64, n: u64, p: u64, e: u32) -> Option<u64> {
350 let br = BarrettReduction::<u128>::new(n as u128);
351 let mul = |x: u64, y: u64| br.rem(x as u128 * y as u128) as u64;
352 let block_size = (p as f64).sqrt().ceil() as u64;
353 let mut baby = HashMap::<u64, u64>::new();
354 let g = pow(a, p.pow(e - 1), &br);
355 let mut xj = 1;
356 for j in 0..block_size {
357 baby.entry(xj).or_insert(j);
358 xj = mul(xj, g);
359 }
360 let xi = modinv(xj, n);
361 let mut t = 0u64;
362 for k in 0..e {
363 let mut h = pow(mul(modinv(pow(a, t, &br), n), b), p.pow(e - 1 - k), &br);
364 let mut ok = false;
365 for i in (0..block_size * block_size).step_by(block_size as usize) {
366 if let Some(j) = baby.get(&h) {
367 t += (i + j) * p.pow(k);
368 ok = true;
369 break;
370 }
371 h = mul(h, xi);
372 }
373 if !ok {
374 return None;
375 }
376 }
377 Some(t)
378}crates/competitive/src/math/arbitrary_mod_binomial.rs (line 8)
3fn pow64(x: u64, mut y: u64, br: &BarrettReduction<u128>) -> u64 {
4 let mut x = x as u128;
5 let mut z: u128 = 1;
6 while y > 0 {
7 if y & 1 == 1 {
8 z = br.rem(z * x);
9 }
10 x = br.rem(x * x);
11 y >>= 1;
12 }
13 z as u64
14}
15
16fn pow32(x: u32, mut y: u64, br: &BarrettReduction<u64>) -> u32 {
17 let mut x = x as u64;
18 let mut z: u64 = 1;
19 while y > 0 {
20 if y & 1 == 1 {
21 z = br.rem(z * x);
22 }
23 x = br.rem(x * x);
24 y >>= 1;
25 }
26 z as u32
27}
28
29#[derive(Debug)]
30struct PrimePowerBinomial {
31 p: u64,
32 e: u32,
33 m: u64,
34 size: usize,
35 fact: Vec<u64>,
36 inv_fact: Vec<u64>,
37 delta: u64,
38 bp: BarrettReduction<u64>,
39 bm: BarrettReduction<u64>,
40 bm128: BarrettReduction<u128>,
41}
42
43impl PrimePowerBinomial {
44 fn new(p: u64, e: u32, max_n: u64) -> Self {
45 let m = p.checked_pow(e).expect("prime power overflow");
46 let bp = BarrettReduction::new(p);
47 let bm = BarrettReduction::new(m);
48 let bm128 = BarrettReduction::new(m as u128);
49 let size = max_n.min(m - 1);
50 assert!(size < usize::MAX as u64);
51 let size = size as usize;
52 let mut fact = vec![1u64; size + 1];
53 let mut inv_fact = vec![1u64; size + 1];
54 if m < 1 << 31 {
55 for i in 2..=size {
56 fact[i] = if bp.rem(i as u64) == 0 {
57 fact[i - 1]
58 } else {
59 bm.rem(fact[i - 1] * i as u64)
60 };
61 }
62 inv_fact[size] = fact[size].mod_inv(m);
63 for i in (3..=size).rev() {
64 inv_fact[i - 1] = if bp.rem(i as u64) == 0 {
65 inv_fact[i]
66 } else {
67 bm.rem(inv_fact[i] * i as u64)
68 };
69 }
70 } else {
71 for i in 2..=size {
72 fact[i] = if bp.rem(i as u64) == 0 {
73 fact[i - 1]
74 } else {
75 bm128.rem(fact[i - 1] as u128 * i as u128) as u64
76 };
77 }
78 inv_fact[size] = fact[size].mod_inv(m);
79 for i in (3..=size).rev() {
80 inv_fact[i - 1] = if bp.rem(i as u64) == 0 {
81 inv_fact[i]
82 } else {
83 bm128.rem(inv_fact[i] as u128 * i as u128) as u64
84 };
85 }
86 }
87 let delta = if p == 2 && e >= 3 { 1 } else { m - 1 };
88 Self {
89 p,
90 e,
91 m,
92 size,
93 fact,
94 inv_fact,
95 delta,
96 bp,
97 bm,
98 bm128,
99 }
100 }
101
102 fn combination(&self, mut n: u64, mut k: u64) -> u64 {
103 if k > n {
104 return 0;
105 }
106 assert!(self.size as u64 >= n.min(self.m - 1));
107 if self.m < 1 << 31 {
108 let mut res = 1u64;
109 if self.e == 1 {
110 while n > 0 {
111 let (nn, n0) = self.bp.div_rem(n);
112 let (nk, k0) = self.bp.div_rem(k);
113 if n0 < k0 {
114 return 0;
115 }
116 res = self.bm.rem(res * self.fact[n0 as usize]);
117 res = self.bm.rem(res * self.inv_fact[k0 as usize]);
118 res = self.bm.rem(res * self.inv_fact[(n0 - k0) as usize]);
119 n = nn;
120 k = nk;
121 }
122 } else {
123 let mut r = n - k;
124 let mut e0 = 0;
125 let mut eq = 0;
126 let mut i = 0;
127 while n > 0 {
128 res = self.bm.rem(res * self.fact[self.bm.rem(n) as usize]);
129 res = self.bm.rem(res * self.inv_fact[self.bm.rem(k) as usize]);
130 res = self.bm.rem(res * self.inv_fact[self.bm.rem(r) as usize]);
131 n = self.bp.div(n);
132 k = self.bp.div(k);
133 r = self.bp.div(r);
134 let eps = n - k - r;
135 e0 += eps;
136 if e0 >= self.e as u64 {
137 return 0;
138 }
139 i += 1;
140 if i >= self.e {
141 eq ^= eps & 1;
142 }
143 }
144 if eq == 1 {
145 res = self.bm.rem(res * self.delta);
146 }
147 res = self
148 .bm
149 .rem(res * pow32(self.p as _, e0 as _, &self.bm) as u64);
150 }
151 res
152 } else {
153 let mut res = 1u128;
154 if self.e == 1 {
155 while n > 0 {
156 let (nn, n0) = self.bp.div_rem(n);
157 let (nk, k0) = self.bp.div_rem(k);
158 if n0 < k0 {
159 return 0;
160 }
161 res = self.bm128.rem(res * self.fact[n0 as usize] as u128);
162 res = self.bm128.rem(res * self.inv_fact[k0 as usize] as u128);
163 res = self
164 .bm128
165 .rem(res * self.inv_fact[(n0 - k0) as usize] as u128);
166 n = nn;
167 k = nk;
168 }
169 } else {
170 let mut r = n - k;
171 let mut e0 = 0;
172 let mut eq = 0;
173 let mut i = 0;
174 while n > 0 {
175 res = self
176 .bm128
177 .rem(res * self.fact[self.bm.rem(n) as usize] as u128);
178 res = self
179 .bm128
180 .rem(res * self.inv_fact[self.bm.rem(k) as usize] as u128);
181 res = self
182 .bm128
183 .rem(res * self.inv_fact[self.bm.rem(r) as usize] as u128);
184 n = self.bp.div(n);
185 k = self.bp.div(k);
186 r = self.bp.div(r);
187 let eps = n - k - r;
188 e0 += eps;
189 if e0 >= self.e as u64 {
190 return 0;
191 }
192 i += 1;
193 if i >= self.e {
194 eq ^= eps & 1;
195 }
196 }
197 if eq == 1 {
198 res = self.bm128.rem(res * self.delta as u128);
199 }
200 res = self.bm128.rem(res * pow64(self.p, e0, &self.bm128) as u128);
201 }
202 res as u64
203 }
204 }crates/competitive/src/math/primitive_root.rs (line 36)
24pub fn check_primitive_root(
25 g: u64,
26 phi: u64,
27 br: &BarrettReduction<u128>,
28 pf: &[(u64, u32)],
29) -> bool {
30 pf.iter().all(|&(q, _)| {
31 let mut g = g as u128;
32 let mut k = phi / q;
33 let mut r: u128 = 1;
34 while k > 0 {
35 if k & 1 == 1 {
36 r = br.rem(r * g);
37 }
38 g = br.rem(g * g);
39 k >>= 1;
40 }
41 r > 1
42 })
43}Trait Implementations§
Source§impl<T: Clone> Clone for BarrettReduction<T>
impl<T: Clone> Clone for BarrettReduction<T>
impl<T: Copy> Copy for BarrettReduction<T>
Auto Trait Implementations§
impl<T> Freeze for BarrettReduction<T>where
T: Freeze,
impl<T> RefUnwindSafe for BarrettReduction<T>where
T: RefUnwindSafe,
impl<T> Send for BarrettReduction<T>where
T: Send,
impl<T> Sync for BarrettReduction<T>where
T: Sync,
impl<T> Unpin for BarrettReduction<T>where
T: Unpin,
impl<T> UnsafeUnpin for BarrettReduction<T>where
T: UnsafeUnpin,
impl<T> UnwindSafe for BarrettReduction<T>where
T: UnwindSafe,
Blanket Implementations§
Source§impl<T> BorrowMut<T> for Twhere
T: ?Sized,
impl<T> BorrowMut<T> for Twhere
T: ?Sized,
Source§fn borrow_mut(&mut self) -> &mut T
fn borrow_mut(&mut self) -> &mut T
Mutably borrows from an owned value. Read more