# Find Ramanujan-Hardy numbers smaller than N

Problem statement: Write a program that takes a command-line argument n and prints all integers less than n that can be expressed as the sum of two cubes in two different ways i.e. distinct positive integers $$\a\$$, $$\b\$$, $$\c\$$, and $$\d\$$ such that $$\a^3 + b^3 = c^3 + d^3\$$.

This is one of my self-imposed challenges in Rust to become better at it. The problem was taken from Sedgewick Exercise 1.3.34.

Here is my code:

use clap::Parser;

#[derive(Debug, Parser)]
struct Arguments {
#[arg(index = 1)]
n: usize,
}

fn main() {
let arguments = Arguments::parse();
let mut ramanujan_numbers = Vec::new();

for a in 1..arguments.n {
if a.pow(3) > arguments.n {
break;
}

for b in (a + 1)..arguments.n {
if b.pow(3) > arguments.n {
break;
}

for c in (b + 1)..arguments.n {
if c.pow(3) > arguments.n {
break;
}

for d in (c + 1)..arguments.n {
if d.pow(3) > arguments.n {
break;
}

let first_condition = a.pow(3) + b.pow(3) == c.pow(3) + d.pow(3);
let second_condition = a.pow(3) + c.pow(3) == b.pow(3) + d.pow(3);
let third_condition = a.pow(3) + d.pow(3) == b.pow(3) + c.pow(3);

if first_condition {
let ramanujan_number = a.pow(3) + b.pow(3);

if ramanujan_number < arguments.n {
ramanujan_numbers.push(ramanujan_number);
}
} else if second_condition {
let ramanujan_number = a.pow(3) + c.pow(3);

if ramanujan_number < arguments.n {
ramanujan_numbers.push(ramanujan_number);
}
} else if third_condition {
let ramanujan_number = a.pow(3) + d.pow(3);

if ramanujan_number < arguments.n {
ramanujan_numbers.push(ramanujan_number);
}
}
}
}
}
}

ramanujan_numbers.sort();

match ramanujan_numbers.len() {
0 => println!(
"No Ramanujan number smaller than {} was found.",
arguments.n
),
1 => println!(
"The Ramanujan number smaller than {} is {:?}.",
arguments.n, ramanujan_numbers
),
_ => println!(
"The Ramanujan numbers smaller than {} are {:?}.",
arguments.n, ramanujan_numbers
),
}
}


Is there any way that I can improve my code?

## Coding Style

### Factor the Algorithm into Helper Functions

Right now, you’ve got two levels of nested if, several with else blocks, nested inside five levels of loops.

This could really benefit from some helper functions. Let’s say I’m giving you the code review and I say, “Oh, we have a function like this already. Maybe you could use that in your algorithm.”

/* Finds all pairs (a,b) such that a**3 + b**3 = sum.
*/
fn taxicab_pairs(sum: usize) -> Vec<(usize, usize)>


The code would get a whole lot simpler now, right? For example, your loop might now look something like

for i in 1..arguments.n {
let cube_pairs = taxicab_pairs(i);
if cube_pairs.len() > 1 {
// Print the pairs out.
}
}


Or if you just want a list of Ramanujan-Hardy numbers like you have now, you might do something like that with iterators.

### Prefer Iterator Expressions to Mutable Objects

You currently initialize

let mut ramanujan_numbers = Vec::new();


and push numbers onto it, one at a time.

But you could also build this list of numbers at once, with a bit of functional programming, and make it immutable. That not only prevents entire categories of bugs, the compiler converts everything into static single assignments during its optimization passes, so you’re making the code easier to optimize, too.

That might look something like

let ramunujan_numbers: Vec<usize> = (1..arguments.n)
.filter(is_taxicab_number)
.collect();


Or if we have taxicab_pairs, maybe we’d rather keep the list of (a, b) and (c, d) pairs we calculate with that, instead of sending it to a print loop.

let sets_of_pairs: Vec<Vec<(usize, usize)>> = (1..arguments.n)
.map(taxicab_pairs) // Gives us a sequence of Vec<(usize, usize)>
.filter(move |v| v.len() > 1) // Removes all with fewer than two pairs
.collect();


Then we can keep the code at the end to handle the case where there aren’t any taxicab numbers in range.

You might also choose to remove the .collect() lines, and let the compiler infer the type, to obtain lazily-evaluated iterators instead of a strictly-evaluated Vec. This lets you generate the sequence as you calculate, and print each output immediately, rather than waiting for all the output to begin printing anything. You would then need to use a different method to see if there are any taxicab numbers in range (such as making it Peekable and peeking ahead to see whether the iterator immediately returns None).

    let mut sets_of_pairs = (1..n)
.map(taxicab_pairs) // Gives us a sequence of Vec<(usize, usize)>
.filter(move |v| v.len() > 1) // Removes all with fewer than two pairs
.peekable();

if sets_of_pairs.peek().is_none() {
println!("No taxicab numbers less than {}.", n);
}
// No need for an else; the loop will just fall through in that case.


### Make the Type of Each Variable Obvious

Right now, you declare a list like

let mut ramanujan_numbers = Vec::new();


One way this could be improved is that we can’t tell what’s in the vector. The compiler tries to infer it from elsewhere, but that’s very fragile to refactoring. Maybe the line that sets this tells the compiler the type moves out of scope, or maybe it also tries to infer the type out of context, or maybe there are multiple places this gets updated, and not all of them agree.

Even if the compiler can tell what the type of ramanujan_numbers is, it’s needlessly complicated for a human maintaining the code. So, annotate it when it’s declared, or (if you aren’t replacing this with a static single assignment) call Vec::<usize>::new() with the turbofish operator, or something else that makes it clear.

### Simplify the Conditionals

The first_condition, second_condition and third_comdition block is unwieldy, although it can be factored out entirely.

### The Command-Line Arguments

You might have chosen the way you did intentionally, but there is support for this in the standard library. You could use that to check for runtime errors and choose which message to print.

## Performance

### Prune the Search Space with Some Easy (Really!) Number Theory

Okay, let’s say we have two whole numbers, a and b. One thing we know about them already is that either a is bigger, b is bigger, or they’re equal. You already noticed that you can throw out all the cases where b > a, and cut your search time nearly in half. That’s a great idea, so let’s stick with it.

Here, we’re interested in an x and a y that add up to some number z. If x and y are equal, they’re both half of z. They can’t both be bigger than half of z, or they wouldn’t add up to z. So one is bigger than z/2 and one is smaller. If they aren’t equal to each other, and we just said that x isn’t bigger than y either, x must be the one that’s smaller than z/2 and y the one that’s bigger.

But, for this problem, we actually want x and y to be perfect cubes! So, x = z/2 and y = z/2. The most interesting part of that is z/2. Since a is positive, we can take the cube root of both sides and solve for a. In other words, we only need to search up to the cube root of half of the target number for a, and we’ll be sure we’ve found every possible smaller cube that adds up with another to the right number.

Now, if we have a, we have x = , and we started with z, we can solve x + y = z for y = z - x. If y is also a perfect cube, we’re done. There are a few ways we chan check this, but one is to convert to floating-point, find the cube root, round off to an integer, and see if its cube is actually the right number.

Since the program specified that we only need to find out whether a number is the sum of two different pairs of positive cubes, not what all those pairs are, we could optimize further. (For example, we can prove that a+b is congruent to z, modulo 3, reducing the search space further.) But I’ll leave it there.

## Putting It All Together

use std::{env, process};

/* Finds all pairs (a,b) such that a**3 + b**3 = sum.
*/
fn taxicab_pairs(sum: usize) -> Vec<(usize, usize)> {
/* If positive integers a**3 + b**3 = i, and WLOG a <= b,
* a <= pivot and b >= pivot.
*/
let pivot: usize = f64::floor(f64::powf(sum as f64 / 2.0, 1.0 / 3.0) + f64::EPSILON) as usize;

(1..=pivot)
.map(move |a| (a, sum - a * a * a))
.map(move |(a, b_cubed)| {
(
a,
f64::round(f64::powf(b_cubed as f64, 1.0 / 3.0)) as usize,
b_cubed,
)
})
.filter(move |&(_, b, b_cubed)| b * b * b == b_cubed)
.map(move |(a, b, _)| (a, b))
.collect()
}

pub fn main() {
const USAGE_MSG: &str = "Call this program with one command-line argument, N.\n";
let args: Vec<String> = env::args().collect();

let n: usize = if args.len() == 2 {
match args[1].parse::<usize>() {
Ok(n) => n,
_ => {
eprintln!("{}", USAGE_MSG);
process::exit(1)
}
}
} else {
eprintln!("{}", USAGE_MSG);
process::exit(1)
};

let sets_of_pairs: Vec<Vec<(usize, usize)>> = (1..n)
.map(taxicab_pairs) // Gives us a sequence of Vec<(usize, usize)>
.filter(move |v| v.len() > 1) // Removes all with fewer than two pairs
.collect();

if sets_of_pairs.is_empty() {
println!("No taxicab numbers less than {}.", n);
}
// No need for an else; the loop will just fall through in that case.

for cube_pairs in sets_of_pairs {
let (a, b) = cube_pairs[0]; // Safe, because we only added sets of at least two pairs.
let i = a * a * a + b * b * b;
let residue = i % 3;

print!("{}", i);
for (a, b) in cube_pairs {
assert_eq!(i, a * a * a + b * b * b);
assert_eq!(residue, (a + b) % 3); // You might prove this congruence.
print!(" = {}³ + {}³", a, b);
}
println!(".");
}
}


You might prefer the use of a helper function to a chain of map calls:

/* Finds all pairs (a,b) such that a**3 + b**3 = sum.
*/
fn taxicab_pairs(sum: usize) -> Vec<(usize, usize)> {
/* Given a and sum, finds b such that a**3 + b**3 == sum, if b exists.
* Returns either Some((a,b)) or None.
*/
fn find_dual(sum: usize, a: usize) -> Option<(usize, usize)> {
let b_cubed = sum - a*a*a;
let b = f64::round(f64::powf(b_cubed as f64, 1.0 / 3.0)) as usize;
if b*b*b == b_cubed {
Some((a, b))
} else {
None
}
}

/* If positive integers a**3 + b**3 = i, and WLOG a <= b,
* a <= pivot and b >= pivot.
*/
let pivot: usize = f64::floor(f64::powf(sum as f64 / 2.0, 1.0 / 3.0) + f64::EPSILON) as usize;

(1..=pivot)
.filter_map(move|a| find_dual(sum, a))
.collect()
}


### Use sort_unstable()

When sorting a list of numbers, stability doesn't matter. So you should prefer sort_unstable() over sort(), as it is faster (and non-allocating).

### Extract arguments.n to a variable

This is subjective, but I would prefer let n = arguments.n; for shorter variable name. This also allows you to use embedded formatting in println!() (println!("{n}")).

### Use format_args_capture (embedded formatting)

Replace:

println!(
"The Ramanujan number smaller than {} is {:?}.",
arguments.n, ramanujan_numbers
)


With:

println!(
"The Ramanujan number smaller than {} is {ramanujan_numbers:?}.",
arguments.n,
)


Not everyone agrees with that: some people prefer all formatting specifiers to be the same, either out-of-line or embedded.

### Use trailing comma for println!() when spanning multiple lines

Usually rustfmt adds the trailing comma, but it can't handle macros so it doesn't do that here.

### Do not rely on Debug formatting

...and generally do not print Debug formatting to the user.

It is not stable, and can change without warning. It is intended to print information to the programmer, and not to the user.

### Unify the two last println!()

Since the difference is only one s and is/are, you can use a nice trick:

if ramanujan_numbers.is_empty() {
println!("No Ramanujan number smaller than {} was found.", arguments.n);
} else {
println!(
"The Ramanujan number{} smaller than {} {} {:?}.",
if ramanujan_numbers.len() > 1 { "s" } else { "" },
arguments.n,
if ramanujan_numbers.len() > 1 { "are" } else { "is" },
ramanujan_numbers
);
}


This also avoids the duplicated code in case of more involved formatting.

### Final code

use clap::Parser;

#[derive(Debug, Parser)]
struct Arguments {
#[arg(index = 1)]
n: usize,
}

fn main() {
let Arguments { n } = Arguments::parse();
let mut ramanujan_numbers = Vec::new();

for a in 1..n {
if a.pow(3) > n {
break;
}

for b in (a + 1)..n {
if b.pow(3) > n {
break;
}

for c in (b + 1)..n {
if c.pow(3) > n {
break;
}

for d in (c + 1)..n {
if d.pow(3) > n {
break;
}

let first_condition = a.pow(3) + b.pow(3) == c.pow(3) + d.pow(3);
let second_condition = a.pow(3) + c.pow(3) == b.pow(3) + d.pow(3);
let third_condition = a.pow(3) + d.pow(3) == b.pow(3) + c.pow(3);

if first_condition {
let ramanujan_number = a.pow(3) + b.pow(3);

if ramanujan_number < n {
ramanujan_numbers.push(ramanujan_number);
}
} else if second_condition {
let ramanujan_number = a.pow(3) + c.pow(3);

if ramanujan_number < n {
ramanujan_numbers.push(ramanujan_number);
}
} else if third_condition {
let ramanujan_number = a.pow(3) + d.pow(3);

if ramanujan_number < n {
ramanujan_numbers.push(ramanujan_number);
}
}
}
}
}
}

ramanujan_numbers.sort_unstable();

if ramanujan_numbers.is_empty() {
println!("No Ramanujan number smaller than {n} was found.");
} else {
println!(
"The Ramanujan number{} smaller than {n} {} {}.",
if ramanujan_numbers.len() > 1 { "s" } else { "" },
if ramanujan_numbers.len() > 1 { "are" } else { "is" },
ramanujan_numbers
.iter()
.map(|num| num.to_string())
.collect::<Vec<_>>()
.join(", "),
);
}
}
$$$$


# Performance

Your code is currently O(n⁴). You can reduce that to O(n²) using the Cartesian product on the candidates. The library itertools comes in handy here.

# Divide and conquer

Your current function main() is pretty long. Consider breaking the problem up into different parts.

# Suggested

Cargo.toml

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

# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html

[dependencies]
clap = { version= "4.3.5", features=["derive"] }
itertools = "0.11.0"


src/lib.rs

use itertools::Itertools;
use std::collections::HashSet;

pub fn ramanujan(limit: u64) -> impl Iterator<Item = (u64, HashSet<(u64, u64)>)> {
(0..=limit)
.map(|n| (n, cubes(n)))
.filter(|(_, cubes)| is_ramanujan(cubes))
}

fn is_ramanujan(cubes: &HashSet<(u64, u64)>) -> bool {
cubes.len() >= 2
}

fn cubes(n: u64) -> HashSet<(u64, u64)> {
let candidates = candidates(n);
candidates
.iter()
.cartesian_product(candidates.iter())
.filter(|(a, b)| a.pow(3) + b.pow(3) == n)
.map(|(a, b)| if a < b { (*a, *b) } else { (*b, *a) })
.collect()
}

fn candidates(limit: u64) -> Vec<u64> {
(0..=limit).take_while(|n| n.pow(3) <= limit).collect()
}


src/main.rs

use clap::Parser;
use ramanujan::ramanujan;

#[derive(Debug, Parser)]
struct Arguments {
#[arg(index = 1)]
n: u64,
}

fn main() {
let args = Arguments::parse();
for (n, cubes) in ramanujan(args.n) {
println!("{n} is a Ramanujan number that can be written as:");

for (a, b) in cubes {
println!("  * {n} = {a}³ + {b}³");
}
}
}
`

The above code is still not fully optimized, but it showcases possible improvements.

• Great suggestions. A lot like what I came up with. The bound can be brought down to at most O(N**(4/3)), with a pruning lemma. Commented Jul 12, 2023 at 7:36
• Yeah, I just read your review. You seem to know a lot more number theory than me. In my review I focused more on structural improvements. Commented Jul 12, 2023 at 7:42