			#include <cassert>
			#include <cstdio>
			#include <cmath>
			#include <iostream>
			#include <chrono>

			#include <fstream>
			#include <sstream>
			#include <string>

			#include <mpirxx.h>
			//#include <gmp-impl.h>

			#include <mpi.h>

			#define DEFAULT_PRECISION 512

			#define SAVE_OUTPUT_FILES false
			#define COMPUTE_CONTRACTION true
			#define MESSAGE_AS_DOUBLE false

			#define DEBUG_TRACE false
			#define DEBUG_FULL_CONTRACTION false

			#if !DEBUG_TRACE
			// This optimization allows to compute only the
			// 2n+1 independent components but it biases the
			// trace check since trace will trivially be forced
			// to 0.
			#define EP0_COMPUTE_OPTIM
			#endif

			constexpr mp_bitcnt_t default_prec = DEFAULT_PRECISION;

			int mpi_size, mpi_rank;
			mpf_class* all_ep0_p;
			char* mpi_buffer_attached;

			//// MPI Datatypes ////
			#ifdef __GMP_SHORT_LIMB
			#define MPI_MPF_LIMB MPI_UNSIGNED
			#define MPI_MPF_SIGNED_LIMB MPI_INT
			#else
			#ifdef _LONG_LONG_LIMB
			#define MPI_MPF_LIMB MPI_UNSIGNED_LONG_LONG
			#define MPI_MPF_SIGNED_LIMB MPI_LONG_LONG
			#else
			#define MPI_MPF_LIMB MPI_UNSIGNED_LONG
			#define MPI_MPF_SIGNED_LIMB MPI_LONG
			#endif
			#endif

			#define MPI_MPF_EXP MPI_LONG

			#include <stdint.h>
			#include <limits.h>

			#if SIZE_MAX == UCHAR_MAX
			#define MPI_SIZE_T MPI_UNSIGNED_CHAR
			#elif SIZE_MAX == USHRT_MAX
			#define MPI_SIZE_T MPI_UNSIGNED_SHORT
			#elif SIZE_MAX == UINT_MAX
			#define MPI_SIZE_T MPI_UNSIGNED
			#elif SIZE_MAX == ULONG_MAX
			#define MPI_SIZE_T MPI_UNSIGNED_LONG
			#elif SIZE_MAX == ULLONG_MAX
			#define MPI_SIZE_T MPI_UNSIGNED_LONG_LONG
			#endif

			MPI_Op MPI_MPF_SUM;
			MPI_Datatype MPI_MPF_CLASS;

			//// End MPI Datatypes ////

			struct Chrono {

			std::chrono::time_point<std::chrono::high_resolution_clock> m_start;
			std::chrono::time_point<std::chrono::high_resolution_clock> m_end;
			std::chrono::duration<double, std::milli> m_ms_timing;

			void start() {
			m_start = std::chrono::high_resolution_clock::now();
			}

			void stop() {
			m_end = std::chrono::high_resolution_clock::now();
			m_ms_timing = m_end - m_start;
			}

			double timing() {
			return m_ms_timing.count();
			}

			};

			inline size_t get_ep0_size(size_t p)
			{
			size_t size = (p+1)*(p+2)/2;
			return size;
			}

			inline size_t get_ep0_idx(size_t p1, size_t p2, size_t p3)
			{
			size_t p = p1+p2+p3;
			size_t idx = (p3(2p+3-p3))/2+p2;
			return idx;
			}

			inline void get_ep0_coords(size_t p, size_t idx, size_t* p1, size_t* p2, size_t* p3)
			{
			size_t ep0_size = get_ep0_size(p);
			size_t reverse_idx = ep0_size - idx - 1;
			size_t A = (size_t)std::ceil( (-1.0 + std::sqrt(1 + 8*(reverse_idx+1)))/2.0 );
			*p3 = p + 1 - A;
			p2 = idx - (ep0_size - A(A+1)/2 );
			p1 = p - p3 - *p2;
			}

			inline mpf_class coeff_Np0(size_t p)
			{
			// (2p)!/(p!)
			mpf_class p2fac_over_pfac = 1.0;
			for(size_t fac = 2*p; fac > p; --fac)
			p2fac_over_pfac *= fac;

			// p!
			mpf_class pfac = 1.0;
			for(size_t fac = 2; fac <= p; ++fac)
			pfac *= fac;

			// 1/2^p
			mpf_class one_over_2p = std::pow(2.0, -(int)p);

			// sqrt( (2p)! / 2^p / (p!)^2 )
			mpf_class result = sqrt(one_over_2p * p2fac_over_pfac / pfac);

			return result;
			}

			inline mpf_class arrangement_factor(size_t n, size_t m)
			{
			assert(m <= n/2);

			size_t m2 = n - 2*m;

			// Note: m <= n/2, thus m <= m2
			mpf_class nfac_over_m2fac = 1.0;
			for(size_t fac = n; fac > m2; --fac)
			nfac_over_m2fac *= fac;

			mpf_class mfac = 1.0;
			for(size_t fac = 2; fac <= m; ++fac)
			mfac *= fac;

			mpf_class one_over_2m = std::pow(2.0, -(int)m);

			mpf_class result = one_over_2m * nfac_over_m2fac / mfac;

			return result;
			}

			mpf_class compute_ep0(size_t p1, size_t p2, size_t p3)
			{
			// Trivial case
			if((p1%2) + (p2%2) > 0)
			return 0.0;

			size_t p = p1+p2+p3;
			size_t k1 = p1/2;
			size_t k2 = p2/2;
			size_t k3 = p3/2;
			mpf_class result = 0.0;

			mpf_class sign = (k1+k2) % 2 == 0 ? -1.0 : 1.0;
			for(size_t m3 = 0; m3 <= k3; ++m3)
			{
			sign *= -1.0;

			mpf_class p1fac_over_k1fac = 1.0;
			for(size_t fac = p1; fac > k1; --fac)
			p1fac_over_k1fac *= fac;

			mpf_class p2fac_over_k2fac = 1.0;
			for(size_t fac = p2; fac > k2; --fac)
			p2fac_over_k2fac *= fac;

			assert( (p + p3 -2m3 -1) % 2 == (2p-1) % 2 );
			// NOTE: Here the denominator is greater than numerator
			mpf_class one_over_double_fac_ratio = 1.0;
			for(size_t fac = 2p-1; fac > p + p3 -2m3 -1; fac -= 2 )
			one_over_double_fac_ratio *= fac;

			mpf_class one_over_2k1k2 = std::pow(2.0, -(int)(k1+k2));

			result += sign * one_over_2k1k2 / one_over_double_fac_ratio
			* p1fac_over_k1fac * p2fac_over_k2fac * arrangement_factor(p3, m3);
			}

			result *= coeff_Np0(p);

			return result;
			}

			#define max(a,b) (((a)<(b))?(b):(a))

			void check_trace_of_ep0(size_t p, mpf_class* ep0_p)
			{
			std::cout << "=== Begin Check Trace (p: " << p << ") ===\n";

			if(p < 2)
			{
			std::cout << " - Trace is only valide for tansor of rank p >= 2.\n";
			return;
			}

			mpf_class max_error = 0.0;

			size_t n = p-2;

			size_t n1 = n;
			size_t n2 = 0;
			size_t n3 = 0;
			do {
			size_t idx1 = get_ep0_idx(n1+2,n2,n3);
			size_t idx2 = get_ep0_idx(n1,n2+2,n3);
			size_t idx3 = get_ep0_idx(n1,n2,n3+2);

			mpf_class trace = ep0_p[idx1]+ep0_p[idx2]+ep0_p[idx3];
			//std::printf("Trace: %g (%g, %g, %g)\n", trace, ep0_p[idx1],ep0_p[idx2],ep0_p[idx3]);

			max_error = max(abs(trace), max_error);

			n2 = (n2+1)%(n+1-n3);
			if(n2 == 0)
			++n3;
			n1 = n - n2 - n3;
			} while(n3 <= n);

			std::cout << " - Max Error: " << max_error << "\n";
			}

			inline size_t get_ep0_offset_from_p(size_t p)
			{
			size_t offset = p(p+1)(p+2)/6;
			return offset;
			}

			inline size_t get_ep0_p_from_offseted_idx(size_t idx)
			{
			// NOTE: Exact solution is not precise enough due to
			// fp errors. This is why I choosed to go with Newton's root finding

			mpf_class x = 0.0;
			mpf_class offset_f = -1.0;
			while(abs(x-offset_f)>1e-4) // magic number, have been tested to be enough for p -> 1'000'000
			{
			offset_f = x;
			x = (2xxx + 3xx + 6idx) / (3xx + 6*x + 2);
			//std::cout << "x = " << x << "\n";
			}

			size_t offset = (size_t)(offset_f.get_d());
			return offset;
			}

			inline size_t get_ep0_offset_from_offseted_idx(size_t idx)
			{
			return get_ep0_offset_from_p(get_ep0_p_from_offseted_idx(idx));
			}

			void fill_ep0(size_t p, mpf_class* ep0_p)
			{
			size_t p1 = p;
			size_t p2 = 0;
			size_t p3 = 0;
			do {
			size_t idx = get_ep0_idx(p1,p2,p3);

			#ifdef EP0_COMPUTE_OPTIM
			if(idx <= 2*p) // <=> idx+1 < 2p+1
			#endif
			{
			ep0_p[idx] = compute_ep0(p1, p2, p3);
			}
			#ifdef EP0_COMPUTE_OPTIM
			else
			#else
			if(false)
			#endif
			{
			assert(p3 >= 2);
			size_t idx1 = get_ep0_idx(p1+2,p2,p3-2);
			size_t idx2 = get_ep0_idx(p1,p2+2,p3-2);
			ep0_p[idx] = - ep0_p[idx1] - ep0_p[idx2];
			}

			//if(ep0_p[idx] != 0.0)
			{
			//std::printf("- (%zu, %zu, %zu) is located at pos %zu", p1, p2, p3, idx);
			//std::cout << " -> value is " << ep0_p[idx] << "\n";
			}

			p2 = (p2+1)%(p+1-p3);
			if(p2 == 0)
			++p3;
			p1 = p - p2 - p3;
			} while(p3 <= p);

			}

			void fill_all_ep0(size_t p_max, mpf_class* all_ep0_p)
			{
			for(size_t p = 0; p <= p_max; ++p)
			{
			mpf_class* ep0_p = all_ep0_p + get_ep0_offset_from_p(p);
			fill_ep0(p, ep0_p);
			}
			}

			void check_trace_of_all_ep0(size_t p_max, mpf_class* all_ep0_p)
			{
			for(size_t p = 0; p <= p_max; ++p)
			check_trace_of_ep0(p, all_ep0_p+get_ep0_offset_from_p(p));
			}

			void dump_mpf_to_buffer(mpf_class* value_p, int* prec_p, int* size_p, mp_exp_t* exp_p, mp_limb_t* limbs_p)
			{
			// NOTE: This procedure is only valid beacause mpf is guarantied to be of fixed precision
			*prec_p = value_p->get_mpf_t()->_mp_prec;
			assert(value_p->get_prec() == default_prec);
			*size_p = value_p->get_mpf_t()->_mp_size;
			*exp_p = value_p->get_mpf_t()->_mp_exp;
			for(size_t i = 0; i < (size_t)abs(*size_p); ++i)
			limbs_p[i] = *(value_p->get_mpf_t()->_mp_d+i);
			}

			void undump_mpf_from_buffer(mpf_class* value_p, int* prec_p, int* size_p, mp_exp_t* exp_p, mp_limb_t* limbs_p)
			{
			// NOTE: This procedure is only valid beacause mpf is guarantied to be of fixed precision
			value_p->get_mpf_t()->_mp_prec = *prec_p;
			value_p->get_mpf_t()->_mp_size = *size_p;
			value_p->get_mpf_t()->_mp_exp = *exp_p;
			for(size_t i = 0; i < (size_t)abs(*size_p); ++i)
			*(value_p->get_mpf_t()->_mp_d+i) = limbs_p[i];
			}

			struct mpf_fixed_class {
			// NOTE: This procedure is only valid beacause mpf is guarantied to be of fixed precision
			int size;
			int prec;
			mp_exp_t exp;
			mp_limb_t limbs[DEFAULT_PRECISION+1];
			mpf_fixed_class & operator=(mpf_class const& other) {
			size = other.get_mpf_t()->_mp_size;
			prec = other.get_mpf_t()->_mp_prec;
			exp = other.get_mpf_t()->_mp_exp;
			prec = other.get_mpf_t()->_mp_prec;
			for(size_t i = 0; i < (size_t)abs(size); ++i)
			limbs[i] = *(other.get_mpf_t()->_mp_d+i);
			return *this;
			}

			void dump_to_mpf_class(mpf_class& other) const {
			other.get_mpf_t()->_mp_size = size;
			other.get_mpf_t()->_mp_prec = prec;
			other.get_mpf_t()->_mp_exp = exp;
			for(size_t i = 0; i < (size_t)abs(size); ++i)
			*(other.get_mpf_t()->_mp_d+i) = limbs[i];
			}

			operator mpf_class() const {
			mpf_class out;
			dump_to_mpf_class(out);
			return out;
			}

			mpf_fixed_class& operator+=(mpf_fixed_class const& other) {
			mpf_class me = (mpf_class)(*this);
			mpf_class ot = (mpf_class)(other);
			me += ot;
			(*this) = me;
			return *this;
			}

			mpf_fixed_class operator+(mpf_fixed_class const& other) const {
			mpf_class me = (mpf_class)(*this);
			mpf_class ot = (mpf_class)(other);
			mpf_fixed_class result;
			result = (me+ot);
			return result;
			}

			};

			void mpi_op_func_sum_mpf_fixed(void void_in, void void_inout, int len_p, MPI_Datatype datatype_p)
			{
			mpf_fixed_class in = (mpf_fixed_class)void_in;
			mpf_fixed_class inout = (mpf_fixed_class)void_inout;
			mpf_fixed_class value;

			for(size_t i = 0; i < *len_p; ++i)
			{
			value = (inout) + (in);
			*inout = value;
			++in; ++inout;
			}
			}

			struct mpf_content_buffer {
			size_t size;
			// ---
			int* size_p;
			int* prec_p;
			mp_exp_t* exp_p;
			mp_limb_t* limbs_p;
			};

			void init_mpf_content_buffer(mpf_content_buffer* buffer, size_t size)
			{
			buffer->size = size;
			buffer->size_p = new int[size];
			buffer->prec_p = new int[size];
			buffer->exp_p = new mp_exp_t[size];
			buffer->limbs_p = new mp_limb_t[size*(default_prec + 1)];
			}

			void clear_mpf_content_buffer(mpf_content_buffer* buffer)
			{
			delete[] buffer->limbs_p;
			delete[] buffer->exp_p;
			delete[] buffer->prec_p;
			delete[] buffer->size_p;
			buffer->size = 0;
			}

			void DEBUG_print_mpf_content_buffer(mpf_content_buffer* buffer, bool all_mpi_tasks = false)
			{
			if(all_mpi_tasks \|\| mpi_rank == 0)
			{
			std::cout << "--- Buffer (size:" << buffer->size << ") ---\n";
			if(all_mpi_tasks) std::cout << "Rank: " << mpi_rank << "\n";
			std::cout << "size:\n "; for(size_t i = 0; i < buffer->size; ++i){ std::cout << buffer->size_p[i] << "\t"; } std::cout << "\n";
			std::cout << "prec:\n "; for(size_t i = 0; i < buffer->size; ++i){ std::cout << buffer->prec_p[i] << "\t"; } std::cout << "\n";
			std::cout << "exp:\n "; for(size_t i = 0; i < buffer->size; ++i){ std::cout << buffer->exp_p[i] << "\t"; } std::cout << "\n";
			std::cout << "limbs:\n "; for(size_t i = 0; i < buffer->size*(default_prec + 1); ++i){ std::cout << buffer->limbs_p[i] << "\t"; } std::cout << "\n";
			}
			}

			void dump_array_of_mpf_to_buffer(mpf_class* value_p, mpf_content_buffer* buffer, size_t size)
			{
			for(size_t i = 0; i < size; ++i)
			{
			dump_mpf_to_buffer(value_p+i, buffer->prec_p+i,
			buffer->size_p+i, buffer->exp_p+i,
			buffer->limbs_p + i*(default_prec+1));
			}
			}

			void undump_array_of_mpf_from_buffer(mpf_class* value_p, mpf_content_buffer* buffer, size_t size)
			{
			for(size_t i = 0; i < size; ++i)
			{
			undump_mpf_from_buffer(value_p+i, buffer->prec_p+i,
			buffer->size_p+i, buffer->exp_p+i,
			buffer->limbs_p + i*(default_prec+1));
			}
			}
			// Thanks to Martin Broadhurst (http://www.martinbroadhurst.com/combinations-of-a-multiset.html)
			unsigned int next_multiset_combination(const unsigned int multiset, unsigned int ar,
			size_t n, size_t k)
			{
			bool finished = false;
			bool changed = false;
			size_t i;

			for (i = k - 1; !finished && !changed; i--) {
			if (ar[i] < multiset[i + (n - k)]) {
			/* Find the successor */
			size_t j;
			for (j = 0; multiset[j] <= ar[i]; j++);
			/* Replace this element with it */
			ar[i] = multiset[j];
			if (i < k - 1) {
			/* Make the elements after it the same as this part of the multiset */
			size_t l;
			for (l = i + 1, j = j + 1; l < k; l++, j++) {
			ar[l] = multiset[j];
			}
			}
			changed = true;
			}
			finished = (i == 0);
			}
			if (!changed) {
			/* Reset to first combination */
			for (i = 0; i < k; i++) {
			ar[i] = multiset[i];
			}
			}
			return changed;
			}

			unsigned int* create_multiset(size_t n1, size_t n2, size_t n3, unsigned int* buffer)
			{
			// Buffer is supposed to be big enough!

			unsigned int* multiset = buffer;

			size_t i = 0;
			for(;i < n1; ++i)
			multiset[i] = 1;
			for(;i < n1+n2; ++i)
			multiset[i] = 2;
			for(;i < n1+n2+n3; ++i)
			multiset[i] = 3;

			return multiset;
			}

			unsigned int* create_first_combination_from_multiset(unsigned int* multiset, size_t k, unsigned int* buffer)
			{
			unsigned int* combination = buffer;
			for(size_t i = 0; i < k; ++i)
			combination[i] = multiset[i];
			return combination;
			}

			void read_combination_of_multiset(unsigned int* combination, size_t k,
			size_t p1, size_t p2, size_t p3,
			size_t* L1, size_t* L2, size_t* L3,
			size_t* R1, size_t* R2, size_t* R3)
			{
			size_t i = 0;
			L1 = L2 = *L3 = 0;
			for(;i<k && combination[i]==1;++i){ ++(L1); } R1 = p1 - *L1;
			for(;i<k && combination[i]==2;++i){ ++(L2); } R2 = p2 - *L2;
			for(;i<k && combination[i]==3;++i){ ++(L3); } R3 = p3 - *L3;
			}

			inline mpf_class binomial(size_t n, size_t p)
			{
			size_t p1 = (p > n-p) ? p : n-p;
			size_t p2 = n-p1;

			mpf_class p2fac = 1.0;
			for(size_t fac = 1; fac <= p2; ++fac) p2fac *= fac;

			mpf_class result = 1.0;
			for(size_t fac = n; fac > p1; --fac) result *= fac;

			result /= p2fac;

			return result;
			}

			inline mpf_class trinomial(size_t n1, size_t n2, size_t n3)
			{
			// This ensure that n1_ is the greatest of the three
			size_t n1_, n2_, n3_;
			if(n1 >= n2 && n1 >= n3) { n1_ = n1; n2_ = n2; n3_ = n3; }
			else if(n2 >= n1 && n2 >= n3) { n1_ = n2; n2_ = n1; n3_ = n3; }
			else { n1_ = n3; n2_ = n1; n3_ = n2; }

			mpf_class n2fac = 1.0;
			mpf_class n3fac = 1.0;
			for(size_t fac = 1; fac <= n2_; ++fac) n2fac *= fac;
			for(size_t fac = 1; fac <= n3_; ++fac) n3fac *= fac;

			mpf_class multinomial = 1.0;
			for(size_t fac = n1+n2+n3; fac > n1_; --fac)
			multinomial *= fac;
			multinomial /= (n2fac*n3fac);

			return multinomial;
			}

			mpf_class compute_EP0S(size_t p, size_t n, size_t pp,
			size_t p1, size_t p2, size_t p3)
			{
			assert(p1+p2+p3 == p + pp - 2*n);

			// Compute the symmetric part of (ep0 .n ep'0)
			mpf_class result = 0.0;

			// Sum over k-combinations:
			size_t k = p-n;

			if(k == 0) return 1.0;

			// NOTE(Sam): This allocation is suboptimal but I had issue with the crapy allocator
			// I implemented eralier and it works like that. I'll try to optimize it if needed
			// in the future.
			unsigned int* multiset = new unsigned int[p1+p2+p3];
			unsigned int* combination = new unsigned int[k];
			multiset = create_multiset(p1,p2,p3, multiset);
			combination = create_first_combination_from_multiset(multiset, k, combination);

			/*
			std::cout << "p1+p2+p3="<<p1+p2+p3<<", k="<<k<<"\n";
			for(size_t i = 0; i < p1+p2+p3; ++i)
			std::cout << multiset[i] << " ";
			std::cout << "\n";
			for(size_t i = 0; i < k; ++i)
			std::cout << combination[i] << " ";
			std::cout << "\n";
			std::cout << "---" "\n";
			//*/

			//std::cout << "---\n";
			do {
			/*
			for(size_t i = 0; i < k; ++i)
			std::cout << combination[i] << " ";
			std::cout << "\n";
			//*/

			size_t L1, L2, L3, R1, R2, R3;
			read_combination_of_multiset(combination, k, p1, p2, p3, &L1, &L2, &L3, &R1, &R2, &R3);
			//std::cout << "(" << L1 << ", " << L2 << ", " << L3 << ") - (" << R1 << ", " << R2 << ", " << R3 << ")\n";

			mpf_class binome_1 = binomial(p1, L1);
			mpf_class binome_2 = binomial(p2, L2);
			mpf_class binome_3 = binomial(p3, L3);

			// Sum with n1+n2+n3=n
			size_t n1 = n;
			size_t n2 = 0;
			size_t n3 = 0;
			do {
			mpf_class multinomial = trinomial(n1,n2,n3);

			size_t idx1 = get_ep0_offset_from_p(p) + get_ep0_idx(L1+n1, L2+n2, L3+n3);
			size_t idx2 = get_ep0_offset_from_p(pp) + get_ep0_idx(R1+n1, R2+n2, R3+n3);

			result += binome_1binome_2binome_3multinomialall_ep0_p[idx1]*all_ep0_p[idx2];

			n2 = (n2+1)%(n+1-n3);
			if(n2 == 0)
			++n3;
			n1 = n - n2 - n3;
			} while(n3 <= n);
			// End sum with n1+n2+n3=n

			} while (next_multiset_combination(multiset, combination, p1+p2+p3, k));
			// End sum over k-combinations

			delete[] combination;
			delete[] multiset;

			// Compute prefactor (p-n)!(pp-n)!/(p+p-2n)!
			size_t upper_left, upper_right;
			size_t lower = p+pp-2*n;
			if(p-n > pp-n)
			{
			upper_left = p-n;
			upper_right = pp-n;
			}
			else
			{
			upper_left = pp-n;
			upper_right = p-n;
			}

			mpf_class lower_fac_over_upper_left_fac = 1.0;
			for(size_t fac = lower; fac > upper_left; --fac)
			lower_fac_over_upper_left_fac *= fac;

			mpf_class upper_right_fac = 1.0;
			for(size_t fac = 1; fac <= upper_right; ++fac)
			upper_right_fac *= fac;

			result *= upper_right_fac / lower_fac_over_upper_left_fac;

			return result;
			}

			mpf_class compute_EP0(size_t p, size_t n, size_t pp,
			size_t p1, size_t p2, size_t p3)
			{
			// Compute the traceless and symmetric part of (ep0 .n ep'0)
			mpf_class result = 0.0;

			size_t k1 = p1/2;
			size_t k2 = p2/2;
			size_t k3 = p3/2;

			for(size_t m1 = 0; m1 <= k1; ++m1)
			for(size_t m2 = 0; m2 <= k2; ++m2)
			for(size_t m3 = 0; m3 <= k3; ++m3)
			{
			mpf_class sign = ( (m1+m2+m3)%2 == 0 ? 1.0 : -1.0 );

			mpf_class arr_factor_1 = arrangement_factor(p1, m1);
			mpf_class arr_factor_2 = arrangement_factor(p2, m2);
			mpf_class arr_factor_3 = arrangement_factor(p3, m3);

			mpf_class prefactor = sign * arr_factor_1 * arr_factor_2 * arr_factor_3;

			if( (2(p+pp-2n) - 2(m1+m2+m3)-1) % 2 == (2(p+pp-2*n)-1) % 2)
			{
			mpf_class one_over_double_fac = 1.0;
			for(size_t fac = (2(p+pp-2n)-1); fac > (2(p+pp-2n)-2*(m1+m2+m3)-1); --fac)
			one_over_double_fac *= fac;
			prefactor /= one_over_double_fac;
			}
			else
			{
			mpf_class numerator = 1.0;
			for(size_t fac = 1; fac <= (2(p+pp-2n)-2(m1+m2+m3)-1); ++fac) numerator = fac;

			mpf_class denominator = 1.0;
			for(size_t fac = 1; fac <= (2(p+pp-2n)-1); ++fac) denominator *= fac;

			prefactor *= numerator/denominator;
			}

			// Sum with h1+h2+h3=m
			size_t h1 = m1+m2+m3;
			size_t h2 = 0;
			size_t h3 = 0;
			do {

			mpf_class multinomial = trinomial(h1,h2,h3);

			result += prefactor * multinomial
			* compute_EP0S(p,n,pp, p1 - 2(m1-h1), p2 - 2(m2-h2), p3 - 2*(m3-h3) );

			h2 = (h2+1)%(m1+m2+m3+1-h3);
			if(h2 == 0)
			++h3;
			h1 = m1+m2+m3 - h2 - h3;
			} while(h3 <= m1+m2+m3);
			// End sum with h1+h2+h3=m
			}

			return result;
			}

			/*
			double avg_counter = 0.0;
			size_t nb_counter = 0;
			//*/

			mpf_class compute_full_tensor_contraction(size_t p, size_t n, size_t pp, size_t mpi_start, size_t mpi_count)
			{
			size_t P = p+pp-2*n;
			mpf_class result = 0.0;

			// NOTE(Sam): This is the loop I want to parallelize since it is the one that is executed the most often
			// void get_ep0_coords(size_t p, size_t idx, size_t* p1, size_t* p2, size_t* p3)
			// size_t get_ep0_size(size_t p)

			size_t nb_iterations = get_ep0_size(P);
			size_t nb_iterations_per_task = (size_t)std::ceil((double)nb_iterations / (double)mpi_count);
			//if(mpi_rank == 0) std::cout << " - Each process will compute " << nb_iterations_per_task << " / " << nb_iterations << " iterations. (P=" << P << ")" "\n";

			if((mpi_rank-mpi_start)*nb_iterations_per_task < nb_iterations)
			{
			// Sum with p1+p2+p3=p+pp-2n
			size_t iterations_counter = 0;
			size_t p1 = P;
			size_t p2 = 0;
			size_t p3 = 0;
			get_ep0_coords(P, (mpi_rank-mpi_start)*nb_iterations_per_task, &p1, &p2, &p3);
			//std::cout << "[Rank " << mpi_rank << "] starts from: (" << p1 << ", " << p2 << ", " << p3 << ")\n";
			do {

			++iterations_counter;

			mpf_class multinomial = trinomial(p1,p2,p3);

			size_t idx1 = get_ep0_offset_from_p(P) + get_ep0_idx(p1,p2,p3);

			result += multinomialall_ep0_p[idx1]compute_EP0(p,n,pp, p1,p2,p3);

			p2 = (p2+1)%(P+1-p3);
			if(p2 == 0)
			++p3;
			p1 = P - p2 - p3;
			} while(iterations_counter < nb_iterations_per_task && p3 <= P);

			//std::cout << "[Rank " << mpi_rank << "] finished with " << iterations_counter << "/" << nb_iterations_per_task << "\n";

			// End sum with p1+p2+p3=p+pp-2n
			}
			else
			{
			//std::cout << "[Rank " << mpi_rank << "] Not enough iterations...\n";
			}

			#if MESSAGE_AS_DOUBLE
			double send_buffer = result.get_d();
			double recv_buffer;

			// This is not possible anymore, we would have to create a custom communicator but
			// we would have a lot of different communicator or we would have to create one on the fly
			//MPI_Reduce(&send_buffer, &recv_buffer, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
			//result = recv_buffer;

			MPI_Bsend(&send_buffer, 1, MPI_DOUBLE, mpi_start, 0, MPI_COMM_WORLD);

			if(mpi_rank == mpi_start)
			{
			result = 0.0;
			for(size_t i = mpi_start; i < mpi_start+mpi_count; ++i)
			{
			MPI_Recv(&recv_buffer, 1, MPI_DOUBLE, i, 0,
			MPI_COMM_WORLD, MPI_STATUS_IGNORE);
			result += recv_buffer;
			}
			}

			#else // Not MESSAGE_AS_DOUBLE

			#if 1

			mpf_fixed_class send_buffer, recv_buffer;
			send_buffer = result;

			// This is not possible anymore, we would have to create a custom communicator but
			// we would have a lot of different communicator or we would have to create one on the fly
			//MPI_Reduce(&send_buffer, &recv_buffer, 1, MPI_MPF_CLASS, MPI_MPF_SUM, 0, MPI_COMM_WORLD);
			//if(mpi_rank == 0) result = (mpf_class)(recv_buffer);

			// Note: We attached a buffer that may contain only 1 mpf_static_class, see MPI_custom_init()
			MPI_Bsend(&send_buffer, 1, MPI_MPF_CLASS, mpi_start, 0, MPI_COMM_WORLD);

			if(mpi_rank == mpi_start)
			{
			result = 0.0;
			for(size_t i = mpi_start; i < mpi_start+mpi_count; ++i)
			{
			MPI_Recv(&recv_buffer, 1, MPI_MPF_CLASS, i, 0,
			MPI_COMM_WORLD, MPI_STATUS_IGNORE);
			result += (mpf_class)(recv_buffer);
			}
			}

			#else

			mpf_content_buffer send_buffer, recv_buffer;
			init_mpf_content_buffer(&send_buffer, 1);
			if(mpi_rank == 0) init_mpf_content_buffer(&recv_buffer, mpi_size);

			//std::cout << "[Rank " << mpi_rank << "] (P:" << P << ") result = " << result << "\n";

			dump_array_of_mpf_to_buffer(&result, &send_buffer, 1);

			//DEBUG_print_mpf_content_buffer(&send_buffer);
			//DEBUG_print_mpf_content_buffer(&recv_buffer);

			MPI_Gather(send_buffer.size_p, 1, MPI_INT, recv_buffer.size_p, 1, MPI_INT, 0, MPI_COMM_WORLD);
			MPI_Gather(send_buffer.prec_p, 1, MPI_INT, recv_buffer.prec_p, 1, MPI_INT, 0, MPI_COMM_WORLD);
			MPI_Gather(send_buffer.exp_p, 1, MPI_MPF_EXP, recv_buffer.exp_p, 1, MPI_MPF_EXP, 0, MPI_COMM_WORLD);
			MPI_Gather(send_buffer.limbs_p, (default_prec + 1), MPI_MPF_LIMB, recv_buffer.limbs_p, (default_prec + 1), MPI_MPF_LIMB, 0, MPI_COMM_WORLD);

			//DEBUG_print_mpf_content_buffer(&recv_buffer);

			if(mpi_rank == 0)
			{
			mpf_class* all_results = new mpf_class[mpi_size];
			for(size_t i = 0; i < mpi_size; ++i) all_results[i] = 0.0;

			undump_array_of_mpf_from_buffer(all_results, &recv_buffer, mpi_size);

			result = 0.0;
			for(size_t i = 0; i < mpi_size; ++i)
			{
			result += all_results[i];
			//std::cout << " >>> [" << i << "] " << all_results[i] << "\n";
			}
			//std::cout << "[Rank " << mpi_rank << "] result final = " << result << "\n";

			delete[] all_results;
			}


			if(mpi_rank == 0) clear_mpf_content_buffer(&recv_buffer);
			clear_mpf_content_buffer(&send_buffer);
			#endif

			#endif

			return result;
			}

			void DEBUG_check_full_contraction(size_t p_max);
			void DEBUG_clear_all_ep0(size_t p_max, mpf_class* all_ep0_p);
			void DEBUG_dump_all_ep0(size_t p_max, mpf_class* all_ep0_p);

			void MPI_custom_init()
			{
			/*
			struct mpf_fixed_class {
			int size;
			int prec;
			mp_exp_t exp;
			mp_limb_t limbs[DEFAULT_PRECISION+1];
			};
			*/

			mpf_fixed_class construction_instances[2];

			constexpr size_t nb_fields = 4;
			int block_length[nb_fields] = {1, 1, 1, DEFAULT_PRECISION+1};
			MPI_Aint displacements[nb_fields];
			MPI_Datatype old_types[nb_fields] = {MPI_INT, MPI_INT, MPI_MPF_EXP, MPI_MPF_LIMB};

			MPI_Aint addr0, addr1;
			MPI_Get_address(&construction_instances[0], &addr0);
			MPI_Get_address(&construction_instances[0].size, &displacements[0]);
			MPI_Get_address(&construction_instances[0].prec, &displacements[1]);
			MPI_Get_address(&construction_instances[0].exp, &displacements[2]);
			MPI_Get_address(&construction_instances[0].limbs[0], &displacements[3]);
			for(size_t i = 0; i < nb_fields; ++i)
			displacements[i] = MPI_Aint_diff(displacements[i], addr0);

			MPI_Datatype MPI_MPF_construct;
			MPI_Type_create_struct(nb_fields, block_length, displacements, old_types, &MPI_MPF_construct);

			MPI_Get_address(&construction_instances[1], &addr1);
			addr1 = MPI_Aint_diff(addr1, addr0);

			// MPI_MPF_CLASS is global
			MPI_Type_create_resized(MPI_MPF_construct, 0, addr1, &MPI_MPF_CLASS);
			MPI_Type_commit(&MPI_MPF_CLASS);

			// Sum of mpf_fixed_class for reduction
			MPI_Op_create( mpi_op_func_sum_mpf_fixed, true, &MPI_MPF_SUM );


			// Attach buffer for Bsend
			int buffer_attached_size = MPI_BSEND_OVERHEAD + 10*sizeof(mpf_fixed_class);
			mpi_buffer_attached = new char[buffer_attached_size];
			MPI_Buffer_attach(mpi_buffer_attached, buffer_attached_size);
			}

			void MPI_custom_free()
			{
			delete[] mpi_buffer_attached;
			MPI_Op_free(&MPI_MPF_SUM);
			MPI_Type_free(&MPI_MPF_CLASS);
			}

			int main(int argc, char **argv)
			{
			MPI_Init(&argc, &argv);
			MPI_custom_init();

			MPI_Comm_size(MPI_COMM_WORLD, &mpi_size);
			MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank);


			if(mpi_rank == 0) std::cout << " ------------------------ \n"
			<< " CONFIG \n"
			<< " ------------------------ \n";

			if(mpi_rank == 0)
			{
			std::cout << "MPI: " << mpi_size << " tasks\n";
			#if MESSAGE_AS_DOUBLE
			std::cout << "- [WARNING] MPI message are passing as double which is less precise!" "\n";
			#endif
			}

			mpf_set_default_prec(default_prec);
			if(mpi_rank == 0)
			std::cout << "Default precision: " << mpf_get_default_prec() << "\n";

			size_t target_p_max = 5;
			if(argc <= 1)
			{
			if(mpi_rank == 0) std::cout << "[WARNING] No argument target_p_max, default is " << target_p_max << "\n";
			}
			else
			{
			std::istringstream iss(argv[1]);
			if(! (iss >> target_p_max ) ) {
			if(mpi_rank == 0) std::cout << "[ERROR] Invalid argument \n";
			MPI_Abort(MPI_COMM_WORLD, 1);
			}
			}
			if(mpi_rank == 0) std::cout << "Target p max: " << target_p_max << "\n";

			size_t ideal_nb_iterations_per_task = 10;
			if(argc <= 2)
			{
			//if(mpi_rank == 0) std::cout << "[WARNING] No argument target_p_max, default is " << target_p_max << "\n";
			}
			else
			{
			std::istringstream iss(argv[2]);
			if(! (iss >> ideal_nb_iterations_per_task ) ) {
			if(mpi_rank == 0) std::cout << "[ERROR] Invalid argument \n";
			MPI_Abort(MPI_COMM_WORLD, 1);
			}
			}
			if(mpi_rank == 0) std::cout << "Ideal number of iterations per tasks: " << ideal_nb_iterations_per_task << "\n";

			if(mpi_rank == 0) std::cout << "Do save output? " << ( (SAVE_OUTPUT_FILES) ? "yes" : "no" ) << "\n";

			if(mpi_rank == 0) std::cout << " ------------------------ \n";

			Chrono outer_chrono;
			outer_chrono.start();
			{
			size_t p_max = 2*target_p_max;
			if(mpi_rank == 0) std::cout << "- We need to precompute " << p_max << " tensors\n";
			size_t all_ep0_size = get_ep0_offset_from_p(p_max+1);
			all_ep0_p = new mpf_class[all_ep0_size];
			//std::printf("For all ep0 up to p=%zu we need %zu independent variables.\n", p_max, all_ep0_size);

			//DEBUG_clear_all_ep0(p_max, all_ep0_p);

			fill_all_ep0(p_max, all_ep0_p);

			}
			outer_chrono.stop();

			if(mpi_rank == 0)
			std::cout << "- ep0 precomputation... " << outer_chrono.timing() << "ms" "\n";

			{
			#if DEBUG_TRACE
			size_t p_max = 2*target_p_max;
			check_trace_of_all_ep0(p_max, all_ep0_p);
			#endif
			}

			#if COMPUTE_CONTRACTION
			size_t total_nb_terms = 0;
			{
			size_t p_max = target_p_max;
			for(size_t p = 0; p < p_max; ++p)
			for(size_t pp = 0; pp <= p; ++pp)
			for(size_t n = 0; n <= std::min(p,pp); ++n)
			{
			++total_nb_terms;
			}
			}

			size_t next_available_task = 0;
			size_t prev_p, prev_pp, prev_n;

			outer_chrono.start();
			{
			size_t nb_terms = 0;
			if(mpi_rank == 0) std::cout << "- We need to compute " << total_nb_terms << " terms\n";
			size_t p_max = target_p_max;
			for(size_t p = 0; p < p_max; ++p)
			for(size_t pp = 0; pp <= p; ++pp)
			for(size_t n = 0; n <= std::min(p,pp); ++n)
			{
			//std::cout << "(p,n,pp) = (" << p <<", " << n << ", " << pp << ")\n";

			size_t P = p+pp-2*n;
			size_t inside_nb_iterations = get_ep0_size(P);

			size_t ideal_nb_of_tasks = (size_t)std::ceil((double)inside_nb_iterations / (double)ideal_nb_iterations_per_task);
			size_t nb_assigned_tasks = ideal_nb_of_tasks < mpi_size ? ideal_nb_of_tasks : mpi_size;

			size_t remaining_available_tasks = mpi_size - next_available_task;
			// TODO(Sam): In order to optimize a bit more, we should distribute remaining tasks if enough

			//if(mpi_rank == 0) std::cout << " -> " << inside_nb_iterations << " would require idealy " << ideal_nb_of_tasks << " tasks (" << remaining_available_tasks << " remaining).\n";
			if(nb_assigned_tasks <= remaining_available_tasks)
			{
			// Launch the computation for 'nb_assigned_tasks' starting from 'next_available_task'
			//if(mpi_rank == 0) std::cout << next_available_task << " + " << nb_assigned_tasks << "\n";
			if(mpi_rank >= next_available_task && mpi_rank < next_available_task+nb_assigned_tasks)
			{
			// COMPUTATION...
			{
			#if 1
			Chrono chrono;
			chrono.start();
			mpf_class value = compute_full_tensor_contraction(p,n,pp, next_available_task, nb_assigned_tasks);
			chrono.stop();

			double timing = chrono.timing();

			++nb_terms;

			#if SAVE_OUTPUT_FILES
			if(mpi_rank == next_available_task)
			{
			std::string to_print;
			mp_exp_t exp;
			to_print = "0."+value.get_str(exp);

			std::stringstream ss;
			ss << "output/" << p << "." << n << "." << pp << ".txt";

			std::ofstream output_file(ss.str());
			output_file << "p=" << p << "\n";
			output_file << "n=" << n << "\n";
			output_file << "pp=" << pp << "\n";
			output_file << "value=" << to_print << "e"<< exp << "\n";
			output_file << "timing=" << timing/1000.0 << "\n";
			output_file.close();
			}
			#endif
			#endif
			}

			// Sync everything
			next_available_task = 0;
			MPI_Bcast(&next_available_task, 1, MPI_SIZE_T, mpi_size-1, MPI_COMM_WORLD);
			MPI_Bcast(&p, 1, MPI_SIZE_T, mpi_size-1, MPI_COMM_WORLD);
			MPI_Bcast(&n, 1, MPI_SIZE_T, mpi_size-1, MPI_COMM_WORLD);
			MPI_Bcast(&pp, 1, MPI_SIZE_T, mpi_size-1, MPI_COMM_WORLD);

			//std::cout << " >>> [" << mpi_rank << "] restarting with (" << p << ", " << n << ", " << pp << ")" "\n";
			//if(mpi_rank == mpi_size -1 ) std::cout << "----- resync -----" "\n";
			//else std::cout << " -- [" << mpi_rank << "] --" "\n";
			}
			else
			{
			if(mpi_rank < next_available_task)
			MPI_Abort(MPI_COMM_WORLD, 1);

			//std::cout << " xxx [" << mpi_rank << "] (" << p << ", " << n << ", " << pp << ") -> " << nb_assigned_tasks << " tasks have been assigned.\n";
			next_available_task += nb_assigned_tasks;

			//std::cout << " NEXT[" << mpi_rank << "]: " << next_available_task << "\n";
			}
			}
			else
			{

			// If we are here that means thate there are no more work to
			// be done for the remaining tasks and we should start
			// another batch (bcast 0)
			next_available_task = 0;

			// We have to recompute the current iteration
			p = prev_p;
			pp = prev_pp;
			n = prev_n;

			MPI_Bcast(&next_available_task, 1, MPI_SIZE_T, mpi_size-1, MPI_COMM_WORLD);
			MPI_Bcast(&p, 1, MPI_SIZE_T, mpi_size-1, MPI_COMM_WORLD);
			MPI_Bcast(&n, 1, MPI_SIZE_T, mpi_size-1, MPI_COMM_WORLD);
			MPI_Bcast(&pp, 1, MPI_SIZE_T, mpi_size-1, MPI_COMM_WORLD);
			//std::cout << " >>> [" << mpi_rank << "] re-restarting with (" << p << ", " << n << ", " << pp << ")" "\n";
			//if(mpi_rank == mpi_size -1 ) std::cout << "----- resync -----" "\n";
			//else std::cout << " -- [" << mpi_rank << "] --" "\n";
			}

			prev_p = p;
			prev_pp = pp;
			prev_n = n;

			// TODO(Sam): Debug
			//std::cout.flush();


			}
			//std::cout << "\n";

			if(next_available_task > 0)
			{
			// If we are a task that is left with no work and if there is no more
			// work to do we will end up here and we must sync with the currently working tasks
			size_t p = p_max; size_t pp = p_max; size_t n = p_max;
			next_available_task = 0;
			MPI_Bcast(&next_available_task, 1, MPI_SIZE_T, mpi_size-1, MPI_COMM_WORLD);
			MPI_Bcast(&p, 1, MPI_SIZE_T, mpi_size-1, MPI_COMM_WORLD);
			MPI_Bcast(&n, 1, MPI_SIZE_T, mpi_size-1, MPI_COMM_WORLD);
			MPI_Bcast(&pp, 1, MPI_SIZE_T, mpi_size-1, MPI_COMM_WORLD);
			}


			}
			outer_chrono.stop();


			if(mpi_rank == 0)
			std::cout << "\r- full tensor contraction... " << outer_chrono.timing() << "ms" "\n";
			#endif


			#if DEBUG_FULL_CONTRACTION
			DEBUG_check_full_contraction(target_p_max);
			#endif



			// NOTE(Sam): This will be indispensable for MPI communication
			/*
			mpf_content_buffer buffer;
			init_mpf_content_buffer(&buffer, all_ep0_size);

			dump_array_of_mpf_to_buffer(all_ep0_p, &buffer, all_ep0_size);
			DEBUG_dump_all_ep0(p_max, all_ep0_p);

			DEBUG_clear_all_ep0(p_max, all_ep0_p);
			DEBUG_dump_all_ep0(p_max, all_ep0_p);

			undump_array_of_mpf_from_buffer(all_ep0_p, &buffer, all_ep0_size);
			DEBUG_dump_all_ep0(p_max, all_ep0_p);

			clear_mpf_content_buffer(&buffer);
			//*/

			//multiset_buffer = new unsigned int[2*p_max];
			//combination_buffer = new unsigned int[2*p_max];




			// TODO(Sam): This could go into tests
			/*
			size_t current_p = 0;
			for(size_t offset = 0; current_p <= 1000000; offset += current_pcurrent_pcurrent_p/1000+1)
			{
			current_p = get_ep0_p_from_offseted_idx(offset);
			size_t offmin = get_ep0_offset_from_p(current_p);
			size_t offmax = get_ep0_offset_from_p(current_p+1);
			std::cout << "Verifying " << offmin << " <= " << offset << " < " << offmax << " (p="<<current_p<<")\n";
			assert(offmin <= offset && offset < offmax);
			}
			//*/

			// TODO(Sam): This could go into tests
			/*
			size_t p1 = p;
			size_t p2 = 0;
			size_t p3 = 0;
			do {
			std::printf("- (%zu, %zu, %zu) is located at pos %zu\n", p1, p2, p3, get_ep0_idx(p1,p2,p3));

			p2 = (++p2)%(p+1-p3);
			if(p2 == 0)
			++p3;
			p1 = p - p2 -p3;
			} while(p3 <= p);

			std::cout << "---------\n";

			for(size_t idx = 0; idx < get_ep0_size(p); ++idx)
			{
			get_ep0_coords(p, idx, &p1, &p2, &p3);
			std::printf("- idx:%zu corresponds to (%zu, %zu, %zu)\n", idx, p1, p2, p3);
			}
			// End of tests
			//*/

			/*
			for(size_t i = 0; i < 20; ++i)
			std::cout << "Coeff for p=" << i << " is " << coeff_Np0(i) << "\n";
			//*/

			/*
			for(size_t n = 0; n < 150; n += 10)
			for(size_t m = 0; m <= n/2; m += 5)
			std::cout << "[" << n << ", " << m << "] = " << arrangement_factor(n, m) << "\n";
			//*/

			//delete[] combination_buffer;
			//delete[] multiset_buffer;

			//delete[] ep0_p;
			delete[] all_ep0_p;

			MPI_custom_free();
			MPI_Finalize();

			return 0;
			}
			void DEBUG_check_full_contraction(size_t p_max)
			{
			std::ifstream mathematica_csv("tests/contracted_tensors_mathematica.csv");
			std::string line;

			mpf_class max_error = 0.0;
			double avg_timing = 0.0;
			double avg_speedup = 0.0;
			size_t count = 0;

			while(std::getline(mathematica_csv, line))
			{
			std::istringstream iss(line);
			size_t p, n, pp;
			char ignore;
			double value, timing;

			if(! (iss >> p >> ignore >> n >> ignore >> pp >> ignore >> value >> ignore >> timing ) ) { break; }

			if(p > p_max) break;

			auto t1 = std::chrono::high_resolution_clock::now();

			mpf_class computed = compute_full_tensor_contraction(p,n,pp, 0, mpi_size);

			auto t2 = std::chrono::high_resolution_clock::now();

			std::chrono::duration<double, std::milli> ms_timing = t2 - t1;

			if(mpi_rank == 0)
			{
			std::cout << "("<< p << ", " << n << ", "<< pp <<") - Error: " << abs(computed - value)
			<< "\t\t\t Time: " << ms_timing.count() << "ms (old: " << timing1000 << "ms) x"<< timing1000/ms_timing.count() << "\n";
			std::cout.flush();
			}

			mpf_class error = abs(computed - value);

			max_error = max_error > error ? max_error : error;
			avg_speedup += timing*1000/ms_timing.count();
			avg_timing += ms_timing.count();
			++count;
			}
			mathematica_csv.close();

			avg_timing /= (double)count;
			avg_speedup /= (double)count;

			if(mpi_rank == 0)
			std::cout << " === Comp. with Mathematica === \n"
			<< " - Max error: " << max_error << "\n"
			<< " - Avg timing: " << avg_timing << "ms \n"
			<< " - Avg speedup: x" << avg_speedup << "\n";
			}

			void DEBUG_clear_all_ep0(size_t p_max, mpf_class* all_ep0_p)
			{
			size_t upper_bound = get_ep0_offset_from_p(p_max+1);
			for(size_t i = 0; i < upper_bound; ++i)
			all_ep0_p[i] = -100.0;
			}

			void DEBUG_dump_all_ep0(size_t p_max, mpf_class* all_ep0_p)
			{
			if(mpi_rank == 0) std::cout << " === DUMP ALL ep0 === \n";
			for(size_t p = 0; p <= p_max; ++p)
			{
			size_t line_counter = 0;
			for(size_t i = 0; i < get_ep0_size(p); ++i)
			{
			std::cout << all_ep0_p[get_ep0_offset_from_p(p)+i] << " ";
			++line_counter;
			if(line_counter % 12 == 0)
			{
			if(mpi_rank == 0) std::cout << "\n";
			}
			}
			if(p < p_max)
			if(mpi_rank == 0) std::cout << "\n -------- ";
			if(mpi_rank == 0) std::cout << "\n";
			}
			}

TensCont/main.cpp
55115a530536master

main.cpp
View Options

TensCont/main.cpp55115a530536master

main.cppView Options

TensCont/main.cpp
55115a530536master

main.cpp
View Options