// 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 findMultipleOfPrimitiveTripel4_OneLoop(int64_t c, int64_t m, TripelPrinter& printer)
{
    int64_t cDivm = c / m;
    if ((cDivm & 1) == 0)
        return;

    int64_t sqrtc = static_cast<int64_t>(std::sqrt(cDivm));

    int64_t v = 1;
    int64_t u = sqrtc;
    if ((u & 1) != 0)
        ++u; // ensure that u + v is not even
    int64_t v2u2_c = (v * v + u * u - cDivm) / 2;

    while (v < u)
    {
        if (v2u2_c == 0)
        {
            if (ggt(u, v) == 1)
            {
                printer.found_uv(u, v, m);
            }
        }

        v2u2_c += v + 1;
        ++v;

        if (v2u2_c < 0)
        {
            v2u2_c += v;
            ++v;
        }
        else
        {
            v2u2_c -= u;
            --u;
        }
    }
}

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


void addPrimfactors(std::vector<int64_t>& divisors, int64_t p, int n)
{
	auto numDivisors = divisors.size();
	int64_t f = p;
	for (int i = 0; i < n; ++i, f *= p)
	{
		for (size_t j = 0; j < numDivisors; ++j)
		{
			divisors.push_back(divisors[j] * f);
		}
	}
}


std::vector<int64_t> findDivisors(int64_t c)
{
	std::vector<int64_t> divisors;
	divisors.reserve(16);
	divisors.push_back(1);

	{
		int n = 0;
		while ((c & 1) == 0)
		{
			c = c / 2;
			++n;
		}
		addPrimfactors(divisors, 2, n);
	}

	int64_t csqrt = static_cast<int64_t>(std::sqrt(c));
	for (int64_t p = 3; p <= csqrt; p += 2)
	{
		if (c % p != 0)
			continue;
		int n = 0;
		do
		{
			c /= p;
			++n;
		} while (c % p == 0);
		addPrimfactors(divisors, p, n);

		csqrt = static_cast<int64_t>(std::sqrt(c));
	}

	if (c != 1)
		addPrimfactors(divisors, c, 1);

	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);
}