3
\$\begingroup\$

I'm trying to efficiently generate all pairs of factors of a number n.

For example, when n=27, the factor pairs are (1, 27), (3, 9), (9, 3), (27, 1)

The order in which the pairs are found is not important for my use case.

I've put together a solution using the crates primal and itertools:

  • The prime sieve yields the prime factorisation of n in (prime, exponent) tuples (e.g. for n=900: (2, 2), (3, 2), (5, 2) )
  • Using the cartesian product of iterators over prime^0..prime^exponent, finding all divisors d by calculating the product of each Vec, and pairing them with n/d.

This is the code I have now:

#[macro_use]
extern crate lazy_static;

use itertools::Itertools;

const LIMIT: usize = 50_000_000;

lazy_static! {
    static ref SIEVE: primal::Sieve = primal::Sieve::new(LIMIT);
}

fn main() {
    for n in 2..LIMIT {
        for (d1, d2) in SIEVE
            .factor(n)
            .unwrap()
            .into_iter()
            .map(|(base, exp)| (0..=(exp as u32)).map(move |e| base.pow(e)))
            .multi_cartesian_product()
            .map(|v| {
                let d = v.into_iter().product::<usize>();
                (d, n / d)
            })
        {
            println!("({d1}, {d2})");
        }
    }
}

It works, and it's much faster than the most naïve approach, but it's still a lot slower than I would like.

Can anyone suggest either improvements to my code or a different approach that runs significantly faster?

\$\endgroup\$
6
  • 1
    \$\begingroup\$ Is the loop of n from 2 to LIMIT just for a benchmark? If you wanted to do that for real, then that structure can be exploited to save some redundant work. \$\endgroup\$
    – user555045
    Mar 2, 2023 at 1:39
  • 1
    \$\begingroup\$ Assign d = max(d1, d2). As soon as d * d > n you can stop. It's impossible for one of n's factors to exceed sqrt(n). \$\endgroup\$
    – J_H
    Mar 2, 2023 at 2:18
  • 1
    \$\begingroup\$ Not strictly true, @J_H - about half of factors are greater than √n, including n itself. But they can all be found by considering the smaller factors, e.g. by producing {d1, d2} and {d2, d1} together, except when they are equal (d1 == d2 == √n). \$\endgroup\$ Mar 2, 2023 at 8:53
  • \$\begingroup\$ The biggest change you can do is remove the println. This suggest that you should rethink your benchmarking strategy. \$\endgroup\$ Mar 2, 2023 at 22:02
  • \$\begingroup\$ Thanks for the comments! I'm using the println just as a placeholder for the logic that actually goes there, and I do actually need to go up to 50m. @harold good point, I've actually managed to eliminate a few cases that don't need to be checked :) TobySpeight that could work in theory, but unfortunately the multi_cartesian_product combinator does not yield them in ascending order. If I had another way to generate them in order, that might help, but it then also requires another .flatten() when yielding 2 tuples on each iteration. \$\endgroup\$
    – Bram
    Mar 3, 2023 at 14:50

1 Answer 1

2
\$\begingroup\$

My solution with a custom iterator runs almost 2x faster. Your main bottleneck is the multi_cartesian_product iterator that produces results on the heap (Vec<usize>). There is no inherent need for getting the factors through the heap.

Speedup for LIMIT=500_000:

test tests::bench_factors  ... bench: 250,582,325 ns/iter (+/- 1,694,499)
test tests::bench_factors2 ... bench: 445,134,232 ns/iter (+/- 2,520,121)

Solution with the custom CombineProducts:

#![feature(test)]

extern crate test;

#[macro_use]
extern crate lazy_static;

use itertools::Itertools;

const LIMIT: usize = 500_000;

lazy_static! {
    static ref SIEVE: primal::Sieve = primal::Sieve::new(LIMIT);
}

struct CombineProducts {
    factors: Vec<Factor>,
    current_product: usize,
    exhausted: bool,
}

struct Factor {
    base: usize,
    exp: usize,
    current_exp: usize,
}

impl CombineProducts {
    fn new(factors: Vec<(usize, usize)>) -> CombineProducts {
        CombineProducts {
            factors: factors.into_iter().map(|(base, exp)| Factor { base, exp, current_exp: 0, }).collect(),
            current_product: 1,
            exhausted: false,
        }
    }
}

impl Iterator for CombineProducts {
    type Item = usize;

    fn next(&mut self) -> Option<Self::Item> {
        if self.exhausted {
            return None;
        }
        let result = self.current_product;
        for factor in &mut self.factors {
            if factor.current_exp == factor.exp {
                factor.current_exp = 0;
                self.current_product /= factor.base.pow(factor.exp as u32);
            } else {
                factor.current_exp += 1;
                self.current_product *= factor.base;
                return Some(result);
            }
        }
        self.exhausted = true;
        Some(result)
    }
}

fn run() -> Vec<(usize, usize)> {
    let mut result = vec![];
    for n in 2..LIMIT {
        result.extend(
            CombineProducts::new(
                SIEVE
                    .factor(n)
                    .unwrap()
            ).map(|d| (d, n / d))
        );
    }
    result
}

fn run_print() -> Vec<(usize, usize)> {
    let mut result = vec![];
    for n in 2..LIMIT {
        result.extend(
            CombineProducts::new(
                SIEVE
                    .factor(n)
                    .unwrap()
            ).map(|d| (d, n / d))
        );
    }
    println!("{:?}", result.len());
    result
}

fn run2() -> Vec<(usize, usize)> {
    let mut result = vec![];
    for n in 2..LIMIT {
        result.extend(
            SIEVE
                .factor(n)
                .unwrap()
                .into_iter()
                .map(|(base, exp)| (0..=(exp as u32)).map(move |e| base.pow(e)))
                .multi_cartesian_product()
                .map(|v| {
                    let d = v.into_iter().product::<usize>();
                    (d, n / d)
                })
        );
    }
    result
}

fn run2_print() -> Vec<(usize, usize)> {
    let mut result = vec![];
    for n in 2..LIMIT {
        result.extend(
            SIEVE
                .factor(n)
                .unwrap()
                .into_iter()
                .map(|(base, exp)| (0..=(exp as u32)).map(move |e| base.pow(e)))
                .multi_cartesian_product()
                .map(|v| {
                    let d = v.into_iter().product::<usize>();
                    (d, n / d)
                })
        );
    }
    println!("{:?}", result.len());
    result
}

#[cfg(test)]
mod tests {
    use super::*;
    use test::bench::Bencher;

    #[test]
    fn test_run() {
        run_print();
    }
    #[test]
    fn test_run2() {
        run2_print();
    }

    #[bench]
    fn bench_factors(b: &mut Bencher) {
        b.iter(|| {
            run()
        });
    }

    #[bench]
    fn bench_factors2(b: &mut Bencher) {
        b.iter(|| {
            run2()
        });
    }
}
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.