Hashing arrays

[tition]

Hi all,

I am writing at the moment template classes for polynomials. (Rewriting actually). I ran into some computational problems with the code I have below. I have tried to make the code as self-explanatory as possible.

//monomials are contained in Polynomial::TheMonomials

template <class Element>
void Polynomial<Element>::MultiplyBy(Polynomial<Element> &p, Polynomial<Element> &output) 
{
  if (p.size==0)
  { this->size=0;
    return;
  }	
  Polynomial<Element> Accum;
  Monomial<Element> tempM;
  Accum.size=0; 
  for (int i=0;i<p.size;i++)
  { for (int j=0;j<this->size;j++)
    { this->TheMonomials[j].MultiplyBy(p.TheMonomials[i],tempM);//tempM contains the result of the multiplication
      Accum.AddMonomial(tempM);
    }
  }
  output->CopyFrom(Accum); 
}

template <class Element>
void Polynomial<Element>::AddMonomial (Monomial<Element>& m) 
{ int j= this->HasSameExponentMonomial(m);
  if (j==-1)
  { this->AddObjectOnTop(m);  
  }
  else
  { this->TheMonomials[j].Coefficient.Add(m.Coefficient);
    if (this->TheMonomials[j].Coefficient.IsEqualTo(*Element::TheRingZero))
    { this->PopIndexSwapWithLast(j); 
    }
  }
}

The (partial) reason for my computational troubles is the following. Let the two polynomials I am multiplying have approximately n monomials. The two for cycles in the Polynomial::MultiplyBy take approximately n^2 operations - you can't get away with that, this is how you multiply polynomials (and matrices by the way, and that is no coincidence).

The trouble lies in Polynomial::AddMonomial: before I add, I must check whether a monomial is already present in my polynomial. Now, if I just cycle through all monomials that my polynomial has, that is n^3 operations. Since I need to apply this to polynomials of 10th degree in 5 variables, which can have max 15 choose 5 = 3003 monomials, that could mean I need 3003^3 operations per multiplication. That is not ok.

Now if I can make the time for the function Polynomial::HasSameExponentMonomial
approximately constant, I would be in great shape. So, my question to you is: what would you do? I have in mind something with a hash function, but only am in the planning for now. How would resolve the problem?

[jsmith]

I've seen this kind of homework assignment before on this board, and I actually implemented a full-fledged polynomial class. My implementation uses a vector<> as the underlying storage, where the nth index in the vector corresponds to the order-n term in the polynomial. [I believe your implementation might be trying to be memory efficient by only storing the orders that have non-zero coefficients, and this is leading to your problem.]
My solution relies on all orders 0 through N to exist with possibly zero coefficient. In this way lookup is constant.

[tition]

If I do what you say, then, to multiply two polynomials with 10^th degree 5 variables, I always need 3003^2 operations. Also, the coefficients of my monomials are not very well-mannered fellows - each is a linear combination of formal exponents of vectors with rational coefficients, so my coefficient adding time is not to be ignored!

[tition]

It just occured to me that

where the nth index in the vector corresponds to the order-n term in the polynomial.

corresponds to having "the hashing function" assign each number to itself :)

However, jsmith, I am dealing with polynomials in *many* variables - how many is user input (for practical purposes up to 6 but it's not hard-coded). The example I need is in 5 variables.

[jsmith]

Can you post an example polynomial?

I'm not sure I'm understanding you.

[turbozedd]

You can reduce multiply to n lg n time if you use FFT (and this is for polynomials assumed to have no zero terms). n lg n is of course for nth degree and one variable.

If you use std::map<order, monomial> you will be able to jump the zero monomials, but map uses some extra memory (you could use hashmap too). You can look at my Matrix class for details on multiplying http://www.cplusplus.com/forum/general/8352/

I would still recommend using the FFT with a vector or an array since most polys are not sparse. Designing something tricky will probably use too much memory for your polynomials.

[tition]

@ jsmith: Computing the polynomials is actually the objective of the program. If my previous posts got you interested, I must compute a polynomial for each of the 10890 "combinatorial chambers" I finally computed with my program. The polynomials in questions are in five variables - a linear one would look like a_0+ a_1 x_1 +a_2 x_2+...+a_5x_5.

Each of these polynomials is computed in a consecutive fashion by

1) Linear substitution adding one extra variable (with respect to which we "integrate"). Call the extra variable t.
2) regrouping with respect to the extra variable.
3) Substituting each power t^k back with a polynomial in 5 variables and k+1-st degree (in fact first I substitute t^k back with something called "Bernoulli polynomials", http://en.wikipedia.org/wiki/Bernoulli_polynomials , and in the new thing I substitute t back with a linear function in x_1,...,x_5).

4) we go back to step 1, 10 times total, starting with a constant polynomial.

I will post an example in the low-dimensional cases (3 variables or less) when I get my program back to running.

@turbozedd (28) Thanks for pointing out FFT, I skimmed over Wikipedia for the FFT in more than one variable but got little out of it. Do you have a simpler link? Also, I actually don't know how to apply the FFT for multiplying polynomials - where can I read a clear explanation of that?

I prefer a hash function to maps, since I also like to make "variable increase", which would mean I would have to reserve memory for 6-variable polynomials. Now, a 10-th degree polynomial in 5 variables can have 3003 monomials, while 10th degree polynomial in 6 variables can have 8008 monomials (more than double the RAM use!)

The ultimate goal of my endeavor is to compute all 10890 polynomials of 10th degree in 5 variables. Thanks a lot for the interesting comments!

[turbozedd]

There are lots of places to read other than wiki... but I think you will find they all say the same thing in the same way.

Since you are doing ALL 5th degree polys I'm pretty sure you will want a vector or an array of monomials. The fancy data structures are too expensive.

FFTW is a package for computing FFTs I think it has multidimensional support.

http://www.fftw.org/

http://www.fftw.org/doc/Multi_002ddimensional-Transforms.html

http://www.fftw.org/doc/Multi_002dDimensional-DFTs-of-Real-Data.html#Multi_002dDimensional-DFTs-of-Real-Data

Are your polys over the reals? will you extend to complex? or maybe restrict to Q or finite fields?

What do you meant by compute? do you just want a big list of polynomials?

[tition]

Thanks a lot for the references!

The polynomials are of 10th degree, not 5th - they are in 5 variables. The number 10890 comes from some geometrical partition of 5 dimensional space - this is the number of different "partition functions" (I am skipping the definition) I need to compute.

For the time being, yes, I want a big list of the polynomials. I will later need to do linear combinations of them - no multiplication, only addition, so no big worries.

In the case I ultimately need, the polynomials are over the rationals extended with formal exponents of vectors with rational coordinates (but there is an upper bound on the height of the rationals that can appear in the exponent). I use a complicated data structure to describe those "quasinumbers".

[turbozedd]

For the FFT you can always try "Numerical Recipes".

And I'm not sure, but you might have to just iterate over your monomial array to do it. So it really would be n lg n . In other words I think all you have to do is:

// do convolve p1 * p2 which is a term by term multiplication in fourier space

fsp1=fft(poly1);
fsp2=fft(poly2);

for(int i=0;i<n;i++){
    fspmult[i]=fsp1[i]*fsp2[i];
}

prodpoly=ifft(fspmult);

You have to be sure that each element of your poly array corresponds to the same term as the other poly array.

!! implementation detail: the coefficients in fourier space are complex !!

read about discrete convolution

[tition]

Now that my program is back to running, here is a typical example of the polynomials I have in mind:

(17/16+1/16e^{2i \pi(1/2X+1/2Y+1/4Z)})+(-1/8-1/8e^{2i \pi(1/2X+1/2Y+1/4Z)})c+(5/4+1/4e^{2i \pi(1/2X+1/2Y+1/4Z)})b-a^2+ab


i stands, as usual, for the imaginary unit. \pi stands for, well, pi :). e stands for, well, e :)

This is a polynomial in three variables - a, b and c. The capital letters X, Y, Z stand for integer parameters, which must be kept as such.


An example of working with the coefficients:

e^{2i\pi (X/2+Y/2)} * e^{2i\pi (X/2)} = e^{2i\pi (X+Y/2)} = e^{2i\pi (Y/2)}


The last equality comes from the usual e^{2i \pi }=1 and the fact that X is an integer parameter.

Other than the difficulty of adding and multiplying the coefficients, operating with the variables a, b, c is just like the usual polynomial operations (hence the template class).

Note: the polynomial above is an actual print out of my program in a low-dimensional case.