二次方程式の解は公式「（－b±√(b*b－4ac)）／2a」を使えば簡単に計算できる。とは言うものの実際にソースコードとして実装しようと思うと結構大変。
このプログラムではヒトが計算するのと同じように平方根の展開や分数の約分を行って解を求めた。
#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で素因数分解を１回行った後は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; //√の中がゼロはおかしいので１にしておく
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("ｉ");
}
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;
}
