Problem 563: Robot Welders

Problem 563 source code

C++

#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <limits>
#include <string>
#include <unordered_map>
#include <vector>

namespace {

using u64 = std::uint64_t;
using u128 = unsigned __int128;

static std::string to_string_u128(u128 value) {
    if (value == 0) return "0";
    std::string s;
    while (value > 0) {
        const int digit = static_cast<int>(value % 10);
        s.push_back(static_cast<char>('0' + digit));
        value /= 10;
    }
    std::reverse(s.begin(), s.end());
    return s;
}

static u64 isqrt_u64(u64 x) {
    u64 r = static_cast<u64>(std::sqrt(static_cast<long double>(x)));
    while (r * r > x) --r;
    while ((r + 1) > r && (r + 1) * (r + 1) <= x) ++r;
    return r;
}

static void gen_smooth_rec(const std::vector<u64>& primes, std::size_t idx, u64 cur, u64 limit, std::vector<u64>& out) {
    if (idx == primes.size()) {
        out.push_back(cur);
        return;
    }
    const u64 p = primes[idx];
    u64 v = cur;
    while (true) {
        gen_smooth_rec(primes, idx + 1, v, limit, out);
        if (v > limit / p) break;
        v *= p;
    }
}

static std::vector<u64> gen_smooth(u64 limit) {
    // Prime factors allowed by multiplying dimensions by any k<=25: primes <= 23.
    const std::vector<u64> primes = {2, 3, 5, 7, 11, 13, 17, 19, 23};
    std::vector<u64> nums;
    nums.reserve(200000);
    gen_smooth_rec(primes, 0, 1, limit, nums);
    std::sort(nums.begin(), nums.end());
    nums.erase(std::unique(nums.begin(), nums.end()), nums.end());
    return nums;
}

struct MResult {
    std::vector<u64> M;  // M[k] = minimal area with exactly k variants (k up to max_k), or INF if missing.
    u64 max_filled = 0;
};

static MResult compute_M(u64 side_limit, int max_k) {
    const auto smooth = gen_smooth(side_limit);

    // Count, for each area A, how many factor pairs (x<=y) satisfy y*10 <= 11*x and x*y=A.
    // For 23-smooth A, every divisor is 23-smooth, so "manufacturable" matches "has such factor pair".
    std::unordered_map<u64, std::uint16_t> cnt;
    cnt.reserve(smooth.size() * 4);

    for (std::size_t i = 0; i < smooth.size(); ++i) {
        const u64 x = smooth[i];
        const u64 ymax = (11 * x) / 10;  // floor(1.1*x)
        auto it_end = std::upper_bound(smooth.begin() + static_cast<std::ptrdiff_t>(i), smooth.end(), ymax);
        const std::size_t j_end = static_cast<std::size_t>(it_end - smooth.begin());
        for (std::size_t j = i; j < j_end; ++j) {
            const u64 y = smooth[j];
            const u128 area128 = static_cast<u128>(x) * static_cast<u128>(y);
            if (area128 > std::numeric_limits<u64>::max()) continue;
            const u64 area = static_cast<u64>(area128);
            auto& v = cnt[area];
            if (v <= static_cast<std::uint16_t>(max_k)) ++v;  // clamp above max_k+1
        }
    }

    const u64 INF = std::numeric_limits<u64>::max();
    std::vector<u64> M(static_cast<std::size_t>(max_k + 1), INF);
    for (const auto& [area, k] : cnt) {
        if (k < 2 || k > max_k) continue;
        if (area < M[k]) M[k] = area;
    }

    u64 max_filled = 0;
    for (int k = 2; k <= max_k; ++k) {
        if (M[k] != INF && M[k] > max_filled) max_filled = M[k];
    }
    return {std::move(M), max_filled};
}

static u128 solve_sum_M(int max_k) {
    u64 side_limit = 1'000'000;  // start moderately; we will raise this if needed.

    while (true) {
        const auto res = compute_M(side_limit, max_k);
        const u64 INF = std::numeric_limits<u64>::max();

        bool all_found = true;
        for (int k = 2; k <= max_k; ++k) {
            if (res.M[k] == INF) {
                all_found = false;
                break;
            }
        }
        if (!all_found) {
            side_limit *= 10;
            continue;
        }

        // If all M(k) are found, ensure the side limit is large enough that every area <= maxM
        // has all its near-square factor pairs represented.
        const u64 maxM = res.max_filled;
        u64 s = isqrt_u64(maxM);
        if (s * s < maxM) ++s;  // ceil(sqrt(maxM))
        const u64 needed = (11 * s + 9) / 10;  // ceil(1.1*s)
        if (side_limit < needed) {
            side_limit = needed;
            continue;
        }

        // Statement check: M(3)=889200.
        assert(res.M[3] == 889200ULL);

        u128 sum = 0;
        for (int k = 2; k <= max_k; ++k) sum += static_cast<u128>(res.M[k]);
        return sum;
    }
}

}  // namespace

int main() {
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);

    std::cout << to_string_u128(solve_sum_M(100)) << '\n';
    return 0;
}

Python

def solve():
    max_k = 100
    primes = [2,3,5,7,11,13,17,19,23]

    def gen_smooth(limit):
        result = []
        def rec(idx, cur):
            if idx == len(primes):
                result.append(cur); return
            p = primes[idx]; v = cur
            while v <= limit:
                rec(idx+1, v)
                if v > limit // p: break
                v *= p
        rec(0, 1)
        return sorted(set(result))

    side_limit = 1000000
    while True:
        smooth = gen_smooth(side_limit)
        cnt = {}
        for i, x in enumerate(smooth):
            ymax = 11 * x // 10
            for j in range(i, len(smooth)):
                y = smooth[j]
                if y > ymax: break
                area = x * y
                cnt[area] = cnt.get(area, 0) + 1

        INF = float('inf')
        M = [INF] * (max_k + 1)
        for area, k in cnt.items():
            if 2 <= k <= max_k and area < M[k]: M[k] = area

        if all(M[k] < INF for k in range(2, max_k+1)):
            maxM = max(M[k] for k in range(2, max_k+1))
            import math
            s = math.isqrt(maxM)
            if s*s < maxM: s += 1
            needed = (11*s + 9) // 10
            if side_limit >= needed:
                return str(sum(M[k] for k in range(2, max_k+1)))
            side_limit = needed
        else:
            side_limit *= 10

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

Java

import java.util.*;

public class Euler563 {
    static void genSmoothRec(long[] primes, int idx, long cur, long limit, List<Long> out) {
        if (idx == primes.length) {
            out.add(cur);
            return;
        }
        long p = primes[idx];
        long v = cur;
        while (true) {
            genSmoothRec(primes, idx + 1, v, limit, out);
            if (v > limit / p)
                break;
            v *= p;
        }
    }

    static List<Long> genSmooth(long limit) {
        long[] primes = { 2, 3, 5, 7, 11, 13, 17, 19, 23 };
        List<Long> nums = new ArrayList<>(200000);
        genSmoothRec(primes, 0, 1, limit, nums);
        Set<Long> set = new HashSet<>(nums);
        List<Long> uniqueNums = new ArrayList<>(set);
        Collections.sort(uniqueNums);
        return uniqueNums;
    }

    static class MResult {
        long[] M;
        long maxFilled;
    }

    static int upperBound(List<Long> list, int from, long val) {
        int left = from;
        int right = list.size();
        while (left < right) {
            int mid = left + (right - left) / 2;
            if (list.get(mid) > val) {
                right = mid;
            } else {
                left = mid + 1;
            }
        }
        return left;
    }

    static MResult computeM(long sideLimit, int maxK) {
        List<Long> smooth = genSmooth(sideLimit);
        Map<Long, Integer> cnt = new HashMap<>(smooth.size() * 4);

        for (int i = 0; i < smooth.size(); i++) {
            long x = smooth.get(i);
            long ymax = (11 * x) / 10;
            int jEnd = upperBound(smooth, i, ymax);
            for (int j = i; j < jEnd; j++) {
                long y = smooth.get(j);
                long area = x * y;
                Integer currentCount = cnt.get(area);
                if (currentCount == null) {
                    cnt.put(area, 1);
                } else if (currentCount <= maxK) {
                    cnt.put(area, currentCount + 1);
                }
            }
        }

        long INF = Long.MAX_VALUE;
        long[] M = new long[maxK + 1];
        Arrays.fill(M, INF);

        for (Map.Entry<Long, Integer> entry : cnt.entrySet()) {
            long area = entry.getKey();
            int k = entry.getValue();
            if (k >= 2 && k <= maxK) {
                if (area < M[k]) {
                    M[k] = area;
                }
            }
        }

        long maxFilled = 0;
        for (int k = 2; k <= maxK; k++) {
            if (M[k] != INF && M[k] > maxFilled) {
                maxFilled = M[k];
            }
        }

        MResult res = new MResult();
        res.M = M;
        res.maxFilled = maxFilled;
        return res;
    }

    static long isqrt(long x) {
        if (x < 0)
            return 0;
        long r = (long) Math.sqrt(x);
        while (r * r > x)
            r--;
        while ((r + 1) * (r + 1) <= x && r + 1 > 0)
            r++;
        return r;
    }

    static long solveSumM(int maxK) {
        long sideLimit = 1000000;

        while (true) {
            MResult res = computeM(sideLimit, maxK);
            long INF = Long.MAX_VALUE;
            boolean allFound = true;

            for (int k = 2; k <= maxK; k++) {
                if (res.M[k] == INF) {
                    allFound = false;
                    break;
                }
            }

            if (!allFound) {
                sideLimit *= 10;
                continue;
            }

            long maxM = res.maxFilled;
            long s = isqrt(maxM);
            if (s * s < maxM)
                s++;

            long needed = (11 * s + 9) / 10;
            if (sideLimit < needed) {
                sideLimit = needed;
                continue;
            }

            long sum = 0;
            for (int k = 2; k <= maxK; k++) {
                sum += res.M[k];
            }
            return sum;
        }
    }

    public static String solve() {
        return Long.toString(solveSumM(100));
    }

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