First, the new edition of NR is really great. The authors are especially to be congratulated for making the absolutely right decision to not build their own elaborate vector/matrix classes but instead to stay focused on the algorithms, so that they can be incorporated into any class architecture. This is a great step forward toward creating a library of generic numerical algorithms, in the spirit of the Standard Template Library.

I would very much like to see NR move toward a more STL-like approach, and I've been making some changes to support this. One change is almost trivial but enhances the re-usability enormously. The algorithms can be wrapped in templates, with template arguments for classes such as DP and Vec_I_DP. Here's an example of the signature for rk4:

template <class DP, class Vec_I_DP, class Vec_O_DP>
void rk4(Vec_I_DP &y, Vec_I_DP &dydx, const DP x, const DP h,
Vec_O_DP &yout, void derivs(const DP, Vec_I_DP &, Vec_O_DP &));

This immediately allows rk4 to be used with any library routine that supports the basic operations of DP, Vec_I_DP, and Vec_O_DP. An example of this is double and the standard library's vector<double> class. So here's the code for a simple falling body, using a generic NR::rk4 algorithm (of course, normally you wouldn't directly call rk4 -- this is just for illustrative purposes):

#include "nr/rk4.h"

#include <iostream>
#include <vector>

void derivs(const double x, vector<double>& y, vector<double>& dydx)
{
dydx[0] = y[1];
dydx[1] = -9.8;
}

int main()
{
vector<double> y(2);
vector<double> dydx(2);
double x = 0.0;
double h = 1.0;
vector<double> yout(2);
derivs(x, y, dydx);
NR::rk4(y, dydx, x, h, yout, derivs);

cout << yout[0] << "\n";
cout << yout[1] << "\n";
}

There are two very nice things about the generic version of the NR code: It can be used exactly as the current NR in C++ code (the signature is the same), and it can automatically adapt to a different class library, so long as it implements the required functions. This makes a lot of pages 29-30 in the book unnecessary, as shown above -- I was able to directly use the standard vector<double> class without making any changes to the generic NR algorithm. As a trivial example, I can shift to vector<float> in my code and not have to make any changes. Of course, if a library doesn't provide exactly the interfaces required, you can always write an adapter class that does so.

Another change to make which would not affect the ability to use the generic versions of NR in exactly the same places as now is to also provide a template argument for function arguments, which would allow any function-like thing (including 'function objects' or 'functors', objects with an overload parentheses operator) to be passed, and not just function pointers. The resulting interface to rk4 would be better expressed as

template <class DP, class Vec_I_DP, class Vec_O_DP, class Function>
void rk4(Vec_I_DP &y, Vec_I_DP &dydx, const DP x, const DP h,
Vec_O_DP &yout, Function derivs)

I really NEED this feature -- my normal practice is to have a derivs 'function' which maintains state information.

That said, there is one other change which might be good to make, but which does not maintain compatibility with the current interfaces, and that is to pass iterators as arguments, rather than the containers themselves. One of the problems with NR in C++ is that they can't be used with the built-in types, such as double[]. It would be nice to be able to use them in these situations, as well. By passing iterators, rather than containers, the algorithms would work for any container, including the built-in types, that implements iterators correctly. This would mean changing the above call to NR for the built-in vector class to something more like:

NR::rk4(y.begin(), y.end(), dydx.begin(), x, h, yout.begin(), derivs);

or, if you were using double y[3], double dydx[3], double yout[3], it would look like this:

NR::rk4(y, y+2, dydx, x, h, yout, derivs)

I hope this is of some help.

Thanks for the kind words, and for a really useful posting!

Once we realized that we weren't going to be the ones to write the be-all-and-end-all of class libraries, it became clear that we should provide every possible syntactic "hook" for smart folks to tap into the underlying (mostly-) generic algorithms at various levels. Your posting provides some great ideas for hooking to the DP, VEC_I_DP, etc., level. The book mostly discusses how to hook to the NRVec, NRMat level. It's interesting to compare and contrast the two approaches.

Cheers,
Bill P.

This is my first post, and I hope this topic has not already been covered .. I would like some views on a matter of style .. as NR moves to genericity, should it use [] or () for member access?

There are evidently some people with very set views on this matter:
http://www.parashift.com/c++-faq-lite/operator-overloading.html#faq-13.9

Personally, I use NR in C most of the time, and am still debating whether moving to NR in C++ makes sense :confused:

I personally don't have strong feelings on [] vs. (). It is hard to be religiously dogmatic in this area because you can make C++ do some remarkable things by combining overloading and automatic type casting, without even getting into the incredible things you can do with templates. You can look, for example, at section 23.4.6 of Stroustrup's 3rd edition, where he provides an implementation of Matrix that provides both notations, and then some.

But for simple, single-indexed things that match well to the concept of random-access iterators, I always use [] and a corresponding * operator, so that x[i] is equivalent to *(x+i), where x is a random-access iterator. This covers most of the cases in NR and all the ones I had worried about up until now, and lets you use arrays of the built-in types as well as all the STL types, etc., in any algorithms.