find_square_sum.cpp

#include<iostream>
#include<vector>
#include<tuple>
#include<cmath>

template<typename Integer>
constexpr Integer sqrt_floor(const Integer& n){
    if(n <= 1)return n;
    Integer r(sqrt(n));
    do r = (r & n / r) + (r ^ n / r) / 2; while(r > n / r);
    return r;
}
template<typename Integer>
constexpr bool is_inside(const Integer& x, const Integer& y, const Integer& n){
    return (!x || x <= n / x) && (!y || y <= (n - x * x) / y);
}
template<typename Integer>
constexpr bool is_cross(const Integer& x, const Integer& y, const Integer& a, const Integer& b, const Integer& n){
    return (a * x - b * y) * (a * x - b * y) / (a * a + b * b) >= x * x + y * y - n;
}
template<typename Integer>
constexpr Integer next_in(const Integer& x, const Integer& y, const Integer& a, const Integer& b, const Integer& n){
    return (b * y - a * x - sqrt_floor((a * a + b * b) * n - (a * y + b * x) * (a * y + b * x)) + a * a + b * b - 1) / (a * a + b * b);
}
template<typename Integer>
constexpr Integer next_out(const Integer& x, const Integer& y, const Integer& a, const Integer& b, const Integer& n){
    return (b * y - a * x + sqrt_floor((a * a + b * b) * n - (a * y + b * x) * (a * y + b * x))) / (a * a + b * b);
}

int main(){
    using namespace std;
    using integer_type = unsigned long;

    integer_type N;
    cin >> N;

    using subtree_type = std::tuple<integer_type, integer_type, integer_type, integer_type>;
    vector<subtree_type> stern_brocot_tree{{0, 1, 1, 0}, {1, 0, 0, 0}};

    using point_type = std::pair<integer_type, integer_type>;
    vector<point_type> result{};

    point_type now_point{0, sqrt_floor(N)};

    if(now_point.second * now_point.second != N)result.emplace_back(now_point);

    while(get<0>(stern_brocot_tree.back())){
        auto& [x, y]{now_point};

        // subtree between a/b and c/d
        auto [a, b, c, d]{stern_brocot_tree.back()};
        stern_brocot_tree.pop_back();

        // next vertex of convex hull
        const auto n_upper{next_out(x, y, a, b, N)};
        x += n_upper * a;
        y -= n_upper * b;
        result.emplace_back(now_point);

        // backtrack
        while(!empty(stern_brocot_tree)){
            tie(a, b, c, d) = stern_brocot_tree.back();
            if(is_inside(x + a, y - b, N))break;
            stern_brocot_tree.pop_back();
        }

        if(empty(stern_brocot_tree))break;

        // search forward
        while(y >= d){
            if(y >= b + d && is_inside(x + a + c, y - b - d, N)){
                const auto n_lower{next_out(x + a + c, y - b - d, c, d, N)};
                tie(a, b, c, d) = make_tuple(a + c * (n_lower + 1), b + d * (n_lower + 1), a + c * (n_lower + 2), b + d * (n_lower + 2));
            }else if(is_cross(x + c, y - d, a, b, N)){
                const auto n_lower{next_in(x + c, y - d, a, b, N)};
                if(y < b * n_lower + d || !is_inside(x + a * n_lower + c, y - b * n_lower - d, N))break;
                tie(a, b, c, d) = make_tuple(a * n_lower + c, b * n_lower + d, a * (n_lower - 1) + c, b * (n_lower - 1) + d);
            }else break;
            stern_brocot_tree.emplace_back(a, b, c, d);
        }
    }

    if(now_point.second)result.emplace_back(now_point.first, 0);

    for(const auto& [x, y] : result)if(x * x + y * y == N)cout << x << " " << y << endl;

    return 0;
}