/*
								+--------------------------------+
								|                                |
								|   ***  Exact arithmetic  ***   |
								|                                |
								|  Copyright © -tHE SWINe- 2014  |
								|                                |
								|            Exact.h             |
								|                                |
								+--------------------------------+
*/

#pragma once
#ifndef __EXACT_MULTIPRECISION_FLOATING_POINT_ARITHMETIC_INCLUDED
#define __EXACT_MULTIPRECISION_FLOATING_POINT_ARITHMETIC_INCLUDED

/**
 *	@file Exact.h
 *	@author -tHE SWINe-
 *	@date 2014
 *	@brief exact multi-precision floating-point arithmetic
 *
 *	@note In msvc 6.0, internal compiler error (compiler file 'E:\9044\vc98\p2\src\P2\main.c',
 *		line 181) was observed. The workaround is to disable inline function expansion.
 *		Installing Visual C++ 6.0 Service Pack 6 solves this issue (although it removes enhanced
 *		instruction set support if a processor pack is installed).
 *	@todo See if there is a better way to fix the error.
 *
 *	@todo Write documentation comments, uphold naming conventions.
 *	@todo Consider adding restrict on pointers in the CExactBase, it could improve speed somewhat.
 *
 *	@todo Need to limit "reasonable" size of the expansions. For types like float and double,
 *		the maximum compressed size is a good measure, but it is not so practical for extended
 *		and quad formats due to large exponent range (but will still mostly avoid excessive
 *		memory thrashing).
 *	@todo Will need to implement base operations that calculate only a partial result when the
 *		result is to be trimmed (approximate mode, like 4 wide expansions which do not produce
 *		wider results).
 */

#include "Integer.h"
#include "Unused.h"
#include "low/HexFloat.h"

/**
 *	@def __EXPANSION_WARN_UNNECESSARILY_TOO_LONG_EXPANSIONS
 *	@brief if defined, warns about assignments to too long expansions
 */
#define __EXPANSION_WARN_UNNECESSARILY_TOO_LONG_EXPANSIONS

/**
 *	@brief exact arithmetic base functions
 *
 *	The code in this class is heavily inspired by Jonathan Richard
 *	Shewchuk's paper "Adaptive Precision Floating Point Arithmetic
 *	and Fast Robust Geometric Predicates", and the accompanying
 *	source code he gave (http://www.cs.cmu.edu/~quake/robust.html).
 *
 *	There are a few differences though. The arguments of the functions
 *	are ordered (output, input) and whenever a value is decomposed to
 *	multiple scalars, they are ordered starting with least-significant,
 *	going to the most-significant. The implementation is also robust,
 *	making sure the caller allocated sufficient size for the result.
 *	There are new algorithms for product and square of an expansion.
 *
 *	@tparam T is scalar data type
 *
 *	@note This implementation does not contain any bit hacks or type puning,
 *		except for the calculation of the epsilon and splitter constants which
 *		could be done at compile-time that way, and except for b_ExpansionCheck()
 *		which needs to check mantissas for overlapping bits.
 */
template <class T>
class CExactBase {
public:
	typedef T TType; /**< @brief scalar type */
	typedef T TInexactType; /**< @brief scalar type, holding an inexact value (only semantics) */

public:
	/**
	 *	@brief prints the expansion to stdout, as a sum of hexadecimal float components
	 *
	 *	@param[in] n_e_size is size of the expansion
	 *	@param[in] p_e is pointer to the first expansion element
	 */
	static void DumpExpansion(unsigned int n_e_size, const TType *p_e);

	/**
	 *	@brief prints the expansion to stdout, as a single long hexadecimal float number
	 *
	 *	@param[in] n_e_size is size of the expansion
	 *	@param[in] p_e is pointer to the first expansion element
	 *
	 *	@note This function throws std::bad_alloc.
	 *	@note Special numbers are not handled the usual way (e.g. zero is "0x0.0p+0").
	 *	@note This is mostly untested, especially with denormal numbers.
	 */
	static void DumpExpansion_Long(unsigned int n_e_size, const TType *p_e); // throw(std::bad_alloc)

	/**
	 *	@brief calculates width of mantissas in the expansion, in bits
	 *
	 *	@param[in] n_e_size is size of the expansion
	 *	@param[in] p_e is pointer to the first expansion element
	 *
	 *	@return Returns the width of mantissas in the expansion,
	 *		in bits (as when concatenated into a single number).
	 */
	static int n_Mantissa_Width(unsigned int n_e_size, const TType *p_e);

	/**
	 *	@brief sums two scalars precisely
	 *
	 *	@param[out] r_f_low is the least significant part of the output expansion
	 *	@param[out] r_f_high is the most significant part of the output expansion
	 *	@param[in] f_a is the first scalar
	 *	@param[in] f_b is the second scalar
	 */
	static inline void Two_Sum(TType &r_f_low, TType &r_f_high, const TType f_a, const TType f_b);

	/**
	 *	@brief subtracts two scalars precisely
	 *
	 *	@param[out] r_f_low is the least significant part of the output expansion
	 *	@param[out] r_f_high is the most significant part of the output expansion
	 *	@param[in] f_a is the first scalar
	 *	@param[in] f_b is the second scalar
	 */
	static inline void Two_Diff(TType &r_f_low, TType &r_f_high, const TType f_a, const TType f_b);

	/**
	 *	@brief sums two ordered scalars precisely
	 *
	 *	@param[out] r_f_low is the least significant part of the output expansion
	 *	@param[out] r_f_high is the most significant part of the output expansion
	 *	@param[in] f_a is the first scalar, must be of greater magnitude than f_b, or must be zero
	 *	@param[in] f_b is the second scalar
	 */
	static inline void Fast_Two_Sum(TType &r_f_low, TType &r_f_high, const TType f_a, const TType f_b);

	/**
	 *	@brief subtracts two ordered scalars precisely
	 *
	 *	@param[out] r_f_low is the least significant part of the output expansion
	 *	@param[out] r_f_high is the most significant part of the output expansion
	 *	@param[in] f_a is the first scalar, must be of greater magnitude than f_b, or must be zero
	 *	@param[in] f_b is the second scalar
	 */
	static inline void Fast_Two_Diff(TType &r_f_low, TType &r_f_high, const TType f_a, const TType f_b);

	/**
	 *	@brief gets epsilon
	 *	@return Returns epsilon, which equals \f$\frac{ULP(1)}{2}\f$.
	 *	@note Matlab's epsilon is twice this value (EPS, with no arguments,
	 *		is the distance from 1.0 to the next larger double), which means one ULP.
	 */
	static inline const TType f_Epsilon();

	/**
	 *	@brief gets multiplication splitter
	 *	@return Returns a multiplication splitter, which equals \f$\frac{1}{\sqrt{epsilon}} + 1\f$.
	 */
	static inline const TType f_Splitter();

	/**
	 *	@brief splits a scalar to a sum of two scalars while distributing
	 *		the significant bits of the mantissa
	 *
	 *	@param[out] f_a_low is the least significant part of the output expansion
	 *	@param[out] f_a_high is the most significant part of the output expansion
	 *	@param[in] f_a is the scalar to be split
	 */
	static inline void Split(TType &f_a_low, TType &f_a_high, const TType f_a);

	/**
	 *	@brief calculates exact square of a scalar
	 *
	 *	@param[out] f_a_low is the least significant part of the output expansion
	 *	@param[out] f_a_high is the most significant part of the output expansion
	 *	@param[in] f_a is the scalar to be squared
	 */
	static inline void Square(TType &r_f_low, TType &r_f_high, const TType f_a);

	/**
	 *	@brief multiplies two scalars precisely
	 *
	 *	@param[out] r_f_low is the least significant part of the output expansion
	 *	@param[out] r_f_high is the most significant part of the output expansion
	 *	@param[in] f_a is the first scalar
	 *	@param[in] f_b is the second scalar
	 */
	static inline void Two_Product(TType &r_f_low, TType &r_f_high, const TType f_a, const TType f_b);

	/**
	 *	@brief multiplies two scalars precisely
	 *
	 *	@param[out] r_f_low is the least significant part of the output expansion
	 *	@param[out] r_f_high is the most significant part of the output expansion
	 *	@param[in] f_a is the first scalar
	 *	@param[in] f_a_low is the least significant part of the split of f_a (obtained by Split())
	 *	@param[in] f_a_high is the most significant part of the split of f_a (obtained by Split())
	 *	@param[in] f_b is the second scalar
	 */
	static inline void Two_Product_Presplit(TType &r_f_low, TType &r_f_high, const TType f_a,
		const TType f_a_low, const TType f_a_high, const TType f_b);

	/**
	 *	@brief multiplies two scalars precisely
	 *
	 *	@param[out] r_f_low is the least significant part of the output expansion
	 *	@param[out] r_f_high is the most significant part of the output expansion
	 *	@param[in] f_a is the first scalar
	 *	@param[in] f_a_low is the least significant part of the split of f_a (obtained by Split())
	 *	@param[in] f_a_high is the most significant part of the split of f_a (obtained by Split())
	 *	@param[in] f_b is the second scalar
	 *	@param[in] f_b_low is the least significant part of the split of f_b (obtained by Split())
	 *	@param[in] f_b_high is the most significant part of the split of f_b (obtained by Split())
	 */
	static inline void Two_Product_Presplit2(TType &r_f_low, TType &r_f_high, const TType f_a,
		const TType f_a_low, const TType f_a_high, const TType f_b, const TType f_b_low, const TType f_b_high);

	/**
	 *	@brief calculates exact sum of an expansion and a scalar
	 *
	 *	@param[in] n_max_result_size is maximum size of the result (must be equal or greater to n_e_size + 1)
	 *	@param[out] p_result is pointer to the destination expansion (must be allocated to n_max_result_size)
	 *	@param[in] n_e_size is size of the expansion
	 *	@param[in] p_e is pointer to the first expansion element (may be the same as p_result)
	 *	@param[in] f_b is the scalar to be added
	 *
	 *	@return Returns the number of coefficients of the resulting expansion (n_e_size + 1).
	 */
	static inline unsigned int n_Grow_Expansion(unsigned int UNUSED(n_max_result_size), TType *p_result,
		unsigned int n_e_size, const TType *p_e, TType f_b); // p_e and p_result can be the same

	/**
	 *	@brief calculates exact sum of an expansion and a scalar with zero elimination
	 *
	 *	@param[in] n_max_result_size is maximum size of the result (must be equal or greater to n_e_size + 1)
	 *	@param[out] p_result is pointer to the destination expansion (must be allocated to n_max_result_size)
	 *	@param[in] n_e_size is size of the expansion
	 *	@param[in] p_e is pointer to the first expansion element (may be the same as p_result)
	 *	@param[in] f_b is the scalar to be added
	 *
	 *	@return Returns the number of coefficients of the resulting expansion (up to n_e_size + 1).
	 */
	static inline unsigned int n_Grow_Expansion_EliminateZeroes(unsigned int UNUSED(n_max_result_size), TType *p_result,
		unsigned int n_e_size, const TType *p_e, TType f_b); // p_e and p_result can be the same

	/**
	 *	@brief calculates exact sum of two expansions
	 *
	 *	@param[in] n_max_result_size is maximum size of the result (must be equal or greater to n_e_size + n_f_size)
	 *	@param[out] p_result is pointer to the destination expansion (must be allocated to n_max_result_size)
	 *	@param[in] n_e_size is size of the first expansion
	 *	@param[in] p_e is pointer to the first element of the first expansion (may be the same as p_result)
	 *	@param[in] n_f_size is size of the second expansion
	 *	@param[in] p_f is pointer to the first element of the second expansion
	 *
	 *	@return Returns the number of coefficients of the resulting expansion (n_e_size + n_f_size).
	 */
	static inline unsigned int n_Expansion_Sum(unsigned int UNUSED(n_max_result_size), TType *p_result,
		unsigned int n_e_size, const TType *p_e, unsigned int n_f_size, const TType *p_f); // p_e and p_result can be the same expansion

	/**
	 *	@brief calculates exact difference of two expansions
	 *
	 *	@param[in] n_max_result_size is maximum size of the result (must be equal or greater to n_e_size + n_f_size)
	 *	@param[out] p_result is pointer to the destination expansion (must be allocated to n_max_result_size)
	 *	@param[in] n_f_size is size of the second expansion
	 *	@param[in] p_f is pointer to the first element of the second expansion
	 *	@param[in] n_e_size is size of the first expansion
	 *	@param[in] p_e is pointer to the first element of the first expansion (may be the same as p_result)
	 *
	 *	@note This calculates \f$f - e\f$.
	 *
	 *	@return Returns the number of coefficients of the resulting expansion (n_e_size + n_f_size).
	 */
	static inline unsigned int n_Expansion_Difference(unsigned int UNUSED(n_max_result_size), TType *p_result,
		unsigned int n_f_size, const TType *p_f, unsigned int n_e_size, const TType *p_e); // p_e and p_result can be the same expansion

	/**
	 *	@brief calculates exact sum of two expansions, with zero elimination
	 *
	 *	@param[in] n_max_result_size is maximum size of the result (must be equal or greater to n_e_size + n_f_size)
	 *	@param[out] p_result is pointer to the destination expansion (must be allocated to n_max_result_size)
	 *	@param[in] n_e_size is size of the first expansion
	 *	@param[in] p_e is pointer to the first element of the first expansion (may be the same as p_result)
	 *	@param[in] n_f_size is size of the second expansion
	 *	@param[in] p_f is pointer to the first element of the second expansion
	 *
	 *	@note This eliminates zeroes at the end of the calculation, in a separate run.
	 *
	 *	@return Returns the number of coefficients of the resulting expansion (up to n_e_size + n_f_size).
	 */
	static inline unsigned int n_Expansion_Sum_EliminateZeroes1(unsigned int UNUSED(n_max_result_size), TType *p_result,
		unsigned int n_e_size, const TType *p_e, unsigned int n_f_size, const TType *p_f); // p_e and p_result can be the same expansion

	/**
	 *	@brief calculates exact difference of two expansions, with zero elimination
	 *
	 *	@param[in] n_max_result_size is maximum size of the result (must be equal or greater to n_e_size + n_f_size)
	 *	@param[out] p_result is pointer to the destination expansion (must be allocated to n_max_result_size)
	 *	@param[in] n_f_size is size of the second expansion
	 *	@param[in] p_f is pointer to the first element of the second expansion
	 *	@param[in] n_e_size is size of the first expansion
	 *	@param[in] p_e is pointer to the first element of the first expansion (may be the same as p_result)
	 *
	 *	@note This calculates \f$f - e\f$.
	 *	@note This eliminates zeroes at the end of the calculation, in a separate run.
	 *
	 *	@return Returns the number of coefficients of the resulting expansion (up to n_e_size + n_f_size).
	 */
	static inline unsigned int n_Expansion_Difference_EliminateZeroes1(unsigned int UNUSED(n_max_result_size), TType *p_result,
		unsigned int n_f_size, const TType *p_f, unsigned int n_e_size, const TType *p_e); // p_e and p_result can be the same expansion

	/**
	 *	@brief calculates exact sum of two expansions, with zero elimination
	 *
	 *	@param[in] n_max_result_size is maximum size of the result (must be equal or greater to n_e_size + n_f_size)
	 *	@param[out] p_result is pointer to the destination expansion (must be allocated to n_max_result_size)
	 *	@param[in] n_e_size is size of the first expansion
	 *	@param[in] p_e is pointer to the first element of the first expansion (may be the same as p_result)
	 *	@param[in] n_f_size is size of the second expansion
	 *	@param[in] p_f is pointer to the first element of the second expansion
	 *
	 *	@note This eliminates zeroes in the progress of the calculation.
	 *
	 *	@return Returns the number of coefficients of the resulting expansion (up to n_e_size + n_f_size).
	 */
	static inline unsigned int n_Expansion_Sum_EliminateZeroes2(unsigned int UNUSED(n_max_result_size), TType *p_result,
		unsigned int n_e_size, const TType *p_e, unsigned int n_f_size, const TType *p_f); // p_e and p_result can be the same expansion

	/**
	 *	@brief calculates exact difference of two expansions, with zero elimination
	 *
	 *	@param[in] n_max_result_size is maximum size of the result (must be equal or greater to n_e_size + n_f_size)
	 *	@param[out] p_result is pointer to the destination expansion (must be allocated to n_max_result_size)
	 *	@param[in] n_f_size is size of the second expansion
	 *	@param[in] p_f is pointer to the first element of the second expansion
	 *	@param[in] n_e_size is size of the first expansion
	 *	@param[in] p_e is pointer to the first element of the first expansion (may be the same as p_result)
	 *
	 *	@note This calculates \f$f - e\f$.
	 *	@note This eliminates zeroes in the progress of the calculation.
	 *
	 *	@return Returns the number of coefficients of the resulting expansion (up to n_e_size + n_f_size).
	 */
	static inline unsigned int n_Expansion_Difference_EliminateZeroes2(unsigned int UNUSED(n_max_result_size), TType *p_result,
		unsigned int n_f_size, const TType *p_f, unsigned int n_e_size, const TType *p_e); // p_e and p_result can be the same expansion

	/**
	 *	@brief calculates exact sum of two expansions, in linear time
	 *
	 *	@param[in] n_max_result_size is maximum size of the result (must be equal or greater to n_e_size + n_f_size)
	 *	@param[out] p_result is pointer to the destination expansion (must be allocated to n_max_result_size)
	 *	@param[in] n_e_size is size of the first expansion
	 *	@param[in] p_e is pointer to the first element of the first expansion (may be the same as p_result)
	 *	@param[in] n_f_size is size of the second expansion
	 *	@param[in] p_f is pointer to the first element of the second expansion
	 *
	 *	@note The linear time is amortized and this is only faster for long expansions, otherwise n_Expansion_Sum() is faster.
	 *
	 *	@return Returns the number of coefficients of the resulting expansion (n_e_size + n_f_size).
	 */
	static inline unsigned int n_Linear_Expansion_Sum(unsigned int UNUSED(n_max_result_size),
		const TType *p_result, unsigned int n_e_size, const TType *p_e, int n_f_size, const TType *p_f); // should be used for long expansions, otherwise n_Expansion_Sum is faster

	/**
	 *	@brief calculates exact sum of two expansions, in linear time and with zero elimination
	 *
	 *	@param[in] n_max_result_size is maximum size of the result (must be equal or greater to n_e_size + n_f_size)
	 *	@param[out] p_result is pointer to the destination expansion (must be allocated to n_max_result_size)
	 *	@param[in] n_e_size is size of the first expansion
	 *	@param[in] p_e is pointer to the first element of the first expansion (may be the same as p_result)
	 *	@param[in] n_f_size is size of the second expansion
	 *	@param[in] p_f is pointer to the first element of the second expansion
	 *
	 *	@note The linear time is amortized and this is only faster for long expansions, otherwise n_Expansion_Sum() is faster.
	 *
	 *	@return Returns the number of coefficients of the resulting expansion (up to n_e_size + n_f_size).
	 */
	static inline unsigned int n_Linear_Expansion_Sum_EliminateZeroes(unsigned int UNUSED(n_max_result_size),
		const TType *p_result, unsigned int n_e_size, const TType *p_e, int n_f_size, const TType *p_f); // should be used for long expansions, otherwise n_Expansion_Sum is faster

	/**
	 *	@brief calculates exact sum of two expansions, faster but without preserving the nonoverlapping property
	 *
	 *	@param[in] n_max_result_size is maximum size of the result (must be equal or greater to n_e_size + n_f_size)
	 *	@param[out] p_result is pointer to the destination expansion (must be allocated to n_max_result_size)
	 *	@param[in] n_e_size is size of the first expansion
	 *	@param[in] p_e is pointer to the first element of the first expansion (may be the same as p_result)
	 *	@param[in] n_f_size is size of the second expansion
	 *	@param[in] p_f is pointer to the first element of the second expansion
	 *
	 *	@note This should be only used at the end of evaluation (e.g. before approximation) and not mixed with
	 *		standard sum or multiplications.
	 *	@note This is only used when called explicitly (never in CExpansion::operator+()).
	 *
	 *	@return Returns the number of coefficients of the resulting expansion (n_e_size + n_f_size).
	 */
	static inline unsigned int n_Fast_Expansion_Sum(unsigned int UNUSED(n_max_result_size),
		TType *p_result, unsigned int n_e_size, const TType *p_e, unsigned int n_f_size, const TType *p_f); // should be used only at the end and not mixed with standard sum or multiplications, only when user requests it (not in operator +)

	/**
	 *	@brief calculates exact sum of two expansions, faster but without preserving the nonoverlapping property, with zero elimination
	 *
	 *	@param[in] n_max_result_size is maximum size of the result (must be equal or greater to n_e_size + n_f_size)
	 *	@param[out] p_result is pointer to the destination expansion (must be allocated to n_max_result_size)
	 *	@param[in] n_e_size is size of the first expansion
	 *	@param[in] p_e is pointer to the first element of the first expansion (may be the same as p_result)
	 *	@param[in] n_f_size is size of the second expansion
	 *	@param[in] p_f is pointer to the first element of the second expansion
	 *
	 *	@note This should be only used at the end of evaluation (e.g. before approximation) and not mixed with
	 *		standard sum or multiplications.
	 *	@note This is only used when called explicitly (never in CExpansion::operator+()).
	 *
	 *	@return Returns the number of coefficients of the resulting expansion (up to n_e_size + n_f_size).
	 */
	static inline unsigned int n_Fast_Expansion_Sum_EliminateZeroes(unsigned int UNUSED(n_max_result_size),
		TType *p_result, unsigned int n_e_size, const TType *p_e, unsigned int n_f_size, const TType *p_f); // should be used only at the end and not mixed with standard sum or multiplications, only when user requests it (not in operator +)

	/**
	 *	@brief calculates exact product of an expansion and a scalar
	 *
	 *	@param[in] n_max_result_size is maximum size of the result (must be equal or greater to n_e_size * 2)
	 *	@param[out] p_result is pointer to the destination expansion (must be allocated to n_max_result_size)
	 *	@param[in] n_e_size is size of the expansion
	 *	@param[in] p_e is pointer to the first expansion element
	 *	@param[in] f_scale is the scale
	 *
	 *	@return Returns the number of coefficients of the resulting expansion (n_e_size * 2).
	 */
	static inline unsigned int n_Scale_Expansion(unsigned int UNUSED(n_max_result_size),
		TType *p_result, unsigned int n_e_size, const TType *p_e, TType f_scale);

	/**
	 *	@brief calculates exact product of an expansion and a scalar, with zero elimination
	 *
	 *	@param[in] n_max_result_size is maximum size of the result (must be equal or greater to n_e_size * 2)
	 *	@param[out] p_result is pointer to the destination expansion (must be allocated to n_max_result_size)
	 *	@param[in] n_e_size is size of the expansion
	 *	@param[in] p_e is pointer to the first expansion element
	 *	@param[in] f_scale is the scale
	 *
	 *	@return Returns the number of coefficients of the resulting expansion (up to n_e_size * 2).
	 */
	static inline unsigned int n_Scale_Expansion_EliminateZeroes(unsigned int UNUSED(n_max_result_size),
		TType *p_result, unsigned int n_e_size, const TType *p_e, TType f_scale);

	/**
	 *	@brief calculates exact product of an expansion and a scalar and store splits of the expansion
	 *
	 *	@param[in] n_max_result_size is maximum size of the result (must be equal or greater to n_e_size * 2)
	 *	@param[out] p_result is pointer to the destination expansion (must be allocated to n_max_result_size)
	 *	@param[in] n_e_size is size of the expansion
	 *	@param[out] p_e_low is pointer to the least significant parts of the splits (must be allocated to n_e_size)
	 *	@param[out] p_e_high is pointer to the most significant parts of the splits (must be allocated to n_e_size)
	 *	@param[in] p_e is pointer to the first expansion element
	 *	@param[in] f_scale is the scale
	 *	@param[in] f_scale_low is the least significant part of the scale
	 *	@param[in] f_scale_high is the most significant part of the scale
	 *
	 *	@return Returns the number of coefficients of the resulting expansion (n_e_size * 2).
	 */
	static inline unsigned int n_Scale_Expansion_Split(unsigned int UNUSED(n_max_result_size),
		TType *p_result, unsigned int n_e_size, TType *p_e_low, TType *p_e_high, const TType *p_e,
		TType f_scale, TType f_scale_low, TType f_scale_high);

	/**
	 *	@brief calculates exact product of an expansion and a scalar,
	 *		while utilizing previously calculated splits of the expansion
	 *
	 *	@param[in] n_max_result_size is maximum size of the result (must be equal or greater to n_e_size * 2)
	 *	@param[out] p_result is pointer to the destination expansion (must be allocated to n_max_result_size)
	 *	@param[in] n_e_size is size of the expansion
	 *	@param[in] p_e is pointer to the first expansion element
	 *	@param[in] p_e_low is pointer to the least significant parts of the splits
	 *	@param[in] p_e_high is pointer to the most significant parts of the splits
	 *	@param[in] f_scale is the scale
	 *	@param[in] f_scale_low is the least significant part of the scale
	 *	@param[in] f_scale_high is the most significant part of the scale
	 *
	 *	@return Returns the number of coefficients of the resulting expansion (n_e_size * 2).
	 */
	static inline unsigned int n_Scale_Expansion_Presplit(unsigned int UNUSED(n_max_result_size),
		TType *p_result, unsigned int n_e_size, const TType *p_e, const TType *p_e_low, const TType *p_e_high,
		TType f_scale, TType f_scale_low, TType f_scale_high);

	/**
	 *	@brief calculates exact product of two expansions
	 *
	 *	@param[in] n_max_result_size is maximum size of the result (must be equal or greater to n_e_size * n_f_size * 2)
	 *	@param[out] p_result is pointer to the destination expansion (must be allocated to n_max_result_size)
	 *	@param[in] n_e_size is size of the first expansion
	 *	@param[in] p_e is pointer to the first element of the first expansion
	 *	@param[out] p_e_low is pointer to intermediate storage, allocated to n_e_size elements
	 *	@param[out] p_e_high is pointer to intermediate storage, allocated to n_e_size elements
	 *	@param[in] n_max_e_scaled_size is at least n_e_size * 2
	 *	@param[out] p_e_scaled is pointer to intermediate storage, allocated to n_max_e_scaled_size elements
	 *	@param[in] n_f_size is size of the second expansion
	 *	@param[in] p_f is pointer to the first element of the second expansion
	 *
	 *	@return Returns the number of coefficients of the resulting expansion (n_e_size * n_f_size * 2).
	 */
	static inline unsigned int n_Expansion_Product(unsigned int n_max_result_size, TType *p_result,
		unsigned int n_e_size, const TType *p_e, TType *p_e_low, TType *p_e_high,
		unsigned int n_max_e_scaled_size, TType *p_e_scaled, unsigned int n_f_size, const TType *p_f);

	/**
	 *	@brief calculates exact product of two expansions, with zero elimination
	 *
	 *	@param[in] n_max_result_size is maximum size of the result (must be equal or greater to n_e_size * n_f_size * 2)
	 *	@param[out] p_result is pointer to the destination expansion (must be allocated to n_max_result_size)
	 *	@param[in] n_e_size is size of the first expansion
	 *	@param[in] p_e is pointer to the first element of the first expansion
	 *	@param[out] p_e_low is pointer to intermediate storage, allocated to n_e_size elements
	 *	@param[out] p_e_high is pointer to intermediate storage, allocated to n_e_size elements
	 *	@param[in] n_max_e_scaled_size is at least n_e_size * 2
	 *	@param[out] p_e_scaled is pointer to intermediate storage, allocated to n_max_e_scaled_size elements
	 *	@param[in] n_f_size is size of the second expansion
	 *	@param[in] p_f is pointer to the first element of the second expansion
	 *
	 *	@return Returns the number of coefficients of the resulting expansion (up to n_e_size * n_f_size * 2).
	 */
	static inline unsigned int n_Expansion_Product_EliminateZeroes(unsigned int n_max_result_size, TType *p_result,
		unsigned int n_e_size, const TType *p_e, TType *p_e_low, TType *p_e_high,
		unsigned int n_max_e_scaled_size, TType *p_e_scaled, unsigned int n_f_size, const TType *p_f);

	/**
	 *	@brief calculates exact square of an expansion
	 *
	 *	@param[in] n_max_result_size is maximum size of the result (must be equal or greater to n_e_size * n_e_size * 2)
	 *	@param[out] p_result is pointer to the destination expansion (must be allocated to n_max_result_size)
	 *	@param[in] n_e_size is size of the first expansion
	 *	@param[in] p_e is pointer to the first element of the first expansion
	 *	@param[out] p_e_low is pointer to intermediate storage, allocated to n_e_size elements
	 *	@param[out] p_e_high is pointer to intermediate storage, allocated to n_e_size elements
	 *	@param[in] n_max_e_scaled_size is at least n_e_size * 2
	 *	@param[out] p_e_scaled is pointer to intermediate storage, allocated to n_max_e_scaled_size elements
	 *
	 *	@return Returns the number of coefficients of the resulting expansion (n_e_size * n_e_size * 2).
	 */
	static inline unsigned int n_Expansion_Square(unsigned int n_max_result_size, TType *p_result,
		unsigned int n_e_size, const TType *p_e, TType *p_e_low, TType *p_e_high,
		unsigned int n_max_e_scaled_size, TType *p_e_scaled);

	/**
	 *	@brief calculates exact square of an expansion, with zero elimination
	 *
	 *	@param[in] n_max_result_size is maximum size of the result (must be equal or greater to n_e_size * n_e_size * 2)
	 *	@param[out] p_result is pointer to the destination expansion (must be allocated to n_max_result_size)
	 *	@param[in] n_e_size is size of the first expansion
	 *	@param[in] p_e is pointer to the first element of the first expansion
	 *	@param[out] p_e_low is pointer to intermediate storage, allocated to n_e_size elements
	 *	@param[out] p_e_high is pointer to intermediate storage, allocated to n_e_size elements
	 *	@param[in] n_max_e_scaled_size is at least n_e_size * 2
	 *	@param[out] p_e_scaled is pointer to intermediate storage, allocated to n_max_e_scaled_size elements
	 *
	 *	@return Returns the number of coefficients of the resulting expansion (up to n_e_size * n_e_size * 2).
	 */
	static inline unsigned int n_Expansion_Square_EliminateZeroes(unsigned int n_max_result_size, TType *p_result,
		unsigned int n_e_size, const TType *p_e, TType *p_e_low, TType *p_e_high,
		unsigned int n_max_e_scaled_size, TType *p_e_scaled);

	/**
	 *	@brief eliminates spurious zeroes in an expansion, inplace
	 *
	 *	@param[in] n_e_size is size of the first expansion
	 *	@param[in,out] p_e is pointer to the first element of the first expansion
	 *
	 *	@return Returns the number of coefficients of the resulting expansion (up to n_e_size).
	 */
	static inline unsigned int n_EliminateZeroes(unsigned int n_e_size, TType *p_e);

	/**
	 *	@brief compresses an expansion
	 *
	 *	This calculates the shortest expansion with the same value as input.
	 *
	 *	@param[in] n_max_result_size is maximum size of the result (must be equal or greater to
	 *		min(n_e_size, CMaxExpansionSize<TType>::n_longest_expansion_size))
	 *	@param[out] p_result is pointer to the destination expansion (must be allocated to n_max_result_size)
	 *	@param[in] n_e_size is size of the expansion
	 *	@param[in] p_e is pointer to the first expansion element (may be the same as p_result)
	 *
	 *	@return Returns the number of coefficients of the resulting expansion (up to n_e_size).
	 *
	 *	@note The maximum possible length of an expansion is limited by TType exponent range, see CMaxExpansionSize.
	 */
	static inline unsigned int n_Compress(unsigned int UNUSED(n_max_result_size),
		TType *p_result, unsigned int n_e_size, const TType *p_e); // p_e and p_result may be the same

	/**
	 *	@brief calculates an approximate value of an expansion
	 *
	 *	@param[in] n_e_size is size of the expansion
	 *	@param[in] p_e is pointer to the first expansion element
	 *
	 *	@return Returns an approximate value of the input expansion.
	 */
	static inline TInexactType f_ApproximateValue(unsigned int n_e_size, const TType *p_e);

	/**
	 *	@brief performs debug checking on an expansion (only if _DEBUG is defined)
	 *
	 *	This makes sure that the expansion is sorted and nonoverlapping.
	 *
	 *	@param[in] n_e_size is size of the expansion
	 *	@param[in] p_e is pointer to the first expansion element
	 *	@param[in] b_allow_NaNs is nonsignalling NaN flag (if set, NaNs are quietly ignored;
	 *		if not set, NaNs trigger an assertion)
	 */
	static void ExpansionCheck(unsigned int n_e_size, const TType *p_e, bool b_allow_NaNs = false); // this only performs the check in debug

	/**
	 *	@brief performs sanity checking on an expansion
	 *
	 *	This makes sure that the expansion is sorted and nonoverlapping.
	 *
	 *	@param[in] n_e_size is size of the expansion
	 *	@param[in] p_e is pointer to the first expansion element
	 *	@param[in] b_allow_NaNs is nonsignalling NaN flag (if set, NaNs are quietly ignored;
	 *		if not set, NaNs trigger an assertion)
	 */
	static bool b_ExpansionCheck(unsigned int n_e_size, const TType *p_e, bool b_allow_NaNs = false);
};

#include "Exact.inl"

/**
 *	@brief namespace with tags for expansion constructor disambiguation
 */
namespace expansion {

extern const expansion_internal::TFromScalar_Tag from_scalar; /**< @brief create expansion from a single scalar */
extern const expansion_internal::TFromScalar_ZeroPad_Tag from_scalar_zero_pad; /**< @brief create expansion from a single scalar and pad with zeros (same as from_scalar, but does not produce the "EXPANSION_UNNECESSARILY_TOO_LONG" warnings when initializing expansion longer than 1 element) */
extern const expansion_internal::TFromArray_Tag from_array; /**< @brief create expansion from an array of scalars, the array size must match expansion size */
extern const expansion_internal::TFromArray_ZeroPad_Tag from_array_zero_pad; /**< @brief create expansion from an array of scalars and pad with zeros (same as from_array, but long expansion may be initialized by a shorter array; initializing a short expansion by a longer array is not permitted) */
extern const expansion_internal::TFromTwoSum_Tag from_sum_of_scalars; /**< @brief create expansion from a sum of two scalars */
extern const expansion_internal::TFromTwoDiff_Tag from_difference_of_scalars; /**< @brief create expansion from a difference of two scalars */
extern const expansion_internal::TFromTwoProd_Tag from_product_of_scalars; /**< @brief create expansion from a product of two scalars */
extern const expansion_internal::TFromScalarSquare_Tag from_square_of_scalar; /**< @brief create expansion from a square of a single scalar */

} // ~expansion

/**
 *	@brief calculates expansion size
 *
 *	A limitation of the arbitrary precision expansions is the limited exponent
 *	range, the same as that of double (or float, respectively). Hence the maximum
 *	precision is slightly less than 2048 + 53 bits in double (40 components),
 *	and 256 + 24 bits in float (12 components). Also note that the maximum precision
 *	can only be used if the first component is near infinity and the last component
 *	is a denormal number.
 *
 *	For 80-bit extended, it is 32768 + 64 bits, up to 513 components.
 *	For quadruple, it is 32768 + 113 bits, up to 300 components.
 *
 *	This, however, does not always mean that longer expansions do not make sense:
 *	in the worst possible case, each component of the expansion may carry only a single
 *	significant bit, allowing for expansions of thousands of elements. On the other
 *	hand, compressed expansions are guaranteed to not be longer than the value
 *	of n_longest_expansion_size calculated here.
 *
 *	@tparam CScalar is floating-point scalar data type (e.g. double)
 */
template <class CScalar>
class CMaxExpansionSize {
public:
	typedef CFloatRepresentationTraits<CScalar> _TyScalarTraits; /**< scalar traits */

	/**
	 *	@brief results, stored as enum
	 */
	enum {
		n_exponent_range = 1 << _TyScalarTraits::n_exponent_bit_num, /**< @brief exponent range */
		n_mantissa_bits = _TyScalarTraits::n_fraction_bit_num + 1, /**< @brief number of mantissa bits */
		n_longest_expansion_bits = n_exponent_range + n_mantissa_bits - 1, /**< @brief number of significant bits of the longest possible finite-value expansion */
		n_longest_expansion_size = (n_longest_expansion_bits + n_mantissa_bits - 1) / n_mantissa_bits /**< @brief size of the longest possible finite-value expansion */
	};
};

/**
 *	@brief expansion class; represents a multiple-precision floating point number
 *
 *	@tparam T is scalar data type (e.g. double)
 *	@tparam _n_max_coeff_num is maximum expansion size, in scalar elements
 *	@tparam _b_allow_dynamic is dynamic expansion size flag (default not set; if set, the expansion
 *		can be shorter than the specified expansion size but not longer and the memory footprint does
 *		not change with the size)
 *
 *	@note See also CMaxExpansionSize, CApproximate, t_Expansion(), f_Approximate() and t_Trimmed().
 */
template <class T, const unsigned int _n_max_coeff_num, const bool _b_allow_dynamic /*= false*/>
class CExpansion : public expansion_internal::CExpansionStorage<T, _n_max_coeff_num, _b_allow_dynamic> {
#if !defined(_MSC_VER) || defined(__MWERKS__) || _MSC_VER > 1200
	template <class T_, const unsigned int _n_max_coeff_num_, const bool _b_allow_dynamic_> friend class CExpansion; // different specializations of expansions want to access each other's members
	template <class CSrcExpansion, class COldScalarType, class CTargetScalarType> friend class expansion_internal::CTypeCastResult; // type casts want to create uninitialized expansions, which is otherwise prohibited
	template <class CSrcExpansion, class CTargetExpansion> friend class expansion_internal::CDynamicCastResult; // dynamic casts want to create uninitialized expansions, which is otherwise prohibited
	// msvc 6.0 can't handle this

	//template <const bool _b_needs_compact, class CResultType> friend expansion_internal::CCompactionAlgorithm<_b_needs_compact, CResultType>::CAlgorithm; // compaction algorithm requires access to uninitialized construction
	// internal compiler error in vs 2008

	/*friend expansion_internal::CCompactionAlgorithm<n_max_coeff_num / (n_max_fully_significant_coeff_num + 1) != 0,
		typename expansion_internal::CChooseType<CCompactSelf, const CSelf&,
		n_max_coeff_num / (n_max_fully_significant_coeff_num + 1) != 0>::CResult>::CAlgorithm;*/
	// does not work, this is not the correct expansion: the CCompactSelf needs to be the friend

	friend expansion_internal::CCompactionAlgorithm<true, CExpansion<T, _n_max_coeff_num, _b_allow_dynamic> >::CAlgorithm;
	// this is CCompactSelf to some other CExpansion specialization (of many other specializations poitentially)
#else // !_MSC_VER || __MWERKS__ || _MSC_VER > 1200
	friend expansion_internal::CExpansionTraits<CExpansion<T, _n_max_coeff_num, _b_allow_dynamic> >;
	// allow CExpansionTraits to access private in msvc 6.0 (can't do it in CExpansionStorage,
	// as in there the other CExpansion template arguments are not known yet and it is not possible
	// to befrient all the possible types (in msvc 6.0)
#endif // !_MSC_VER || __MWERKS__ || _MSC_VER > 1200

public:
	/**
	 *	@brief intermediates stored as an enum
	 */
	enum {
		n_max_fully_significant_coeff_num = (_n_max_coeff_num <
			CMaxExpansionSize<T>::n_longest_expansion_size)? _n_max_coeff_num :
			CMaxExpansionSize<T>::n_longest_expansion_size /**< @brief maximum size of an expansion where all the coefficients' mantissas are fully significant (unlikely to occur unless working with numbers in the 1e+300 magnitude range) */
	};

	typedef CExpansion<T, _n_max_coeff_num, _b_allow_dynamic> CSelf; /**< @brief type of self */
	typedef CExpansion<T, _n_max_coeff_num, true> CDynamicSelf; /**< @brief type of self with the dynamic flag set */
	typedef CExpansion<T, n_max_fully_significant_coeff_num, _b_allow_dynamic> CCompactSelf; /**< @brief type of self, limited to the minimum size required to represent any number, representable by CSelf */
	typedef CExpansion<T, n_max_fully_significant_coeff_num, true> CDynamicCompactSelf; /**< @brief type of self with the dynamic flag set, and the size limited to the minimum size required to represent any number, representable by CSelf */

	typedef CExactBase<TType> CBase; /**< @brief exact operations implementation */

	typedef void (CSelf::*bool_type)(void) const; /**< @brief safe bool idiom; such pointer cannot be deleted and is not indexable, seems to work in msvc 6.0 as well */

	/*template <class _CType, const int _size, const bool _varsize>
	struct CSpecialize {
		typedef CExpansion<_CType, _size, _varsize> CResult;
	};*/ // not needed anymore

public:
#if !defined(_MSC_VER) || defined(__MWERKS__) || _MSC_VER > 1200

	/**
	 *	@brief copy-constructor
	 *
	 *	@tparam T2 is scalar type of the assigned expansion
	 *	@tparam _n_max_coeff_num2 is maximum number of the coeffs of the assigned expansion
	 *	@tparam _b_allow_dynamic2 is dynamic expansion size flag of the assigned expansion
	 *
	 *	@param[in] r_other is the assigned expansion
	 *
	 *	@note This does not allow assignment of longer expansions to shorter ones.
	 *		This is either detected at compile time if the assigned expansion is static,
	 *		or at runtime (debug only) if it is dynamic.
	 *	@note This emits warnings if this expansion is longer than the maximum size
	 *		of the assigned expansion (independent of the dynamic size).
	 */
	template <class T2, const unsigned int _n_max_coeff_num2, const bool _b_allow_dynamic2>
	inline CExpansion(const CExpansion<T2, _n_max_coeff_num2, _b_allow_dynamic2> &r_other)
	{
		typedef CExpansion<T2, _n_max_coeff_num2, _b_allow_dynamic2> CExpansion2;
		enum {
			n_required_size = (_b_allow_dynamic2)?
				/*1*/ // would trigger size warning when assigning dynamic to static
				/*n_max_coeff_num*/ // does not trigger warning, but ignores assignments of dynamic expansions that are *always* shorter
				((n_max_coeff_num > _n_max_coeff_num2)? _n_max_coeff_num2 : n_max_coeff_num) // if the assigned expansion is shorter, do warn. if it is longer, leave it to the runtime check.
				: _n_max_coeff_num2 // not dynamic
		};
		typedef typename expansion_internal::CCompatibleTypeCheck<CSelf, CExpansion2>::CAssert CTypeCheck; // without CAssert it will not check
		/*typedef*/ typename CCheckSize<n_required_size>::CAssert /*CSizeCheck*/(); // without CAssert it will not check
		_ASSERTE(r_other.m_n_size <= n_max_coeff_num); // check at runtime
		for(unsigned int i = 0, n = r_other.m_n_size; i < n; ++ i)
			m_p_value[i] = r_other.m_p_value[i];
		Uninitialized_Contract(r_other.m_n_size);
	}

#else // !_MSC_VER || __MWERKS__ || _MSC_VER > 1200

private:
	inline CExpansion(long double f); // initialization using scalar not allowed, use CExpansion(T,TFromScalar_Tag)
	inline CExpansion(double f); // initialization using scalar not allowed, use CExpansion(T,TFromScalar_Tag)
	inline CExpansion(float f); // initialization using scalar not allowed, use CExpansion(T,TFromScalar_Tag)

public:
	/**
	 *	@brief copy-constructor (msvc 6.0)
	 *
	 *	@tparam CExpansion2 is type of the assigned expansion (must be CExpansion specialization)
	 *
	 *	@param[in] r_other is the assigned expansion
	 *
	 *	@note This does not allow assignment of longer expansions to shorter ones.
	 *		This is either detected at compile time if the assigned expansion is static,
	 *		or at runtime (debug only) if it is dynamic.
	 *	@note This emits warnings if this expansion is longer than the maximum size
	 *		of the assigned expansion (independent of the dynamic size).
	 */
	template <class CExpansion2>
	inline CExpansion(const CExpansion2 &r_other) // msvc 6.0 can only do this
	{
		typedef typename expansion_internal::CCompatibleTypeCheck<CSelf, CExpansion2>::CAssert CTypeCheck; // without CAssert it will not check
		enum {
			//n_required_size = (CExpansion2::b_allow_dynamic)? /*1*/n_max_coeff_num : CExpansion2::n_max_coeff_num
			_b_allow_dynamic2 = CExpansion2::b_allow_dynamic,
			_n_max_coeff_num2 = CExpansion2::n_max_coeff_num,
			n_required_size = (_b_allow_dynamic2)?
				/*1*/ // would trigger size warning when assigning dynamic to static
				/*n_max_coeff_num*/ // does not trigger warning, but ignores assignments of dynamic expansions that are *always* shorter
				((n_max_coeff_num > _n_max_coeff_num2)? _n_max_coeff_num2 : n_max_coeff_num) // if the assigned expansion is shorter, do warn. if it is longer, leave it to the runtime check.
				: _n_max_coeff_num2 // not dynamic
		};
		/*typedef*/ typename CCheckSize<n_required_size>::CAssert /*CSizeCheck*/(); // without CAssert it will not check
		_ASSERTE(r_other.n_Size() <= n_max_coeff_num); // check at runtime
		for(unsigned int i = 0, n = r_other.n_Size(); i < n; ++ i)
			m_p_value[i] = r_other[i];
		Uninitialized_Contract(r_other.n_Size());
	}

#endif // !_MSC_VER || __MWERKS__ || _MSC_VER > 1200

	/**
	 *	@brief constructor; initializes an expansion from a scalar
	 *
	 *	@param[in] f_value is scalar to initialize with
	 *	@param[in] t_tag is type scheduling tag (unused at runtime)
	 *
	 *	@note This warns if initializing a static expansion of size greater
	 *		than one (see also CExpansion(TType,TFromScalar_ZeroPad_Tag)).
	 */
	inline CExpansion(TType f_value, expansion_internal::TFromScalar_Tag UNUSED(t_tag));

	/**
	 *	@brief constructor; initializes an expansion from a scalar, does not warn about implicit zero padding
	 *
	 *	@param[in] f_value is scalar to initialize with
	 *	@param[in] t_tag is type scheduling tag (unused at runtime)
	 */
	inline CExpansion(TType f_value, expansion_internal::TFromScalar_ZeroPad_Tag UNUSED(t_tag));

	/**
	 *	@brief constructor; initializes an expansion from raw data
	 *
	 *	@param[in] n_size is size of the input expansion
	 *	@param[in] p_value is pointer to the first element of the input expansion
	 *	@param[in] t_tag is type scheduling tag (unused at runtime)
	 *
	 *	@note This warns if initializing a static expansion of size greater
	 *		than n_size (see also CExpansion(TType,TFromArray_ZeroPad_Tag)).
	 */
	inline CExpansion(unsigned int n_size, const TType *p_value, expansion_internal::TFromArray_Tag UNUSED(t_tag));

	/**
	 *	@brief constructor; initializes an expansion from raw data, does not warn about implicit zero padding
	 *
	 *	@param[in] n_size is size of the input expansion
	 *	@param[in] p_value is pointer to the first element of the input expansion
	 *	@param[in] t_tag is type scheduling tag (unused at runtime)
	 */
	inline CExpansion(unsigned int n_size, const TType *p_value, expansion_internal::TFromArray_ZeroPad_Tag UNUSED(t_tag));

	/**
	 *	@brief constructor; initializes an expansion from a sum of two scalars
	 *
	 *	@param[in] f_left is scalar to initialize with
	 *	@param[in] f_right is scalar to initialize with
	 *	@param[in] t_tag is type scheduling tag (unused at runtime)
	 *
	 *	@note This warns if initializing a static expansion of size greater than two.
	 *	@note This is the same as writing CExpansion<TType, 2>(t_Expansion(f_left) + f_right).
	 */
	inline CExpansion(TType f_left, TType f_right, expansion_internal::TFromTwoSum_Tag UNUSED(t_tag));

	/**
	 *	@brief constructor; initializes an expansion from a difference of two scalars
	 *
	 *	@param[in] f_left is scalar to initialize with
	 *	@param[in] f_right is scalar to initialize with
	 *	@param[in] t_tag is type scheduling tag (unused at runtime)
	 *
	 *	@note This warns if initializing a static expansion of size greater than two.
	 *	@note This is the same as writing CExpansion<TType, 2>(t_Expansion(f_left) - f_right).
	 */
	inline CExpansion(TType f_left, TType f_right, expansion_internal::TFromTwoDiff_Tag UNUSED(t_tag));

	/**
	 *	@brief constructor; initializes an expansion from a product of two scalars
	 *
	 *	@param[in] f_left is scalar to initialize with
	 *	@param[in] f_right is scalar to initialize with
	 *	@param[in] t_tag is type scheduling tag (unused at runtime)
	 *
	 *	@note This warns if initializing a static expansion of size greater than two.
	 *	@note This is the same as writing CExpansion<TType, 2>(t_Expansion(f_left) * f_right).
	 */
	inline CExpansion(TType f_left, TType f_right, expansion_internal::TFromTwoProd_Tag UNUSED(t_tag));

	/**
	 *	@brief constructor; initializes an expansion from a square of a scalar
	 *
	 *	@param[in] f_scalar is scalar to initialize with
	 *	@param[in] t_tag is type scheduling tag (unused at runtime)
	 *
	 *	@note This warns if initializing a static expansion of size greater than two.
	 *	@note This is the same as writing CExpansion<TType, 2>(t_Expansion(f_scalar).t_Square()).
	 */
	inline CExpansion(TType f_scalar, expansion_internal::TFromScalarSquare_Tag UNUSED(t_tag));

	/**
	 *	@brief destructor; checks the integrity of the expansion
	 *	@note This is empty in release.
	 */
	~CExpansion(); // empty in release

	/**
	 *	@brief expansion element accessor
	 *
	 *	@param[in] n_index is zero-based index of the queried element
	 *
	 *	@return Returns the element with the specified index.
	 *
	 *	@note The elements of expansions are sorted in ascending order, with possible interspersed zeroes.
	 *	@note The last element of an expansion may not reflect its value accurately; use f_Approximate().
	 *	@note Writing the elements of expansions is not permitted intentionally, as it is error-prone.
	 *		To do that, one can either use CExpansion(unsigned int,TType*,TFromArray_Tag) or
	 *		expansion_internal::CExpansionTraits::p_Data().
	 */
	inline TType operator [](unsigned int n_index) const;

	/**
	 *	@brief prints the expansion to stdout, as a sum of hexadecimal float components
	 */
	inline void Dump() const;

	/**
	 *	@brief prints the expansion to stdout, as a single long hexadecimal float number
	 *
	 *	@note This function throws std::bad_alloc.
	 *	@note Special numbers are not handled the usual way (e.g. zero is "0x0.0p0").
	 *	@note Denormals are not handled correctly.
	 */
	inline void Dump_Long() const; // throw(std::bad_alloc)

	/**
	 *	@brief calculates the width of the joint mantissa of the expansion, in bits
	 *	@return Returns the width of the joint mantissa in the expansion,
	 *		in bits (as when concatenated into a single number).
	 */
	inline unsigned int n_Mantissa_Bit_Num() const;

	/**
	 *	@brief calculates the sum of widths of mantissas in the expansion, in bits
	 *	@return Returns the sum of widths of mantissas in the expansion,
	 *		in bits (as when concatenated into a single number).
	 */
	inline unsigned int n_Mantissa_RepresentedBit_Num() const;

	/**
	 *	@brief gets the size of this expansion
	 *	@return Returns the size of this expansion (equal to _n_max_coeff_num unless this expansion is dynamic).
	 */
	inline unsigned int n_Size() const;

	/**
	 *	@brief casts this expansion to dynamic, implying zero elimination in operations with it
	 *	@return Returns this expansion, cast as a dynamic one.
	 *	@note Note that calling this on a dynamic expansion is a no-op (provided the return value optimization is enabled).
	 */
	inline typename expansion_internal::CDynamicCastResult<CSelf, CDynamicSelf>::CResult t_AsDynamic() const;

	/**
	 *	@brief conditionally compresses the expansion inpace, so that it is shorter than n_max_fully_significant_coeff_num
	 *	@note If this expansion is static, this consolidates zeroes at the end
	 *		of the expansion (provided that a shorter representation exists).
	 */
	inline void Compact();

	/**
	 *	@brief conditionally compresses the expansion, so that it is shorter than n_max_fully_significant_coeff_num
	 *
	 *	@return Returns an expansion of the same value, with its size bound by n_max_fully_significant_coeff_num.
	 *
	 *	@note If this expansion is static then CCompactSelf is also static.
	 *	@note If this expansion is static, this consolidates zeroes at the end
	 *		of the expansion (provided that a shorter representation exists).
	 */
	inline CCompactSelf t_Compact() const;

	// t_odo - think about void DropZeros() const and CDynamicSelf t_DropZeros() const as a cheaper variant to t_Compressed() // nah, the bounds on that are too lose (if exponent range is 2048 bits and mantissa is 53 bits, then the expansion with zeroes dropped can be 2048+52 long in the worst case)

	/**
	 *	@brief compresses the expansion inpace
	 *
	 *	@note If this expansion is static, this consolidates zeroes at the end
	 *		of the expansion (provided that a shorter representation exists).
	 *	@note After the compression, any expansion is only up to n_max_fully_significant_coeff_num long.
	 */
	inline void Compress();

	/**
	 *	@brief compresses the expansion
	 *
	 *	@return Returns the shortest expansion of the same value as this.
	 *
	 *	@note The result of this operation is always a dynamic expansion.
	 *	@note After the compression, any expansion is only up to n_max_fully_significant_coeff_num long.
	 */
	inline CDynamicCompactSelf t_Compressed() const; // probably the only function which does not increase size of the expansion

	/**
	 *	@brief cast to an expansion of a different type
	 *
	 *	To cast to double, one can do any of:
	 *	@code
	 *	this->t_Cast(.0);
	 *	this->t_Cast(double());
	 *	double x = ...; // not const double!
	 *	this->t_Cast(x);
	 *	@endcode
	 *
	 *	@param[in] t_tag is new scalar type tag (msvc 6.0 requires a tag; value unused at runtime)
	 *
	 *	@return Returns expansion equal to this, using the selected type.
	 *
	 *	@note The casts must be approached with great care, and in general the only seemingly
	 *		safe cast is from a smaller IEEE floating point type to a larger IEEE floating point type.
	 *	@note Note that casting to an expansion of the same type as this is a no-op
	 *		(provided the return value optimization is enabled).
	 */
	template <class TNewScalar>
	inline CExpansion<TNewScalar, _n_max_coeff_num, _b_allow_dynamic> t_Cast(TNewScalar UNUSED(t_tag)) const // msvc 6.0 requires a tag
	{
		_ASSERTE(TNewScalar(0.5) > 0); // make sure that thit is a floating-point type
		_ASSERTE(sizeof(TNewScalar) >= sizeof(T));
		// only casts from float to double can be done this way, casting the other way would
		// both lose precision and also likely violate the non-overlapping property. that would
		// need to be done using distillation

		typedef expansion_internal::CTypeCastResult<CSelf, T, TNewScalar>::CImplementation CAlgorithm;
		return CAlgorithm::t_Run(*this);
		// run the cast algorithm

		/*CExpansion<TNewScalar, _n_max_coeff_num, _b_allow_dynamic> other;
		for(unsigned int i = 0, n = other.m_n_size = m_n_size; i < n; ++ i)
			other.m_p_value[i] = TNewScalar(m_p_value[i]);
		return other;*/
		// old version of the algorithm
	}

	/**
	 *	@brief calculates an approximate value
	 *	@return Returns an approximate value of this.
	 */
	inline T f_Approximate() const;

	// t_odo - add fabs()
	// t_odo - add comparison operators (most likely needs to be implemented using subtraction, as there are many arbitrary ways to represent the same value using expansions)

	/**
	 *	This can be now used to calculate long sums as follows:
	 *	@code
	 *	const double numbers[] = {M_PI, 2 * M_PI, 3 * M_PI, 4 * M_PI,
	 *		5 * M_PI, 6 * M_PI, 7 * M_PI, 8 * M_PI, 9 * M_PI, 10 * M_PI}; // some numbers to sum up
	 *
	 *	CExpansion<double, 2> accum(0, expansion::from_scalar_zero_pad);
	 *	for(int i = 0; i < 10; ++ i) {
	 *		CExpansion<double, 3, true> accum_new((accum + numbers[i]).t_Compressed());
	 *		// accumulate and compress
	 *
	 *		accum_new.Trim(min(accum_new.n_Size(), 2U));
	 *		// trim to two elements (the caller acknowledges a possible loss of precision)
	 *
	 *		accum = accum_new;
	 *		// assign back, the assignment succeeds as the new expansion is never longer than two elements
	 *	}
	 *	@endcode
	 *	or like this:
	 *	@code
	 *	CExpansion<double, 2, true> accum(0, expansion::from_scalar);
	 *	for(int i = 0; i < 10; ++ i) {
	 *		CExpansion<double, 3, true> accum_new2(accum2 + numbers[i]);
	 *		// accumulate with zero elimination
	 *
	 *		accum_new2.Trim(min(accum_new2.n_Size(), 2U));
	 *		// trim to two elements (the caller acknowledges a possible loss of precision)
	 *
	 *		accum2 = accum_new2;
	 *		// assign back, the assignment succeeds as the new expansion is never longer than two elements
	 *	}
	 *	@endcode
	 *
	 *	@note Also see ::t_Trimmed().
	 */
	void Trim(unsigned int n_size); // clear entries n_size till the last one to zero; if dynamic size is allowed, it is simply changed to n_size

	// t_odo - write .Trim() which sets the expansion size (expects compression, or maybe as an argument), which again could be implemented using CTrimExpansionExpr and making a CExpansionTraits for it
	// now i'm not sure what would that be used for

	// t_odo - write .From(), or a constructor tag from_expansion, which implements expansion conversion (dynamic to non-dynamic)
	// did not write that finally, just wrote a copy-constructor with conversions, without a tag

	/**
	 *	@brief unary minus operator
	 *	@return Returns an expansion with the opposite sign.
	 */
	CSelf operator -() const;

	/**
	 *	@brief absolute value
	 *	@return Returns an expansion with the absolute value.
	 */
	CSelf t_Abs() const;

	/**
	 *	@brief mutliplies the expansion by a power of two factor
	 *
	 *	@param[in] n_power_of_two_factor is the factor, must be a power of two, can be negative
	 *
	 *	@return Returns *this * n_power_of_two_factor.
	 *
	 *	@note Multiplication by a power of two amounts to merely changing the exponent,
	 *		which is used here. Over/underflow to infinity is not handled in any way.
	 *	@note As opposed to *this * TType(n_power_of_two_factor), the size of the result expansion is not increased.
	 */
	CSelf t_Multiply_by_POT(int n_power_of_two_factor) const;

	/**
	 *	@brief mutliplies the expansion by a power of two factor inplace
	 *
	 *	@param[in] n_power_of_two_factor is the factor, must be a power of two, can be negative
	 *
	 *	@note Multiplication by a power of two amounts to merely changing the exponent,
	 *		which is used here. Over/underflow to infinity is not handled in any way.
	 *	@note As opposed to *this * TType(n_power_of_two_factor), the size of the result expansion is not increased.
	 */
	void Multiply_by_POT(int n_power_of_two_factor); // can work in-place

	/**
	 *	@brief divides the expansion by a power of two factor
	 *
	 *	@param[in] n_power_of_two_factor is the factor, must be a power of two, can be negative
	 *
	 *	@return Returns *this / n_power_of_two_factor.
	 *
	 *	@note Division by a power of two amounts to merely changing the exponent,
	 *		which is used here. Over/underflow to infinity is not handled in any way.
	 *	@note Division by all other values cannot be done exactly, and is implemented in CApproximate::t_Divide().
	 */
	CSelf t_Divide_by_POT(int n_power_of_two_factor) const;

	/**
	 *	@brief divides the expansion by a power of two factor inplace
	 *
	 *	@param[in] n_power_of_two_factor is the factor, must be a power of two, can be negative
	 *
	 *	@note Division by a power of two amounts to merely changing the exponent,
	 *		which is used here. Over/underflow to infinity is not handled in any way.
	 *	@note Division by all other values cannot be done exactly, and is implemented in CApproximate::t_Divide().
	 */
	void Divide_by_POT(int n_power_of_two_factor); // can work in-place

	/**
	 *	@brief calculates square of the expansion
	 *
	 *	@return Returns square of this.
	 *
	 *	@note This is slightly faster than doing *this * *this, as there are some redundant operations.
	 *	@note For expansions of size 2, there is an optimized algorithm returning an expansion of size 6 (instead of 8).
	 */
	CExpansion<T, expansion_internal::CSquareCoeffNum<_n_max_coeff_num>::n_result, _b_allow_dynamic> t_Square() const;
	// t_odo - make a template that calculates the size of the expansion to make this shorter
	/*typename expansion_internal::CMultiplicationResult<CSelf, CSelf>::CResult*/ // msvc 6.0 won't parse this
	/*typename expansion_internal::CSquareResult<CSelf>::CResult*/ // msvc 6.0 won't parse this either

	/*inline bool operator !() const // this is safer than the more generic operator bool(), and this is what is mostly needed, anyway
	{
		return f_Approximate() == T(0);
	}*/

	/**
	 *	@brief conversion to bool
	 *	@return Returns nonzero (not 1) if this does not equal to zero, otherwise returns null.
	 *	@note This uses the safe bool idiom to avoid mixing expansions in unsafe arithmetic expressions.
	 */
	inline operator bool_type() const
	{
		return (f_Approximate() != T(0))? &CSelf::dummy_bool_type_value : (bool_type)0;
	}

	//inline operator bool() const;

	/**
	 *	@brief less-than comparison
	 *	@tparam CRightExpansion is type on the right (may be another expansion, or TType)
	 *	@param[in] r_right is value of the operand on the right
	 *	@return Returns true if the value of this < the value of r_right.
	 *	@note This uses \f$6mn\f$ FLOPs where \f$m\f$ is size of this and \f$n\f$ is size of r_right.
	 */
	template <class CRightExpansion>
	inline bool operator <(const CRightExpansion &r_right) const
	{
		return (*this - r_right).f_Approximate() < 0; // probably not the fastest way that works
	}

	/**
	 *	@brief greater-than comparison
	 *	@tparam CRightExpansion is type on the right (may be another expansion, or TType)
	 *	@param[in] r_right is value of the operand on the right
	 *	@return Returns true if the value of this > the value of r_right.
	 *	@note This uses \f$6mn\f$ FLOPs where \f$m\f$ is size of this and \f$n\f$ is size of r_right.
	 */
	template <class CRightExpansion>
	inline bool operator >(const CRightExpansion &r_right) const
	{
		return (*this - r_right).f_Approximate() > 0; // probably not the fastest way that works
		//return (r_right < *this);// wont work in case r_right is a scalar
	}

	/**
	 *	@brief equality comparison
	 *	@tparam CRightExpansion is type on the right (may be another expansion, or TType)
	 *	@param[in] r_right is value of the operand on the right
	 *	@return Returns true if the value of this == the value of r_right.
	 *	@note This uses \f$6mn\f$ FLOPs where \f$m\f$ is size of this and \f$n\f$ is size of r_right.
	 */
	template <class CRightExpansion>
	inline bool operator ==(const CRightExpansion &r_right) const
	{
		return (*this - r_right).f_Approximate() == 0; // probably not the fastest way that works
	}

	/**
	 *	@brief less than or equal comparison
	 *	@tparam CRightExpansion is type on the right (may be another expansion, or TType)
	 *	@param[in] r_right is value of the operand on the right
	 *	@return Returns true if the value of this <= the value of r_right.
	 *	@note This uses \f$6mn\f$ FLOPs where \f$m\f$ is size of this and \f$n\f$ is size of r_right.
	 */
	template <class CRightExpansion>
	inline bool operator <=(const CRightExpansion &r_right) const
	{
		return !(*this > r_right);
	}

	/**
	 *	@brief greater than or equal comparison
	 *	@tparam CRightExpansion is type on the right (may be another expansion, or TType)
	 *	@param[in] r_right is value of the operand on the right
	 *	@return Returns true if the value of this >= the value of r_right.
	 *	@note This uses \f$6mn\f$ FLOPs where \f$m\f$ is size of this and \f$n\f$ is size of r_right.
	 */
	template <class CRightExpansion>
	inline bool operator >=(const CRightExpansion &r_right) const
	{
		return !(*this < r_right);
	}

	/**
	 *	@brief inequality comparison
	 *	@tparam CRightExpansion is type on the right (may be another expansion, or TType)
	 *	@param[in] r_right is value of the operand on the right
	 *	@return Returns true if the value of this != the value of r_right.
	 *	@note This uses \f$6mn\f$ FLOPs where \f$m\f$ is size of this and \f$n\f$ is size of r_right.
	 */
	template <class CRightExpansion>
	inline bool operator !=(const CRightExpansion &r_right) const
	{
		return !(*this == r_right);
	}

	/**
	 *	@brief expansion addition
	 *	@tparam CRightExpansion is type on the right (may be another expansion, or TType)
	 *	@param[in] r_right is value of the operand on the right
	 *	@return Returns value of the sum.
	 *	@note This uses \f$6mn\f$ FLOPs where \f$m\f$ is size of this and \f$n\f$ is size
	 *		of r_right, the size of the output is \f$m + n\f$,
	 */
	template <class CRightExpansion>
	inline typename expansion_internal::CAdditionResult<CSelf, CRightExpansion/*, (b_allow_dynamic ||
		expansion_internal::CExpansionTraits<CRightExpansion>::b_allow_dynamic)*/>::CResult // t_odo - see if msvc 6.0 compiles this // no, problems ensue
		operator +(const CRightExpansion &r_right) const
	{
		typedef typename expansion_internal::CAdditionResult<CSelf, CRightExpansion/*, (b_allow_dynamic ||
			expansion_internal::CExpansionTraits<CRightExpansion>::b_allow_dynamic)*/> CAlgorithm;
		typedef typename CAlgorithm::CResult CResult;
		CResult result(expansion_internal::uninitialized_construct);

		unsigned int n_result_size = CAlgorithm::CImplementation::n_Run(
			expansion_internal::CExpansionTraits<CResult>::p_Data(result), m_n_size, m_p_value,
			expansion_internal::CExpansionTraits<CRightExpansion>::n_Size(r_right), // neet to go through the traits, as CRightExpansion::b_allow_dynamic is not available if CRightExpansion is TType, and SFINAE guarantees no comperhensible error message is output in such case
			expansion_internal::CExpansionTraits<CRightExpansion>::p_Data(r_right));
		result.Uninitialized_Contract(n_result_size);
		// add two expansions

		return result;
	}

	/**
	 *	@brief expansion subtraction
	 *	@tparam CRightExpansion is type on the right (may be another expansion, or TType)
	 *	@param[in] r_right is value of the operand on the right
	 *	@return Returns value of the difference.
	 *	@note This uses \f$6mn\f$ FLOPs where \f$m\f$ is size of this and \f$n\f$ is size
	 *		of r_right, the size of the output is \f$m + n\f$,
	 */
	template <class CRightExpansion>
	inline typename expansion_internal::CSubtractionResult<CSelf, CRightExpansion>::CResult
		operator -(const CRightExpansion &r_right) const
	{
		typedef typename expansion_internal::CSubtractionResult<CSelf, CRightExpansion> CAlgorithm;
		typedef typename CAlgorithm::CResult CResult;
		CResult result(expansion_internal::uninitialized_construct);

		unsigned int n_result_size = CAlgorithm::CImplementation::n_Run(
			expansion_internal::CExpansionTraits<CResult>::p_Data(result), m_n_size, m_p_value,
			expansion_internal::CExpansionTraits<CRightExpansion>::n_Size(r_right),
			expansion_internal::CExpansionTraits<CRightExpansion>::p_Data(r_right));
		result.Uninitialized_Contract(n_result_size);
		// subtract two expansions

		return result;
	}

	/**
	 *	@brief expansion multiplication
	 *	@tparam CRightExpansion is type on the right (may be another expansion, or TType)
	 *	@param[in] r_right is value of the operand on the right
	 *	@return Returns value of the product.
	 *	@note The size of the output is \f$m * n\f$, where \f$m\f$ is size of this and
	 *		\f$n\f$ is size of r_right,
	 */
	template <class CRightExpansion>
	inline typename expansion_internal::CMultiplicationResult<CSelf, CRightExpansion>::CResult
		operator *(const CRightExpansion &r_right) const
	{
		typedef typename expansion_internal::CMultiplicationResult<CSelf, CRightExpansion> CAlgorithm;
		typedef typename CAlgorithm::CResult CResult;
		CResult result(expansion_internal::uninitialized_construct);

		/*int n_rdyn = expansion_internal::CExpansionTraits<CRightExpansion>::b_allow_dynamic;
		int ld = CAlgorithm::b_left_dynamic, rd = CAlgorithm::b_right_dynamic;
		int n_dyn = CAlgorithm::b_left_or_right_allows_dynamic,
			n_sw = CAlgorithm::b_need_swap, n_l = CAlgorithm::n_left_coeff_num,
			n_r = CAlgorithm::n_right_coeff_num;*/
		// debug

		unsigned int n_result_size = CAlgorithm::CImplementation::n_Run(
			expansion_internal::CExpansionTraits<CResult>::p_Data(result), m_n_size, m_p_value,
			expansion_internal::CExpansionTraits<CRightExpansion>::n_Size(r_right),
			expansion_internal::CExpansionTraits<CRightExpansion>::p_Data(r_right));
		result.Uninitialized_Contract(n_result_size);
		// multiply two expansions

		return result;
	}

protected:
	/**
	 *	@brief safe bool idiom; has no effect
	 */
	void dummy_bool_type_value() const {}

#if defined(_MSC_VER) && !defined(__MWERKS__) && _MSC_VER <= 1200
private:
	// element accessor for use by CExpansionTraits in msvc 6.0
	// note that this can't be operator [](), as it would shadow the const one and disallow users to read expansions
	inline TType &f_At(unsigned int n_index);

public: // msvc 6.0 unfortunately can't befriend all template specializations and the code below needs to be left public
#else // _MSC_VER && !__MWERKS__ && _MSC_VER <= 1200
protected: // other compilers can friend and this can be protected
#endif // _MSC_VER && !__MWERKS__ && _MSC_VER <= 1200
	/**
	 *	@brief default constructor
	 *	@note This is not allowed to be used from the outside, as it leaves the expansion
	 *		uninitialized, impending an assertion failure from its destructor.
	 */
	inline CExpansion(expansion_internal::TUninitializedConstruct_Tag UNUSED(t_tag));

	/**
	 *	@brief contracts the expansion to a given size
	 *
	 *	Clear entries n_size till the last one to zero. If dynamic size
	 *	is allowed, it is simply changed to n_size.
	 *
	 *	@param[in] n_size is required target size of the expansion
	 *	@note This potentially throws away the most significant value of the expansion,
	 *		and can therefore only be done safely from the functions of CExpanasion.
	 */
	void Uninitialized_Contract(unsigned int n_size);

};

#if defined(_MSC_VER) && !defined(__MWERKS__) && _MSC_VER <= 1200

/**
 *	@brief converts an expansion to approximate value
 *	@tparam TExpansion is specialization of CExpansion
 *	@param[in] r_t_expansion is expansion to be approximated
 *	@return Returns an approximate value of the input.
 */
template <class TExpansion>
inline typename expansion_internal::CExpansionTraits<TExpansion>::TType
	f_Approximate(const TExpansion &r_t_expansion)
{
	return r_t_expansion.f_Approximate();
}

/**
 *	@brief converts a scalar to approximate value (no-op)
 *	@param[in] f_value is approximate value to be approximated
 *	@return Returns value of the input.
 */
inline long double f_Approximate(long double f_value)
{
	return f_value;
}

/**
 *	@brief converts a scalar to approximate value (no-op)
 *	@param[in] f_value is approximate value to be approximated
 *	@return Returns value of the input.
 */
inline double f_Approximate(double f_value)
{
	return f_value;
}

/**
 *	@brief converts a scalar to approximate value (no-op)
 *	@param[in] f_value is approximate value to be approximated
 *	@return Returns value of the input.
 */
inline float f_Approximate(float f_value)
{
	return f_value;
}

/**
 *	@brief converts an expansion to an expansion (no-op)
 *	@tparam TExpansion is specialization of CExpansion
 *	@param[in] r_t_expansion is expansion to be approximated
 *	@return Returns value of the input.
 */
template <class TExpansion>
inline TExpansion t_Expansion(const TExpansion &r_t_expansion)
{
	return r_t_expansion;
}

#if 0 // can't use long double, don't have CFloatRepresentationTraits<long double>

/**
 *	@brief converts a scalar to an expansion
 *	@param[in] f_scalar is value to be converted
 *	@return Returns value of the input as size 1 expansion.
 */
inline CExpansion<long double, 1> t_Expansion(long double f_scalar)
{
	return CExpansion<long double, 1>(f_scalar, expansion::from_scalar);
}

#endif // 0

/**
 *	@brief converts a scalar to an expansion
 *	@param[in] f_scalar is value to be converted
 *	@return Returns value of the input as size 1 expansion.
 */
inline CExpansion<double, 1> t_Expansion(double f_scalar)
{
	return CExpansion<double, 1>(f_scalar, expansion::from_scalar);
}

/**
 *	@brief converts a scalar to an expansion
 *	@param[in] f_scalar is value to be converted
 *	@return Returns value of the input as size 1 expansion.
 */
inline CExpansion<float, 1> t_Expansion(float f_scalar)
{
	return CExpansion<float, 1>(f_scalar, expansion::from_scalar);
}

#else // _MSC_VER && !__MWERKS__ && _MSC_VER <= 1200

/**
 *	@brief converts an expansion to approximate value
 *
 *	@tparam T is scalar type used by the expansion
 *	@tparam _n_max_coeff_num is maximum size of the expansion
 *	@tparam _b_allow_dynamic is dynamic flag of the expansion
 *
 *	@param[in] r_t_expansion is expansion to be approximated
 *
 *	@return Returns an approximate value of the input.
 */
template <class T, const unsigned int _n_max_coeff_num, const bool _b_allow_dynamic>
inline T f_Approximate(const CExpansion<T, _n_max_coeff_num, _b_allow_dynamic> &r_t_expansion) // msvc 6.0 can't parse this
{
	return r_t_expansion.f_Approximate();
}

/**
 *	@brief converts a scalar to approximate value (no-op)
 *	@tparam TScalar is type of a scalar
 *	@param[in] f_value is approximate value to be approximated
 *	@return Returns value of the input.
 */
template <class TScalar>
inline TScalar f_Approximate(TScalar f_value)
{
	return f_value;
}

/**
 *	@brief converts an expansion to an expansion (no-op)
 *
 *	@tparam T is scalar type used by the expansion
 *	@tparam _n_max_coeff_num is maximum size of the expansion
 *	@tparam _b_allow_dynamic is dynamic flag of the expansion
 *
 *	@param[in] r_t_expansion is expansion to be approximated
 *
 *	@return Returns value of the input.
 */
template <class T, const unsigned int _n_max_coeff_num, const bool _b_allow_dynamic>
inline CExpansion<T, _n_max_coeff_num, _b_allow_dynamic> t_Expansion(const CExpansion<T, _n_max_coeff_num, _b_allow_dynamic> &r_t_expansion) // msvc 6.0 can't parse this
{
	return r_t_expansion;
}

/**
 *	@brief converts a scalar to an expansion
 *	@tparam TScalar is type of a scalar
 *	@param[in] f_scalar is value to be converted
 *	@return Returns value of the input as size 1 expansion.
 */
template <class TScalar>
inline CExpansion<TScalar, 1> t_Expansion(TScalar f_scalar)
{
	return CExpansion<TScalar, 1>(f_scalar, expansion::from_scalar);
}

//template <class TScalar>
//inline CExpansion<TScalar, 1> t_Precise(TScalar f_scalar) // just a different name
//{
//	return t_Expansion(f_scalar);
//}

#endif // _MSC_VER && !__MWERKS__ && _MSC_VER <= 1200

/**
 *	@brief trims an expansions to a given size
 *
 *	Can use this as follows:
 *	@code
 *	const double numbers[] = {M_PI, 2 * M_PI, 3 * M_PI, 4 * M_PI,
 *		5 * M_PI, 6 * M_PI, 7 * M_PI, 8 * M_PI, 9 * M_PI, 10 * M_PI}; // some numbers to sum up
 *
 *	CExpansion<double, 2> accum(0, expansion::from_scalar);
 *	for(int i = 0; i < 10; ++ i) {
 *		accum = t_Trimmed<2, false>(accum.t_AsDynamic() + numbers[i]); // false = do not warn about the loss of precision
 *		// accumulate with zero elimination and trim
 *	}
 *	@endcode
 *	or as:
 *	@code
 *	CExpansion<double, 2> accum(0, expansion::from_scalar_zero_pad);
 *	for(int i = 0; i < 10; ++ i) {
 *		accum = t_Trimmed<2, false>((accum + numbers[i]).t_Compressed()); // false = do not warn about the loss of precision
 *		// accumulate, compress and trim
 *	}
 *	@endcode
 *	or as:
 *	@code
 *	CExpansion<double, 2, true> accum(0, expansion::from_scalar);
 *	for(int i = 0; i < 10; ++ i) {
 *		accum = t_Trimmed<2, false>(accum + numbers[i]); // false = do not warn about the loss of precision
 *		// accumulate with zero elimination and trim
 *	}
 *	@endcode
 *	or to simplify this:
 *	@code
 *	CExpansion<double, 20> long_exp(d * a * a * -27.0 + b * c * a * 9.0 - b * b * 2.0); // a long complex expansion
 *	CExpansion<double, 20 * 20 * 2> squared(long_exp.t_Square()); // becomes untractable
 *	@endcode
 *	into this:
 *	@code
 *	CExpansion<double, 4> long_exp(t_Trimmed<4, true>(d * a * a * -27.0 + b * c * a * 9.0 - b * b * 2.0));
 *	// will rarely need a longer expansion
 *
 *	CExpansion<double, 8> squared(t_Trimmed<8, true>(long_exp.t_Square())); // simple
 *	@endcode
 *
 *	@tparam n_max_size is maximum size of the result (can be less if the input is shorter)
 *	@tparam b_warn_losing_precision is precision loss warning flag (applies in debug; if set, expansions size being smaller than the given size is asserted)
 *	@tparam TExpansion is an expansion (specialization of CExpansion)
 *
 *	@param[in] t_exp is input expansion
 *
 *	@return Returns expansion with maximum compile-time size n_max_size.
 *
 *	@note See also t_Trimmed_CC().
 */
template <const unsigned int n_max_size, const bool b_warn_losing_precision, class TExpansion>
inline CExpansion<typename TExpansion::TType, n_max_size,
	expansion_internal::CExpansionTraits<TExpansion>::b_allow_dynamic> t_Trimmed(const TExpansion &r_t_exp)
{
	enum {
		n_src_max_coeff_num = expansion_internal::CExpansionTraits<TExpansion>::n_max_coeff_num,
		b_useful_trim = n_max_size < n_src_max_coeff_num // to allow zero-trim (n_max_size == n_src_max_coeff_num, we would need a specialized implementation which just returns the input)
	};
	typedef expansion_internal::CStaticAssert<b_useful_trim>::USELESS_TRIM_PERFORMED CAssert0; // make sure that the trim is useful
	// if this triggers, then the selected size is too big and the trim has no effect
	// (other than using more memory and making further computation slower; possibly a typo in trim size)

	CExpansion<MSVC6_OMMIT(typename) TExpansion::TType, n_src_max_coeff_num, true> t_exp(r_t_exp);
	// needs to be dynamic, otherwise the Trim() will not have any effect
	// note that the copy is needed either way, as the expansion is modified in here

	_ASSERTE(!b_warn_losing_precision || t_exp.n_Size() <= n_max_size); // otherwise losing precision
	t_exp.Trim(min(t_exp.n_Size(), n_max_size));
	return t_exp; // the copy-constructor takes care of the rest
}

/**
 *	@brief trims an expansions to a given size with conditional compression
 *
 *	@tparam n_max_size is maximum size of the result (can be less if the input is shorter)
 *	@tparam b_warn_losing_precision is precision loss warning flag (applies in debug; if set, expansions size being smaller than the given size is asserted)
 *	@tparam TExpansion is an expansion (specialization of CExpansion)
 *
 *	@param[in] t_exp is input expansion
 *
 *	@return Returns expansion with maximum compile-time size n_max_size,
 *		which is compressed if the input expansion would not fit the limit.
 *
 *	@note For usage samples, see t_Trimmed().
 *	@note The conditional compression tends to be slower if applied
 *		to expansions which are used for further computation. Use with care.
 */
template <const unsigned int n_max_size, const bool b_warn_losing_precision, class TExpansion>
inline CExpansion<typename TExpansion::TType, n_max_size,
	expansion_internal::CExpansionTraits<TExpansion>::b_allow_dynamic> t_Trimmed_CC(const TExpansion &r_t_exp)
{
	enum {
		n_src_max_coeff_num = expansion_internal::CExpansionTraits<TExpansion>::n_max_coeff_num,
		b_useful_trim = n_max_size < n_src_max_coeff_num // to allow zero-trim (n_max_size == n_src_max_coeff_num, we would need a specialized implementation which just returns the input)
	};
	typedef expansion_internal::CStaticAssert<b_useful_trim>::USELESS_TRIM_PERFORMED CAssert0; // make sure that the trim is useful
	// if this triggers, then the selected size is too big and the trim has no effect
	// (other than using more memory and making further computation slower; possibly a typo in trim size)

	CExpansion<MSVC6_OMMIT(typename) TExpansion::TType, n_src_max_coeff_num, true> t_exp(r_t_exp);
	// needs to be dynamic, otherwise the Trim() will not have any effect
	// note that the copy is needed either way, as the expansion is modified in here

	if(t_exp.n_Size() > n_max_size)
		t_exp.Compress();
	// conditionally compress to make it fit

	_ASSERTE(!b_warn_losing_precision || t_exp.n_Size() <= n_max_size); // otherwise losing precision
	t_exp.Trim(min(t_exp.n_Size(), n_max_size));
	return t_exp; // the copy-constructor takes care of the rest
}

typedef CExpansion<double, 1, false> CExpansion1d; /**< size 1 static double expansion */
typedef CExpansion<double, 2, false> CExpansion2d; /**< size 2 static double expansion */
typedef CExpansion<double, 3, false> CExpansion3d; /**< size 3 static double expansion */
typedef CExpansion<double, 4, false> CExpansion4d; /**< size 4 static double expansion */
typedef CExpansion<double, 5, false> CExpansion5d; /**< size 5 static double expansion */
typedef CExpansion<double, 6, false> CExpansion6d; /**< size 6 static double expansion */
typedef CExpansion<double, 7, false> CExpansion7d; /**< size 7 static double expansion */
typedef CExpansion<double, 8, false> CExpansion8d; /**< size 8 static double expansion */

/**
 *	@brief approximate operations on expansions
 *	@tparam n_result_precision is desired result precision
 */
template <const unsigned int n_result_precision>
class CApproximate {
public:
	/**
	 *	@brief calculates ratio of two numbers
	 *
	 *	@tparam L is type of numerator; may be scalar or CExpansion specialization
	 *	@tparam R is type of denominator; may be scalar or CExpansion specialization
	 *
	 *	@param[in] r_t_numerator is value of numerator
	 *	@param[in] r_t_denominator is value of denominator
	 *
	 *	@return Returns value of the fraction, evaluated to the desired precision.
	 *
	 *	@note Although it can be exact in some cases, division is approximate: the result may not
	 *		be representable by finite number of digits. To see if the division was approximate,
	 *		multiply back (which is exact unles over/underflow occurs), subtract and compare to zero.
	 *	@note This can be used to evaluate precise ratios of scalars (L and R may both be scalars).
	 */
	template <class L, class R>
	static CExpansion<typename expansion_internal::CExpansionTraits<L>::TType, n_result_precision,
		expansion_internal::CExpansionTraits<L>::b_allow_dynamic ||
		expansion_internal::CExpansionTraits<R>::b_allow_dynamic>
		t_Divide(const L &r_t_numerator, const R &r_t_denominator)
	{
		typedef typename expansion_internal::CCompatibleTypeCheck<L, R>::CAssert CTypeCheck; // without CAssert it will not check
		typedef typename expansion_internal::CExpansionTraits<L>::TType TScalar;
		enum {
			b_left_or_right_allows_dynamic = expansion_internal::CExpansionTraits<L>::b_allow_dynamic ||
				expansion_internal::CExpansionTraits<R>::b_allow_dynamic,
			n_denom_precision = expansion_internal::CExpansionTraits<R>::n_max_coeff_num,
			n_max_compressed = CMaxExpansionSize<TScalar>::n_longest_expansion_size,
			n_corection_precision_ = (n_result_precision + 1) * n_denom_precision * 2,
			n_corection_precision = (n_corection_precision_ < n_max_compressed)? n_corection_precision_ : n_max_compressed
		};
		//typedef CExpansion<TScalar, n_result_precision, b_left_or_right_allows_dynamic> CResult; // not required
		typedef CExpansion<TScalar, n_result_precision + 1, true> CIntermediate; // todo - make this dynamic only if b_left_or_right_allows_dynamic, to follow the traditional semantic

		const TScalar f_num = f_Approximate(r_t_numerator), f_den = f_Approximate(r_t_denominator);
		// division done approximately

		CIntermediate t_result(f_num / f_den, expansion::from_scalar);
		// calculate the most significant element of the result expansion

		for(int i = 0; i < n_result_precision; ++ i) {
			// up to thousands, but only affects the stack size. could lead to integral constant overflow in its implementation of t_Square() but that cannot be avoided
			// square size = ((n_result_precision + 1) * n_denom_precision * 2)^2 * 2 =
			//				 ((n_result_precision + 1) * n_denom_precision)^2 * 8 =
			//				 (n_result_precision + 1)^2 * n_denom_precision^2 * 8
			// n_result_precision and n_denom_precision of size >127 will cause 32-bit integer overflow in square size calculation
			// 127:			2081157128
			// int32_max:	2147483647
			// 128:			2147483648
			// could have a specialized branch for very long expansions to first compress them and then divide. the caller can do this, however.

			CExpansion<TScalar, n_corection_precision, true> t_correction((t_result * r_t_denominator).t_Compact()); // msvc 6.0 requires this
			TScalar f_correction = (t_correction - r_t_numerator).f_Approximate();
			TScalar f_result_i = f_correction / f_den;
			if(b_left_or_right_allows_dynamic && !f_result_i)
				break;
			t_result = t_result - f_result_i; // must always fit, trim not required
			//t_result = t_Trimmed<n_result_precision + 1, true>(t_result - f_result_i); // must always fit
		}
		// calculate n_result_precision corrections (one more than the result in order to provide correct rounding)

		t_result.Compress(); // does not reduce the size of the expansion but possibly shifts the bits
		// (depending on how the results of the division are rounded), slightly improving the precision. could be omitted

		t_result.Trim(min(t_result.n_Size(), n_result_precision)); // trims almost always
		// compress to the desired precision

		return t_result;
	}

	/**
	 *	@brief reciprocal square root
	 *
	 *	@tparam E is type of the argument; may be scalar or CExpansion specialization
	 *
	 *	@param[in] r_t_x is value of the argument
	 *	@param[in] n_max_iteration_num is number of Newton iterations to perform
	 *		(three are usually "good enough" for reasonably short expansions)
	 *
	 *	@return Returns one over the value of the square root, evaluated to the desired precision.
	 *
	 *	@note This uses sqrt() as an initial guess, which is then made more precise by Newton iteration.
	 */
	template <class E>
	static CExpansion<typename expansion_internal::CExpansionTraits<E>::TType, n_result_precision,
		expansion_internal::CExpansionTraits<E>::b_allow_dynamic>
		t_RSqrt(const E &r_t_x, const unsigned int n_max_iteration_num = 3)
	{
		typedef typename expansion_internal::CExpansionTraits<E>::TType TScalar;
		enum {
			n_input_precision = expansion_internal::CExpansionTraits<E>::n_max_coeff_num,
			n_max_compressed = CMaxExpansionSize<TScalar>::n_longest_expansion_size,
			b_input_allows_dynamic = expansion_internal::CExpansionTraits<E>::b_allow_dynamic,
			n_square_result_precision_ = expansion_internal::CSquareCoeffNum<n_result_precision>::n_result, // for 4, this is 32
			n_square_result_precision = (n_square_result_precision_ < n_max_compressed)? n_square_result_precision_ : n_max_compressed,
			n_intermediate0_precision_ = n_input_precision * n_square_result_precision * 2 + 1, // this is large even for modest inputs ... 128, 256, or larger are common
			n_intermediate0_precision = (n_intermediate0_precision_ < n_max_compressed)? n_intermediate0_precision_ : n_max_compressed,
			n_intermediate0_precision_1 = (n_intermediate0_precision_ - 1 < n_max_compressed)? n_intermediate0_precision_ - 1 : n_max_compressed, // !!
			n_intermediate1_precision_ = n_intermediate0_precision * n_result_precision * 2, // thousands
			n_intermediate1_precision = (n_intermediate1_precision_ < n_max_compressed)? n_intermediate1_precision_ : n_max_compressed,
			n_intermediate2_precision_ = n_result_precision + n_intermediate1_precision,
			n_intermediate2_precision = (n_intermediate2_precision_ < n_max_compressed)? n_intermediate2_precision_ : n_max_compressed
		};
		//typedef CExpansion<TScalar, n_result_precision, b_input_allows_dynamic> CResult; // not required
		typedef CExpansion<TScalar, n_result_precision, true> CIntermediate; // todo - make this dynamic only if b_left_or_right_allows_dynamic, to follow the traditional semantic

		CIntermediate xi(1 / (TScalar)sqrt(t_Expansion(r_t_x).t_Cast(double()).f_Approximate()), expansion::from_scalar);
		// makes an initial guess, perform the square root in double precision

		CExpansion<TScalar, n_input_precision, b_input_allows_dynamic> t_x_half(t_Expansion(r_t_x).t_Divide_by_POT(2));
		// reuse calculation in Newton iteration

		for(unsigned int i = 0; i < n_max_iteration_num; ++ i) {
			CExpansion<TScalar, n_square_result_precision, true> xi_square(xi.t_Square().t_Compact());
			CExpansion<TScalar, n_intermediate0_precision_1, true> t_x_xi_square_half((t_x_half * xi_square).t_Compact());
			CExpansion<TScalar, n_intermediate0_precision, true> t_x_xi_square_half_minus_half((t_x_xi_square_half - TScalar(.5)).t_Compact());
			CExpansion<TScalar, n_intermediate1_precision, true> t_xi_delta((xi * t_x_xi_square_half_minus_half).t_Compact());
			CExpansion<TScalar, n_intermediate2_precision, true> t_xiplus1((xi - t_xi_delta).t_Compact());
			t_xiplus1.Compress();
			t_xiplus1.Trim(min(t_xiplus1.n_Size(), n_result_precision));
			xi = t_xiplus1;
			// todo - reverse error analysis, reduce storage
		}
		// iterate Newton; note that it is hard to tell how many iterations are required

		return xi;
	}

	/**
	 *	@brief square root
	 *
	 *	@tparam E is type of the argument; may be scalar or CExpansion specialization
	 *
	 *	@param[in] r_t_x is value of the argument
	 *	@param[in] n_max_iteration_num is number of Newton iterations to perform
	 *		(three are usually "good enough" for reasonably short expansions)
	 *
	 *	@return Returns value of the square root, evaluated to the desired precision.
	 *
	 *	@note This uses sqrt() as an initial guess, which is then made more precise by Newton iteration.
	 */
	template <class E>
	static CExpansion<typename expansion_internal::CExpansionTraits<E>::TType, n_result_precision,
		expansion_internal::CExpansionTraits<E>::b_allow_dynamic>
		t_Sqrt(const E &r_t_x, const unsigned int n_max_iteration_num = 3)
	{
		typedef typename expansion_internal::CExpansionTraits<E>::TType TScalar;
		typedef CExpansion<MSVC6_OMMIT(typename) expansion_internal::CExpansionTraits<E>::TType, n_result_precision,
			expansion_internal::CExpansionTraits<E>::b_allow_dynamic> CResult;
		if(!r_t_x)
			return CResult(TScalar(0), expansion::from_scalar_zero_pad); // otherwise calculates a NaN by taking a product 0 * inf
		return t_Trimmed<n_result_precision, false>((t_RSqrt(r_t_x, n_max_iteration_num) * r_t_x).t_Compressed()); // x^(1/2) = x * x^(-1/2)
		// not sure how to guarantee sufficient precision here, I think if r_t_x is one digit, then
		// rsqrt must be expansion_internal digits to produce a result precise to expansion_internal digits
		// multiplication doubles the number of digits, but that might be redundant in some cases,
		// so expansion_internal / 2 digits would be too little. seems to be correct this way, or more precise than needed.
	}

	/**
	 *	@brief reciprocal cube root
	 *
	 *	@tparam E is type of the argument; may be scalar or CExpansion specialization
	 *
	 *	@param[in] r_t_x is value of the argument
	 *	@param[in] n_max_iteration_num is number of Newton iterations to perform
	 *		(three are usually "good enough" for reasonably short expansions)
	 *
	 *	@return Returns one over the value of the cube root, evaluated to the desired precision.
	 *
	 *	@note This uses pow(r_t_x, 1.0 / 3) as an initial guess, which is then made more precise by Newton iteration.
	 */
	template <class E>
	static CExpansion<typename expansion_internal::CExpansionTraits<E>::TType, n_result_precision,
		expansion_internal::CExpansionTraits<E>::b_allow_dynamic>
		t_RCurt(const E &r_t_x, const unsigned int n_max_iteration_num = 3)
	{
		typedef typename expansion_internal::CExpansionTraits<E>::TType TScalar;
		enum {
			n_input_precision = expansion_internal::CExpansionTraits<E>::n_max_coeff_num,
			b_input_allows_dynamic = expansion_internal::CExpansionTraits<E>::b_allow_dynamic,
			n_max_compressed = CMaxExpansionSize<TScalar>::n_longest_expansion_size,

			n_square_result_precision_ = expansion_internal::CSquareCoeffNum<n_result_precision>::n_result, // for 4, this is 32
			n_square_result_precision = (n_square_result_precision_ < n_max_compressed)? n_square_result_precision_ : n_max_compressed,
			n_cube_result_precision_ = n_square_result_precision * n_result_precision * 2, // for 4, this is 256
			n_cube_result_precision = (n_cube_result_precision_ < n_max_compressed)? n_cube_result_precision_ : n_max_compressed,
			n_intermediate0_precision_ = n_input_precision * n_cube_result_precision * 2 + 1, // this is very large even for modest inputs ... 2048, 4096, or larger are common
			n_intermediate0_precision = (n_intermediate0_precision_ < n_max_compressed)? n_intermediate0_precision_ : n_max_compressed,
			n_intermediate0_precision_1 = (n_intermediate0_precision_ - 1 < n_max_compressed)? n_intermediate0_precision_ - 1 : n_max_compressed, // !!
			n_intermediate1_precision_ = n_intermediate0_precision * n_result_precision * 2, // tens of thousands 16384 and so on
			n_intermediate1_precision = (n_intermediate1_precision_ < n_max_compressed)? n_intermediate1_precision_ : n_max_compressed
		};
		//typedef CExpansion<TScalar, n_result_precision, b_input_allows_dynamic> CResult; // not required
		typedef CExpansion<TScalar, n_result_precision, true> CIntermediate; // todo - make this dynamic only if b_left_or_right_allows_dynamic, to follow the traditional semantic

		CIntermediate xi((TScalar)pow(r_t_x.t_Cast(double()).f_Approximate(), -1.0 / 3), expansion::from_scalar);
		// makes an initial guess, perform the cube root in double precision

		/*if(!r_t_x)
			return xi;*/ // seems to be never used
		// infinity

		for(unsigned int i = 0; i < n_max_iteration_num; ++ i) {
			CExpansion<TScalar, n_square_result_precision, true> xi_square(xi.t_Square().t_Compact());
			CExpansion<TScalar, n_cube_result_precision, true> xi_cube((xi * xi_square).t_Compact());
			CExpansion<TScalar, n_intermediate0_precision_1, true> t_x_xi_cube((r_t_x * xi_cube).t_Compact());
			CExpansion<TScalar, n_intermediate0_precision, true> t_x_xi_cube_minus_one((t_x_xi_cube - TScalar(1)).t_Compact());
			CExpansion<TScalar, n_intermediate1_precision, true> t_xi_3delta((xi * t_x_xi_cube_minus_one).t_Compact());
			CExpansion<TScalar, n_result_precision> t_xi_delta(t_Divide(t_xi_3delta, TScalar(3)));
			CExpansion<TScalar, n_result_precision * 2, true> t_xiplus1(xi - t_xi_delta);
			t_xiplus1.Compress();
			t_xiplus1.Trim(min(t_xiplus1.n_Size(), n_result_precision));
			xi = t_xiplus1;
			// todo - reverse error analysis, reduce storage
		}
		// iterate Newton; note that it is hard to tell how many iterations are required

		return xi;
	}

	/**
	 *	@brief cube root
	 *
	 *	@tparam E is type of the argument; may be scalar or CExpansion specialization
	 *
	 *	@param[in] r_t_x is value of the argument
	 *	@param[in] n_max_iteration_num is number of Newton iterations to perform
	 *		(three are usually "good enough" for reasonably short expansions)
	 *
	 *	@return Returns value of the cube root, evaluated to the desired precision.
	 *
	 *	@note This uses pow(r_t_x, 1.0 / 3) as an initial guess, which is then made more precise by Newton iteration.
	 */
	template <class E>
	static CExpansion<typename expansion_internal::CExpansionTraits<E>::TType, n_result_precision,
		expansion_internal::CExpansionTraits<E>::b_allow_dynamic>
		t_Curt(const E &r_t_x, const unsigned int n_max_iteration_num = 3)
	{
		typedef typename expansion_internal::CExpansionTraits<E>::TType TScalar;
		typedef CExpansion<MSVC6_OMMIT(typename) expansion_internal::CExpansionTraits<E>::TType, n_result_precision,
			expansion_internal::CExpansionTraits<E>::b_allow_dynamic> CResult;
		if(!r_t_x)
			return CResult(TScalar(0), expansion::from_scalar_zero_pad); // otherwise calculates NaN by taking a product 0 * inf
		CResult t_reciprocal(t_RCurt(r_t_x, n_max_iteration_num));
		return t_Divide(TScalar(1), t_reciprocal);
	}

	/**
	 *	@brief adds three scalars approximately, recovering three most significant terms
	 *
	 *	@param[out] r_f_low is the third most significant part of the output expansion
	 *	@param[out] r_f_mid is the second most significant part of the output expansion
	 *	@param[out] r_f_high is the most significant part of the output expansion
	 *	@param[in] f_a is the first scalar
	 *	@param[in] f_b is the second scalar
	 *	@param[in] f_c is the third scalar
	 */
	template <class TType>
	static inline void Three_Sum_Approx3(TType &r_f_low, TType &r_f_mid, TType &r_f_high,
		const TType a, const TType b, const TType c)
	{
		TType t1, t2, t3;
		CExactBase<TType>::Two_Sum(t2, t1, a, b);
		CExactBase<TType>::Two_Sum(t3, r_f_high, c, t1);
		CExactBase<TType>::Two_Sum(r_f_low, r_f_mid, t2, t3);
	}

	/**
	 *	@brief adds three scalars approximately, recovering two most significant terms
	 *
	 *	@param[out] r_f_low is the second most significant part of the output expansion
	 *	@param[out] r_f_high is the most significant part of the output expansion
	 *	@param[in] f_a is the first scalar
	 *	@param[in] f_b is the second scalar
	 *	@param[in] f_c is the third scalar
	 */
	template <class TType>
	static inline void Three_Sum_Approx2(TType &r_f_low, TType &r_f_high,
		const TType f_a, const TType f_b, const TType f_c)
	{
		TType t1, t2;
		CExactBase<TType>::Two_Sum(t2, t1, f_a, f_b);
		CExactBase<TType>::Two_Sum(r_f_low, r_f_high, f_c, t1);
		r_f_low += t2;
	}

	// t_odo - t_CubeRoot(const E &r_t_x)
	// todo - CExpansion::t_Power() with (positive, signed for checking) integer exponent
	// todo - t_RNthRoot(const E &r_t_x, int n_power)
	// todo - exp, log, sin, cos, asin, acos, atan
};

#include "Exact2.inl"

#endif // !__EXACT_MULTIPRECISION_FLOATING_POINT_ARITHMETIC_INCLUDED
