二次方程式の解を求める

test113.gif
二次方程式の解は公式「(-b±√(b*b-4ac))/2a」を使えば簡単に計算できる。とは言うものの実際にソースコードとして実装しようと思うと結構大変。
このプログラムではヒトが計算するのと同じように平方根の展開や分数の約分を行って解を求めた。

依存環境:ATL
#include "math.h"
#include "atlcoll.h"
#include "atlstr.h"



class	CElementInfo
{
public:
	long	_nElement;
	ULONG	_nCount;

	CElementInfo()
	{
		_nElement = 0;
		_nCount = 0;
	}

	CElementInfo(long nElement,ULONG nCount)
	{
		_nElement = nElement;
		_nCount = nCount;
	}

	CElementInfo&	operator++()
	{
		_nCount++;
		return	*this;
	}
	CElementInfo&	operator++(int)
	{
		_nCount++;
		return	*this;
	}

	CElementInfo&	operator--()
	{
		if(_nCount)
			_nCount--;
		return	*this;
	}
	CElementInfo&	operator--(int)
	{
		if(_nCount)
			_nCount--;
		return	*this;
	}
};



//
//	素因数分解
//
bool	Digest(long nData,CAtlArray<CElementInfo>* pacElements)
{
	long	nElement;
	ULONG	nCount;
	long	nLoopMax;

	if(pacElements == NULL)
		return	false;
	pacElements->RemoveAll();

	//負の数なら素因数として-1を与える
	if(nData < 0)
	{
		pacElements->Add(CElementInfo(-1,1));
		nData *= -1;
	}
	if(nData == 0 || nData == 1)
		return	true;

	nLoopMax = (long)sqrt((double)nData) + 1;
	for(nElement = 2; nElement <= nLoopMax; nElement += 2)
	{
		if(nElement == 4)	//2で素因数分解を1回行った後は3,5,7...2n+1で処理を
			nElement--;	//行なうように4のときに1を減算

		nCount = 0;
		while(1)
		{
			if(nData % nElement)
				break;

			nCount++;
			nData /= nElement;
		}
		if(nCount > 0)
			pacElements->Add(CElementInfo(nElement,nCount));
		if(nData == 1)
			break;
	}
	if(nData > 1)
		pacElements->Add(CElementInfo(nData,1));

	return	true;
}





//
//	平方根
//
//pacInnerElementsはルートの中側
//pacOuterElementsはルートの外側
//
bool	SquareRoot(long nData,CAtlArray<CElementInfo>* pacInnerElements,CAtlArray<CElementInfo>* pacOuterElements)
{
	bool	ret;
	CAtlArray<CElementInfo>	acDigest;

	if(pacInnerElements)
		pacInnerElements->RemoveAll();
	if(pacOuterElements)
		pacOuterElements->RemoveAll();
	if(pacInnerElements == NULL || pacOuterElements == NULL)
		return	false;

	if(nData == 0)
		return	true;

	pacInnerElements->Add(CElementInfo((nData < 0) ? -1 : 1,1));
	pacOuterElements->Add(CElementInfo(1,1));
	if(nData == 1 || nData == -1)
		return	true;

	if(nData < 0)
		nData *= -1;

	//素因数分解
	ret = Digest(nData,&acDigest);
	if(ret == false)
		return	false;

	size_t	i;
	size_t	nSize;

	nSize = acDigest.GetCount();
	for(i = 0; i < nSize; i++)
	{
		if(acDigest[i]._nCount == 1)
			pacInnerElements->Add(CElementInfo(acDigest[i]._nElement,1));
		else if(acDigest[i]._nCount % 2 == 0)
			pacOuterElements->Add(CElementInfo(acDigest[i]._nElement,acDigest[i]._nCount / 2));
		else
		{
			if(acDigest[i]._nCount != 0)
			{
				pacInnerElements->Add(CElementInfo(acDigest[i]._nElement,1));
				pacOuterElements->Add(CElementInfo(acDigest[i]._nElement,acDigest[i]._nCount / 2));
			}
		}
	}

	return	true;
}





long	GetValue(const CAtlArray<CElementInfo>& acElements)
{
	long	ret;
	size_t	i;
	size_t	j;
	size_t	nSizeI;
	size_t	nSizeJ;

	nSizeI = acElements.GetCount();
	if(nSizeI == 0)
		return	0;

	ret = 1;
	for(i = 0; i < nSizeI; i++)
	{
		nSizeJ = acElements[i]._nCount;
		for(j = 0; j < nSizeJ; j++)
			ret *= acElements[i]._nElement;
	}

	return	ret;
}



//
//	平方根
//
//pacInnerElementsはルートの中側
//pacOuterElementsはルートの外側
//
bool	SquareRoot(long nData,long* pnInner,long* pnOuter)
{
	bool	ret;
	CAtlArray<CElementInfo>	acInnerElements;
	CAtlArray<CElementInfo>	acOuterElements;

	if(pnInner)
		*pnInner = 0;
	if(pnOuter)
		*pnOuter = 0;
	if(pnInner == NULL || pnOuter == NULL)
		return	false;
	if(nData == 0)
		return	true;

	////////////////////////////////////
	//平方根の計算
	//
	ret = SquareRoot(nData,&acInnerElements,&acOuterElements);
	if(ret == false)
		return	false;


	////////////////////////////////////
	//結果値の計算
	//

	*pnInner = GetValue(acInnerElements);
	*pnOuter = GetValue(acOuterElements);

	if(*pnInner == 0)
		*pnInner = 1;

	return	true;
}




//
//	約分
//
//nDenominator	:分母
//nNumerator	:分子
//
bool	Reduction(long nDenominator,long nNumerator,CAtlArray<CElementInfo>* pacDenominatorElements,CAtlArray<CElementInfo>* pacNumeratorElements)
{
	CAtlArray<CElementInfo>*	pA;
	CAtlArray<CElementInfo>*	pB;

	if(pacDenominatorElements)
		pacDenominatorElements->RemoveAll();
	if(pacNumeratorElements)
		pacNumeratorElements->RemoveAll();
	if(pacNumeratorElements == NULL || pacDenominatorElements == NULL)
		return	false;

	if(nDenominator == 0 || nNumerator == 0)
		return	true;
	if(nDenominator == -1)
	{
		//分母が-1だったら、分子に-1をかけて終わり
		nDenominator = 1;
		nNumerator *= -1;
		return	true;
	}
	if(nDenominator == 1 || nNumerator == 1)
		return	true;


	///////////////////////////////
	//約分前処理
	//

	if(nDenominator < 0)
	{
		//分母がマイナスだったら分子にマイナス記号を移す
		nDenominator *= -1;
		nNumerator *= -1;
	}

	//分母・分子ともに素因数分解
	Digest(nDenominator,pacDenominatorElements);
	Digest(nNumerator,pacNumeratorElements);

	//素因数の数が少ない方をpA、多い方をpBにする
	if(pacDenominatorElements->GetCount() < pacNumeratorElements->GetCount())
	{
		pA = pacDenominatorElements;
		pB = pacNumeratorElements;
	}
	else
	{
		pA = pacNumeratorElements;
		pB = pacDenominatorElements;
	}


	///////////////////////////////
	//約分処理
	//
	size_t	i;
	size_t	j;
	size_t	nSizeA;
	size_t	nSizeB;

	i = 0;
	nSizeA = pA->GetCount();
	nSizeB = pB->GetCount();
	for(j = 0; j < nSizeB; j++)
	{
		for(; i < nSizeA; i++)
		{
			if((*pA)[i]._nElement < (*pB)[j]._nElement)
				continue;
			if((*pA)[i]._nElement > (*pB)[j]._nElement)
				break;

			//共通する素因数の数だけnCountを減らす
			if((*pA)[i]._nCount < (*pB)[j]._nCount)
			{
				(*pB)[j]._nCount -= (*pA)[i]._nCount;
				(*pA)[i]._nCount = 0;
			}
			else if((*pA)[i]._nCount > (*pB)[j]._nCount)
			{
				(*pA)[i]._nCount -= (*pB)[j]._nCount;
				(*pB)[j]._nCount = 0;
			}
			else
			{
				(*pA)[i]._nCount = 0;
				(*pB)[j]._nCount = 0;
			}
			break;
		}
	}

	return	true;
}


//
//	約分
//
//nDenominator	:分母
//nNumerator	:分子
//
bool	Reduction(long& nDenominator,long& nNumerator)
{
	bool	ret;
	CAtlArray<CElementInfo>		acElementsA;
	CAtlArray<CElementInfo>		acElementsB;

	///////////////////////////////
	//約分
	//
	ret = Reduction(nDenominator,nNumerator,&acElementsA,&acElementsB);
	if(ret == false)
		return	false;


	///////////////////////////////
	//約分結果の数値計算
	//
	nDenominator = GetValue(acElementsA);
	nNumerator = GetValue(acElementsB);

	return	true;
}





//
//	二次方程式の解
//
//
//入力:	ax^2 + bx + c = 0
//出力:	(nA1 + nB1√nC1)/nD1	、	(nA2 + nB2√nC2)/nD2
//
bool	QuadricEquation(long a,long b,long c,long& nA1,long& nB1,long& nC1,long& nD1,long& nA2,long& nB2,long& nC2,long& nD2)
{
	bool	ret;
	long	b24ac;
	CAtlArray<CElementInfo>		acBElements;
	CAtlArray<CElementInfo>		ac2AElements;
	CAtlArray<CElementInfo>		acInnerElements;
	CAtlArray<CElementInfo>		acOuterElements;

	if(a == 0)
		return	false;

	b24ac = b * b - 4 * a * c;		//√内を計算
	if(b24ac == 0)
	{
		//重解
		//√内がゼロなら重解

		a *= 2;
		b *= -1;
		ret = Reduction(a,b);
		if(ret == false)
			return	false;

		//ATLTRACE("重解:%d/%d",b,a);

		nA1 = b;
		nB1 = 0;
		nC1 = 1;
		nD1 = a;

		//重解なので値を同じにする
		nA2 = nA1;
		nB2 = nB1;
		nC2 = nC1;
		nD2 = nD1;

		return	true;
	}

	//√を計算して素因数分解
	ret = SquareRoot(b24ac,&acInnerElements,&acOuterElements);
	if(ret == false)
		return	false;

	//(-b)を素因数分解
	ret = Digest(-1 * b,&acBElements);
	if(ret == false)
		return	false;

	//(2a)を素因数分解
	ret = Digest(2 * a,&ac2AElements);
	if(ret == false)
		return	false;

	bool	bFound;
	size_t	i;
	size_t	j;
	size_t	k;
	size_t	nSizeI;
	size_t	nSizeJ;
	size_t	nSizeK;

	//素因数分解した値を元に解の公式に当てはめて約分
	nSizeI = ac2AElements.GetCount();
	for(i = 0; i < nSizeI; i++)
	{
		//分母分子に含まれる共通因子を探す

		bFound = false;
		nSizeJ = acBElements.GetCount();
		for(j = 0; j < nSizeJ; j++)
		{
			if(acBElements[j]._nElement != ac2AElements[i]._nElement)
				continue;
			bFound = true;
			break;
		}
		if(bFound == false || acBElements[j]._nCount == 0)
			continue;

		bFound = false;
		nSizeK = acOuterElements.GetCount();
		for(k = 0; k < nSizeK; k++)
		{
			if(acOuterElements[k]._nElement != ac2AElements[i]._nElement)
				continue;
			bFound = true;
			break;
		}
		if(bFound == false || acOuterElements[k]._nCount == 0)
			continue;


		//共通因子が見つかった!

		ULONG	nMin;

		//約分できる共通因子の数を求める
		nMin = min(ac2AElements[i]._nCount,min(acBElements[j]._nCount,acOuterElements[k]._nCount));

		//見つかった共通因子で約分
		ac2AElements[i]._nCount		-= nMin;
		acBElements[j]._nCount		-= nMin;
		acOuterElements[k]._nCount	-= nMin;
	}


	////////////////////////////
	//求まった解の整形
	//
	//(nA ± nB√nC)/nD という形の解を√内が1のときはnA±nBを計算して整形する
	//

	long	nA;
	long	nB;
	long	nC;
	long	nD;

	nA = GetValue(acBElements);
	nB = GetValue(acOuterElements);
	nC = GetValue(acInnerElements);
	nD = GetValue(ac2AElements);
	if(nC == 0 && nB)
		nC = 1;			//√の中がゼロはおかしいので1にしておく

	if(nC != 1)
	{
		nA1 = nA;
		nB1 = nB;
		nC1 = nC;
		nD1 = nD;

		nA2 = nA;
		nB2 = -1 * nB;
		nC2 = nC;
		nD2 = nD;

		return	true;
	}

	nA1 = nA + nB;
	nB1 = 0;
	nC1 = 1;
	nD1 = nD;
	ret = Reduction(nD1,nA1);	//約分
	if(ret == false)
		return	false;

	nA2 = nA - nB;
	nB2 = 0;
	nC2 = 1;
	nD2 = nD;
	ret = Reduction(nD2,nA2);	//約分
	if(ret == false)
		return	false;

	return	true;
}





bool	Test()
{
	CAtlString	strBuff;
	CAtlString	strMessage;

	bool	ret;
	long	a;
	long	b;
	long	c;
	long	nA1;
	long	nB1;
	long	nC1;
	long	nD1;
	long	nA2;
	long	nB2;
	long	nC2;
	long	nD2;

	//////////////////////////
	//二次方程式の解を計算
	//

	a = 2;
	b = 6;
	c = 2;

	ret = QuadricEquation(a,b,c,nA1,nB1,nC1,nD1,nA2,nB2,nC2,nD2);
	if(ret == false)
		return	false;


	//////////////////////////
	//結果を整形して表示
	//

	strMessage.Format(_T("二次方程式:(%d)x*x + (%d)x + (%d)=0の解は\n\n"),a,b,c);

	if(nA1 == nA2 && nB1 == nB2 && nC1 == nC2 && nD1 == nD2)
	{
		//重解

		if(nD1 == 1)
			strBuff.Format(_T("%d"),nA1);
		else
			strBuff.Format(_T("%d/%d"),nA1,nD1);

		//重解なのでnBとnCは無視してOK

		strMessage += strBuff;
		::MessageBox(NULL,strMessage,_T(""),MB_OK);
		return	true;
	}

	if(nA1 == nA2 && nB1 == -1 * nB2 && nC1 == nC2 && nD1 == nD2)
	{
		//通常の±形式の解

		if(nD1 != 1)
		{
			if(nB1)
				strMessage += _T("(");
		}

		if(nA1)
		{
			strBuff.Format(_T("%d"),nA1);
			strMessage += strBuff;
		}

		if(nB1)
		{
			strMessage += _T("±");
			if(nB1 < 0)
				nB1 *= -1;
			if(nB1 != 1)
			{
				strBuff.Format(_T("%d"),nB1);
				strMessage += strBuff;
			}

			if(nC1 != 1 && nC1 != -1)
			{
				strBuff.Format(_T("√%d"),(nC1 > 0) ? nC1 : -1 * nC1);
				strMessage += strBuff;
			}
			if(nC1 < 0)
				strMessage += _T("i");
		}

		if(nD1 != 1)
		{
			if(nB1)
				strBuff.Format(_T("%s)/%d"),strMessage,nD1);
			else
				strBuff.Format(_T("%s/%d"),strMessage,nD1);
			strMessage = strBuff;
		}

		::MessageBox(NULL,strMessage,_T(""),MB_OK);
		return	true;
	}

	//ここに来る場合は以下のATLASSERT条件を満たすはず
	ATLASSERT(nB1 == 0);
	ATLASSERT(nC1 == 1);
	ATLASSERT(nB2 == 0);
	ATLASSERT(nC2 == 1);


	if(nD1 == 1)
		strBuff.Format(_T("%d"),nA1);
	else
		strBuff.Format(_T("%d/%d"),nA1,nD1);
	strMessage += strBuff;

	strMessage += _T("\n");

	if(nD2 == 1)
		strBuff.Format(_T("%d"),nA2);
	else
		strBuff.Format(_T("%d/%d"),nA2,nD2);
	strMessage += strBuff;

	::MessageBox(NULL,strMessage,_T(""),MB_OK);
	return	true;
}

プロジェクトファイルをダウンロード


カテゴリー「数学」 のエントリー