/** To be compiled with
* g++ -O2 -fopenmp -o tsumcube tsumcube.cxx -lgmpxx -lgmp
* For the usage see main().
*/

#include <unistd.h>

#include <iostream>
/* gnu multiprecision library of https://gmplib.org/
* used to deal with products and arithmetics of numbers that
* may be larger than the typical long long int.
* Use s.th. like "zypper in gnump-devel" in openSUSE or equivalent
* in Linuxes to intall the header files.
*/
#include <gmpxx.h>

/*!*************
* @brief Generating a table of c such that the perfect cubes c^3
*  can be written as sums of distinct tetrahedral numbers.

* Overview of the algorithm:
* @author Richard J. Mathar
* @since 2023-03-20
*/
class TSumCube {
public:
	/** modulus for the lookup tables
	*/
	int m;

	/** a (6m)*(6m) table with -1 entries if
	* T(i)+T(j) cannot be a perfect cube based on
	* using mod m arithmetic for the product and the squares
	*/
	int ** tablmod ;

	TSumCube(int mod=32) ;
	~TSumCube() ;

	void run(mpz_class maxz, bool verb,int nThrd) const ;
	static mpz_class icqrt(mpz_class n) ;
	static mpz_class Tetr(int n) ;
	static mpz_class invTetr(mpz_class n) ;
protected:
private:
} ;

#include <time.h>
#include <fstream>
#include <vector>
#include <omp.h>

using namespace std ;

/*!**********************************
* @brief ctor given the modulus for the lookup tables
* @param m modulus >=2.
*/
TSumCube::TSumCube(int mod) : m(mod)
{
	/* the values of n^3 (mod m)
	*/
	vector<int> sqrsmodm ;
	for(int n=0 ; n < m ; n++)
	{
		int res = (n*n*n) % m ;
		bool known = false ;
		for(int i=0 ; i < sqrsmodm.size() && ! known; i++)
		{
			if ( sqrsmodm[i] == res)
				known = true ;
		}
		if ( ! known)
			sqrsmodm.push_back(res) ;
	}

#if 0
	/* debugging: showing how sparse the positive values
	* of the moduli are in the square array, i.e., how efficient
	* the modulus m is.
	*/
	cout << " res mod " << m << endl ;
	for(int i=0 ; i < sqrsmodm.size() ; i++)
		cout << " " << sqrsmodm[i] ;
	cout << endl ;
#endif

	tablmod = new int*[6*m] ;
	for(int i=0 ; i < 6*m ; i++)
		tablmod[i] = new int[6*m] ;

	for(int i=0 ; i < 6*m ; i++)
	{
		for(int j=0 ; j < 6*m ; j++)
		{
			int res = ( (i*(i+1)*(i+2)+j*(j+1)*(j+2))/6 ) % m;
			bool known = false ;
			for(int k=0 ; k < sqrsmodm.size() && ! known; k++)
			{
				if ( sqrsmodm[k] == res)
					known = true ;
			}
			tablmod[i][j] = (known ? res : -1 ) ;
			//cout << " " << tablmod[i][j] ;
		}
		//cout << endl ;
	}
} /* ctor */



/*!**********************************
* @brief dtor 
*  De-allocate the square array of the moduli created in the ctor.
*/
TSumCube::~TSumCube()
{
	for(int i=0 ; i < 6*m ; i++)
		delete tablmod[i] ;
	delete [] tablmod ;
} /* dtor */

/*!*********************
* @brief write snapshots of b-files to b054208.txt.
* @param maxc maximum number tabulated for T(i)*T(j)=c^3, c <= maxc
* @param verb If true known values (not generated necessarily in order)
*   are printed to stdout to give some earlier idea that the
*   values are generated.
*/
void TSumCube::run(const mpz_class maxc, const bool verb,const int nThrd) const
{
	/* sorted list of all known entries so far
	*/
	vector<mpz_class> tabc ;

	/* include 0 as a product of T(0) with any other tetrahedral number.
	tabc.push_back(0) ;
	*/

	const mpz_class maxc3 = maxc*maxc*maxc ;

	/* T(i)=1 starting at i=1 */
	mpz_class Ti(1) ;

	int imodm = 1 ;
	for(long i=1 ; ; i++)
	{
		if ( nThrd <= 1)
		{
			// mpz_class Tj = Ti+(i+1)*(i+2)/2 ;
			mpz_class Tj = Ti ;
			/* innerloop over j, T(j), considering j>=i
			*/
			int jmodm = imodm ;
			for(long j=i ; ; j++)
			{
				/* candidate for a perfect cube is Ti+Tj
				*/
				const mpz_class Tij = Ti+Tj ;
	
				if ( Ti+Tj > maxc3)
				{
					//cout << "i " << i << " j " << j << " Ti+Tj " << Tij << " exc " << maxc3 << endl ;
					break ;
				}
	
				/* speedup by checking that the product might
				* produce a perfect squaure
				*/
				//if (tablmod[i% (6*m)][j% (6*m)] >= 0 )
				if (tablmod[imodm][jmodm] >= 0 )
				{
					const mpz_class c2 = icqrt(Tij) ;
					/* check that this is a perfect cube; attach to results, keep sorted
					*/
					if ( c2*c2*c2 == Tij)
					{
						if ( tabc.empty() )
						{
							tabc.push_back(c2) ;
							if ( verb) 
								cout << i << " " << j << " " << Ti << " " << Tj << " " << c2 << endl ;
						}
						else
						{
							vector<mpz_class>::iterator k = tabc.begin() ;
							bool known = false ;
							for( ; k != tabc.end() && !known ; k++)
							{
								if ( *k == c2)
								{
									known = true ;
									break ;
								}
								else if ( *k > c2)
								{
									if ( verb) 
										cout << i << " " << j << " " << Ti << " " << Tj << " " << c2 << endl ;
									tabc.insert(k,c2) ;
									known = true ;
								}
							}
							if ( ! known)
							{
								tabc.push_back(c2) ;
								if ( verb) 
									cout << i << " " << j << " " << Ti << " " << Tj << " " << c2 << endl ;
							}
						}
					}
				}
				/* faster than T(j) = j*(j+1)*(j+2)/6 by just updating Tj
				Tj += (j+1)*(j+2)/2 ;
				*/
				Tj = Tetr(j+1) ;
				jmodm = (jmodm+1)%(6*m) ;
			}
		}
		else
		{
			/* parallelized variant , nThrd >=2
			* upper j index such that Ti+Tj <= maxc3
			* Tj <= maxc3-Ti
			*/
			// mpz_class maxj = invTetr(maxc3-Ti) ;
			const int maxj = invTetr(maxc3-Ti).get_si() ;
#pragma omp parallel for shared(tabc)
			for(int j=i ; j < maxj; j++)
			{
				const mpz_class Tj = Tetr(j) ;
				/* candidate for a perfect cube is Ti+Tj
				*/
				const mpz_class Tij = Ti+Tj ;
	
				/* speedup by checking that the product might
				* produce a perfect squaure
				*/
				if (tablmod[i% (6*m)][j% (6*m)] >= 0 )
				{
					const mpz_class c2 = icqrt(Tij) ;
					/* check that this is a perfect cube; attach to results, keep sorted
					*/
					if ( c2*c2*c2 == Tij)
					{
#pragma omp cricital
						if ( tabc.empty() )
						{
							tabc.push_back(c2) ;
							if ( verb) 
								cout << i << " " << j << " " << Ti << " " << Tj << " " << c2 << endl ;
						}
						else
						{
							vector<mpz_class>::iterator k = tabc.begin() ;
							bool known = false ;
							for( ; k != tabc.end() && !known ; k++)
							{
								if ( *k == c2)
								{
									known = true ;
									break ;
								}
								else if ( *k > c2)
								{
									if ( verb) 
										cout << i << " " << j << " " << Ti << " " << Tj << " " << c2 << endl ;
									tabc.insert(k,c2) ;
									known = true ;
								}
							}
							if ( ! known)
							{
								tabc.push_back(c2) ;
								if ( verb) 
									cout << i << " " << j << " " << Ti << " " << Tj << " " << c2 << endl ;
							}
						}
					}
				}
			} /* end j-loop parallized */
		}

		/* write an intermediate version into b054208.txt
		*/
		time_t now = time(NULL) ;
		ofstream bf("b054208.txt") ;
		bf << "# >= " << tabc.size() << " values <= " << maxc << " " << ctime(&now)  ;
		for(int k=0 ; k < tabc.size() ; k++)
		{
			/* trusted correct order: we have j>=i so all
			* values c^3 <= 2*Ti are complete in the list.
			* Add question mark if there may be missing smaller values.
			*/
			bf << (k+1) << " "  << tabc[k]  ;
			if ( tabc[k]*tabc[k]*tabc[k] > 2*Ti)
				bf << " ?" ;
			bf << endl ;
		}
		bf.flush() ;
		bf.close() ;

		/* T(i+1) = (i+1)*(i+2)*(i+3)/6 = T(i)+(i+1)(i+2)/2
		Ti += (i+1)*(i+2)/2 ;
		*/
		Ti = Tetr(i+1) ;
		imodm = (imodm+1) %(6*m) ;

		/* since j>i and T(j) > T(i) and we're increaseing T(j)
		* in the outer loop, all results scanned if 2*T(i) > maxc in
		* the next loop.
		*/
		if ( 2*Ti > maxc3 )
		{
			//cout << "Ti " << Ti << " exc " << maxc3 << endl ;
			break ;
		}
	} /* outer loop over increasing i, T(i) */
} /* run */


/*!**********************************
* @brief Integer square root of a positive number n.
* @param n The argument of the function
* @return The largest integer whose square is not larger than n.
*  Concrete: If the return value is x, then x^2 <=n and (x+1)^2>n.
*  -1 if n was smaller than zero and the result is undefined.
*/
mpz_class TSumCube::icqrt(mpz_class n) 
{
	mpz_srcptr nagain ;
	nagain = n.get_mpz_t() ;
	mpz_t result ;
	mpz_init (result) ;
	mpz_root(result,nagain,3) ;
	mpz_class x(result) ;
	mpz_clear(result) ;
	return x ;
} /* icqrt */

/*!**********************************
* @brief Inverse lookup of tetrahedral number index.
* @param n The argument of the function
* @return The largest index x such that x*(x+1)*(x+2)/6 <= n.
*/
mpz_class TSumCube::invTetr(const mpz_class n) 
{
	/* solve x*(x+1)*(x+2) <= 6n with a newton method
	*/
	if ( n < 0 )
		return mpz_class(-1) ;

	mpz_class x=TSumCube::icqrt(6*n) ;
	while ( true )
	{
		x -= (x*(x+1)*(x+2)-6*n)/(3*x*x+6*x+2) ;
		mpz_class Tx = x*(x+1)*(x+2)/6 ;
		if ( Tx == n )
			return x ;
		else
		{
			mpz_class Tx1 = (x+1)*(x+2)*(x+3)/6 ;
			if ( Tx1 > n && Tx <= n)
				return x ;
			Tx1 = (x-1)*x*(x+1)/6 ;
			if ( Tx > n && Tx1 <= n)
				return x-1 ;
		}
	}
} /* invTetr */

/*!**********************************
* @brief tetrahedral number n*(n+1)*(n+2)/6 ;
* @param n The argument of the function
* @return The triple product without overflows
*/
mpz_class TSumCube::Tetr(const int n) 
{
	mpz_class x(n) ;
	x *= n+1 ;
	x *= n+2 ;
	x /= 6 ;
	return x ;
} /* Tetr */

using namespace std ;
/**
* Produces b054208.txt by the algorithm in TSumCube.cxx.
* Usage:
*   tsumcube [-b baseformod] [-j] maxentry
* Example:
*   tsumcube -b 7 2000000
* where maxentry in A053208  is 167986 as of 2023-03-20
* If the option -b is not used a default of 5 is inserted.
* @author Richard J. Mathar
* @since 2023-03-20 
*/
int main(int argc, char *argv[])
{
	/* modular sieve for early recognition of
	* impossible indices of triangular numbers in the sum.
	*/
	int modb = 5 ;
	bool verb = false ;

	/* number of threads for the generator
	*/
	int nThrd =1 ;

	char oc;
	while ( (oc=getopt(argc,argv,"vb:j")) != -1)
	{
		switch(oc)
		{
		case 'v' :
			verb = true ;
			break ;
		case 'b' :
                        modb=atoi(optarg) ;
                        break ;
		case 'j' :
                        // nThrd=atoi(optarg) ;
			// we're using opemp so nThrd  is a proxy 
                        nThrd= 2 ;
                        break ;
                case '?' : 
                        cerr << "Invalid command line option " << optopt <<endl ;
                        return 1 ;
                }
        }

	/* layout of the lookout table for
	* the square array which removes some mixes of
	* products T(i) +T(j) right away from the candidates of
	* squares, if j (mod modb) and j (mod modb) indicate that
	* no cube would be created by the sum anyway.
	*/ 
	TSumCube t(modb) ;

	/* convert the string argument and command line argument
	* to a (huge) integer number which supplies the main
	* algorithm with a maximum of integer indices of
	* triangular numbers that need to be considered such that
	* their sum stays below the cube of that argument.
	*/
	mpz_class maxz(argv[argc-1],10) ;

	/* start the scan of products and overwrite b175497.txt.
	*/
	t.run(maxz,verb,nThrd) ;
} /* main */