Problem 738: Counting Ordered Factorisations

Problem 738 source code

C++

#include <cassert>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <unordered_map>
#include <vector>

namespace {

using u64 = std::uint64_t;
using u128 = __uint128_t;

constexpr u64 kMod = 1'000'000'007ULL;

bool power_leq(const u64 base, const int exp, const u64 limit) {
    u128 p = 1;
    for (int i = 0; i < exp; ++i) {
        p *= static_cast<u128>(base);
        if (p > static_cast<u128>(limit)) {
            return false;
        }
    }
    return true;
}

u64 int_root(const u64 n, const int k) {
    long double guess = std::powl(static_cast<long double>(n), 1.0L / static_cast<long double>(k));
    u64 r = static_cast<u64>(guess);
    if (r < 1ULL) {
        r = 1ULL;
    }
    while (power_leq(r + 1ULL, k, n)) {
        ++r;
    }
    while (!power_leq(r, k, n)) {
        --r;
    }
    return r;
}

struct Key {
    u64 limit;
    u64 min_factor;
    std::uint16_t parts;
    bool operator==(const Key& o) const {
        return limit == o.limit && min_factor == o.min_factor && parts == o.parts;
    }
};

struct KeyHash {
    std::size_t operator()(const Key& k) const {
        std::size_t h1 = std::hash<u64>{}(k.limit);
        std::size_t h2 = std::hash<u64>{}(k.min_factor);
        std::size_t h3 = std::hash<std::uint16_t>{}(k.parts);
        return h1 ^ (h2 + 0x9e3779b97f4a7c15ULL + (h1 << 6U) + (h1 >> 2U)) ^
               (h3 + 0x9e3779b97f4a7c15ULL + (h2 << 6U) + (h2 >> 2U));
    }
};

u64 count_nondecreasing(const u64 limit,
                        const u64 min_factor,
                        const int parts,
                        std::unordered_map<Key, u64, KeyHash>& memo) {
    if (parts == 0) {
        return 1ULL;
    }
    if (parts == 1) {
        if (limit < min_factor) {
            return 0ULL;
        }
        return limit - min_factor + 1ULL;
    }

    const Key key{limit, min_factor, static_cast<std::uint16_t>(parts)};
    const auto it = memo.find(key);
    if (it != memo.end()) {
        return it->second;
    }

    const u64 r = int_root(limit, parts);
    if (r < min_factor) {
        memo.emplace(key, 0ULL);
        return 0ULL;
    }

    u64 total = 0ULL;
    for (u64 x = min_factor; x <= r; ++x) {
        total += count_nondecreasing(limit / x, x, parts - 1, memo);
    }
    memo.emplace(key, total);
    return total;
}

u64 D(const u64 N, const u64 K) {
    int max_parts = 0;
    for (u64 p = 1ULL; p <= N; p <<= 1U) {
        ++max_parts;
    }
    --max_parts;

    std::vector<u64> counts(static_cast<std::size_t>(max_parts + 1), 0ULL);
    counts[0] = 1ULL;
    std::unordered_map<Key, u64, KeyHash> memo;
    memo.reserve(1 << 20);
    for (int m = 1; m <= max_parts; ++m) {
        counts[static_cast<std::size_t>(m)] = count_nondecreasing(N, 2ULL, m, memo);
    }

    u64 answer = K % kMod;
    for (int m = 1; m <= max_parts; ++m) {
        if (static_cast<u64>(m) > K) {
            break;
        }
        const u64 ways = counts[static_cast<std::size_t>(m)] % kMod;
        const u64 multiplicity = (K - static_cast<u64>(m) + 1ULL) % kMod;
        answer = (answer + static_cast<u64>((static_cast<u128>(ways) * multiplicity) % kMod)) % kMod;
    }

    return answer;
}

}  // namespace

int main() {
    assert(D(10ULL, 10ULL) == 153ULL);
    assert(D(100ULL, 100ULL) == 35'384ULL);

    std::cout << D(10'000'000'000ULL, 10'000'000'000ULL) << '\n';
    return 0;
}

Python

import math

def solve():
    MOD = 1000000007
    N = K = 10000000000

    def int_root(n, k):
        r = max(1, int(n ** (1.0/k)))
        while True:
            try:
                if (r+1)**k <= n: r += 1
                else: break
            except: break
        while r**k > n: r -= 1
        return r

    max_parts = 0
    p = 1
    while p <= N: max_parts += 1; p <<= 1
    max_parts -= 1

    memo = {}
    def count_nd(limit, min_f, parts):
        if parts == 0: return 1
        if parts == 1:
            return max(0, limit - min_f + 1) if limit >= min_f else 0
        key = (limit, min_f, parts)
        if key in memo: return memo[key]
        r = int_root(limit, parts)
        if r < min_f: memo[key] = 0; return 0
        total = 0
        for x in range(min_f, r + 1):
            total += count_nd(limit // x, x, parts - 1)
        memo[key] = total; return total

    counts = [0] * (max_parts + 1); counts[0] = 1
    for m in range(1, max_parts + 1):
        counts[m] = count_nd(N, 2, m)

    answer = K % MOD
    for m in range(1, max_parts + 1):
        if m > K: break
        ways = counts[m] % MOD
        mult = (K - m + 1) % MOD
        answer = (answer + ways * mult) % MOD

    return str(answer)

if __name__ == '__main__':
    print(solve())

Java

import java.util.HashMap;
import java.util.Map;
import java.util.Objects;

public class Euler738 {
    static final long kMod = 1000000007L;

    static boolean powerLeq(long base, int exp, long limit) {
        long p = 1;
        for (int i = 0; i < exp; ++i) {
            if (Long.MAX_VALUE / base < p)
                return false;
            p *= base;
            if (p > limit)
                return false;
        }
        return true;
    }

    static long intRoot(long n, int k) {
        if (n == 0)
            return 0;
        double guess = Math.pow((double) n, 1.0 / k);
        long r = (long) guess;
        if (r < 1L)
            r = 1L;

        while (powerLeq(r + 1, k, n)) {
            r++;
        }
        while (!powerLeq(r, k, n)) {
            r--;
        }
        return r;
    }

    static class Key {
        long limit, minFactor;
        int parts;

        Key(long limit, long minFactor, int parts) {
            this.limit = limit;
            this.minFactor = minFactor;
            this.parts = parts;
        }

        @Override
        public boolean equals(Object o) {
            if (this == o)
                return true;
            if (o == null || getClass() != o.getClass())
                return false;
            Key key = (Key) o;
            return limit == key.limit && minFactor == key.minFactor && parts == key.parts;
        }

        @Override
        public int hashCode() {
            return Objects.hash(limit, minFactor, parts);
        }
    }

    static long countNondecreasing(long limit, long minFactor, int parts, Map<Key, Long> memo) {
        if (parts == 0)
            return 1L;
        if (parts == 1) {
            if (limit < minFactor)
                return 0L;
            return limit - minFactor + 1L;
        }

        Key key = new Key(limit, minFactor, parts);
        Long cached = memo.get(key);
        if (cached != null)
            return cached;

        long r = intRoot(limit, parts);
        if (r < minFactor) {
            memo.put(key, 0L);
            return 0L;
        }

        long total = 0L;
        for (long x = minFactor; x <= r; ++x) {
            total += countNondecreasing(limit / x, x, parts - 1, memo);
        }

        memo.put(key, total);
        return total;
    }

    public static String solve() {
        long N = 10000000000L;
        long K = 10000000000L;

        int maxParts = 0;
        for (long p = 1; p <= N; p <<= 1) {
            maxParts++;
        }
        maxParts--;

        long[] counts = new long[maxParts + 1];
        counts[0] = 1L;

        Map<Key, Long> memo = new HashMap<>(1 << 20);

        for (int m = 1; m <= maxParts; ++m) {
            counts[m] = countNondecreasing(N, 2L, m, memo);
        }

        long answer = K % kMod;
        for (int m = 1; m <= maxParts; ++m) {
            if (m > K)
                break;
            long ways = counts[m] % kMod;
            long multiplicity = (K - m + 1L) % kMod;
            answer = (answer + ways * multiplicity) % kMod;
        }

        return Long.toString(answer);
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}