determinant algorithms

Blake · #1 27-Feb-2006, 09:00

I'm running an algorithm that involves taking the determinants of 4 x 4 matrices with polynomial entries. The polynomials are in three variables and their inverses, over the finite field F2. I need to make it as fast as possible. The problem is that polynomials form a ring, not a field; In other words there is no division. Thus I can't use the standard algorithm of row reducing first and then taking the product of the diagonal. So far, the method of minors (breaking the matrix down into submatrices and taking their determinants) is the fastest thing I have tried. I'm wondering if anyone knows of a faster algorithm.

I've considered implementing the field of rational functions (ratios of polynomials), but then I would have to constantly get them in lowest terms to avoid winding up with huge denominators. I'm also aware that row reducing without division is possible, but much slower.

Does anyone have any suggestions? I have been poking around on the internet for awhile, and I haven't found anything useful.

QED · #2 01-Mar-2006, 08:28

Interesting problem. Are you implementing this in C or C++? (That is, procedural or object-oriented?)

If the matrix size is fixed (4x4) and you do not need this solution to be extensible, why not hard-code the determinant formula into your algorithm?

I.e., in pseudocode

Code:

det(A) = a11 * a22 * a33 * a44
         - a11 * a22 * a34 * a43
         + a11 * a23 * a34 * a42
         - a11 * a23 * a32 * a44
         + a11 * a24 * a32 * a43
         - a11 * a24 * a33 * a42
         - a12 * a23 * a34 * a41
         + a12 * a23 * a31 * a44
         - a12 * a24 * a31 * a43
         + a12 * a24 * a33 * a41
         - a12 * a21 * a33 * a44
         + a12 * a21 * a34 * a43
         + a13 * a24 * a31 * a42
         - a13 * a24 * a32 * a41
         + a13 * a21 * a32 * a44
         - a13 * a21 * a34 * a42
         + a13 * a22 * a34 * a41
         - a13 * a22 * a31 * a44
         - a14 * a21 * a32 * a43
         + a14 * a21 * a33 * a42
         - a14 * a22 * a33 * a41
         + a14 * a22 * a31 * a43
         - a14 * a23 * a31 * a42
         + a14 * a23 * a32 * a41

If this seems an ugly and brute-force approach, I imageine the same effeiciency could probably be achieved by carefully constructed functions expanded inline.

Matthew

Blake · #3 01-Mar-2006, 08:45

That brings up another interesting issue that I forgot to mention. I multiplied out the determinant, and then exploited the relations between the polynomials (they're not just any polynomials) to cut the number of multiplications down to 23 by factoring and conjugating things. Then I wrote an algorithm to test the new method that counts the number of calls to all the relevant functions. The result was that everything (copy constructor calls, addition/subtraction calls, multiplication calls, conjugation calls, etc.) was significantly less. However, the number of iterations of the inner loop on the multiplication operator was much larger for all but the smallest polynomials. As a result, that method, even with reduced number of multiplications, was slower. That is another mystery that I have yet to solve.

I'm working in C++. I wrote a polynomial class with overloaded +, -, and * operators. It uses an stl map with the list of exponents as a key (sorted lexicographically) and the coefficient as a value. Thus inserting runs in logarithmic time.

Thanks for the help.

QED · #4 01-Mar-2006, 09:22

You mean the inner loop of the polynomial multiplication operator, right? Because the method I suggested above requires no loop (other than that).

So in any case, by hard-coding the expansion you should still get some (small) speedup over actually implementing the method of minors, independent of whether you try to optimize the poly multiplications or not.

The slowdown you mention in the poly multiplication is interesting... but rather hard to grasp without some more detail. Can you post code? or is your work restricted/proprietary?

Matthew

Blake · #5 01-Mar-2006, 10:23

Yes, I mean the inner loop of the * operator.

I have done what you suggested and hard coded the determinant. That was actually slower than the method of minors. I'm not entirely sure why.

The new method I mentioned involved first multiplying out the determinant, and then factoring it to reduce the number of multiplications. Here's how it works:

All 16 entries in the matrix are generated from 4 polynomials A, B, C, and D in variables a, b, c and their inverses. The formula for the determinant is:

AAxAyA + BBxByBz + CCxCyCz + DDxDyDz - AxAyBBzc^-1 - AAzBxByc -AAxCyCza - AyAzCCxa^-1 - AAyDxDzb - AxAzDDyb^-1 - BByCxCzb
- BxBzCCyb^-1 - BBxDyDza^-1 -ByBzDDxa - CxCyDDzc - CCzDxDyc^-1 + AxBCzDyc^-1 + AyBzCDxc^-1 + ABxCyDzc + AzByCxDc + AByCzDxab + AzBxCDya^-1b^-1 + AxBzCyDab^-1 + AyBCxDza^-1b

For a polynomial P in a,b,c,a^-1,b^-1,c^-1, Px means invert the exponents on b and c but leave a, Py means invert the exponents on a and c but leave b, and Pz means invert the exponents on a and b but leave c.

It is easy to varify that (Px)y = (Py)x = Pz,
(Pz)x = (Px)z = Py,
P(y)z = (Pz)y = Px,
and (Px)x = (Py)y = (Pz)z = P

Using this I can write the determinant as the sum of conjugates of 12 polynomials P0 - P11, defined as follows:

AAxAyAz = P0
BBxByBz = P1
CCxCyCz = P2
DDxDyDz = P3

-AxAyBBzc-1 = P4x
-AAzBxByc = P4

-AAxCyCza = P5y
-AyAzCCxa-1 = P5

-AAyDxDzb = P6
-AxAzDDyb-1 = P6x

-BByCxCzb = P7x
-BxBzCCyb-1 = P7

-BBxDyDza-1 = P8
-ByBzDDxa = P8y

-CxCyDDzc = P9x
-CCzDxDyc-1 = P9

AxBCzDyc-1 = P10z
AyBzCDxc-1 = P10
ABxCyDzc = P10y
AzByCxDc = P10x

AByCzDxab = P11z
AzBxCDya-1b-1 = P11
AxBzCyDab-1 = P11y
AyBCxDza-1b = P11x

Thus the determinant is
P0 + P1 + P2 + P3 + P4 + P4x + P5 + P5y + P6 + P6x + P7 + P7x + P8 + P8y + P9 + P9x + P10 + P10x + P10y + P10z + P11 + P11x + P11y + P11z

Now, there is also redundancy in the formulas for the Pi's that I can exploit. Define the following 8 polynomials:

T0 = AAx
T1 = BBx
T2 = CCx
T3 = DDx
T4 = AzBx
T5 = CDy
T6 = AyDx
T7 = BzC

Then I can write P0 through P11 as

P0 = T0T0y
P1 = T1T1y
P2 = T2T2y
P3 = T3T3y
P4 = -T4T4zc
P5 = -T0yT2a-1
P6 = -T6T6yb
P7 = -T7T7yb-1
P8 = -T1T3ya-1
P9 = -T5T5zc-1
P10 = T6T7c-1
P11 = T4T5a-1b-1

Note that each Ti is used exactly twice.

Finally, there is still more redundancy in the addition. Let

Q0 = P10 + P11
Q1 = P4 + P6 + P7 + P9 + Q0
Q2 = P5 + P8

Then the determinant becomes

P0 + P1 + P2 + P3 + Q1 + Q2 + Q1x + (Q0 + Q2)y + Q0z.

That is how the "fast" method works. I could post the code if you want to see it, but I think that explanation will be more helpfull.

QED · #6 01-Mar-2006, 11:14

That does seem odd. Although, I assume that you are compiling with all speed optimizations on, in which case some of the redundancy that you manually eliminated should be likewise automatically eliminated by the compiler.

For example, if I write

CPP / C++ / C Code:

x = a * b + a * b;

the compiler-optimized code will not perform the multiplication a * b twice. I.e., writing

CPP / C++ / C Code:

c = a * b;
x = c + c;

will not be any faster.

As for your other mystery, I can't say that I have any insight. If you can post code, I'd be interested in browsing through it. Mostly for my own curiosity.

Matthew

Blake · #7 01-Mar-2006, 12:16

The kind of optomization I'm doing is (for example) if I have

CPP / C++ / C Code:

X = a*a*b*b;
Y = c*c*d*d;
Z = a*b*c*d;

then I calculate

CPP / C++ / C Code:

P = a*b;
Q = c*d;

and then use

CPP / C++ / C Code:

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:

CPP / C++ / C Code:

// 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:

CPP / C++ / C Code:


#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;

}

QED · #8 07-Mar-2006, 08:47

Hi, blake. Sorry for the delay. I think I can help you boost the speed of your code. Using g++ v4, I observed a 1.7 times speedup with a relatively few changes. I also have a reworked version of the code that runs about a little more than 2 times faster on my setup, but it includes some bigger changes.

My tests used the following polynomials:
A = 2x + 1
B = 3y + 4z
C = 5xy + 6yz
D = 7xz + 8xyz

There are three basic things I did to speed up the process. First, I eliminated the construction of temporary pair object wherever possible. For example, instead of

CPP / C++ / C Code:

terms.insert(make_pair(E, C));

I used

CPP / C++ / C Code:

terms[E] = C;

In the former, we make a temporary pair object, which means we copy the expression E; then this pair is again copied into the map, which means we make a second copy of E. Two call to the copy constructor for class expList. In the latter, we make just one call to the copy constructor. This operator should be used with care however, because if E is not already a key in the map, it is inserted along with value C, but if it does exist the value is overwritten. It is useful in your code in places where we know we are inserting something for the first time.

Another important change is the use of operator+= instead of operator+, or operator*= instead of operator*. For example, instead of

CPP / C++ / C Code:

x = x * b;

we can use

CPP / C++ / C Code:

x *= b;

For complex structures like expList and Polynomial, this can save some processing. The former requires a temporary Polynomial to be created inside the operator body which is then copied when the function returns. By using the assignment version of the operator, we save a call to the copy constructor. I might have expected a good optimizing compiler to be able to convert the former statement to the latter automatically, but g++ v4 does not. Upon reflection, I suppose it may not be practical for the compiler to do such an optimization, since it may be that not both operators are defined, or they are defined in such a way that the two statements are not identical.

The third change involves using C++ references wherever possible in the determinant function. That is, instead of

CPP / C++ / C Code:

Polynomial<COEFF> Q0(P10 + P11);

we can write

CPP / C++ / C Code:

Polynomial<COEFF>& Q0 = P10;
Q0 += P11;

The idea here is to avoid an unecessary call to the copy constructor, since Q0 is now jsut an alias for P10. This can be done if we know that we no longer need P10 beyond this point in the code.

I'll post you full code with my changes at the end. For the >2 times speedup mentioned above, I change methods like invert() to return void and operator directly on the calling object, and created getInverse() which implements the existing behavior of making a temporary object. I also replaced std::map with std::set and made the expList class become a Term class which includes the coefficient in it.

Here are some diagnostic outputs that illustrate where the speedup was achieved.
Original:

Code:

Determinant computed in 0.001879s
  term_multiplies: 222
  term_arr_constructs: 19
  term_cpy_constructs: 3869
  term_siz_constructs: 162

Modified(minor changes):

Code:

Determinant computed in 0.001099s
  term_multiplies: 222
  term_arr_constructs: 16
  term_cpy_constructs: 1333
  term_siz_constructs: 381

Modified(bigger changes):

Code:

Determinant computed in 0.000831s
  term_multiplies: 222
  term_arr_constructs: 16
  term_cpy_constructs: 1150
  term_siz_constructs: 0

Of course, with a different compiler on a different platform, these figures may differ slightly.

Here is the changed code:

CPP / C++ / C Code:

// polynomial.h
// Written by Blake Foster
// Last updated February 7, 2006

#ifndef POLYNOMIAL_H
#define POLYNOMIAL_H

#include <vector>
#include <set>
#include <map>
#include <cassert>
#include <iostream>
#include <string>
#include <algorithm>

// macro used for naming variables
#define alphabet "abcdefghijklmnopqrstuvwxyz"

using namespace std;

// 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:
  vector<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:
  typedef typename map<expList,COEFF>::iterator PolynomialIterator;
  typedef typename map<expList,COEFF>::iterator ConstPolynomialIterator;

  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>&);
  Polynomial<COEFF> operator-(const Polynomial<COEFF>&) const;
  Polynomial<COEFF> operator-() const;
  Polynomial<COEFF> operator*(const Polynomial<COEFF>&) const;
  Polynomial<COEFF>& operator*=(const Polynomial<COEFF>&);
  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 begin(); // return an iterator pointing to the
  // first term
  PolynomialIterator 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&, 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;
};

// expList class

inline expList::expList(unsigned int n, int exps[])
  : exponents(n)
{
  ++term_arr_constructs;
  for (unsigned int i=0; i<n; ++i) exponents[i] = exps[i];
}

inline expList::expList(const vector<int>& V)
  : exponents(V)
{
  ++term_vec_constructs;
}

inline expList::expList(const expList& E)
  : exponents(E.exponents)
{
  ++term_cpy_constructs;
}

inline expList::expList(const int& N)
  : exponents(N)
{
  ++term_siz_constructs;
}

inline expList::expList()
  : exponents()
{
  ++term_def_constructs;
}

inline expList::~expList()
{
}

inline int& expList::operator[](const int& i)
{
  return exponents[i];
}

inline bool expList::operator<(const expList& L) const
{
  assert (size()==L.size());
  for (unsigned int i=0; i<size(); ++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 (size()==L.size());
  for (unsigned int i=0; i<size(); ++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<size(); ++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<size(); ++i)
  {
    if (L.exponents[i]!=exponents[i]) return true;
  }
  return false;

}

inline void expList::operator=(const expList& E)
{
  exponents = E.exponents;
}

inline expList expList::operator+(const expList& L) const
{
  ++multiply_term_count;
  assert(size()==L.size());
  expList sum(size());
  transform(exponents.begin(), exponents.end(), L.exponents.begin(),
            sum.exponents.begin(), plus<int>());
  return sum;
}

inline void expList::operator+=(const expList& L)
{
  ++multiply_term_count;
  assert(size()==L.size());
  for (unsigned int i=0; i<size(); i++) exponents[i] += L.exponents[i];
}

inline expList expList::operator-() const
{
  expList inverse(size());
  for (unsigned int i=0; i<size(); ++i)
  {
    inverse.exponents[i] = -exponents[i];
  }
  return inverse;
}

inline unsigned int expList::size() const
{
  return exponents.size();
}

inline expList expList::conj(const unsigned int& V) const
{
  expList E(size());
  for (unsigned int i=0; i<size(); 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(size());
  for (unsigned int i=0; i<size(); i++)
  {
    E.exponents[i]=-exponents[i];
  }
  return E;
}

inline bool expList::constant() const
{
  for (unsigned int i=0; i<size(); i++) if (exponents[i]!=0) return false;
  return true;
}

inline void expList::printV(ostream& ost) const
{
  assert(size() > 0);
  ost<<'(';
  ost<<exponents[0];
  for (unsigned int i=1; i<size(); ++i)
  {
    ost<<", "<<exponents[i];
  }
  ost<<')';
}

inline void expList::print(ostream& ost) const
{
  for (unsigned int i=0; i<size(); ++i)
  {
    if (exponents[i]!=0)
    {
      if (size()<=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<size(); ++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> >()),
    numVars(E.size())
{
  //terms.insert(make_pair(E,C));
  terms[E] = C;
}

template <typename COEFF>
inline Polynomial<COEFF>::Polynomial(const int& n)
  : terms(map<expList,COEFF,less<expList> >(less<expList>())),
    numVars(n)
{}

template <typename COEFF>
inline Polynomial<COEFF>::Polynomial()
  : terms(map<expList,COEFF,less<expList> >(less<expList>())),
    numVars(0)
{
  ++default_construct_count;
}

template <typename COEFF>
inline void Polynomial<COEFF>::operator=(const Polynomial& P)
{
  ++assignment_count;
  terms=P.terms;
  numVars=P.numVars;
}

template <typename COEFF>
inline Polynomial<COEFF> Polynomial<COEFF>::operator+(const Polynomial& P) const
{
  if (terms.size() > P.terms.size())
  {
    Polynomial<COEFF> sum(*this);
    for (ConstPolynomialIterator iter = P.terms.begin();
         iter != P.terms.end(); ++iter)
    {
      sum.insert(iter->second,iter->first);
    }
    return sum;
  }
  else
  {
    Polynomial<COEFF> sum(P);
    for (ConstPolynomialIterator 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)
{
  for (ConstPolynomialIterator iter = P.terms.begin();
       iter != P.terms.end(); ++iter)
  {
    insert(iter->second,iter->first);
  }
  return *this;
}

template <typename COEFF>
inline Polynomial<COEFF> Polynomial<COEFF>::operator-(const Polynomial& P) const
{
  ++subtract_count;
  if (terms.size() > P.terms.size())
  {
    Polynomial<COEFF> diff(*this);
    for (ConstPolynomialIterator 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 (ConstPolynomialIterator 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 (PolynomialIterator 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 (ConstPolynomialIterator iter1 = terms.begin();
       iter1 != terms.end(); ++iter1)
  {
    for (ConstPolynomialIterator 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 Polynomial<COEFF>& Polynomial<COEFF>::operator*=(const Polynomial& P)
{
  assert(numVars==P.numVars);
  map<expList,COEFF,less<expList> > tmp(terms);
  terms.clear();
  for (ConstPolynomialIterator iter1 = tmp.begin();
       iter1 != tmp.end(); ++iter1)
  {
    for (ConstPolynomialIterator iter2 = P.terms.rbegin();
         iter2 != P.terms.rend(); ++iter2)
    {
      insert((iter1->second)*(iter2->second),(iter1->first+iter2->first));
    }
  }
  return *this;
}

template <typename COEFF>
inline void Polynomial<COEFF>::insert(const COEFF& C, const expList& E)
{
  assert(E.size() == numVars);
  if (C != ZERO)
  {
    PolynomialIterator pos = terms.find(E);
    if (pos == terms.end())
    {
      terms[E] = C;
    }
    else
    {
      pos->second += C;
      if (pos->second == 0)
      {
        ++reduce_terms_count;
        terms.erase(pos);
      }
    }
    //pair<typename 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);
    //  }
    //}
  }
}

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 (ConstPolynomialIterator iter = terms.begin();
         iter != terms.end(); ++iter)
    {
      //P.terms.insert(make_pair(iter->first.conj(V), iter->second));
      P.terms[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 (ConstPolynomialIterator iter = terms.begin();
         iter != terms.end(); ++iter)
    {
      //P.terms.insert(P.terms.begin(), make_pair(-iter->first, iter->second));
      P.terms[-iter->first] = iter->second;
    }
    return P;
  }
}

template <typename COEFF>
void Polynomial<COEFF>::printV() const
{
  if (terms.empty()) cout << "0";
  else
  {
    ConstPolynomialIterator 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
  {
    ConstPolynomialIterator 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
  {
    ConstPolynomialIterator 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 (ConstPolynomialIterator 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 terms.size();
}

template <typename COEFF>
inline void Polynomial<COEFF>::clear()
{
  terms.clear();
}

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);

template <typename COEFF>
inline PolynomialIterator Polynomial<COEFF>::begin()
{
  return terms.begin();
}

template <typename COEFF>
inline PolynomialIterator Polynomial<COEFF>::end()
{
  return 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 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 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 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 iter = terms.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& 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 (ConstPolynomialIterator iter = terms.begin();
       iter != terms.end(); ++iter)
  {
    //product.terms.insert(product.terms.end(),
    //                     make_pair(iter->first + E, iter->second));
    product.terms[iter->first + E] = iter->second;
  }
  return product;
}

#endif

Determinant function:

CPP / C++ / C Code:

template <typename COEFF_TYPE>
inline void
determinant(Polynomial<COEFF_TYPE> const& A,
            Polynomial<COEFF_TYPE> const& B,
            Polynomial<COEFF_TYPE> const& C,
            Polynomial<COEFF_TYPE> const& D,
            Polynomial<COEFF_TYPE>& R)
{
  const COEFF_TYPE ONE(1);

  Polynomial<COEFF_TYPE> Ax(A.CONJX);
  Polynomial<COEFF_TYPE> Bx(B.CONJX);
  Polynomial<COEFF_TYPE> Cx(C.CONJX);
  Polynomial<COEFF_TYPE> Dx(D.CONJX);
  Polynomial<COEFF_TYPE> By(B.CONJY);
  Polynomial<COEFF_TYPE> Dy(D.CONJY);
  Polynomial<COEFF_TYPE> Bz(B.CONJZ);
  Polynomial<COEFF_TYPE> Dz(D.CONJZ);

  Polynomial<COEFF_TYPE>& T0 = Ax; T0 *= A;
  Polynomial<COEFF_TYPE>& T1 = Bx; T1 *= B;
  Polynomial<COEFF_TYPE>& T2 = Cx; T2 *= C;
  Polynomial<COEFF_TYPE>& T3 = Dx; T3 *= D;
  Polynomial<COEFF_TYPE>& T4 = By; T4 *= A;
  Polynomial<COEFF_TYPE>& T5 = Dy; T5 *= C;
  Polynomial<COEFF_TYPE>& T6 = Dz; T6 *= A;
  Polynomial<COEFF_TYPE>& T7 = Bz; T7 *= C;

  Polynomial<COEFF_TYPE> T0y(T0.CONJY);
  Polynomial<COEFF_TYPE> T1y(T1.CONJY);
  Polynomial<COEFF_TYPE> T2y(T2.CONJY);
  Polynomial<COEFF_TYPE> T3y(T3.CONJY);
  Polynomial<COEFF_TYPE> T4z(T4.CONJZ);
  Polynomial<COEFF_TYPE> T5z(T5.CONJZ);
  Polynomial<COEFF_TYPE> T6y(T6.CONJY);
  Polynomial<COEFF_TYPE> T7y(T7.CONJY);

  Polynomial<COEFF_TYPE>& P0 = T0;  P0 *= T0y;
  Polynomial<COEFF_TYPE>& P1 = T1y; P1 *= T1;
  Polynomial<COEFF_TYPE>& P2 = T2y; P2 *= T2;
  Polynomial<COEFF_TYPE>& P3 = T3;  P3 *= T3y;
  Polynomial<COEFF_TYPE> P4(-T4.multByVars(makeTerm(0,0,1))*T4z);
  Polynomial<COEFF_TYPE>& P5 = T2; P5 *= -T0y.multByVars(makeTerm(-1,0,0));
  Polynomial<COEFF_TYPE> P6(-T6.multByVars(makeTerm(0,1,0))*T6y);
  Polynomial<COEFF_TYPE>& P7 = T7y; P7 *= -T7.multByVars(makeTerm(0,-1,0));
  Polynomial<COEFF_TYPE>& P8 = T3y; P8 *= -T1.multByVars(makeTerm(-1,0,0));
  Polynomial<COEFF_TYPE>& P9 = T5z; P9 *= -T5.multByVars(makeTerm(0,0,-1));
  Polynomial<COEFF_TYPE>& P10 = T7; P10 *= T6y.multByVars(makeTerm(0,0,-1));
  Polynomial<COEFF_TYPE>& P11 = T5; P11 *= T4z.multByVars(makeTerm(-1,-1,0));

  Polynomial<COEFF_TYPE>& Q0 = P10; Q0 += P11;
  Polynomial<COEFF_TYPE>& Q1 = P4;  Q1 += P6; Q1 += P7; Q1 += P9; Q1 += Q0;
  Polynomial<COEFF_TYPE>& Q2 = P5;  Q2 += P8;

  R = P0;
  R += P1;
  R += P2;
  R += P3;
  R += Q1;
  R += Q2;
  R += Q1.CONJX;
  R += (Q0 + Q2).CONJY;
  R += Q0.CONJZ;
}

A final note: I removed the numTerms variable in Polynomial since the terms map already has a size() method. I also removed the PolynomialIterator class and replaced it with a simpler typedef. There are many similar style changes that can be done to make the code more readable and maintainable. Also, some of your for loops can probably be replaced by STL algorithms. The is already an algorithm for lexicographical comparison

Cheers,
Matthew

P.S. Forgot to mention that I also made user of std::vector in your expList class. No sense in managing memory explicitly if you don't have to, I say.

Blake · #9 09-Mar-2006, 10:14

Wow, that's an amazing speed up! Thank you VERY much for the help. It looks to me like it's time for another major overhaul of the polynomial class. This time I'll write a monomial class that stores the exponent, and use an stl set instead of a map as you suggested. I'll also define the * operator on monomials to replace the multByVars function, which assumes the exponent is one. I should probably add a special case in the * operator for polynomials that have one term as well, since that won't change the order (thus it will run in linear time instead of nlog(n)).
I think I'm also going to switch from lex to deg-lex ordering on the monomials. In deg-lex, for terms A and B, we say A < B if the total degree of A is less than the total degree of B. If A and B have the same total degree, then we do a lexicographic compare on the exponents. If I have a data member that keeps track of the total degree, I think that will be a little faster.

Again, thank you so much for the help!

QED · #10 09-Mar-2006, 10:50

My pleasure.

The deg-lex ordering does sound like it could be faster; you'll get the biggest gain with polynomials in many variables.

Also, you will certainly find out by testing, but I suspect a special case for multiplication by a single term will not give much speedup. True that std::map insertion is at most log(n) complexity, but you can ensure that it is merely constant time if the insertion order is the same as the sort order, by using insert(pos, x) where pos is simply end(). I suppose handling it as a special case would save n-1 logical comparisons. Probably worth it if you absolutely have to make this as fast as possible.

A profiler could help pinpoint where the bulk of the processing is done (though you may already know that much). That'll tell you where to spend your efforts trying to improve efficiency.

Another thing you might consider is using std::deque instead of std::set. You lose the automatic ordering and check for duplicates that set gives you, but by manually sorting and checking for duplicates only at those times when you know it needs to be done (e.g., it doesn't need to be done when multiplying by a monomial!) you might be able to save some cycles.

Do you have a specific goal in mind for the execution time? When will you stop trying to optimize?

Generally optimization is best left until everything else is done, if even then. Although, I must admit, I'm a sucker for this sort of thing, too.

Matthew

Thread Tools	Search this Thread
Show Printable Version Email this Page	Search this Thread: Advanced Search
Rate This Thread
Rate This Thread:

determinant algorithms

Re: determinant algorithms

Re: determinant algorithms

Re: determinant algorithms

Re: determinant algorithms

Re: determinant algorithms

Re: determinant algorithms

Re: determinant algorithms

Re: determinant algorithms

Re: determinant algorithms