CODEWITHSUNDEEP
arraysalgorithmmath
ANSI C Mastah
ANSI C Mastah

Reputation: 149

Compact storage coefficients of a multivariate polynomial

The setup

I am writing a code for dealing with polynomials of degree n over d-dimensional variable x and ran into a problem that others have likely faced in the past. Such polynomial can be characterized by coefficients c(alpha) corresponding to x^alpha, where alpha is a length d multi-index specifying the powers the d variables must be raised to.

The dimension and order are completely general, but known at compile time, and could be easily as high as n = 30 and d = 10, though probably not at the same time. The coefficients are dense, in the sense that most coefficients are non-zero.

The number of coefficients required to specify such a polynomial is n + d choose n, which in high dimensions is much less than n^d coefficients that could fill a cube of side length n. As a result, in my situation I have to store the coefficients rather compactly. This is at a price, because retrieving a coefficient for a given multi-index alpha requires knowing its location.

The question

Is there a (straightforward) function mapping a d-dimensional multi-index alpha to a position in an array of length (n + d) choose n?

Upvotes: 3

Views: 566

Answers (2)

Tom De Caluwé
Tom De Caluwé

Reputation: 926

Ordering combinations

A well-known way to order combinations can be found on this wikipedia page. Very briefly you order the combinations lexically so you can easily count the number of lower combinations. An explanation can be found in the sections Ordering combinations and Place of a combination in the ordering.

Precomputing the binomial coefficients will speed up the index calculation.

Associating monomials with combinations

If we can now associate each monomial with a combination we can effectively order them with the method above. Since each coefficient corresponds with such a monomial this would provide the answer you're looking for. Luckily if

alpha = (a[1], a[2], ..., a[d])

then the combination you're looking for is

combination = (a[1] + 0, a[1] + a[2] + 1, ..., a[1] + a[2] + ... + a[d] + d - 1)

The index can then readily be calculated with the formula from the wikipedia page.

Upvotes: 3

duffymo
duffymo

Reputation: 308928

A better, more object oriented solution, would be to create Monomial and Polynomial classes. The Polynomial class would encapsulate a collection of Monomials. That way you can easily model a pathological case like 

y(x) = 1.0 + x^50

using just two terms rather than 51.

Another solution would be a map/dictionary where the key was the exponent and the value is the coefficient. That would only require two entries for my pathological case. You're in business if you have a C/C++ hash map.

Personally, I don't think doing it the naive way with arrays is so terrible, even with a polynomial containing 1000 terms. RAM is cheap; that array won't make or break you.

Upvotes: 1

Related Questions