5
\$\begingroup\$

For fun I wanted to investigate Brocard's Conjecutre. That is, does there exist \$n, j\$ such that \$n! + 1 = j^2\$. The method I use takes inspiration from looking at divisibility rules which I explain in the doc comment at the top of my Rust file, but to reiterate at a high level:

The idea is squares can only take certain values \$\mod m\$ for some value of \$m\$. I basically check a lot of modulo classes to see if \$n! + 1\$ is equal to some number a square is allowed to take. If it takes some value a square is not allowed to take then we can eliminate it as a potential candidate. Note, there can be false positives, but I have not found one yet, any advice for speeding up would be greatly appreciated (of course, I run with cargo run --release)!

Cargo.toml

[package]
name = "brocards_problem"
version = "0.1.0"
edition = "2021"

[dependencies]
rustc-hash = "*"

src/main.rs

/// Brocard's Problem Candidate Finder 
/// [https://en.wikipedia.org/wiki/Brocard's_problem]
///
/// Uses divisibility rules to determine if n! + 1 != j^2 for some n and j.
///
/// That is, compute what possible values for a square value to be mod m and
/// determine if n! + 1 mod m, m+1, ..., m+RULES_UPPER if it lands on one of
/// those values for each modularity class. If it does, it is potentially a
/// square and therefore a candidate solultion to Brocard's problem. 
/// Otherwise, it can be eliminated as a candidate.
///
/// Note that n! = 0 mod m when n > m so only values of m > n are computed.

use rustc_hash::FxHashSet;

/// Computes the factorial of a number but does so mod m as to prevent
/// blowup in size.
fn mod_factorial(mut n: usize, m: usize) -> usize {
    let mut res = 1;
    while n > 0 {
        res *= n;
        res %= m;
        n -= 1;
    }
    res
}

fn main() {
    // Upper bound on the number of divisibility rules to check.
    const RULES_UPPER: usize = 300;
    let mut n = 8;
    'outer: loop {
        println!("Testing {n:?}");
        for i in n+1..RULES_UPPER+n {
            let mut rules: FxHashSet<usize> = FxHashSet::default();
            for j in 0..i {
                rules.insert((j * j) % i);
            }
            if !rules.contains(&((mod_factorial(n, i) + 1) % i)) {
                n += 1;
                continue 'outer;
            }
        }
        println!("Found {n:?}");
        return    
    }
}
\$\endgroup\$

2 Answers 2

3
\$\begingroup\$

Eliminate Redundant Modular Checks

You spend a lot of time doing mod operations, which are slow on most CPUs. This includes recalculating \$n!+1\$ (mod \$m\$) for each \$m\$. But most of these are redundant.

By the Chinese Remainder Theorem, checking any set of co-prime numbers is equivalent to checking their product. Thus, if you’ve checked for an odd number \$m\$, checking for \$2m\$ is redundant (and ¼ of all whole numbers are twice an odd number). It’s equivalent to testing both mod \$m\$ and mod \$2\$. But checking mod \$2\$ tells us nothing about whether the original number might be a square number: those can be even or odd.

Similarly, testing modulo a number also tests for all factors of that number. A special case of this is that testing a prime power tests for all smaller powers of the prime, and testing powers of multiple primes tests all products of powers of those primes (of equal or lower degree).

If you were testing arbitrary numbers to find if they were potentially square, it would be efficient to check mod a product of small primes: checking mod 3 × 5 × 7 × 11 × 13 = 15 015 is enough to filter out more than 97% of numbers, and going up to 23, more than 99% of them, with a single mod operation and a table lookup. However, that won’t work well here, because \$n! + 1 \cong 1\$ modulo any number from 2 to \$n\$, so the numbers we’re interested in always pass that check.

Check Modulo Prime Powers

If \$m \cong p \mod p^2\$ for prime \$p\$, \$p\$ is a prime factor with multiplicity 1, so \$m\$ cannot be a perfect square. (The converse is not true: 8 is evenly divisible by the square of its only prime factor, but is not a perfect square.) More generally, a prime \$p\$ is a factor of \$m\$ with multiplicity \$i\$ iff, for \$j > i\$, \$m \cong p^i \mod p^j\$. And any number that is not a perfect square has at least one prime factor with odd multiplicity.

This check might be worth pulling out separately because you can vectorize the check for every prime in range using very little memory. Or, since all powers of different primes are co-prime to each other, you can use the Chinese Remainder Theorem to test multiple squared primes at once (and this result will fall out naturally). Don’t bother testing any number smaller than half the factorial; the modular residue by its square will be exactly 1.

You could also test modulo large powers of primes less than or equal to \$\frac{n! + 1}{2}\$. At minimum, checking a few of these in constant time tells us that any false positive must have a factor that is not one of the prime powers we checked, narrowing the search space drastically. This will also tell us if a lower power of the prime is a factor—and any number with a prime factor of odd multiplicity cannot be a perfect square. If we additionally want to check the modular residue for a potential square, and not only check for prime factors of odd multiplicity, we can either (given enough gigabytes of program text) pre-calculate the lookup tables for every prime, or generate them lazily as needed and memo-ize them. Speaking of which:

Use Dynamic Programming

A program to test this conjecture would be checking a large number of values. It therefore wants to remember its tables and re-use them on future runs. Lazy static initialization would work well here in Rust, or you could attempt to meta-program, and compile the lookup tables into the program text.

Use Bitfields

A set is probably implemented as a binary search tree. A bitfield has far superior lookup and insertion time, and probably, for most of the numbers you’re crunching, the same or less memory usage.

Handle Powers of 2 Specially

Although I said above that % is the slowest ALU operation, calculating modulo a power of 2 is one of the fastest. Just calculating the factorial as u8 would eliminate 87% of numbers. This is, unfortunately, not much value for this problem, since \$n! + 1 \cong 1\$ (mod 256) for all \$n \ge 10\$. You could, however, get the multiplicity of 2 in the prime factorization very quickly, by counting the trailing zeroes in the binary representation, eliminate the candidate if the multiplicity is odd, and get the odd modular residue of the candidate by right-shifting.

\$\endgroup\$
6
  • \$\begingroup\$ Hi, thanks for the review! One thing I am confused about though: I don't believe I ever check modulo 2. Perhaps I should have commented this better, but I check for m > n (as mentioned in the top comment). m is taken by the role i in my program. Thus, for i in n+1..RULES_UPPER+n precludes mod 2. Unless I am misunderstanding something. \$\endgroup\$
    – Dair
    Commented May 15 at 20:39
  • \$\begingroup\$ @Dair You’re correct. I was trying to introduce the concept that checking by the product of 2 and any odd number already checked is unnecessary, because that is the same as checking bu the odd number and by 2 separately. Also, 10! and up is divisible by 256, not 16!. \$\endgroup\$
    – Davislor
    Commented May 15 at 20:45
  • \$\begingroup\$ Sorry to bother you with another question: How do you use a bit field as a set? Is the idea to encode every remainder as a segment in the bitfield? Also, I believe the FxHashSet does not use a tree (this would be BTreeSet) but a hash function if that matters. \$\endgroup\$
    – Dair
    Commented May 15 at 21:08
  • \$\begingroup\$ @Dair A hash function would work better than a BST, but still not as well as a bitfield. If your set is a bitfield, you find the i-th bit of the field and set or test that bit. If the field is an array of u8, you take the bitwise or of word a[i/8] and the bitmask 1u8 << (i%8). This can also work for other word sizes. \$\endgroup\$
    – Davislor
    Commented May 15 at 21:18
  • \$\begingroup\$ I bet there’s a library for this already, but I’m about to go AFK, so I can’t look it up right now. \$\endgroup\$
    – Davislor
    Commented May 15 at 21:20
6
\$\begingroup\$

Do not use "*" as a dependency version requirement. You probably didn't care what version was picked when starting and practically it doesn't matter for rustc-hash (since there is only one major version and a 2.0 seems unlikely), but it is a bad habit to get into. If there is ever a new breaking version of rustc-hash, cargo may choose to use that and may break your code. Even with a Cargo.lock file preserving versions, you may do a cargo update, you may lose your lock file, or you may have another contributor/use that takes your code without the lock file and may end up having problems. Publishing your package to crates.io will be rejected if you have a "*" dependency.

In general you should specify the version you are actually using, even minor and patch versions, since those may have introduced features you are using that may be missing in a prior version (and possibly break your code by extensions of above reasons). You may want to use a lesser minimal version like 1.0.0 if creating a library for maximum usability (though you should be sure to actually test that your code works with 1.0.0).

So put this:

rustc-hash = "1.1.0"

If you use cargo add rustc-hash the current version will be fully specified by default, or if browsing crates.io it has the TOML dependency version requirement spelled out on the right side you can just copy-paste.

\$\endgroup\$
1
  • \$\begingroup\$ Noted: I will be sure to do this in the future! Also, thanks for the cargo add tip. \$\endgroup\$
    – Dair
    Commented May 14 at 18:03

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.