Need help optimizing....

[JustSomeGuy]

I have a routine that evaluates a polynomial equation that have 3 variables
x,y,z of orders 1,2,3 the coefficients of the polynomial are in an array.
This routine is quite slow and I'd like to optimize it.
Any suggestions?

simply put pts is a vector of
typdef struct
double x;
double y;
double z;
double val;
} vect

xorder, yorder and zorder are the maximum polynomial exponents for an
equation
like val = fn(x,y,z).

void EvaluatePoly(Vector<vect> & pts, int xorder, int yorder, int
zorder)
{
for (int pn=0; pn < pts.length(); ++pn)
{
int cn=0;
for (int k=0; k <= zorder; ++k)
{
float zexp = pow(pts[pn].z, k);
for (int j=0; j <= yorder; ++j)
{
float zyexp = zexp * pow(pts[pn].y, j);
for (int i=0; i <= xorder; ++i)
{
pts[pn].val = c[cn++] * pow(pts[pn].x, i) * zyexp;
}
}
}
}
}

 

[Kevin Goodsell]

JustSomeGuy wrote:
<snip>

Why did you cross-post this to several groups? Massive cross-posting is
very bad. I don't know the topics of all the groups, but I very much
doubt your question is topical in all of them.

If you feel the need to cross post to 2 or 3 groups, stop and think very
carefully about whether your message is topical in all of them. If you
feel the need to cross post to *more* than 2 or 3 groups, don't.

-Kevin

 

[David B. Held]

I have a routine that evaluates a polynomial equation that
have 3 variables x,y,z of orders 1,2,3 the coefficients of the
polynomial are in an array. This routine is quite slow and I'd
like to optimize it. Any suggestions?
[...]

Yes. 1) Don't cross-post to 9 newsgroups. You will almost
certainly receive a number of angry responses, which will
not help you solve your optimization problem (this is a
meta-optimization tip for the next time you have a question).
2) Take a look at Blitz++ and/or the GNU Scientific Library.
While I am not sure either one is totally relevant to your
problem, they might be. But your question is not a C++-
specific question (at least, not the way you framed it).
However, if you try Blitz++, it might become one. My
guess is that comp.programming would have been the
most appropriate group to try first.

Dave

 

[JustSomeGuy]

I have a routine that evaluates a polynomial equation that
have 3 variables x,y,z of orders 1,2,3 the coefficients of the
polynomial are in an array. This routine is quite slow and I'd
like to optimize it. Any suggestions?
[...]

Yes. 1) Don't cross-post to 9 newsgroups. You will almost
certainly receive a number of angry responses, which will
not help you solve your optimization problem (this is a
meta-optimization tip for the next time you have a question).

People who get 'Angry' in regards to news group posting need to
seek help desperatly. But thanks for the info I will take in the
good spirit.

 

[Gianni Mariani]

I have a routine that evaluates a polynomial equation that have 3 variables
x,y,z of orders 1,2,3 the coefficients of the polynomial are in an array.
This routine is quite slow and I'd like to optimize it.
Any suggestions?

Now that you have been significantly chastized for the massive cross
post, comp.lang.c++ is NOT the best NG for numerical stuff.

Not only that, optimization questions not relating to the C++ standard
are off topic in comp.lang.c++.

I'll however give you a few hints.

a) It seems like you may have done SOME profiling but let me state that
YOU REALLY NEED TO PROFILE THE CODE.

b) Let's assume you did some profiling and found that EvaluatePoly is
the culprit. You would have found that pow() was called a large number
of times more than EvaluatePoly.

c) Assuming EvaluatePoly needs to scream, I'd suggest making a small
test harness which allows you to run EvaluatePoly as a stand-alone app
and this will allow you to iterate your changes faster.

simply put pts is a vector of
typdef struct
double x;
double y;
double z;
double val;
} vect

xorder, yorder and zorder are the maximum polynomial exponents for an
equation
like val = fn(x,y,z).

void EvaluatePoly(Vector<vect> & pts, int xorder, int yorder, int
zorder)
{
for (int pn=0; pn < pts.length(); ++pn)
{
int cn=0;
for (int k=0; k <= zorder; ++k)
{
float zexp = pow(pts[pn].z, k);

^^^^^^^^^^^^^^^^^^^^^^^^^ why floats ?
^^^^^^^^^^^^^^^^^^^^^^^^^ why call pow ? These may be computed using
multiplication per iterarion. Significantly fast

for (int j=0; j <= yorder; ++j)
{
float zyexp = zexp * pow(pts[pn].y, j);

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ why floats ?

for (int i=0; i <= xorder; ++i)
{
pts[pn].val = c[cn++] * pow(pts[pn].x, i) * zyexp;

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
what is this ? don't you want this to be a sum of all the terms ?

}
}
}
}
}

So, give this a go. Compared to what is above, it should be
significantly faster and it might actually work.

Without knowing the data set, it's kind of hard to know better what you
need.


#include <vector>

typedef struct {
double x;
double y;
double z;
double val;
} vect;


void EvaluatePoly(
std::vector<vect> & pts,
int xorder,
int yorder,
int zorder,
double * c
)
{

std::vector<vect>::iterator p_pts = pts.begin();
std::vector<vect>::iterator pts_end = pts.end();

for (; p_pts != pts_end; ++ p_pts )
{
int cn=0;

double val_x = p_pts->x;
double val_y = p_pts->y;
double val_z = p_pts->z;

double sum = 0;
double pow_z = 1;
for (int k=0; k <= zorder; ++k)
{
double pow_y = pow_z;
for (int j=0; j <= yorder; ++j)
{
double pow_x = pow_y;
for (int i=0; i <= xorder; ++i)
{
sum += c[cn++] * pow_x;

pow_x *= val_x;
}

pow_y *= val_y;
}

pow_z *= val_z;
}

p_pts->val = sum;
}
}

 

[Robert Vienneau]

It may be helpful to realize:

a*x^3 + b*x^2 + c*x + d = d + x*(c + b*(x + a*x))

and so on.

I think that's called Horner's rule, and it is a standard trick when
optimizing the evalutaion of polynomials.

--
Try http://csf.colorado.edu/pkt/pktauthors/Vienneau.Robert/Bukharin.html
To solve Linear Programs: .../LPSolver.html
r c A game: .../Keynes.html
v s a Whether strength of body or of mind, or wisdom, or
i m p virtue, are found in proportion to the power or wealth
e a e of a man is a question fit perhaps to be discussed by
n e . slaves in the hearing of their masters, but highly
@ r c m unbecoming to reasonable and free men in search of
d o the truth. -- Rousseau

 

[LibraryUser]

I have a routine that evaluates a polynomial equation that
have 3 variables x,y,z of orders 1,2,3 the coefficients of
the polynomial are in an array. This routine is quite slow
and I'd like to optimize it. Any suggestions?
[...]

Yes. 1) Don't cross-post to 9 newsgroups. You will almost
certainly receive a number of angry responses, which will
not help you solve your optimization problem (this is a
meta-optimization tip for the next time you have a question).

People who get 'Angry' in regards to news group posting need to
seek help desperatly. But thanks for the info I will take in the
good spirit.

Inconsiderate people who do such massive crossposting need the
help. This leads to interminable off-topic threads and
continuous annoyance. Once the original post is made the harm is
done.

Usually any originally cross-posted queries should include a
followup to the single most appropriate newsgroup. That allows a
query to be spread over 3 or 4 groups, without the penalties.
Note that many of the better ISPs simply destroy any messages
crossposted to an excessive extent. If not, the better
newsreaders will automatically reject them. So the whole
business is counter-productive.

 

[JustSomeGuy]

It may be helpful to realize:

a*x^3 + b*x^2 + c*x + d = d + x*(c + b*(x + a*x))

and so on.

I think that's called Horner's rule, and it is a standard trick when
optimizing the evalutaion of polynomials.

I think your are right....

horner's rule in c.

n=5, p=c[n];
for (j=n-1; j >=0; j=j-1)
p *= x + c[j];

however I wonder how well, or even if its possible, to do this
with and equation with 3 variables, as in my example?

check this out.
http://www.dspguru.com/comp.dsp/tricks/alg/horner.htm

ps someone was awake enough to see the error
pts[pn].val = c[cn++] * pow(pts[pn].x, i) * zyexp;
which should have been:
pts[pn++].val = c[cn++] * pow(pts[pn].x, i) * zyexp;

So if any one has any suggestions on using this trick with more than one
variable let me know please.

 

[JustSomeGuy]

I have a routine that evaluates a polynomial equation that have 3 variables
x,y,z of orders 1,2,3 the coefficients of the polynomial are in an array.
This routine is quite slow and I'd like to optimize it.
Any suggestions?

Now that you have been significantly chastized for the massive cross
post, comp.lang.c++ is NOT the best NG for numerical stuff.

Not only that, optimization questions not relating to the C++ standard
are off topic in comp.lang.c++.

I'll however give you a few hints.

a) It seems like you may have done SOME profiling but let me state that
YOU REALLY NEED TO PROFILE THE CODE.

b) Let's assume you did some profiling and found that EvaluatePoly is
the culprit. You would have found that pow() was called a large number
of times more than EvaluatePoly.

c) Assuming EvaluatePoly needs to scream, I'd suggest making a small
test harness which allows you to run EvaluatePoly as a stand-alone app
and this will allow you to iterate your changes faster.

simply put pts is a vector of
typdef struct
double x;
double y;
double z;
double val;
} vect

xorder, yorder and zorder are the maximum polynomial exponents for an
equation
like val = fn(x,y,z).

void EvaluatePoly(Vector<vect> & pts, int xorder, int yorder, int
zorder)
{
for (int pn=0; pn < pts.length(); ++pn)
{
int cn=0;
for (int k=0; k <= zorder; ++k)
{
float zexp = pow(pts[pn].z, k);

^^^^^^^^^^^^^^^^^^^^^^^^^ why floats ?
^^^^^^^^^^^^^^^^^^^^^^^^^ why call pow ? These may be computed using
multiplication per iterarion. Significantly fast

for (int j=0; j <= yorder; ++j)
{
float zyexp = zexp * pow(pts[pn].y, j);

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ why floats ?

for (int i=0; i <= xorder; ++i)
{
pts[pn].val = c[cn++] * pow(pts[pn].x, i) *

zyexp;
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
what is this ? don't you want this to be a sum of all the terms ?

}
}
}
}
}

So, give this a go. Compared to what is above, it should be
significantly faster and it might actually work.

Without knowing the data set, it's kind of hard to know better what you
need.

Yes I've gotten rid of the pow function and I've also started looking
into horner's rule, but I'm not sure I can apply it to polynomials with more
than one variable. If anyone knows... please let me know.

#include <vector>

typedef struct {
double x;
double y;
double z;
double val;
} vect;


void EvaluatePoly(
std::vector<vect> & pts,
int xorder,
int yorder,
int zorder,
double * c
)
{

std::vector<vect>::iterator p_pts = pts.begin();
std::vector<vect>::iterator pts_end = pts.end();

for (; p_pts != pts_end; ++ p_pts )
{
int cn=0;

double val_x = p_pts->x;
double val_y = p_pts->y;
double val_z = p_pts->z;

double sum = 0;
double pow_z = 1;
for (int k=0; k <= zorder; ++k)
{
double pow_y = pow_z;
for (int j=0; j <= yorder; ++j)
{
double pow_x = pow_y;
for (int i=0; i <= xorder; ++i)
{
sum += c[cn++] * pow_x;

pow_x *= val_x;
}

pow_y *= val_y;
}

pow_z *= val_z;
}

p_pts->val = sum;
}
}

 

[Gianni Mariani]

Yes I've gotten rid of the pow function and I've also started looking
into horner's rule, but I'm not sure I can apply it to polynomials with more
than one variable. If anyone knows... please let me know.

Horner's rule only effects the inner loop but in this case it does not
save you anything because there are still 2 multiplies and an add per
inner loop iteration.

Here it is,

horner_sum = *(coef--)*pow_y + horner_sum * val_x;

- 2 multiplies and an add

The first code I wrote

sum += c[cn++] * pow_x;
pow_x *= val_x;

- two multiplies and an add !

I would go for simple.


#include <cmath>
#include <vector>

typedef struct {
double x;
double y;
double z;
double val;
} vect;


void EvaluatePoly(
std::vector<vect> & pts,
int xorder,
int yorder,
int zorder,
double * c
)
{

std::vector<vect>::iterator p_pts = pts.begin();
std::vector<vect>::iterator pts_end = pts.end();

for (; p_pts != pts_end; ++ p_pts )
{
int cn=0;

double val_x = p_pts->x;
double val_y = p_pts->y;
double val_z = p_pts->z;

double sum = 0;
double pow_z = 1;
for (int k=0; k <= zorder; ++k)
{
double pow_y = pow_z;
for (int j=0; j <= yorder; ++j)
{
double horner_sum = 0;
double * coef = c + cn + xorder + 1;
for (int i=xorder; i >= 0; i --)
{
horner_sum = *(coef--)*pow_y + horner_sum * val_x;
}

cn += xorder + 1;
sum += horner_sum;

pow_y *= val_y;
}

pow_z *= val_z;
}

p_pts->val = sum;
}
}

 

[Doug Norris]

People who get 'Angry' in regards to news group posting need to
seek help desperatly. But thanks for the info I will take in the
good spirit.

People who read what David wrote to you and infer that he's "angry"
need to either (1) read it again, or (2) read it again. HTH.

Doug

 

[jeffc]

People who read what David wrote to you and infer that he's "angry"
need to either (1) read it again, or (2) read it again. HTH.

I don't think SomeGuy inferred that David was angry any more than David
implied it to begin with. They were both speaking in the third person, no?

 

[Slartibartfast]

People who read what David wrote to you and infer that he's "angry"
need to either (1) read it again, or (2) read it again. HTH.

People who think that "JustSomeGuy" was referring to David himself,
rather than to the senders of the "number of angry responses" to which
David refers, need to either (1) read it again, or (2) read it again. HTH.

 

[Robert Vienneau]

I think your are right....

horner's rule in c.

n=5, p=c[n];
for (j=n-1; j >=0; j=j-1)
p *= x + c[j];

however I wonder how well, or even if its possible, to do this
with and equation with 3 variables, as in my example?

I did not look at your code too closely. But frequently one can
reverse the order of nested loops. Can you do such? And will it
make some such trick as Horner's rule useful?

check this out.
http://www.dspguru.com/comp.dsp/tricks/alg/horner.htm

--
Try http://csf.colorado.edu/pkt/pktauthors/Vienneau.Robert/Bukharin.html
To solve Linear Programs: .../LPSolver.html
r c A game: .../Keynes.html
v s a Whether strength of body or of mind, or wisdom, or
i m p virtue, are found in proportion to the power or wealth
e a e of a man is a question fit perhaps to be discussed by
n e . slaves in the hearing of their masters, but highly
@ r c m unbecoming to reasonable and free men in search of
d o the truth. -- Rousseau

 

[JustSomeGuy]

I don't think SomeGuy inferred that David was angry any more than David
implied it to begin with. They were both speaking in the third person, no?

Correct, I didn't think Dave was angry... Simply stating a fact.

 

[JustSomeGuy]

I think your are right....

horner's rule in c.

n=5, p=c[n];
for (j=n-1; j >=0; j=j-1)
p *= x + c[j];

however I wonder how well, or even if its possible, to do this
with and equation with 3 variables, as in my example?

I did not look at your code too closely. But frequently one can
reverse the order of nested loops. Can you do such? And will it
make some such trick as Horner's rule useful?

Yes I certainly can (And have) reversed the orders of the loops.
It really depends on weather I want the highest order polynomials
on the left or on the right... horner's rule I believe requires the
highest order polynomial on the right..

 

[mjm]

It may be helpful to realize:

a*x^3 + b*x^2 + c*x + d = d + x*(c + b*(x + a*x))

and so on.

I think that's called Horner's rule, and it is a standard trick when
optimizing the evalutaion of polynomials.

I think your are right....

horner's rule in c.

n=5, p=c[n];
for (j=n-1; j >=0; j=j-1)
p *= x + c[j];

however I wonder how well, or even if its possible, to do this
with and equation with 3 variables, as in my example?

check this out.
http://www.dspguru.com/comp.dsp/tricks/alg/horner.htm

ps someone was awake enough to see the error
pts[pn].val = c[cn++] * pow(pts[pn].x, i) * zyexp;
which should have been:
pts[pn++].val = c[cn++] * pow(pts[pn].x, i) * zyexp;

So if any one has any suggestions on using this trick with more than one
variable let me know please.

Try this: a polynomial p(x,y) in 2 variables x,y is a polynomial in y
with coefficients c_j(x) in R[x] (polynomials in x), for example

x^2y^2+3xy^2+2xy+y+x+1 = y^2(x^2+3x) + y(2x+1) +(x+1)

Compute the coefficients c_j(x} using Horner's rule and then with
these
compute p(x,y) again using Horner's rule. The same principle applies
to polynomials in 3 or more variables. You could build this up as
follows:

1. a routine for p(x)
2. a routine for p(x,y) using the routine for p(x) on the coefficients
c_j(x)
3. a routine for p(x,y,z) using the routine for p(x,y) on the
coefficients c_ij(x,y)

and son on.