Numerical Recipes Forum  

Go Back   Numerical Recipes Forum > Obsolete Editions Forum > Methods: Chapters 9 and 10

Reply
 
Thread Tools Display Modes
  #1  
Old 12-16-2011, 01:24 PM
willus willus is offline
Registered User
 
Join Date: Dec 2011
Posts: 1
SQR compiler ambiguity in powell (sequence points)

I apologize--I don't have NR 3rd edition--only NR 1st and 2nd editions. I did a quick search on this issue and didn't come up with anything.

I wanted to point out an interesting change between the 1st and 2nd editions of NR which (apparently) fixes a sequence point issue, but I'm not sure if the authors intended that or just got lucky. In the 1st and 2nd editions, the powell() function, when determining if it needs to calculate a new direction set, uses the SQR() macro twice in one expression:

t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del)-del*SQR(fp-fptt);

In NR 1st edition, SQR is defined like so:
static float sqrarg;
#define SQR(a) (sqrarg=(a),sqrarg*sqrarg)

In NR 2nd edition, SQR is defined as so in nrutil.h:
static float sqrarg;
#define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 : sqrarg*sqrarg)

That change in the SQR definition between editions is enough that gcc 4.x compiles the first edition powell() function incorrectly (or rather, correctly but with an unintended result) but the second edition correctly (intended result). I discovered this because I was using the code from the 1st edition and my direction set was never changing (because t was always >= 0). I submitted this issue to gcc bugzilla as a possible gcc bug, but was told by their experts (since confirmed by some research) that because there is no "sequence point" between the evaluations of SQR() in the t=... expression above, the SQR expressions can legally be ordered by the compiler such that the assignments to sqrarg are both done first, and then the evaluations of sqrarg*sqrarg both done after. I am not sure if the authors just got lucky with the 2nd edition of NR and there is still a possibility that the compiler could legally order the evaluations in an unintended order, or if the ?: operator in the SQR definition forces the correct ordering. In either case, I'd strongly encourage the authors to replace the t=... assignment above with something like this (if it hasn't already been done in the latest edition), as I did:

...
{
double delsq1,delsq2;
delsq1 = SQR(fp-(*fret)-del);
delsq2 = SQR(fp-fptt);
t=2.0*(fp-2.0*(*fret)+fptt)*delsq1-del*delsq2;
...

This would remove any compiler or sequence-point ambiguity and ensure the intended execution.

There may be other expressions in other functions in the book that have a similar issue. I haven't checked.

For more on sequence points, Google "sequence points c programming".

Here is a simple sequence point example to test on your compiler if you want:

#include <stdio.h>
void main(void)
{
int x;
printf("1^2 + 2^2 = %d\n",(x=1,x*x)+(x=2,x*x));
}

Gcc 4.5.2-compiled exe reports this:

1^2 + 2^2 = 8

Compiling with Tiny CC (tinycc.org), on the other hand, reports this:

1^2 + 2^2 = 5

Both results are legally allowed.
Reply With Quote
Reply

Tags
ambiguity, powell, sequence points, sqr

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is Off
HTML code is Off

Forum Jump


All times are GMT -5. The time now is 10:55 PM.


Powered by vBulletin® Version 3.8.6
Copyright ©2000 - 2014, Jelsoft Enterprises Ltd.