Struct PrimeList

Source

pub struct PrimeList { /* private fields */ }

Implementations§

Source §

impl PrimeList

Source

pub fn new(max_n: u64) -> Self

Examples found in repository ?

crates/competitive/src/math/discrete_logarithm.rs (line 62)

60    fn new() -> Self {
61        Self {
62            primes: PrimeList::new(2),
63            br_primes: Default::default(),
64            ic: Default::default(),
65        }
66    }

More examples

Hide additional examples

crates/library_checker/src/math/enumerate_primes.rs (line 10)

6pub fn enumerate_primes(reader: impl Read, mut writer: impl Write) {
7    let s = read_all_unchecked(reader);
8    let mut scanner = Scanner::new(&s);
9    scan!(scanner, n: u64, a, b);
10    let primes = PrimeList::new(n);
11    let iter = primes.primes().iter().skip(b).step_by(a);
12    writeln!(writer, "{} {}", primes.primes().len(), iter.clone().len()).ok();
13    for p in iter {
14        write!(writer, "{} ", p).ok();
15    }
16    writeln!(writer).ok();
17}

Source

pub fn primes(&self) -> &[u64]

Examples found in repository ?

crates/library_checker/src/math/enumerate_primes.rs (line 11)

6pub fn enumerate_primes(reader: impl Read, mut writer: impl Write) {
7    let s = read_all_unchecked(reader);
8    let mut scanner = Scanner::new(&s);
9    scan!(scanner, n: u64, a, b);
10    let primes = PrimeList::new(n);
11    let iter = primes.primes().iter().skip(b).step_by(a);
12    writeln!(writer, "{} {}", primes.primes().len(), iter.clone().len()).ok();
13    for p in iter {
14        write!(writer, "{} ", p).ok();
15    }
16    writeln!(writer).ok();
17}

Source

pub fn primes_lte(&self, n: u64) -> &[u64]

Examples found in repository ?

crates/competitive/src/math/lcm_convolve.rs (line 16)

13    pub fn zeta_transform(f: &mut [M::T]) {
14        let n = f.len().saturating_sub(1) as u64;
15        with_prime_list(n, |pl| {
16            for &p in pl.primes_lte(n).iter() {
17                for (i, j) in (0..f.len()).step_by(p as _).enumerate() {
18                    f[j] = M::operate(&f[j], &f[i]);
19                }
20            }
21        })
22    }
23}
24
25impl<G> LcmConvolve<G>
26where
27    G: Group,
28{
29    /// $$f(m) = \sum_{n \mid m}h(n)$$
30    pub fn mobius_transform(f: &mut [G::T]) {
31        let n = f.len().saturating_sub(1) as u64;
32        with_prime_list(n, |pl| {
33            for &p in pl.primes_lte(n).iter() {
34                for (i, j) in (0..f.len()).step_by(p as _).enumerate().rev() {
35                    f[j] = G::rinv_operate(&f[j], &f[i]);
36                }
37            }
38        })
39    }

More examples

Hide additional examples

crates/competitive/src/math/gcd_convolve.rs (line 16)

13    pub fn zeta_transform(f: &mut [M::T]) {
14        let n = f.len().saturating_sub(1) as u64;
15        with_prime_list(n, |pl| {
16            for &p in pl.primes_lte(n).iter() {
17                for (i, j) in (0..f.len()).step_by(p as _).enumerate().rev() {
18                    f[i] = M::operate(&f[i], &f[j]);
19                }
20            }
21        })
22    }
23}
24
25impl<G> GcdConvolve<G>
26where
27    G: Group,
28{
29    /// $$f(m) = \sum_{n \mid m}h(n)$$
30    pub fn mobius_transform(f: &mut [G::T]) {
31        let n = f.len().saturating_sub(1) as u64;
32        with_prime_list(n, |pl| {
33            for &p in pl.primes_lte(n).iter() {
34                for (i, j) in (0..f.len()).step_by(p as _).enumerate() {
35                    f[i] = G::rinv_operate(&f[i], &f[j]);
36                }
37            }
38        })
39    }

crates/competitive/src/math/quotient_array.rs (line 66)

61    pub fn lucy_dp<G>(mut self, mut mul_p: impl FnMut(T, u64) -> T) -> Self
62    where
63        G: Group<T = T>,
64    {
65        with_prime_list(self.isqrtn, |pl| {
66            for &p in pl.primes_lte(self.isqrtn) {
67                let k = self.quotient_index(p - 1);
68                let p2 = p * p;
69                for (i, q) in Self::index_iter(self.n, self.isqrtn).enumerate() {
70                    if q < p2 {
71                        break;
72                    }
73                    let diff = mul_p(G::rinv_operate(&self[q / p], &self.data[k]), p);
74                    G::rinv_operate_assign(&mut self.data[i], &diff);
75                }
76            }
77        });
78        self
79    }
80
81    /// convert $\sum_{i\leq n, i\text{ is prime}} f(i)$ to $\sum_{i\leq n} f(i)$
82    pub fn min_25_sieve<R>(&self, mut f: impl FnMut(u64, u32) -> T) -> Self
83    where
84        T: Clone + One,
85        R: Ring<T = T>,
86        R::Additive: Invertible,
87    {
88        let mut dp = self.clone();
89        with_prime_list(self.isqrtn, |pl| {
90            for &p in pl.primes_lte(self.isqrtn).iter().rev() {
91                let k = self.quotient_index(p);
92                for (i, q) in Self::index_iter(self.n, self.isqrtn).enumerate() {
93                    let mut pc = p;
94                    if pc * p > q {
95                        break;
96                    }
97                    let mut c = 1;
98                    while q / p >= pc {
99                        let x = R::mul(&f(p, c), &(R::sub(&dp[q / pc], &self.data[k])));
100                        let x = R::add(&x, &f(p, c + 1));
101                        dp.data[i] = R::add(&dp.data[i], &x);
102                        c += 1;
103                        pc *= p;
104                    }
105                }
106            }
107        });
108        for x in &mut dp.data {
109            *x = R::add(x, &T::one());
110        }
111        dp
112    }

crates/competitive/src/math/discrete_logarithm.rs (line 71)

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 u64;
70        self.primes.reserve(lim);
71        let primes = self.primes.primes_lte(lim);
72        while self.br_primes.len() < primes.len() {
73            let br = BarrettReduction::<u64>::new(primes[self.br_primes.len()]);
74            self.br_primes.push(br);
75        }
76        let br_primes = &self.br_primes[..primes.len()];
77        self.ic
78            .entry(p)
79            .or_insert_with(|| IndexCalculusWithPrimitiveRoot::new(p, br_primes))
80            .discrete_logarithm(a, b, br_primes)
81    }

Source

pub fn is_prime(&self, n: u64) -> bool

Source

pub fn trial_division(&self, n: u64) -> PrimeListTrialDivision<'_>

Examples found in repository ?

crates/competitive/src/math/prime_list.rs (line 45)

44    pub fn prime_factors(&self, n: u64) -> Vec<(u64, u32)> {
45        self.trial_division(n).collect()
46    }
47    pub fn count_divisors(&self, n: u64) -> u64 {
48        let mut divisor_cnt = 1u64;
49        for (_, cnt) in self.trial_division(n) {
50            divisor_cnt *= cnt as u64 + 1;
51        }
52        divisor_cnt
53    }
54    pub fn divisors(&self, n: u64) -> Vec<u64> {
55        let mut d = vec![1u64];
56        for (p, c) in self.trial_division(n) {
57            let k = d.len();
58            let mut acc = 1;
59            for _ in 0..c {
60                acc *= p;
61                for i in 0..k {
62                    d.push(d[i] * acc);
63                }
64            }
65        }
66        d.sort_unstable();
67        d
68    }

Source

pub fn reserve(&mut self, max_n: u64)

list primes less than or equal to max_n by segmented sieve

Examples found in repository ?

crates/competitive/src/math/prime_list.rs (line 21)

19    pub fn new(max_n: u64) -> Self {
20        let mut self_: Self = Default::default();
21        self_.reserve(max_n);
22        self_
23    }
24    pub fn primes(&self) -> &[u64] {
25        self.primes.as_slice()
26    }
27    pub fn primes_lte(&self, n: u64) -> &[u64] {
28        assert!(n <= self.max_n, "expected `n={} <= {}`", n, self.max_n);
29        let i = self.primes.partition_point(|&p| p <= n);
30        &self.primes[..i]
31    }
32    pub fn is_prime(&self, n: u64) -> bool {
33        assert!(n <= self.max_n, "expected `n={} <= {}`", n, self.max_n);
34        self.primes.binary_search(&n).is_ok()
35    }
36    pub fn trial_division(&self, n: u64) -> PrimeListTrialDivision<'_> {
37        let bound = self.max_n.saturating_mul(self.max_n);
38        assert!(n <= bound, "expected `n={} <= {}`", n, bound);
39        PrimeListTrialDivision {
40            primes: self.primes.iter(),
41            n,
42        }
43    }
44    pub fn prime_factors(&self, n: u64) -> Vec<(u64, u32)> {
45        self.trial_division(n).collect()
46    }
47    pub fn count_divisors(&self, n: u64) -> u64 {
48        let mut divisor_cnt = 1u64;
49        for (_, cnt) in self.trial_division(n) {
50            divisor_cnt *= cnt as u64 + 1;
51        }
52        divisor_cnt
53    }
54    pub fn divisors(&self, n: u64) -> Vec<u64> {
55        let mut d = vec![1u64];
56        for (p, c) in self.trial_division(n) {
57            let k = d.len();
58            let mut acc = 1;
59            for _ in 0..c {
60                acc *= p;
61                for i in 0..k {
62                    d.push(d[i] * acc);
63                }
64            }
65        }
66        d.sort_unstable();
67        d
68    }
69    /// list primes less than or equal to `max_n` by segmented sieve
70    pub fn reserve(&mut self, max_n: u64) {
71        if max_n <= self.max_n || max_n < 2 {
72            return;
73        }
74
75        if self.primes.is_empty() {
76            self.primes.push(2);
77            self.max_n = 2;
78        }
79        if max_n == 2 {
80            return;
81        }
82
83        let max_n = (max_n + 1) / 2 * 2; // odd
84        let sqrt_n = ((max_n as f64).sqrt() as usize + 1) / 2 * 2; // even
85        let mut table = Vec::with_capacity(sqrt_n >> 1);
86        if self.max_n < sqrt_n as u64 {
87            let start = (self.max_n as usize + 1) | 1; // odd
88            let end = sqrt_n + 1;
89            let sqrt_end = (sqrt_n as f64).sqrt() as usize;
90            let plen = self.primes[1..]
91                .binary_search(&(sqrt_end as u64 + 1))
92                .unwrap_or_else(|x| x);
93            table.resize(end / 2 - start / 2, false);
94            for &p in self.primes.iter().skip(1).take(plen) {
95                let y = p.max((start as u64 + p - 1) / (2 * p) * 2 + 1) * p / 2;
96                (y as usize - start / 2..end / 2 - start / 2)
97                    .step_by(p as usize)
98                    .for_each(|i| table[i] = true);
99            }
100            for i in 0..=(sqrt_end / 2).saturating_sub(start / 2) {
101                if !table[i] {
102                    let p = (i + start / 2) * 2 + 1;
103                    for j in (p * p / 2 - start / 2..sqrt_n / 2 - start / 2).step_by(p) {
104                        table[j] = true;
105                    }
106                }
107            }
108            self.primes
109                .extend(table.iter().cloned().enumerate().filter_map(|(i, b)| {
110                    if !b {
111                        Some((i + start / 2) as u64 * 2 + 1)
112                    } else {
113                        None
114                    }
115                }));
116            self.max_n = sqrt_n as u64;
117        }
118
119        let sqrt_n = sqrt_n as u64;
120        for start in (self.max_n + 1..=max_n).step_by(sqrt_n as usize) {
121            let end = (start + sqrt_n).min(max_n + 1);
122            let sqrt_end = (end as f64).sqrt() as u64;
123            let length = end - start;
124            let plen = self.primes[1..]
125                .binary_search(&(sqrt_end + 1))
126                .unwrap_or_else(|x| x);
127            table.clear();
128            table.resize(length as usize / 2, false);
129            for &p in self.primes.iter().skip(1).take(plen) {
130                let y = p.max((start + p - 1) / (2 * p) * 2 + 1) * p / 2;
131                ((y - start / 2) as usize..length as usize / 2)
132                    .step_by(p as usize)
133                    .for_each(|i| table[i] = true);
134            }
135            self.primes
136                .extend(table.iter().cloned().enumerate().filter_map(|(i, b)| {
137                    if !b {
138                        Some((i as u64 + start / 2) * 2 + 1)
139                    } else {
140                        None
141                    }
142                }));
143        }
144        self.max_n = max_n;
145    }
146}
147
148#[derive(Debug, Clone)]
149pub struct PrimeListTrialDivision<'p> {
150    primes: Iter<'p, u64>,
151    n: u64,
152}
153impl Iterator for PrimeListTrialDivision<'_> {
154    type Item = (u64, u32);
155    fn next(&mut self) -> Option<Self::Item> {
156        if self.n <= 1 {
157            return None;
158        }
159        loop {
160            match self.primes.next() {
161                Some(&p) if p * p <= self.n => {
162                    if self.n % p == 0 {
163                        let mut cnt = 1u32;
164                        self.n /= p;
165                        while self.n % p == 0 {
166                            cnt += 1;
167                            self.n /= p;
168                        }
169                        return Some((p, cnt));
170                    }
171                }
172                _ => break,
173            }
174        }
175        if self.n > 1 {
176            return Some((replace(&mut self.n, 1), 1));
177        }
178        None
179    }
180}
181
182pub fn with_prime_list<F>(max_n: u64, f: F)
183where
184    F: FnOnce(&PrimeList),
185{
186    thread_local!(static PRIME_LIST: UnsafeCell<PrimeList> = Default::default());
187    PRIME_LIST.with(|cell| {
188        unsafe {
189            let pl = &mut *cell.get();
190            pl.reserve(max_n);
191            f(pl);
192        };
193    });
194}

More examples

Hide additional examples

crates/competitive/src/math/discrete_logarithm.rs (line 70)

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 u64;
70        self.primes.reserve(lim);
71        let primes = self.primes.primes_lte(lim);
72        while self.br_primes.len() < primes.len() {
73            let br = BarrettReduction::<u64>::new(primes[self.br_primes.len()]);
74            self.br_primes.push(br);
75        }
76        let br_primes = &self.br_primes[..primes.len()];
77        self.ic
78            .entry(p)
79            .or_insert_with(|| IndexCalculusWithPrimitiveRoot::new(p, br_primes))
80            .discrete_logarithm(a, b, br_primes)
81    }