#include "stdafx.h"
#include "DBLMath.h"	// For almosteq.
//#include "UtilMisc.h"		// For POW2.
#include "Bhaskara.h"	// For FCBhaskara.
#include "Cardano.h"	// This header.

#ifdef _DEBUG
#define new DEBUG_NEW
#undef THIS_FILE
static char THIS_FILE[] = __FILE__;
#endif

inline BOOL figeq (const dblcompex& c1, const dblcompex& c2, int iFig)
{
	return figeq(c1.real(), c2.real(), iFig) && figeq(c1.imag(), c2.imag(), iFig);
} // figeq (const dblcompex& c1, const dblcompex& c2).

/*************************************************************************
 * <_xZkqqHx_iKnGjAjE_xHqqkZx_>	FCCardano::Set
 *
 * <_xZkqqHx_iLeAjEfM_xHqqkZx_>	_xZkqqHx_jFpLjCpGiOkOiCpAjAnNjCoIiClHiCoJiBeC_xHqqkZx_
 *
 * <_xZkqqHx_iIpIjAjE_xHqqkZx_>	a	:_xZkqqHx_iCfCiOjPiNiAiCmMiMfHjAjEiCpAiOhHjCoIiBeC_xHqqkZx_
 *			b	:_xZkqqHx_iCfBiOjPiNiAiCmMiMfHjAjEiCpAiOhHjCoIiBeC_xHqqkZx_
 *			c	:_xZkqqHx_iCfAiOjPiNiAiCmMiMfHjAjEiCpAiOhHjCoIiBeC_xHqqkZx_
 *			d	:_xZkqqHx_iCePiOjPiNiAiCmMiMfHjAjEiCpAiOhHjCoIiBeC_xHqqkZx_
 *
 * <_xZkqqHx_iJpAjAoA_xHqqkZx_>	_xZkqqHx_jFmPjAjE_xHqqkZx_x_xZkqqHx_iCmJjBmOiClHiCoJiCfCiOjPjFpLjCpGiOkO_xHqqkZx_ ax^3+bx^2+cx+d=0 _xZkqqHx_iCmMiKgFiMfHjAjEiCpAjAnNjCoIiClHiCoJiBeC_xHqqkZx_
 *
 * <_xZkqqHx_jHjKjHpA_xHqqkZx_>	05.04.01 Fukushiro M. _xZkqqHx_iNoMjAkM_xHqqkZx_
 *************************************************************************/
void FCCardano::Set (double a, double b, double c, double d)
{
	m_dA = a;
	m_dB = b;
	m_dC = c;
	m_dD = d;
} // FCCardano::Set.

/*************************************************************************
 * <_xZkqqHx_iKnGjAjE_xHqqkZx_>	FCCardano::Solve
 *
 * <_xZkqqHx_iLeAjEfM_xHqqkZx_>	_xZkqqHx_jFpLjCpGiOkOiCmMiJpAiCpAiLiBiCnPiCoJiBeC_xHqqkZx_
 *
 * <_xZkqqHx_jFnEjCgM_xHqqkZx_>	_xZkqqHx_iJpAiCmMjAjEiBeC_xHqqkZx_
 *
 * <_xZkqqHx_iJpAjAoA_xHqqkZx_>	Set_xZkqqHx_iKnGjAjEiCmFjAnNjCoIiClDiCoKiClNjFpLjCpGiOkOiCpAiJpAiCkNiBeC_xHqqkZx_\doc\_xZkqqHx_iDhIiDfHiDeHiLmIjApM_xHqqkZx_\_xZkqqHx_iCmM_xHqqkZx_
 *			_xZkqqHx_iDeKiDiLiDfPiDgNjGeAiKnGjIeBiOjBjHlPiCpAiOfBiPmGiBeC_xHqqkZx_
 *
 * <_xZkqqHx_jHjKjHpA_xHqqkZx_>	05.04.01 Fukushiro M. _xZkqqHx_iNoMjAkM_xHqqkZx_
 *************************************************************************/
long FCCardano::Solve ()
{
	if (almost0(m_dA))
	//----- a=0_xZkqqHx_iCmMiPoKiNiH_xHqqkZx_ -----
	{
		// b*x^2 + c*x + d = 0 _xZkqqHx_iCmGiCkCiCkEiOjPjFpLjCpGiOkOiCmGiCmIiCoJiCmMiCmFiBeB_xHqqkZx_
		// _xZkqqHx_iDgPiBfLiDfIiDeKiDiJiCmMiCfBiOjPjFpLjCpGiOkOiJpAjGeAiCmFiJpAiCkNiBeC_xHqqkZx_
		FCBhaskara bh;
		bh.Set(m_dB, m_dC, m_dD);
		const long lCount = bh.Solve();
		for (long lC = 0; lC != lCount; lC++)
			m_cSol[lC] = bh.m_cSol[lC];
		return lCount;
	}
	// _xZkqqHx_iCfCiOjPiNiAiCpAiCfAiCmGiClFiCmEiBeBiImIiJlKiCmMiMgAiCmJjFmPiMgAiClHiCoJiBeC_xHqqkZx_
	// x^3 + dA * x^2 + dB * x + dC = 0
	const double dA = m_dB / m_dA;
	const double dB = m_dC / m_dA;
	const double dC = m_dD / m_dA;

//----- 05.04.13 Fukushiro M. _xZkqqHx_jFmPiNfIjBeP_xHqqkZx_ ()-----
//	double d = dB - POW2(dA) / 3.0;
//	double e = 2.0 * POW3(dA) / 27.0 - dA * dB / 3.0 + dC;
//	double tmp = POW2(e) / 4 + POW3(d) / 27;
//----- 05.04.13 Fukushiro M. _xZkqqHx_jFmPiNfIiMoD_xHqqkZx_ ()-----
	// _xZkqqHx_iImIiJlKiCmMiPiIjHjNiCmFjBoFiCkLiCmIiMoLiNlHiCkKjEkNjAlGiClFiCmEiBeB_xHqqkZx_bc == xxx _xZkqqHx_iCmMjEoEiKhCiCmFiLlGiCkCiCkKjAlGiClGiClNiBeC_xHqqkZx_
	// _xZkqqHx_iMoLiNlHiCpAiPgPjHiIiCoJiClOiCkPiMlIiCoHiClHiClNiCnPiBeBjHnNiPoGiCmMjBePiCmJiKiEiCmBiCmEiCkIiCkNiCoGiCkEiCmJjFmPiNfIiBeC_xHqqkZx_
	// double d = dB - POW2(dA) / 3.0;
	// double e = 2.0 * POW3(dA) / 27.0 - dA * dB / 3.0 + dC;
	// double tmp = POW2(e) / 4.0 + POW3(d) / 27.0;
	double d = dB - POW2(dA / 1.7320508075688772935274463415059);
	double e = 2.0 * POW3(dA / 3.0) - dA * dB / 3.0 + dC;
	double tmp = POW2(e / 2.0) + POW3(d / 3.0);
//----- 05.04.13 Fukushiro M. _xZkqqHx_jFmPiNfIiPeJ_xHqqkZx_ ()-----
	dblcompex b;
	dblcompex c;
	if (tmp < 0.0)
	{
		tmp = sqrt(-tmp);
		b = pow(dblcompex(-e / 2.0, tmp), 1.0 / 3.0);
		c = pow(dblcompex(-e / 2.0, -tmp), 1.0 / 3.0);
	} else
	{
		tmp = sqrt(tmp);
		double b3 = -e / 2.0 + tmp;
		double c3 = -e / 2.0 - tmp;
		// pow(r, 1/3), pow(complex(r,0), 1/3)_xZkqqHx_iCmGiCoAiCmJiBeB_xHqqkZx_
		// r_xZkqqHx_iCkKjFiJiCmMiOmAjAjEiCmMiPoKiNiHiCmNiDeHiDiJiBfLjCgMiCpAjFnEiClHiBeC_xHqqkZx_
		// _xZkqqHx_iClBiCmMiClNiCnPiBeB_xHqqkZx_r_xZkqqHx_iCkKjFiJiCmMiPoKiNiHiCmNjFiEiNiGiCpAjFoCjAlDiClHiCoJiBeC_xHqqkZx_
		if (b3 < 0.0)
			b = dblcompex(-pow(-b3, 1.0 / 3.0), 0.0);
		else
			b = dblcompex(pow(b3, 1.0 / 3.0), 0.0);
		if (c3 < 0.0)
			c = dblcompex(-pow(-c3, 1.0 / 3.0), 0.0);
		else
			c = dblcompex(pow(c3, 1.0 / 3.0), 0.0);
	}
	dblcompex bc = b * c;
	dblcompex xxx(-d / 3.0, 0.0);

	dblcompex r(-0.5, 0.86602540378443864);
	dblcompex s(-0.5, -0.86602540378443864);

	// _xZkqqHx_jHeMiMpIiMiFjAjEiCpAiCfEiMiFiCmFjEoEiKhCiBeCiMoLiNlHiCkKjBoFiCkLiCkCiClNiCnPiBeBiCfEiMiFiClFiCkJiOoGiCoKiCmIiCkCiBeC_xHqqkZx_
	if (figeq(bc, xxx, 5))
	{
	} else
	if (figeq(bc * r, xxx, 5))
	{
		b *= r;
	} else
	{
		b *= s;
	}
	dblcompex aaa(dA / 3.0, 0.0);
	m_cSol[0] = b + c - aaa;
	m_cSol[1] = b * r + c * s - aaa;
	m_cSol[2] = b * s + c * r - aaa;
	return 3;
} // FCCardano::Solve.

