Hello;

There appears to be a bug in Subroutine DDPOLY ( C, NC, X, PD, ND) as listed in NR, F77, page 138.

For a 3rd deg poly with the following (printed within the subroutine):
NC= 4
C = -4886568.5 1.5202154E+7 -1.5764277E+7 5448923.
X = 0.9595
ND=2

the routine returns the value of the poly at X as:
PD(1) = -1.0
and the 1st derivative at X as:
PD(2) = -2.0

The correct values should be (calculated manually, and by Excel VBA, and also by a computer Algebraic package):
PD(1) = -1.72061
PD(2) = 176.2873926

Has anyone identified the error in the code of DDPOLY() ??

Your help would be greatly appreciated.
Monir

Hello,

I am trying to duplicate the results that you have shown.

I am assuming that your coefficients are are loaded into the coefficient vector in this order

-4886568.5, 1.5202154E+7, -1.5764277E+7, 5448923.

I am not using FORTRAN, instead I am using C++ with the NR3 routines.

First I will use the routine poly.h to evaluate the 3rd order polynomial to see what the output from that routine is in order to try and track down a possible error in the ddpoly().

My output from Poly() with x=0.9595 is 0.015367. I recieve the same output from Excel.

My output from ddpoly() in pd[0] which is equivalent to pd[1] in FORTRAN is 0.015367.

Since the results from Poly() and ddpoly() are the same, I do not think there is a mistake in the two routines in the C++ language. However, that doesn't really help you out all that much since you are using FORTRAN but it does show that something is wrong because we are getting different values from Excel.

My output from ddpoly() in pd[1] (the first derivative) is -0.0641478

I am not able to duplicate your results of

PD(1) = -1.72061
PD(2) = 176.2873926

Is there something that I am not doing correctly? Is my oder of the coefficients incorrect?

Is c[0] = -4886568.5 or equivalently for FORTRAN c[1] = -4886568.5 ?

Can you post the equation in its true form?

Thanks
Matt

Hi Matt;

Thank you kindly for your prompt reply and for taking the time to check it out.

1) The four points are:
x y
0.9595 -1.72061
0.9640 -2.0312
0.9685 -2.5635
0.9730 -0.3379

2) The poly is:
...y =a0 +a1*x +a2*x^2 +a3*x^3
and the poly coeff {input array C(1:NC) in Sub DDPOLY(), and also confirmed by Excel} are:
NC = 4
a0 = -4886554.376442
a1 = 15202246.681715
a2 = -15764705.662669
a3 = 5449251.131605
ND = 2

3) Back substitution at x = 0.9595 using the above poly coeff gives:
poly value = -1.72061 (which agrees with the y-value in item 1 above)
1st derivative = 176.2873926

4) As indicated in my OP, DDPOLY() returns instead the value of the poly at x = 0.9595 as:
PD(1) = -1.0
and the 1st derivative as:
PD(2) = -2.0
(just having both values returned as "integers" by DDPOLY is suspicious!)

I wonder if in your NR3 C++ version of ddpoly(), the poly has the same representation as in item 2 above!

Thank you.
Monir

Monir,

These are the values in your second post.

a0 = -4886554.376442
a1 = 15202246.681715
a2 = -15764705.662669
a3 = 5449251.131605

The values in your first post are;

C = -4886568.5, 1.5202154E+7, -1.5764277E+7, 5448923.

As you can see they are not the same in each post.

With this set of values,

a0 = -4886554.376442
a1 = 15202246.681715
a2 = -15764705.662669
a3 = 5449251.131605

From Excel and Poly() from NR3 the result for the value of x = 0.9595 is -1.72060545

The result for the first derivative is pd[1] = 176.287.

Now we have the same results. Therefore the NR3 routines are working correctly.

I will try to duplicate your results with FORTRAN now.

PS:

Have you checked your code for typo(s)?


Thanks
Matt

Matt;

Glad your NR3 C++ version of DDPOLY() confirms the Excel results.

The F77 DDPOLY() is about 12 lines of code. I'm using "implicit none" to identify any typo, and have just re-checked the code for the 100th time!!

There's a much complex NR F90 version of DDPOLY() in: "Numerical Recipes in Fortran 90, Second Edition (1996), Chapter B5, p. 1071". It requires however a number of subsidiary and nr utility routines in the form of Modules. It would require major (and risky) changes to the code's structure, logic and syntax, and thus I don't think incorporating the F90 version would justify the required effort.

In the meantime, I'm anxious to know if your were able to get the F77 DDPOLY() to return the correct values for the quoted example!!
I'm still of the opinion that the code listing of the F77 DDPOLY() routine in p. 138 is at fault. But I could be wrong!

Regards.
Monir

...have just re-checked the code for the 100th time!! In my book it is on page 168, and there is no typo.

Has anyone identified the error in the code of DDPOLY() ??Maybe the error is in your expectations.
...help would be greatly appreciated...
ddpoly.for uses single precision arithmetic (a maximum of six or seven significant decimal digits).

Evaluation of your example polynomial at your test point involves subtracting nearly equal terms, so loss of significance is inevitable.

With the coefficients of your original post, here's what I got from ddpoly.for using Numerical Recipes Version 2 legacy code from the NR3 CD. The listing of ddpoly.for is on page 168 of Numerical Recipes in Fortran, Second Edition. http://www.nrbook.com/a/bookfpdf/f5-3.pdf

c(1) = -4.8865685E+06
c(2) = 1.5202154E+07
c(3) = -1.5764277E+07
c(4) = 5.4489230E+06
ddpoly( 0.959500): pd(1) = 0.5210320
ddpoly( 0.959500): pd(2) = 0.5582905
ddpoly( 0.959500): pd(3) = -159103.6093750


If the answers from ddpoly are different from yours, I'm guessing it's because you didn't put those exact values into your program where it called ddpoly.

Now, forget ddpoly for a minute and just write a program to evaluate the polynomial. I won't worry about the derivatives for now. My thinking is that, if the polynomial can't be evaluated with any precision, the derivatives are even less likely to be correct.

Here is an example with single precision arithmetic:


PROGRAM EVAL
C
INTEGER NC
PARAMETER(NC=4)
INTEGER i
REAL x,c(NC)
REAL y, ysum
DATA c/-4886568.5, 1.5202154E+7, -1.5764277E+7, 5448923/
C
do 1 i=1,NC
write(*,'(" c(",I1,") = ", 1pe23.15)')i, c(i)
1 continue
write(*,*)
C
x = 0.9595
write(*,'(" Single precision arithmetic with x = ",f10.7/)') x
write(*,'("=================================================="/)')
write(*,'(" Naive (brute-force) evaluation:"/)')
write(*,'(15x,"term", 20x,"sum")')
ysum = 0.0
do 2 i=1,nc
y = c(i)*x**(i)
write(*,'(I3,":",2(1pe23.15))') i, c(i)*x**i, ysum
ysum = ysum + y
write(*,'(27x,1pe23.15/)') ysum
2 continue
write(*,'(" Value = ",f23.15)/') ysum

write(*,'("=================================================="/)')
write(*,'(" With Horner''s nesting scheme:"/)')
write(*,'(15x,"term", 20x,"sum")')
ysum = 0.0
do 3 i = 1,nc
write(*,'(I3,":",2(1pe23.15))') i, c(nc-i+1)-ysum*(1.0-x), ysum
ysum = ysum * x + c(nc-i+1)
write(*,'(27x,1pe23.15/)') ysum
3 continue
write(*,'(" Value = ",f23.15)') ysum

END


I did the evaluation a couple of ways. For each method I printed out the value of the term that was added to the total at each step.
Output

c(1) = -4.886568500000000E+06
c(2) = 1.520215400000000E+07
c(3) = -1.576427700000000E+07
c(4) = 5.448923000000000E+06

Single precision arithmetic with x = 0.9595000

==================================================

Naive (brute-force) evaluation:

term sum
1: -4.688662500000000E+06 0.000000000000000E+00
-4.688662500000000E+06

2: 1.399571500000000E+07 -4.688662500000000E+06
9.307052000000000E+06

3: -1.392544300000000E+07 9.307052000000000E+06
-4.618391000000000E+06

4: 4.618390000000000E+06 -4.618391000000000E+06
-1.000000000000000E+00

Value = -1.000000000000000
==================================================

With Horner's nesting scheme:

term sum
1: 5.448923000000000E+06 0.000000000000000E+00
5.448923000000000E+06

2: -1.598495800000000E+07 5.448923000000000E+06
-1.053603500000000E+07

3: 1.562886300000000E+07 -1.053603500000000E+07
5.092828500000000E+06

4: -5.092828000000000E+06 5.092828500000000E+06
5.210319757461548E-01

Value = 0.521031975746155



A reminder: Single precision arithmetic is only good for about six or seven significant decimal digits. I printed out more just to see how they would be converted, but everything after the seventh is guaranteed garbage, no matter what calculations were involved.

Note that the final value from Horner's scheme in my test program is the same as I got with ddpoly.

Notice that with both methods, the last operation involves subtracting two nearly-equal numbers. That is, the magnitudes of the numbers agree to seven decimal digits, so subtraction results in total loss of significance. Both answers are equally wrong.

When I took that program and make everything double precision, here's the output:

c(1) = -4.886568500000000D+06
c(2) = 1.520215400000000D+07
c(3) = -1.576427700000000D+07
c(4) = 5.448923000000000D+06

Double precision arithmetic with x = 0.9595000

==================================================

Naive (brute-force) evaluation:

term sum
1: -4.688662475750000D+06 0.000000000000000D+00
-4.688662475750000D+06

2: 1.399571485909850D+07 -4.688662475750000D+06
9.307052383348498D+06

3: -1.392544218765611D+07 9.307052383348498D+06
-4.618389804307608D+06

4: 4.618389819052237D+06 -4.618389804307608D+06
1.474462915211916D-02

Value = 0.014744629152119
==================================================

With Horner's nesting scheme:

term sum
1: 5.448923000000000D+06 0.000000000000000D+00
5.448923000000000D+06

2: -1.598495838150000D+07 5.448923000000000D+06
-1.053603538150000D+07

3: 1.562886343295075D+07 -1.053603538150000D+07
5.092828051450750D+06

4: -5.092828036083755D+06 5.092828051450750D+06
1.536699458938529D-02

Value = 0.015366994589385


Now, the magnitude of the operands in the last operation in the Horner's scheme agree to seven significant decimal digits, so at least seven digits of precision are lost, but maybe, just maybe, some significance is left. Note that the value shown is the same Matt saw from NR3 routine (done in double precision). Note also that the two mathematically equivalent routines have different round-off error. Which do you think is the "better" answer. In other words, how many correct significant digits would you expect?

Next-to-Bottom line: One of the most common causes of loss of precision in numerical analysis is subtraction of nearly equal terms.

Bottom line: There is a difference between seven significant digits and seven correct significant digits.

Regards,

Dave

Dave,

Thanks for clearing this one up. Another programmer and I were discussing this same issue.

One of the most common causes of loss of precision in numerical analysis is subtraction of nearly equal terms.


Good post Dave.

Thanks
Matt

Dave and Matt;

I apologize for the delay in getting back to you. I had to do some basic analysis just to make sure.

1) You're absolutely correct!
The use of NR routine DDPOLY() in Single Precision (SP) was the problem.
Declaring both the input poly coeffs c(1:nc) and the returned values pd(1:nd) as double Precision (DB) produce the correct results; consistent with those by Excel and others. No surprise there, since these applications use double precision arithmetic (by default).

2) Interestingly enough, the sample quoted in the OP was selected at random from thousands of analytical data points for debugging purposes, with the hope to identifying a serious computational error in the algorithm.
As I've just found out, routine DDPOLY() in SP (as listed in NR, 1992, on page 138) works fine for almost all the data points and its results are accurate to at least five/six decimal places. The exception being a handful or so data points which require special handling (i.e.; DB) where loss of precision is (apparently!) serious.
It was a pure luck that one of those handful points was picked up for spot-check, otherwise the problem would have not been easily identified and eventually solved.

3) This is clearly a classic case illustrating that using the wrong precision can seriously affect the result of the calculation; producing not only "wrong" and "unexpected" but also "strange" results!
For the example quoted in the OP, the SP DDPOLY() returns the "strange" poly value and 1st deriv of the 3rd deg poly at x = 0.9595 as:
...PD(1) = -1. (which agrees with your Naive evaluation)
...PD(2) = -2.

4) The correct values by the DB DDPOLY() using the poly coeffs provided in the OP are:
...PD(1) = 0.01536699291318655
(compared with yours 0.014744629152119, and Excel's 0.015367.....)
...PD(2) = -0.0664996188133955
(compared with Excel's -0.06414775.....)

5) Using the exact DB poly coeffs computed internally by the program, the DB DDPOLY() gives the correct values at x = 0.9595 as:
...PD(1) = -1.7206060886383057
(compared with y = -1.72061 quoted in the OP)
...PD(2) = 176.28550790436566

PROBLEM SOLVED!!
Thanks again Dave and Matt for your help.

Kind regards.
Monir