// 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 findMultipleOfPrimitiveTripel1_OneLoop(int64_t c, int64_t m, TripelPrinter& printer)
{
    int64_t sqrtc = static_cast<int64_t>(std::sqrt(c / m));

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

void findTripelEuklid(int64_t c)
{
	TripelPrinter printer(c);
	
    auto divisors = findDivisors(c);
    for (int64_t m : divisors)
    {
        findMultipleOfPrimitiveTripel1_OneLoop(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);
}