式を展開せよ

数学 プログラム C++

中学１年生あたりで習うと思うのだが、いわゆる分配法則とか頑張って多項式に落とす作業です。
(x-y)^2 = x^2 + 2xy + y^2
(x-y)(x+y) = x^2 - y^2
等々。
やり方は完全に決まっているので覚えてしまえば誰でも出来るのですが式が大きくなってくると割と大変です。なのでそういう作業はコンピュータにやってもらおうと思って走り書きしてみました。

template< class COEFFICIENT >
class Term
{
	enum{ m_size = 256 };
public:
	Term( char variable, COEFFICIENT coefficient ):m_coefficient(coefficient){
		for( int i = 0; i < m_size; ++i ){
			m_aDegree[i] = 0;
		}
		m_aDegree[variable] = true;
	}
	Term( const Term< COEFFICIENT >& term ):m_coefficient(term.m_coefficient){
		for( int i = 0; i < m_size; ++i ){
			m_aDegree[i] = term.m_aDegree[i];
		}
	}
	virtual ~Term(){}
	
	bool operator==( const Term< COEFFICIENT >& term )const{
		for( int i = 0; i < m_size; ++i ){
			if( term.m_aDegree[i] != m_aDegree[i] ){
				return false;
			}
		}
		return true;
	}
	Term< COEFFICIENT > operator*( const Term< COEFFICIENT >& term )const{
		Term temp;
		temp.m_coefficient = m_coefficient * term.m_coefficient;
		for( int i = 0; i < m_size; ++i ){
			temp.m_aDegree[i] = m_aDegree[i] + term.m_aDegree[i];
		}
		return temp;
	}
	Term< COEFFICIENT > operator-()const{
		Term< COEFFICIENT > temp( *this );
		temp.m_coefficient = -m_coefficient;
		return temp;
	}
private:
	COEFFICIENT m_coefficient;
	int m_aDegree[m_size];


	Term(){}
	friend std::ostream& operator<<( std::ostream& out, const Term< COEFFICIENT >& term ){
		if( term.m_coefficient == -1 ){
			out << '-';
		}else if( term.m_coefficient != 1 ){
			out << term.m_coefficient;
		}
		for( int i = 0; i < m_size; ++i ){
			if( term.m_aDegree[i] != 0 ){
				out << static_cast< char >( i );
				if( term.m_aDegree[i] < 0 ){
					out << "^(" << term.m_aDegree[i] << ')';
				}else if( term.m_aDegree[i] != 1 ){
					out << "^" << term.m_aDegree[i];
				}
			}
		}
		return out;
	}
	template< class COEFFICIENT > friend class Polynomial;
};


template< class COEFFICIENT >
class Polynomial
{
public:
	Polynomial( char variable, COEFFICIENT coefficient = 1 ){
		m_aTerms.push_back( Term< COEFFICIENT >( variable, coefficient ) );
	}
	Polynomial( const Polynomial< COEFFICIENT >& poly ):m_aTerms( poly.m_aTerms ){}
	virtual~ Polynomial(){}
	
	//足し算
	Polynomial< COEFFICIENT > operator+( const Polynomial< COEFFICIENT >& poly )const{
		Polynomial< COEFFICIENT > result( *this );
		for( auto pTerm = poly.m_aTerms.begin(); pTerm != poly.m_aTerms.end(); ++pTerm ){
			result += *pTerm;
		}
		return result;
	}
	//引き算
	Polynomial< COEFFICIENT > operator-( const Polynomial< COEFFICIENT >& poly )const{
		Polynomial< COEFFICIENT > result( *this );
		for( auto pTerm = poly.m_aTerms.begin(); pTerm != poly.m_aTerms.end(); ++pTerm ){
			result += -*pTerm;
		}
		return result;
	}
	//掛け算
	Polynomial< COEFFICIENT > operator*( const Polynomial< COEFFICIENT >& poly )const{
		Polynomial< COEFFICIENT > result;
		for( auto pTerm1 = poly.m_aTerms.begin(); pTerm1 != poly.m_aTerms.end(); ++pTerm1 ){
			for( auto pTerm2 = m_aTerms.begin(); pTerm2 != m_aTerms.end(); ++pTerm2 ){
				result += *pTerm1 * *pTerm2;
			}
		}
		return result;
	}
	Polynomial< COEFFICIENT > operator*( const COEFFICIENT& d )const{
		Polynomial< COEFFICIENT > result( *this );
		for( auto pTerm = result.m_aTerms.begin(); pTerm != result.m_aTerms.end(); ++pTerm ){
			pTerm->m_coefficient *= d;
		}
		return result;
	}
	
	//符号反転
	Polynomial< COEFFICIENT > operator-()const{
		Polynomial< COEFFICIENT > result( *this );
		for( auto pTerm = result.m_aTerms.begin(); pTerm != result.m_aTerms.end(); ++pTerm ){
			pTerm->m_coefficient = -pTerm->m_coefficient;
		}
		return result;
	}

	//項を得る。
	Polynomial< COEFFICIENT > GetTerm( char c, int n ){
		Polynomial< COEFFICIENT > result;
		for( auto pTerm = m_aTerms.begin(); pTerm != m_aTerms.end(); ++pTerm ){
			if( pTerm->m_aDegree[c] == n ){
				result += *pTerm;
			}
		}
		for( auto pTerm = result.m_aTerms.begin(); pTerm != result.m_aTerms.end(); ++pTerm ){
			pTerm->m_aDegree[c] = 0;
		}
		return result;
	}


	//出力
	friend std::ostream& operator<<( std::ostream& out, const Polynomial< COEFFICIENT >& poly ){
		if( poly.m_aTerms.empty() ){
			out << '0';
		}
		for( auto pTerm = poly.m_aTerms.begin(); pTerm != poly.m_aTerms.end(); ++pTerm ){
			out << *pTerm << " + ";
		}
		return out;
	}

private:
	Polynomial< COEFFICIENT >& operator+=( const Term< COEFFICIENT >& term )
	{
		for( auto pThisTerm = m_aTerms.begin(); pThisTerm != m_aTerms.end(); ++pThisTerm ){
			if( *pThisTerm == term ){
				const COEFFICIENT tempCoef = pThisTerm->m_coefficient + term.m_coefficient;
				if( tempCoef != 0 ){
					pThisTerm->m_coefficient += term.m_coefficient;
				}else{
					m_aTerms.erase( pThisTerm );
				}
				return *this;
			}
		}
		m_aTerms.push_back( term );
		return *this;
	}
	Polynomial< COEFFICIENT >& operator*=( const Term< COEFFICIENT >& term )
	{
		for( auto pThisTerm = m_aTerms.begin(); pThisTerm != m_aTerms.end(); ++pThisTerm ){
			if( *pThisTerm == term ){
				const COEFFICIENT tempCoef = pThisTerm->m_coefficient + term.m_coefficient;
				if( tempCoef != 0 ){
					pThisTerm->m_coefficient += term.m_coefficient;
				}else{
					m_aTerms.erase( pthisTern );
				}
				return *this;
			}
		}
		m_aTerms.push_back( term );
		return *this;
	}
	
	Polynomial(){}
private:
	std::vector< Term< COEFFICIENT > > m_aTerms;

};

２時間ぐらいしか時間掛けてないので作りは結構適当です。まあ結果がわかればいいや程度の作りです。使い方はこんな感じ

Polynomial< int > a( 'a' ), b( 'b' );
std::cout << (a+b)*(a+b) << std::endl;
std::cout << (a+b)*(a-b) << std::endl;

こうすれば展開された式が表示される。まあ最後に＋が出たりして表示が若干おかしいが作り込むのめどいので頑張って解釈してください（＾＾；

で、なんでこんなの作ったかというと、３×３の行列を対角化しようと思って３次方程式を展開しようと思ったわけです。普通は対角化といえばHouseholderとかで頑張ると思うのですが３×３程度なら特性方程式を解いた方が早いかもとか思ったわけです。で、試しに{{a,d,e},{d,b,f},{e,f,c}}の対角行列の特性方程式を立てて係数を取り出して３次方程式の解の公式のルートの中身だけ表示してみたんですよ。

Polynomial< int > a( 'a' ), b( 'b' ), c( 'c' ), d( 'd' ), e( 'e' ), f( 'f' ), n( 'n' );
//特性方程式
Polynomial< int > z = (a-n)*(b-n)*(c-n)+d*f*e+d*f*e-(a-n)*f*f-(b-n)*e*e-(c-n)*d*d;
//各項の取り出し
Polynomial< int > A = -z.GetTerm( 'n', 3 );
Polynomial< int > B = -z.GetTerm( 'n', 2 );
Polynomial< int > C = -z.GetTerm( 'n', 1 );
Polynomial< int > D = -z.GetTerm( 'n', 0 );
//３次方程式のルートの中身
Polynomial< int > p = B*B*B*2 - A*B*C*9 + A*A*D*27;
Polynomial< int > q = B*B - A*C*3;
Polynomial< int > pd = p*p - q*q*q*4;
std::cout << "pd = " << pd << std::endl;

そしたら結果は以下のとおり

pd = -27b^2c^4 + 54abc^4 + -27a^2c^4 + 54b^3c^3 + -54ab^2c^3 + -54a^2bc^3 + 54a^3c^3 + 54bc^3f^2 + -54ac^3f^2 + -54bc^3e^2 + 54ac^3e^2 + -108c^4d^2 + 216bc^3d^2 + 216ac^3d^2 + 216c^3def + -27b^4c^2 + -54ab^3c^2 + 162a^2b^2c^2 + -54a^3bc^2 + -216b^2c^2f^2 + 270abc^2f^2 + 270abc^2e^2 + -54b^2c^2d^2 + -540abc^2d^2 + -324bc^2def + -27a^4c^2 + -216a^2c^2e^2 + -54a^2c^2d^2 + -324ac^2def + 54ab^4c + -54a^2b^3c + -54a^3b^2c + 54b^3cf^2 + 270ab^2cf^2 + 216b^3ce^2 + -540ab^2ce^2 + -54b^3cd^2 + 270ab^2cd^2 + -324b^2cdef + 54a^4bc + -540a^2bcf^2 + 270a^2bce^2 + 270a^2bcd^2 + 1296abcdef + 216a^3cf^2 + 54a^3ce^2 + -54a^3cd^2 + -324a^2cdef + -27a^2b^4 + 54a^3b^3 + -54ab^3f^2 + -108b^4e^2 + 216ab^3e^2 + 54ab^3d^2 + 216b^3def + -27a^4b^2 + -54a^2b^2e^2 + -216a^2b^2d^2 + -324ab^2def + 216a^3bf^2 + -54a^3be^2 + 54a^3bd^2 + -324a^2bdef + -108a^4f^2 + 216a^3def + -27c^2f^4 + 270bcf^4 + -216acf^4 + -54c^2e^2f^2 + 54ace^2f^2 + -540c^2d^2f^2 + 54bcd^2f^2 + 1026acd^2f^2 + 972cdef^3 + -27b^2f^4 + -216abf^4 + -540b^2e^2f^2 + 1026abe^2f^2 + -54b^2d^2f^2 + 54abd^2f^2 + 972bdef^3 + -54a^2c^2f^2 + -54a^2b^2f^2 + 216a^2f^4 + -540a^2e^2f^2 + -540a^2d^2f^2 + -1944adef^3 + -27c^2e^4 + -216bce^4 + 270ace^4 + -540c^2d^2e^2 + 1026bcd^2e^2 + 54acd^2e^2 + 972cde^3f + -54b^2c^2e^2 + 54bce^2f^2 + 216b^2e^4 + -216abe^4 + -540b^2d^2e^2 + -1944bde^3f + -27a^2e^4 + -54a^2d^2e^2 + 972ade^3f + 216c^2d^4 + -216bcd^4 + -216acd^4 + -1944cd^3ef + -27b^2d^4 + 270abd^4 + 972bd^3ef + 54abd^2e^2 + -27a^2d^4 + 972ad^3ef + 2268d^2e^2f^2 + -108f^6 + -324e^2f^4 + -324d^2f^4 + -324e^4f^2 + -324d^4f^2 + -108e^6 + -324d^2e^4 + -324d^4e^2 + -108d^6 +

・・・orz

手計算しなくてよかった。

だれか因数分解してくれるプログラム書いてください（ぉ