# Need help optimizing....

Discussion in 'C++' started by JustSomeGuy, Sep 22, 2003.

1. ### JustSomeGuyGuest

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

JustSomeGuy, Sep 22, 2003

2. ### Kevin GoodsellGuest

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
--
My email address is valid, but changes periodically.

Kevin Goodsell, Sep 22, 2003

3. ### David B. HeldGuest

"JustSomeGuy" <> wrote in message
news:MBtbb.8974\$TM4.4494@pd7tw2no...
> 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
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

David B. Held, Sep 22, 2003
4. ### JustSomeGuyGuest

"David B. Held" <> wrote in message
news:bklp3g\$jsh\$...
> "JustSomeGuy" <> wrote in message
> news:MBtbb.8974\$TM4.4494@pd7tw2no...
> > 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
> 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.

> 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.
>

Thanks.
B.

> Dave
>
>

JustSomeGuy, Sep 22, 2003
5. ### Gianni MarianiGuest

Re: Need help optimizing.... [Off topic]

JustSomeGuy wrote:
> 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;
}
}

Gianni Mariani, Sep 22, 2003
6. ### Robert VienneauGuest

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.

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

Robert Vienneau, Sep 22, 2003
7. ### LibraryUserGuest

JustSomeGuy wrote:
> "David B. Held" <> wrote in message
> > "JustSomeGuy" <> wrote in message
> >
> > > 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
> > 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
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

--
Replies should be to the newsgroup
Chuck Falconer, on vacation.

LibraryUser, Sep 22, 2003
8. ### JustSomeGuyGuest

"Robert Vienneau" <> wrote in message
news:...
> 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.
>

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

JustSomeGuy, Sep 23, 2003
9. ### JustSomeGuyGuest

Re: Need help optimizing.... [Off topic]

"Gianni Mariani" <> wrote in message
news:bkm68m\$...
> JustSomeGuy wrote:
> > 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;
> }
> }
>

JustSomeGuy, Sep 23, 2003
10. ### Gianni MarianiGuest

Re: Need help optimizing.... [Off topic]

JustSomeGuy wrote:
> "Gianni Mariani" <> wrote in message
> news:bkm68m\$...
>
>>JustSomeGuy wrote:
>>
>>>I have a routine that evaluates a polynomial equation that have 3

....
>>
>>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.
>

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

Gianni Mariani, Sep 23, 2003
11. ### Doug NorrisGuest

"JustSomeGuy" <> writes:

>"David B. Held" <> wrote in message
>news:bklp3g\$jsh\$...
>>
>> Yes. 1) Don't cross-post to 9 newsgroups. You will almost
>> certainly receive a number of angry responses, which will
>> 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.

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

Doug Norris, Sep 23, 2003
12. ### jeffcGuest

"Doug Norris" <> wrote in message
news:...
> "JustSomeGuy" <> writes:
>
> >"David B. Held" <> wrote in message
> >news:bklp3g\$jsh\$...
> >>
> >> Yes. 1) Don't cross-post to 9 newsgroups. You will almost
> >> certainly receive a number of angry responses, which will
> >> 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.

>
> 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?

jeffc, Sep 23, 2003
13. ### SlartibartfastGuest

"Doug Norris" <> wrote in message news:...
> "JustSomeGuy" <> writes:
>
> >"David B. Held" <> wrote in message
> >news:bklp3g\$jsh\$...
> >>
> >> Yes. 1) Don't cross-post to 9 newsgroups. You will almost
> >> certainly receive a number of angry responses, which will
> >> 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.

>
> 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.

--
#include <stdio.h>
char*f="#include <stdio.h>%cchar*f=%c%s%c;%cint main(void){printf(f,10,34,f,34,10,10);return 0;}%c";
int main(void){printf(f,10,34,f,34,10,10);return 0;}

Slartibartfast, Sep 23, 2003
14. ### Robert VienneauGuest

In article <bko9op\$1348\$>, "JustSomeGuy"
<> wrote:

> "Robert Vienneau" <> wrote in message
> news:...
> > 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?

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

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

Robert Vienneau, Sep 23, 2003
15. ### JustSomeGuyGuest

"jeffc" <> wrote in message
news:...
>
> "Doug Norris" <> wrote in message
> news:...
> > "JustSomeGuy" <> writes:
> >
> > >"David B. Held" <> wrote in message
> > >news:bklp3g\$jsh\$...
> > >>
> > >> Yes. 1) Don't cross-post to 9 newsgroups. You will almost
> > >> certainly receive a number of angry responses, which will
> > >> 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.

> >
> > 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?
>
>

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

JustSomeGuy, Sep 24, 2003
16. ### JustSomeGuyGuest

"Robert Vienneau" <> wrote in message
news:...
> In article <bko9op\$1348\$>, "JustSomeGuy"
> <> wrote:
>
> > "Robert Vienneau" <> wrote in message
> > news:...
> > > 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?

>
> 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..

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

>
> --
> 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, Sep 24, 2003
17. ### mjmGuest

"JustSomeGuy" <> wrote in message news:<bko9op\$1348\$>...
> "Robert Vienneau" <> wrote in message
> news:...
> > 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.

mjm, Sep 24, 2003