Problem 247: Squares Under a Hyperbola

Problem 247 source code

C++

#include <cmath>
#include <cstdint>
#include <functional>
#include <iostream>
#include <queue>
#include <string>

namespace {

using u64 = std::uint64_t;

struct Options {
    int target_left = 3;
    int target_below = 3;
    bool run_checkpoints = true;
};

bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }

    int parsed = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<int>(c - '0');
    }
    value = parsed;
    return true;
}

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_int_after_prefix(arg, "--left=", options.target_left) ||
            parse_int_after_prefix(arg, "--below=", options.target_below)) {
            continue;
        }

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }

    return options.target_left >= 0 && options.target_below >= 0;
}

struct Node {
    int left = 0;
    int below = 0;
    long double x = 1.0L;
    long double y = 0.0L;
    long double side = 0.0L;
};

long double largest_side(const long double x, const long double y) {
    // Solve y + t = 1/(x + t) for t>0.
    return (std::sqrt((x - y) * (x - y) + 4.0L) - (x + y)) / 2.0L;
}

u64 binom(const int n, const int k) {
    if (k < 0 || k > n) {
        return 0;
    }
    if (k == 0 || k == n) {
        return 1;
    }
    u64 numer = 1;
    u64 denom = 1;
    const int kk = std::min(k, n - k);
    for (int i = 1; i <= kk; ++i) {
        numer *= static_cast<u64>(n - kk + i);
        denom *= static_cast<u64>(i);
    }
    return numer / denom;
}

u64 solve(const int target_left, const int target_below) {
    const u64 target_occurrences = binom(target_left + target_below, target_left);

    auto cmp = [](const Node& a, const Node& b) {
        return a.side < b.side;
    };
    std::priority_queue<Node, std::vector<Node>, decltype(cmp)> pq(cmp);

    Node root;
    root.left = 0;
    root.below = 0;
    root.x = 1.0L;
    root.y = 0.0L;
    root.side = largest_side(root.x, root.y);
    pq.push(root);

    u64 index = 0;
    u64 seen = 0;
    u64 answer = 0;

    while (seen < target_occurrences) {
        const Node cur = pq.top();
        pq.pop();

        ++index;
        if (cur.left == target_left && cur.below == target_below) {
            ++seen;
            answer = index;
        }

        Node right;
        right.left = cur.left + 1;
        right.below = cur.below;
        right.x = cur.x + cur.side;
        right.y = cur.y;
        right.side = largest_side(right.x, right.y);
        pq.push(right);

        Node up;
        up.left = cur.left;
        up.below = cur.below + 1;
        up.x = cur.x;
        up.y = cur.y + cur.side;
        up.side = largest_side(up.x, up.y);
        pq.push(up);
    }

    return answer;
}

bool run_checkpoints() {
    if (solve(1, 1) != 50ULL) {
        std::cerr << "Checkpoint failed for index (1,1)" << '\n';
        return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }

    if (options.run_checkpoints && !run_checkpoints()) {
        return 2;
    }

    std::cout << solve(options.target_left, options.target_below) << '\n';
    return 0;
}

Python

import math
import heapq

def largest_side(x, y):
    # Solve y + t = 1/(x + t) for t > 0
    return (math.sqrt((x - y)**2 + 4.0) - (x + y)) / 2.0

def binom(n, k):
    if k < 0 or k > n: return 0
    if k == 0 or k == n: return 1
    numer = 1
    denom = 1
    kk = min(k, n - k)
    for i in range(1, kk + 1):
        numer *= (n - kk + i)
        denom *= i
    return numer // denom

class Node:
    def __init__(self, left, below, x, y, side):
        self.left = left
        self.below = below
        self.x = x
        self.y = y
        self.side = side
        
    def __lt__(self, other):
        # We want a max-heap based on side, so return True if self.side > other.side
        return self.side > other.side

def solve(target_left=3, target_below=3):
    target_occurrences = binom(target_left + target_below, target_left)
    
    pq = []
    
    root_x = 1.0
    root_y = 0.0
    root_side = largest_side(root_x, root_y)
    root = Node(0, 0, root_x, root_y, root_side)
    
    heapq.heappush(pq, root)
    
    index = 0
    seen = 0
    answer = 0
    
    while seen < target_occurrences:
        cur = heapq.heappop(pq)
        index += 1
        
        if cur.left == target_left and cur.below == target_below:
            seen += 1
            answer = index
            
        right_x = cur.x + cur.side
        right_y = cur.y
        right_side = largest_side(right_x, right_y)
        right = Node(cur.left + 1, cur.below, right_x, right_y, right_side)
        heapq.heappush(pq, right)
        
        up_x = cur.x
        up_y = cur.y + cur.side
        up_side = largest_side(up_x, up_y)
        up = Node(cur.left, cur.below + 1, up_x, up_y, up_side)
        heapq.heappush(pq, up)
        
    return str(answer)

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

Java

import java.util.PriorityQueue;

public class Euler247 {
    static class Node implements Comparable<Node> {
        int left, below;
        double x, y, side;

        Node(int left, int below, double x, double y, double side) {
            this.left = left;
            this.below = below;
            this.x = x;
            this.y = y;
            this.side = side;
        }

        @Override
        public int compareTo(Node other) {
            // Max-heap based on side lengths mapping to highest priority
            return Double.compare(other.side, this.side);
        }
    }

    static double largestSide(double x, double y) {
        return (Math.sqrt((x - y) * (x - y) + 4.0) - (x + y)) / 2.0;
    }

    static long binom(int n, int k) {
        if (k < 0 || k > n)
            return 0;
        if (k == 0 || k == n)
            return 1;
        long numer = 1;
        long denom = 1;
        int kk = Math.min(k, n - k);
        for (int i = 1; i <= kk; ++i) {
            numer *= (n - kk + i);
            denom *= i;
        }
        return numer / denom;
    }

    public static String solve() {
        int targetLeft = 3;
        int targetBelow = 3;

        long targetOccurrences = binom(targetLeft + targetBelow, targetLeft);

        PriorityQueue<Node> pq = new PriorityQueue<>();

        double rootX = 1.0;
        double rootY = 0.0;
        double rootSide = largestSide(rootX, rootY);
        pq.add(new Node(0, 0, rootX, rootY, rootSide));

        long index = 0;
        long seen = 0;
        long answer = 0;

        while (seen < targetOccurrences) {
            Node cur = pq.poll();
            index++;

            if (cur.left == targetLeft && cur.below == targetBelow) {
                seen++;
                answer = index;
            }

            double rightX = cur.x + cur.side;
            double rightY = cur.y;
            double rightSide = largestSide(rightX, rightY);
            pq.add(new Node(cur.left + 1, cur.below, rightX, rightY, rightSide));

            double upX = cur.x;
            double upY = cur.y + cur.side;
            double upSide = largestSide(upX, upY);
            pq.add(new Node(cur.left, cur.below + 1, upX, upY, upSide));
        }

        return String.valueOf(answer);
    }

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