|
Re: determinant algorithms
The kind of optomization I'm doing is (for example) if I have
X = a*a*b*b;
Y = c*c*d*d;
Z = a*b*c*d;
then I calculate
and then use
X = P*P;
Y = Q*Q;
Z = P*Q
Notice that the number of * operators is down from 9 to 5. Do you know if the compiler will do that automatically?
Here's the code.
The polynomial class:
// polynomial.h
// Written by Blake Foster
// Last updated February 7, 2006
#ifndef POLYNOMIAL_H
#define POLYNOMIAL_H
#include <vector>
#include <map>
#include <cassert>
#include <iostream>
#include <string>
#include <algorithm>
// macro used for naming variables
#define alphabet "abcdefghijklmnopqrstuvwxyz"
using namespace std;
template <typename COEFF>
class PolynomialIterator;
// data structure for storing the exponents of a single term.
class expList
{
public:
expList(unsigned int,int[]); // construct from array
expList(const vector<int>&); // construct from vector
expList(const expList&); // copy constructor
expList(const int&); // expList of length determined by int.
expList(); // construct an empty expList
~expList(); // destructor
int& operator[](const int&); // return exponent from index
bool operator>(const expList&) const; // lexicographic compare
bool operator<(const expList&) const; // lexicographic compare
bool operator==(const expList&) const;
bool operator!=(const expList&) const;
void operator=(const expList&);
expList operator+(const expList&) const; // vector addition (used to add
// exponents when multiplying terms)
void operator+=(const expList&); // vector addition
expList operator-() const; // negate exponents
unsigned int size() const; // the number of terms
expList conj(const unsigned int&) const; // negate all exponents except the
// the exponent indexed by the parameter
expList conjInv(const unsigned int&) const; // negate exponent indexed by
// the parameter
expList invert() const; // negate all exponents
bool constant() const; // determine if all exponents are 0
void printV(ostream&) const; // write to ostream& as a vector
void print(ostream&) const; // write to ostream& as a monomial (no
// coefficient). Variables are named a,b,c, ... ,z if there are <= 26
// variables, or x0, x1, ... , xn for n > 26 variables.
void print(ostream&, const vector<string>*) const; // write to ostream& with
// variables names determined by vector<string>*.
private:
unsigned int numVars; // the number of variables
int* exponents; // array containing exponents
};
// A class representing a polynomial in arbitrarily many variables and their
// inverses over an arbitrary field. The template parameter can be any class
// representing a field element (+,-,*,/ defined). If the constructor
// COEFF(int) defined. This is used to construct additive and multiplicative
// identities, 0 and 1.
template <typename COEFF>
class Polynomial
{
public:
Polynomial<COEFF>(const Polynomial<COEFF>&); // copy constructor
Polynomial<COEFF>(const COEFF&, const expList&); // construct a monomial
// with coefficient determined by COEFF and exponents determined by expList
Polynomial<COEFF>(const int&); // construct the 0 polynomial. Number of
// variables determined by int.
Polynomial<COEFF>(); // construct the 0 polynomial in 0 variables
void operator=(const Polynomial<COEFF>&);
// Polynomial algebra:
Polynomial<COEFF> operator+(const Polynomial<COEFF>&) const;
Polynomial<COEFF> operator-(const Polynomial<COEFF>&) const;
Polynomial<COEFF> operator-() const;
Polynomial<COEFF> operator*(const Polynomial<COEFF>&) const;
Polynomial<COEFF> operator*(const COEFF&) const;
void insert(const COEFF&, const expList&); // insert a term
Polynomial<COEFF> conj(const unsigned int&) const; // Invert the exponents
// of all variables but the variable indexed by int.
Polynomial<COEFF> invert() const; // Invert all exponents
void printV() const; // print the exponents of each term as a vectors
void print(ostream&) const; // print the polynomial with variables denoted
// a,b,c, ... ,z if there are <= 26 variables, or as x1, x2, ..., xn for
// n > 26 variables.
void print(ostream&, const vector<string>*) const; // print polynomial with
// variable names determined by the vector<string>* (each element represents
// a variables name, in the order.
int countTerms() const; // count the number of terms
int getNumTerms() const; // return numTerms
void clear(); // delete all terms
bool zero(); // return true if THIS is the zero polynomial (no terms, or all
// coefficients are 0)
bool constant(); // return true if this is a constant polynomial (only 1
// term and all exponents 0)
PolynomialIterator<COEFF> begin(); // return an iterator pointing to the
// first term
PolynomialIterator<COEFF> end(); // return an iterator pointing past the
// last term
unsigned int getNumVars() const; // return the number of variables
int maxDeg(int); // return the highest power of the variables indexed by int
int minDeg(int); // return the lowest power of the variable indexed by int
Polynomial<COEFF> plugIn(const int&,const COEFF&); // plug in COEFF for the
// variable indexed by int
COEFF plugIn(const vector<COEFF>&); // evaluate plolynomial with values
// deterined by vector<COEFF>.
Polynomial<COEFF> multByVars(const expList&) const;
private:
COEFF evalTerm(PolynomialIterator<COEFF>&, vector<COEFF>& values); // plug
// in values for a single term (used by plugIn)
map<expList,COEFF,less<expList> > terms; // map for storing terms. The
// exponent list is used as the key (sorted using the < operator defined in
// the expList class) and the coefficient is the value.
unsigned int numVars; // the number of variables
// additive identity of ground field. Defined as a data member since the
// insert member function is called repeatedly by the + and * operator, and
// every time the coefficient must be tested for equality to zero.
static const COEFF ZERO;
int numTerms;
};
// Iterator to iterate over all terms of a Polynomial sorted in lexicographic
// order.
template <typename COEFF>
class PolynomialIterator
{
friend class Polynomial<COEFF>;
public:
void operator++();
pair<expList, COEFF> operator*(); // return the current term as a pair
// containing the coefficient and the exponent list.
bool operator==(const PolynomialIterator<COEFF>&) const;
bool operator!=(const PolynomialIterator<COEFF>&) const;
protected:
// constructors can only be accessed by the Polynomial class. Iterators
// should be constructed using the begin() and end() functions of the
// Polynomial class.
PolynomialIterator(map<expList,COEFF,less<expList> >::iterator&);
PolynomialIterator();
map<expList,COEFF,less<expList> >::iterator value;
};
// expList class
inline expList::expList(unsigned int n, int exps[])
: exponents(new int [n]), numVars(n)
{
for (unsigned int i=0; i<numVars; ++i) exponents[i] = exps[i];
}
inline expList::expList(const vector<int>& V)
: exponents(new int [V.size()]), numVars(V.size())
{
for (unsigned int i=0; i<numVars; ++i) exponents[i] = V[i];
}
inline expList::expList(const expList& E)
: exponents(new int [E.numVars]), numVars(E.numVars)
{
for (unsigned int i=0; i<numVars; ++i) exponents[i] = E.exponents[i];
}
inline expList::expList(const int& N)
: exponents(new int [N]), numVars(N)
{}
inline expList::expList()
: exponents(NULL), numVars(0)
{}
inline expList::~expList()
{
delete [] exponents;
}
inline int& expList::operator[](const int& i)
{
return exponents[i];
}
inline bool expList::operator<(const expList& L) const
{
assert (numVars==L.numVars);
for (unsigned int i=0; i<numVars; ++i)
{
if (exponents[i]<L.exponents[i]) return false;
if (exponents[i]>L.exponents[i]) return true;
}
return false;
}
inline bool expList::operator>(const expList& L) const
{
assert (numVars==L.numVars);
for (unsigned int i=0; i<numVars; ++i)
{
if (exponents[i]<L.exponents[i]) return true;
if (exponents[i]>L.exponents[i]) return false;
}
return false;
}
inline bool expList::operator==(const expList& L) const
{
assert(L.size()==size());
for (unsigned int i=0; i<numVars; ++i)
{
if (L.exponents[i]!=exponents[i]) return false;
}
return true;
}
inline bool expList::operator!=(const expList& L) const
{
assert(L.size()==size());
for (unsigned int i=0; i<numVars; ++i)
{
if (L.exponents[i]!=exponents[i]) return true;
}
return false;
}
inline void expList::operator=(const expList& E)
{
delete exponents;
exponents = new int [E.numVars];
numVars = E.numVars;
for (unsigned int i=0; i<numVars; ++i) exponents[i] = E.exponents[i];
}
inline expList expList::operator+(const expList& L) const
{
assert(numVars==L.numVars);
expList sum(*this);
for (unsigned int i=0; i<numVars; i++) sum.exponents[i]+=L.exponents[i];
return sum;
}
inline void expList::operator+=(const expList& L)
{
assert(numVars==L.numVars);
for (unsigned int i=0; i<numVars; i++) exponents[i] += L.exponents[i];
}
inline expList expList::operator-() const
{
expList inverse(numVars);
for (unsigned int i=0; i<numVars; ++i)
{
inverse.exponents[i] = -exponents[i];
}
return inverse;
}
inline unsigned int expList::size() const
{
return numVars;
}
inline expList expList::conj(const unsigned int& V) const
{
expList E(numVars);
for (unsigned int i=0; i<numVars; i++)
{
if (i!=V) E.exponents[i]=-exponents[i];
else E.exponents[i]=exponents[i];
}
return E;
}
inline expList expList::conjInv(const unsigned int& V) const
{
expList E(*this);
E.exponents[V] = -exponents[V];
return E;
}
inline expList expList::invert() const
{
expList E(numVars);
for (unsigned int i=0; i<numVars; i++)
{
E.exponents[i]=-exponents[i];
}
return E;
}
inline bool expList::constant() const
{
for (unsigned int i=0; i<numVars; i++) if (exponents[i]!=0) return false;
return true;
}
inline void expList::printV(ostream& ost) const
{
assert(numVars > 0);
ost<<'(';
ost<<exponents[0];
for (unsigned int i=1; i<numVars; ++i)
{
ost<<", "<<exponents[i];
}
ost<<')';
}
inline void expList::print(ostream& ost) const
{
for (unsigned int i=0; i<numVars; ++i)
{
if (exponents[i]!=0)
{
if (numVars<=26) ost<<alphabet[i];
else ost<<'x'<<i;
if (exponents[i]!=1) ost<<'^'<<exponents[i];
}
}
}
inline void expList::print(ostream& ost, const vector<string>* variable_names) const
{
for (unsigned int i=0; i<numVars; ++i)
{
if (exponents[i]!=0)
{
ost<<(*variable_names)[i];
if (exponents[i]!=1) ost<<'^'<<exponents[i];
}
}
}
// Polynomial class
template <typename COEFF>
inline Polynomial<COEFF>::Polynomial(const Polynomial<COEFF>& P)
: terms(P.terms), numVars(P.numVars), numTerms(P.numTerms)
{}
template <typename COEFF>
inline Polynomial<COEFF>::Polynomial(const COEFF& C, const expList& E)
: terms(map<expList,COEFF,less<expList> >(less<expList>())), numVars(E.size()),
numTerms(1)
{
terms.insert(make_pair(E,C));
}
template <typename COEFF>
inline Polynomial<COEFF>::Polynomial(const int& n)
: terms(map<expList,COEFF,less<expList> >(less<expList>())), numVars(n),
numTerms(0)
{}
template <typename COEFF>
inline Polynomial<COEFF>::Polynomial()
: terms(map<expList,COEFF,less<expList> >(less<expList>())), numVars(0),
numTerms(0)
{
++default_construct_count;
}
template <typename COEFF>
inline void Polynomial<COEFF>::operator=(const Polynomial& P)
{
++assignment_count;
terms=P.terms;
numVars=P.numVars;
numTerms=P.numTerms;
}
template <typename COEFF>
inline Polynomial<COEFF> Polynomial<COEFF>::operator+(const Polynomial& P) const
{
if (numTerms > P.numTerms)
{
Polynomial<COEFF> sum(*this);
for (map<expList,COEFF,less<expList> >::const_iterator iter=P.terms.begin();
iter!=P.terms.end();
++iter)
{
sum.insert(iter->second,iter->first);
}
return sum;
}
else
{
Polynomial<COEFF> sum(P);
for (map<expList,COEFF,less<expList> >::const_iterator iter=terms.begin();
iter!=terms.end();
++iter)
{
sum.insert(iter->second,iter->first);
}
return sum;
}
}
template <typename COEFF>
inline Polynomial<COEFF> Polynomial<COEFF>::operator-(const Polynomial& P) const
{
++subtract_count;
if (numTerms > P.numTerms)
{
Polynomial<COEFF> diff(*this);
for (map<expList,COEFF,less<expList> >::const_iterator iter=P.terms.begin();
iter!=P.terms.end();
++iter)
{
++insert_calls_by_sub;
diff.insert(-(iter->second),iter->first);
}
return diff;
}
else
{
Polynomial<COEFF> diff(-P);
for (map<expList,COEFF,less<expList> >::const_iterator iter=terms.begin();
iter!=terms.end();
++iter)
{
++insert_calls_by_sub;
diff.insert((iter->second),iter->first);
}
return diff;
}
}
template <typename COEFF>
inline Polynomial<COEFF> Polynomial<COEFF>::operator-() const
{
Polynomial<COEFF> neg(*this);
for (map<expList,COEFF,less<expList> >::iterator iter = neg.terms.begin();
iter!=neg.terms.end();
++iter)
{
iter->second = -(iter->second);
}
return neg;
}
template <typename COEFF>
inline Polynomial<COEFF> Polynomial<COEFF>::operator*(const Polynomial& P) const
{
assert(numVars==P.numVars);
Polynomial<COEFF> product(numVars);
for (map<expList,COEFF,less<expList> >::const_iterator iter1=terms.begin(); iter1!=terms.end(); ++iter1)
{
for (map<expList,COEFF,less<expList> >::const_reverse_iterator iter2=P.terms.rbegin(); iter2!=P.terms.rend(); ++iter2)
{
product.insert((iter1->second)*(iter2->second),(iter1->first+iter2->first));
}
}
return product;
}
template <typename COEFF>
inline void Polynomial<COEFF>::insert(const COEFF& C, const expList& E)
{
assert(E.size()==numVars);
if (C!=ZERO)
{
pair<map<expList,COEFF,less<expList> >::iterator, bool>
P(terms.insert(make_pair(E, C)));
if (!P.second)
{
P.first->second += C;
if (P.first->second == 0)
{
terms.erase(P.first);
--numTerms;
}
}
else
{
++numTerms;
}
}
}
template <typename COEFF>
inline Polynomial<COEFF> Polynomial<COEFF>::conj(const unsigned int& V) const
{
if (terms.empty()) return *this;
else
{
Polynomial<COEFF> P(numVars);
for (map<expList,COEFF,less<expList> >::const_iterator iter=terms.begin();
iter!=terms.end();
++iter)
{
P.terms.insert(make_pair(iter->first.conj(V), iter->second));
}
return P;
}
}
template <typename COEFF>
inline Polynomial<COEFF> Polynomial<COEFF>::invert() const
{
if (terms.empty()) return Polynomial<COEFF>(numVars);
else
{
Polynomial<COEFF> P(numVars);
for (map<expList,COEFF,less<expList> >::const_iterator iter = terms.begin();
iter!=terms.end();
++iter)
{
P.terms.insert(P.terms.begin(), make_pair(-iter->first, iter->second));
}
return P;
}
}
template <typename COEFF>
void Polynomial<COEFF>::printV() const
{
if (terms.empty()) cout<<"0";
else
{
map<expList,COEFF,less<expList> >::const_iterator iter=terms.begin();
cout<<(iter->second);
iter->first.printV();
++iter;
for (; iter!=terms.end(); ++iter)
{
cout<<", "<<(iter->second);
iter->first.printV();
}
}
}
template <typename COEFF>
void Polynomial<COEFF>::print(ostream& ost) const
{
if (terms.empty()) ost<<"0";
else
{
map<expList,COEFF,less<expList> >::const_iterator iter=terms.begin();
COEFF C=iter->second;
if(C!=C*C||iter->first.constant()) ost<<C;
iter->first.print(ost);
++iter;
for (; iter!=terms.end(); ++iter)
{
ost<<'+';
C=iter->second;
if(C!=C*C||iter->first.constant()) ost<<C;
iter->first.print(ost);
}
}
}
template <typename COEFF>
void Polynomial<COEFF>::print(ostream& ost, const vector<string>* variable_names) const
{
if (terms.empty()) ost<<"0";
else
{
map<expList,COEFF,less<expList> >::const_iterator iter=terms.begin();
COEFF C=iter->second;
if(C!=C*C||iter->first.constant()) ost<<C;
iter->first.print(ost, variable_names);
++iter;
for (; iter!=terms.end(); ++iter)
{
ost<<'+';
C=iter->second;
if(C!=C*C||iter->first.constant()) ost<<C;
iter->first.print(ost, variable_names);
}
}
}
template <typename COEFF>
inline int Polynomial<COEFF>::countTerms() const
{
int count = 0;
for (map<expList,COEFF,less<expList> >::const_iterator iter=terms.begin();
iter!=terms.end();
++iter)
{
++count;
assert(iter->second!=(iter->second)-(iter->second));
}
return count;
}
template <typename COEFF>
inline int Polynomial<COEFF>::getNumTerms() const
{
return numTerms;
}
template <typename COEFF>
inline void Polynomial<COEFF>::clear()
{
terms.clear();
numTerms = 0;
}
template <typename COEFF>
inline bool Polynomial<COEFF>::zero()
{
return terms.empty();
}
template <typename COEFF>
inline bool Polynomial<COEFF>::constant()
{
return terms.empty() || (terms.size()==1
&& (terms.begin()->first).constant());
}
template <typename COEFF>
const COEFF Polynomial<COEFF>::ZERO = COEFF(0);
// PolynomialIterator class
template <typename COEFF>
inline PolynomialIterator<COEFF> Polynomial<COEFF>::begin()
{
return PolynomialIterator<COEFF>(terms.begin());
}
template <typename COEFF>
inline PolynomialIterator<COEFF> Polynomial<COEFF>::end()
{
return PolynomialIterator<COEFF>(terms.end());
}
template <typename COEFF>
inline unsigned int Polynomial<COEFF>::getNumVars() const
{
return numVars;
}
template <typename COEFF>
inline int Polynomial<COEFF>::maxDeg(int var)
{
int deg = ((*begin()).first)[var];
for (PolynomialIterator<int> iter = begin(); iter!= end(); ++iter)
{
deg = max(deg, (*iter).first[var]);
}
return deg;
}
template <typename COEFF>
inline int Polynomial<COEFF>::minDeg(int var)
{
int deg = ((*begin()).first)[var];
for (PolynomialIterator<int> iter = begin(); iter!= end(); ++iter)
{
deg = min(deg, (*iter).first[var]);
}
return deg;
}
template <typename COEFF>
inline Polynomial<COEFF> Polynomial<COEFF>::plugIn(const int& var,
const COEFF& value)
{
assert(var<numVars);
Polynomial<COEFF> newPoly(numVars);
for (PolynomialIterator<COEFF> iter = begin(); iter!=end(); ++iter)
{
expList E = (*iter).first;
COEFF C = (*iter).second;
if (E[var] > 0)
{
for (int i=0; i<E[var]; ++i)
{
C = C*value;
}
}
else if (E[var] < 0)
if (E[var] > 0)
{
for (int i=0; i<E[var]; ++i)
{
C = C*value;
}
}
else
{
for (int i=0; i>E[var]; --i)
{
C = C/value;
}
}
E[var] = 0;
newPoly.insert(C,E);
}
return newPoly;
}
template <typename COEFF>
inline COEFF Polynomial<COEFF>::plugIn(const vector<COEFF>& values)
{
assert(values.size() == numVars);
PolynomialIterator<COEFF> iter = begin();
COEFF C = evalTerm(iter, values);
++iter;
for (; iter!=end(); ++iter)
{
C = C + evalTerm(iter, values);
}
return C;
}
template <typename COEFF>
inline COEFF Polynomial<COEFF>::evalTerm(PolynomialIterator<COEFF>& iter, vector<COEFF>& values)
{
assert(values.size() == numVars);
COEFF C = (*iter).second;
for (unsigned int i=0; i<numVars; ++i)
{
if ((*iter).first[i] > 0)
{
for (int j=0; j<(*iter).first[i]; ++j)
{
C = C*values[i];
}
}
else if ((*iter).first[i] < 0)
{
for (int j=0; j>(*iter).first[i]; --j)
{
C = C/values[i];
}
}
}
return C;
}
template <typename COEFF>
inline Polynomial<COEFF> Polynomial<COEFF>::multByVars(const expList& E) const
{
Polynomial<COEFF> product(numVars);
for (map<expList,COEFF,less<expList> >::const_iterator iter
= terms.begin();
iter!=terms.end();
++iter)
{
product.terms.insert(product.terms.end(),
make_pair(iter->first + E, iter->second));
}
return product;
}
// PolynomialIterator class
template <typename COEFF>
inline PolynomialIterator<COEFF>::PolynomialIterator(map<expList,COEFF,less<expList> >::iterator& iter)
: value(iter)
{}
template <typename COEFF>
inline PolynomialIterator<COEFF>::PolynomialIterator()
: value(map<expList,COEFF,less<expList> >::iterator())
{}
template <typename COEFF>
inline void PolynomialIterator<COEFF>::operator++()
{
++value;
}
template <typename COEFF>
inline pair<expList, COEFF> PolynomialIterator<COEFF>::operator*()
{
return *value;
}
template <typename COEFF>
inline bool PolynomialIterator<COEFF>::operator==(const PolynomialIterator<COEFF>& iter) const
{
return value == iter.value;
}
template <typename COEFF>
inline bool PolynomialIterator<COEFF>::operator!=(const PolynomialIterator<COEFF>& iter) const
{
return value != iter.value;
}
#endif
and the determinant function:
#define CONJX conj(0)
#define CONJY conj(1)
#define CONJZ conj(2)
// function to construct a monomial in 3 variables with coefficient coeff and
// exponents a, b, and c.
template <typename COEFF_TYPE>
inline Polynomial<COEFF_TYPE> monomial(const COEFF_TYPE& coeff,
const int& a, const int& b, const int& c)
{
int exps [] = {a, b, c};
expList E(3, exps);
return Polynomial<COEFF_TYPE>(coeff,E);
}
// function to construct an exponent list containing exponents a, b, and c.
inline expList makeTerm(const int& a, const int& b, const int& c)
{
int exps [] = {a, b, c};
expList E(3, exps);
delete [] exps;
return E;
}
template <typename COEFF_TYPE>
inline Polynomial<COEFF_TYPE> determinant(Polynomial<COEFF_TYPE>& A,
Polynomial<COEFF_TYPE>& B,
Polynomial<COEFF_TYPE>& C,
Polynomial<COEFF_TYPE>& D)
{
const COEFF_TYPE ONE(1);
Polynomial<COEFF_TYPE> T0(A*(A.CONJX));
Polynomial<COEFF_TYPE> T1(B*(B.CONJX));
Polynomial<COEFF_TYPE> T2(C*(C.CONJX));
Polynomial<COEFF_TYPE> T3(D*(D.CONJX));
Polynomial<COEFF_TYPE> T4(A*(B.CONJY));
Polynomial<COEFF_TYPE> T5(C*(D.CONJY));
Polynomial<COEFF_TYPE> T6(A*(D.CONJZ));
Polynomial<COEFF_TYPE> T7(C*(B.CONJZ));
Polynomial<COEFF_TYPE> T0y(T0.CONJY);
Polynomial<COEFF_TYPE> T3y(T3.CONJY);
Polynomial<COEFF_TYPE> T4z(T4.CONJZ);
Polynomial<COEFF_TYPE> T6y(T6.CONJY);
Polynomial<COEFF_TYPE> a(monomial(ONE, 1, 0, 0));
Polynomial<COEFF_TYPE> b(monomial(ONE, 0, 1, 0));
Polynomial<COEFF_TYPE> c(monomial(ONE, 0, 0, 1));
Polynomial<COEFF_TYPE> aInv(a.invert());
Polynomial<COEFF_TYPE> bInv(b.invert());
Polynomial<COEFF_TYPE> cInv(c.invert());
Polynomial<COEFF_TYPE> P0(T0*T0y);
Polynomial<COEFF_TYPE> P1(T1*(T1.CONJY));
Polynomial<COEFF_TYPE> P2(T2*(T2.CONJY));
Polynomial<COEFF_TYPE> P3(T3*T3y);
Polynomial<COEFF_TYPE> P4(-T4.multByVars(makeTerm(0,0,1))*T4z);
Polynomial<COEFF_TYPE> P5(-T0y.multByVars(makeTerm(-1,0,0))*T2);
Polynomial<COEFF_TYPE> P6(-T6.multByVars(makeTerm(0,1,0))*(T6.CONJY));
Polynomial<COEFF_TYPE> P7(-T7.multByVars(makeTerm(0,-1,0))*(T7.CONJY));
Polynomial<COEFF_TYPE> P8(-T1.multByVars(makeTerm(-1,0,0))*T3y);
Polynomial<COEFF_TYPE> P9(-T5.multByVars(makeTerm(0,0,-1))*(T5.CONJZ));
Polynomial<COEFF_TYPE> P10(T6y.multByVars(makeTerm(0,0,-1))*T7);
Polynomial<COEFF_TYPE> P11(T4z.multByVars(makeTerm(-1, -1, 0))*T5);
Polynomial<COEFF_TYPE> Q0(P10 + P11);
Polynomial<COEFF_TYPE> Q1(P4 + P6 + P7 + P9 + Q0);
Polynomial<COEFF_TYPE> Q2(P5 + P8);
return P0 + P1 + P2 + P3 + Q1 + Q2 + Q1.CONJX + (Q0 + Q2).CONJY + Q0.CONJZ;
}
|
27-Feb-2006, 10:00
