Project Euler 408.

[Locked]

I am stuck on this.

Problem link

Let's call a lattice point (x, y) inadmissible if x, y and x + y are all positive perfect squares.
For example, (9, 16) is inadmissible, while (0, 4), (3, 1) and (9, 4) are not.

Consider a path from point (x1, y1) to point (x2, y2) using only unit steps north or east.
Let's call such a path admissible if none of its intermediate points are inadmissible.

Let P(n) be the number of admissible paths from (0, 0) to (n, n).
It can be verified that P(5) = 252, P(16) = 596994440 and P(1000) mod 1 000 000 007 = 341920854.

Find P(10 000 000) mod 1 000 000 007.

Let us consider the fact that if there are no inadmissible points, then the function to determine the number of total paths is defined as nCr(2n,n). This is especially true because when you consider P(5), there are 0 inadmissible paths and the function nCr(10,5) == 252, which is the answer to the problem. However, consider the other givens, P(16) has 2 inadmissible nodes, both (16,9) and (9,16). However, finding out nCr(32,16) leads us with 601080390, yet the answer is 596994440 for a difference of 4085950 for just two inadmissible nodes. P(1000) has about 26 nodes and the results are even larger. The arbitrary digit would take up quite an amount of space here.

I figured there had to be a different way of doing it so I thought about it and if we set x = a^2, and y = b^2, then x+y is the same as a^2 + b^2 = c^2, so in order for a point to be inadmissible, it has to be a Pythagorean triple, which might improve the runtime of the code by a ton compared to finding the square root of the number to figure out if it's a perfect square. However as of that point, I'm lost. I have a couple of ideas to throw around and linger but I'm unsure of how to calculate N amount of Pythagorean triples that threshold to that limit. I tried applying Euclid's formula in C++ here:

#include <iostream>
#include <vector>
#include <utility>
#include <cmath>

using ullong = unsigned long long;

constexpr ullong gcd(ullong x, ullong y) {
    return (y == 0) ? x : gcd(y,(x%y));
}

template<typename T>
std::vector<std::pair<T,T>> generateTriples(ullong limit) {
    std::vector<std::pair<T,T>> triples;
    ullong a, b;
    for(ullong n = 1; n < limit; ++n) {
        if(n*n + (n+1)*(n+1) > limit)
            break;
        for(ullong m = n+1; m < limit; m += 2) {
            if(n*n + m*m > limit)
                break;
            if(gcd(n,m) > 1)
                continue;
            a = m*m - n*n;
            b = 2*m*n;
            triples.push_back(std::make_pair(a*a,b*b));
        }
    }
    return triples;
}

int main() {
    std::vector<std::pair<unsigned,unsigned>> triples = generateTriples<unsigned>(100);
    for(auto& i : triples)
        std::cout << "a: " << i.first << " b: " << i.second << '\n';
} 

Output

a: 9 b: 16
a: 225 b: 64
a: 1225 b: 144
a: 3969 b: 256
a: 25 b: 144
a: 441 b: 400
a: 2025 b: 784
a: 5929 b: 1296
a: 49 b: 576
a: 3025 b: 2304
a: 81 b: 1600
a: 1089 b: 3136
a: 4225 b: 5184
a: 121 b: 3600
a: 1521 b: 6400
a: 169 b: 7056

However, I'm still kind of unsure where else to go about with this. What do? :(

My questions:
Is there a better way to solve it?
How do I generate Pythagorean Triples with an arbitrary limit of N?

I will post updates I guess.

Thank you for your help.

Update 1:
I changed the Pythagorean Triple to:

template<typename T>
std::vector<std::pair<T,T>> generateTriples(ullong limit, ullong c) {
    std::vector<std::pair<T,T>> triples;
    ullong a, b;
    for(ullong i = 1; i < c; ++i) {
        for(ullong j = 1; j < i; ++j) {
            a = i*i - j*j;
            b = 2*i*j;
            if(a*a <= limit && b*b <= limit) {
                triples.push_back(std::make_pair(a*a,b*b));
                triples.push_back(std::make_pair(b*b,a*a));
            }
        }
    }
    return triples;
} 

It works for both P(5) and P(16) but fails for P(1000) % 1000000007.