3
\$\begingroup\$

I've been playing with perfect, amicable and social numbers lately, and created a class for investigating these. At present, it has functions to return the perfect and amicable numbers in a specified range (and also the primes, since they are found as a useful side-effect). It could be expanded to explore the implied directed-graph to find social numbers and classify them by their cycle length, and perhaps other interesting analyses of the the sums of proper factors.

I've tried to make it reasonably performant without obscuring the mode of operation too much. I found that I needed to comment more than I normally would, to get it both fast and readable.

Performance benchmarks on a couple of systems I have available (both compiled using GCC 13 with -O3 -march=native, examining numbers up to 50,000,000:

CPU time (s)
Intel i7-3770 @ 3.40GHz 66
Intel i7-6700 @ 3.40GHz 54

One attempted optimisation that didn't work was to return a view of amicable pairs rather than constructing a vector. This actually cost around a second of execution time:

constexpr auto amicable_pairs() const
{
    namespace v = std::views;
    return v::iota(1uz, sums.size())
        | v::transform([this](Number i){
            return std::pair{i, sums[i]};
        })
        | v::filter([this](auto p) {
            auto [a, b] = p;
            return  a < b  &&  b < sums.size()
                && sums[b] == a;
        });
}

The Code

#include <algorithm>
#include <cmath>
#include <concepts>
#include <cstdlib>
#include <functional>
#include <limits>
#include <numeric>
#include <ranges>
#include <utility>
#include <vector>

template<std::unsigned_integral Number>
class aliquot_sums
{
    std::vector<Number> sums;
    std::vector<Number> prime_numbers;

    static constexpr inline Number ipow(Number base, unsigned exp)
    {
        // Binary exponentiation
        Number result = 1;
        for (auto m = base;  exp;  m *= m, exp /= 2) {
            if (exp % 2) { result *= m; }
        }
        return result;
    }

    static constexpr inline Number sum_powers(Number p, unsigned i)
    {
        // Each prime factor p with cardinality i contributes
        //    1 + p + p² + ... + pⁱ to the product
        // That simplifies to
        //    (pⁱ⁺¹ - 1) / (p - 1)
        return (ipow(p, i+1) - 1) / (p - 1);
    };

    static constexpr inline bool greater_than_sqrt(Number a, Number b)
    {
        // Test a > √b efficiently without overflow
        static constexpr const Number max_root = std::sqrt(std::numeric_limits<Number>::max());
        if (a <= max_root) [[likely]] {
            return a * a > b;
        } else {
            return a > b / a;
        }
    }


public:
    constexpr aliquot_sums(Number maxval)
        : sums{ 0 },
          prime_numbers{2}
    {
        sums.reserve(maxval);
        // Estimate number of primes using Gauss/Legendre approximation: π(x) ≅ x / ln(x)
        prime_numbers.reserve(static_cast<std::size_t>(std::ceil(static_cast<double>(maxval) / std::log(maxval))));

        for (auto const number: std::views::iota(Number{1}, maxval)) {
            // From the prime factorisation n = 2ᵃ·3ᵇ·5ᶜ·…·pˣ·…,
            // we can reconstruct the factors as all 2ⁱ·3ʲ·5ᵏ·…·pʸ·…
            // where 0 ≤ i ≤ a, 0 ≤ j ≤ b, etc.
            // The sum of these is thus
            // (2⁰+2¹+2²+…+2ᵃ)·(3⁰+3¹+3²+…+3ᵇ)·…·(p⁰+p¹+p²+…+pˣ)·…
            Number aliquot = 1;
            auto n = number;
            for (auto const p: prime_numbers) {
                unsigned count = 0; // cardinality of this prime
                while (n % p == 0) {
                    ++count;
                    n /= p;
                }
                if (count > 0) {
                    aliquot *= sum_powers(p, count);
                };
                if (greater_than_sqrt(p, n)) {
                    break;          // p > √n  ⇒  n is 1 or prime
                }
            }
            if (n > 1) {
                // We ended on a prime; its count is 1
                aliquot *= n + 1;
                if (n > prime_numbers.back()) {
                    // not already seen
                    prime_numbers.push_back(n);
                }
            }
            // We have summed _all_ the factors, so need to subtract
            // the number itself to get the _proper_ factors.
            sums.push_back(aliquot - number);
        }
    }

    constexpr auto const& primes() const
    {
        return prime_numbers;
    }

    constexpr auto perfect_numbers() const
    {
        std::vector<Number> result;
        for (auto const i: std::views::iota(1u, sums.size())) {
            if (sums[i] == i) {
                result.push_back(i);
            }
        }
        return result;
    }

    constexpr auto amicable_pairs() const
    {
        std::vector<std::pair<Number,Number>> pairs;
        for (auto const a: std::views::iota(1uz, sums.size())) {
            auto const b = sums[a];
            if (a < b && b < sums.size() && sums[b] == a) {
                pairs.emplace_back(a, b);
            }
        }
        return pairs;
    }
};
#include <cstdint>
#include <iostream>
#include <limits>
#include <stdexcept>
#include <string>

int main(int argc, char **argv)
{
    using Number = std::uint_fast32_t;

    Number maxval = 1'000'000;     // default
    switch (argc) {
    case 0:
    case 1:
        break;
    case 2:
        try {
            std::size_t endpos;
            auto n = std::stoul(argv[1], &endpos, 0);
            if (n > std::numeric_limits<Number>::max()) { throw std::out_of_range{argv[1]}; }
            if (argv[1][endpos]) { throw std::invalid_argument{argv[1]}; }
            maxval = static_cast<Number>(n);
        } catch (std::exception& e) {
            std::cerr << "Invalid argument: " << e.what();
            return EXIT_FAILURE;
        }
        break;
    default:
        std::cerr << "Too many arguments";
        return EXIT_FAILURE;
    }

    auto const sums = aliquot_sums{maxval};

    std::cout << "Perfect Numbers:\n";
    for (auto n: sums.perfect_numbers()) {
        std::cout << n << '\n';
    }

    std::cout << "\nAmicable Pairs:\n";
    for (auto [a, b]: sums.amicable_pairs()) {
        std::cout << a << ',' << b << '\n';
    }
}

Review objectives

Is the code clear and maintainable? Are the names, layout and comments all good?

Am I using const and constexpr appropriately?

Have I missed any reasonable optimisation opportunities?

Is there anything else you particularly like or dislike?

\$\endgroup\$

1 Answer 1

2
\$\begingroup\$

Is the code clear and maintainable? Are the names, layout and comments all good?

Mostly yes. I especially like the comments; they clearly explain what methods are used to calculate certain things. You could even include links, but I think everything you wrote is searchable on the Internet.

Some variable names are still very short. Consider for example p and i in sum_powers(). Why not prime and cardinality? You have the comments inside the function explaining what the variables mean, but it's even better if you don't have to. Also consider that in most code editors, if they show a popup with a function's signature, you wouldn't see the explanation, only the variable names.

Am I using const and constexpr appropriately?

In principle yes, but the standard math functions are not constexpr until C++26.

Have I missed any reasonable optimisation opportunities?

Some vectors can be reserve()d better. Also consider that you initialize prime_numbers to {2}, but then in the body of the constructor call reserve(). This will most likely lead to two memory allocations. It might be better to reserve first and then add 2 to prime_numbers.

Is there anything else you particularly like or dislike?

Yes, I'd rather have aliquot_sums not be a class, but rather a function that returns a simple data structure with the results. perfect_numbers() and amicable_pairs() can then be made out-of-class functions. This is better for extensibility.

It would be even better if you had aliquot_sums(maxval) return a range or a generator, so the caller can then use std::ranges::to() to convert it to a container of its own choosing, instead of hardcoding the std::vector. However, since it needs to keep track of sums and prime_numbers while generating numbers, that wouldn't much of an optimization. Maybe instead you can pass it a template parameter for the type of container to use?

Some of your error messages lack a newline. I hope standard library writers will make it possible to use std::println() sooner rather than later.

\$\endgroup\$
3
  • \$\begingroup\$ Thanks for a good review! I had a feeling this one would interest you. I actually started with a simple function, but was then passing an optional pointer to let the caller access the primes by-product. Also, with the class, we can make the result generation more incremental - we could easily pick up where we left off (I was considering that for benchmarking - choose a time limit rather than a numeric limit). \$\endgroup\$ Commented Aug 25, 2023 at 13:18
  • 2
    \$\begingroup\$ So sums and prime_numbers are just like a cache to help produce more perfect, amicable and social numbers. Maybe it should be treated as such then? So have a free function amicable_pairs(state, maxval). You could even make the state a static inline singleton object so it doesn't have to be explicitly passed around. I'm saying all this because I think it's worthwhile to reduce the responsibilities of a class like aliquot_sums, as that will make the code more maintainable and extensible. \$\endgroup\$
    – G. Sliepen
    Commented Aug 25, 2023 at 13:57
  • 1
    \$\begingroup\$ Yes, good suggestion. I'll see how that turns out. Thanks again. \$\endgroup\$ Commented Aug 25, 2023 at 14:09

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.