// floriankirsch.de/c++/tripel.html

#include <chrono>
#include <cmath>
#include <iostream>
#include <map>
#include <vector>

int64_t ggt(int64_t a, int64_t b); // assumes a > b
std::vector<int64_t> findDivisors(int64_t c);

class TripelPrinter {
    int64_t c;
    std::map<int64_t, int64_t> solutions;
    void print(int64_t a, int64_t b, int64_t c);

public:
	TripelPrinter(int64_t argc);
	~TripelPrinter();
	
    void found(int64_t a, int64_t b);
    void found_uv(int64_t u, int64_t v, int64_t m);
};

void findMultipleOfPrimitiveTripel0_TwoNestedLoops(int64_t c, int64_t m, TripelPrinter& printer)
{
    int64_t sqrtc = static_cast<int64_t>(std::sqrt(c / m));

    for (int64_t v = 1; v < sqrtc; ++v)
    {
        for (int64_t u = v + 1; u <= sqrtc; ++u)
        {
            if (v * v * m + u * u * m == c)
            {
                if (ggt(u, v) == 1 &&
                    (v + u & 1) == 1)
                {
                    printer.found_uv(u, v, m);
                }
            }
        }
    }
}

void findTripelEuklid(int64_t c)
{
	TripelPrinter printer(c);
	
    auto divisors = findDivisors(c);
    for (int64_t m : divisors)
    {
        findMultipleOfPrimitiveTripel0_TwoNestedLoops(c, m, printer);
    }
}

int main()
{
	std::cout << "Pythagorean tripels found in 1 millisecond\n";
	auto start = std::chrono::steady_clock::now();
	for (int64_t c=1; ; ++c)
	{
	    findTripelEuklid(c);
		auto duration = std::chrono::duration_cast<std::chrono::microseconds>(std::chrono::steady_clock::now() - start);
		if (duration.count() > 1000)
		    break;
	}
	return 0;
}


int64_t ggt(int64_t a, int64_t b)
{
    if (b == 0)
        return a;
    else
        return ggt(b, a % b);
}

std::vector<int64_t> findDivisors(int64_t c)
{
    std::vector<int64_t> divisors;

    int64_t mLimit = static_cast<int64_t>(std::sqrt(c));
    for (int64_t m1 = 1; m1 <= mLimit; ++m1)
    {
        if (c % m1 != 0)
            continue;
        int64_t m2 = c / m1;

        divisors.push_back(m1);
        if (m2 != m1)
            divisors.push_back(m2);
    }

    return divisors;
}


TripelPrinter::TripelPrinter(int64_t argc) : c(argc) { }

TripelPrinter::~TripelPrinter()
{
	for (auto s : solutions)
        print(s.first, s.second, c);
}

void TripelPrinter::print(int64_t a, int64_t b, int64_t c)
{
    std::cout << a << "² + " << b << "² = " << c << "²\n";
}

void TripelPrinter::found(int64_t a, int64_t b)
{
    if (a * a + b * b != c * c)
        std::cout << c << ": algorithm error!\n";
    else if (solutions.find(a) == solutions.end())
        solutions.emplace(a, b);
}

void TripelPrinter::found_uv(int64_t u, int64_t v, int64_t m)
{
    int64_t a = u * u * m - v * v * m;
    int64_t b = 2 * u * v * m;
    if (a < b)
        found(a, b);
    else
        found(b, a);
}