Skip to content

divisor_table.hpp

SECTIONMath INCLUDEnoya/divisor_table.hpp

Linear-sieve table of divisor count tau(n) and divisor sum sigma(n) for every integer through a limit.

\[ \displaystyle (\tau(n),\sigma(n)) \]

Implementation

View on GitHub

#ifndef NOYA_DIVISOR_TABLE_HPP
#define NOYA_DIVISOR_TABLE_HPP 1

/// @complexity Time: O(n) build and O(1) table lookup.
/// Space: O(n).

#include <cassert>
#include <cstdint>
#include <vector>

namespace noya {

/// @brief Linear-sieve table of divisor count tau(n) and divisor sum sigma(n)
/// for every integer through a limit.
struct divisor_table {
  std::vector<int> primes;
  std::vector<int> least_prime;
  std::vector<int> exponent;
  std::vector<std::uint64_t> prime_power;
  std::vector<std::uint64_t> divisor_count;
  std::vector<std::uint64_t> divisor_sum;

  divisor_table() = default;
  explicit divisor_table(int limit) { build(limit); }

  void build(int limit) {
    assert(limit >= 0);
    primes.clear();
    least_prime.assign(limit + 1, 0);
    exponent.assign(limit + 1, 0);
    prime_power.assign(limit + 1, 1);
    divisor_count.assign(limit + 1, 0);
    divisor_sum.assign(limit + 1, 0);
    if (limit >= 1) {
      divisor_count[1] = divisor_sum[1] = 1;
    }
    for (int value = 2; value <= limit; value++) {
      if (least_prime[value] == 0) {
        least_prime[value] = value;
        exponent[value] = 1;
        prime_power[value] = value;
        divisor_count[value] = 2;
        divisor_sum[value] = value + 1;
        primes.push_back(value);
      }
      for (int prime : primes) {
        if (value > limit / prime) {
          break;
        }
        int product = value * prime;
        least_prime[product] = prime;
        if (value % prime == 0) {
          exponent[product] = exponent[value] + 1;
          prime_power[product] = prime_power[value] * prime;
          int rest = int(value / prime_power[value]);
          divisor_count[product] =
              divisor_count[rest] * std::uint64_t(exponent[product] + 1);
          divisor_sum[product] =
              divisor_sum[rest] *
              ((prime_power[product] * prime - 1) / (prime - 1));
          break;
        }
        exponent[product] = 1;
        prime_power[product] = prime;
        divisor_count[product] = divisor_count[value] * 2;
        divisor_sum[product] = divisor_sum[value] * (prime + 1);
      }
    }
  }
};

} // namespace noya

#endif // NOYA_DIVISOR_TABLE_HPP
#include <cassert>
#include <cstdint>
#include <vector>

/// @complexity Time: O(n) build and O(1) table lookup.
/// Space: O(n).

namespace noya {

/// @brief Linear-sieve table of divisor count tau(n) and divisor sum sigma(n)
/// for every integer through a limit.
struct divisor_table {
  std::vector<int> primes;
  std::vector<int> least_prime;
  std::vector<int> exponent;
  std::vector<std::uint64_t> prime_power;
  std::vector<std::uint64_t> divisor_count;
  std::vector<std::uint64_t> divisor_sum;

  divisor_table() = default;
  explicit divisor_table(int limit) { build(limit); }

  void build(int limit) {
    assert(limit >= 0);
    primes.clear();
    least_prime.assign(limit + 1, 0);
    exponent.assign(limit + 1, 0);
    prime_power.assign(limit + 1, 1);
    divisor_count.assign(limit + 1, 0);
    divisor_sum.assign(limit + 1, 0);
    if (limit >= 1) {
      divisor_count[1] = divisor_sum[1] = 1;
    }
    for (int value = 2; value <= limit; value++) {
      if (least_prime[value] == 0) {
        least_prime[value] = value;
        exponent[value] = 1;
        prime_power[value] = value;
        divisor_count[value] = 2;
        divisor_sum[value] = value + 1;
        primes.push_back(value);
      }
      for (int prime : primes) {
        if (value > limit / prime) {
          break;
        }
        int product = value * prime;
        least_prime[product] = prime;
        if (value % prime == 0) {
          exponent[product] = exponent[value] + 1;
          prime_power[product] = prime_power[value] * prime;
          int rest = int(value / prime_power[value]);
          divisor_count[product] =
              divisor_count[rest] * std::uint64_t(exponent[product] + 1);
          divisor_sum[product] =
              divisor_sum[rest] *
              ((prime_power[product] * prime - 1) / (prime - 1));
          break;
        }
        exponent[product] = 1;
        prime_power[product] = prime;
        divisor_count[product] = divisor_count[value] * 2;
        divisor_sum[product] = divisor_sum[value] * (prime + 1);
      }
    }
  }
};

} // namespace noya