It's been a very long time since I've used NR, but it was pretty terrible code that was a poor competitor for NAG and IMSL.

The discussion of algorithms, being informal and well-illustrated, was the high point of the books. The code, not so much.

I had the Fortran book and so, when the C book came out, I bought that as well. The C code seemed to have been auto-generated by the Fortran code, and it was really quite terrible.

Functions in the fortran code had arguments for what I'd call "work arrays". There was some risk in that, because if the documentation told you to give space for 2N+1 floating point values and you gave 2N, you might get into problems. I don't think that was a real problem in practice, because people were accustomed to this general way of working and so they were careful.

The C code took a different approach. The functions allocated and deallocated space as needed. And, since lots of functions called other functions for very simple things, there was a massive amount of memory juggling going on. This made the functions be wildly slow. So, the point of the boycott was not terribly relevant, because I think pretty much everyone trying to use that C code just rewrote it, anyway.

It's worth putting this book in historical context. It was first published in 1986. Numerical programming like this was very much a black art and the book was an well researched, documented, usable set of algorithms with (mostly) working code. Open source code was a very new idea (FSF was only established a year before, and while there are earlier precedents like BSD source it was all very new.) Source code distributions were awkward; absolutely nothing like github of course, just a few FTP servers which maybe you could get code from if you knew they existed and had Internet access. Typing in code from a book seemed perfectly reasonable.

There was a numerical analyst at JPL back in the day that absolutely savaged NR. Remember, the kind of engineering they do at JPL has to be correct, not just plausible. Also, noncommercial alternatives to things like Matlab (octave, the python interfaces to Netlib libraries, etc) weren't as mainstream then, and science and engineering research had a bigger "get a book and roll your own" ethic. https://web.archive.org/web/19990202153232/http://math.jpl.n...

I've seen it breed horrors - "senior" people instructing interns to implement its matrix inversion routine line-by-line in Matlab. Today it would be Python.

Oh, and here's a list of alternatives, from the heyday of NR: https://web.archive.org/web/19990202141044/http://math.jpl.n...

I remember being disappointed by this book.

I wanted to read about algorithms for choosing a random float in [0,1], and IIRC its only suggestion was randint()/INT_MAX. This is wrong in part because the smallest nonzero value you'll ever choose is 1/INT_MAX, which is much larger than the smallest nonzero float.

Is there a better book out there? I was hoping for something fabulous like Hacker's Delight but for floating point.

There is no date on this but it is pretty old, since I remember thekerfluffle from probably 10+ years ago. I don't know if the NR license was ever adjusted but the software was nothing special anyway. I did learn a lot from the books, from a semi-neophyte perspective. I think serious numerics experts pooh pooh them though.

  Last-Modified: Thu, 10 Aug 2006 17:08:26 GMT

The GSL library is indeed great! It provides a free implementarion of all methods described in NR, and many more. The source code of GSL is somewhat readable, but I find NR still very helpful for the textual descriptions of the methods.

Also, the venerable ccmath library by the late Daniel Atkinson is an extraordinary jewel of readability, simplicity, portability and performance among numerical codes.

I've used Numerical Recipes in C heavily since 1996, but I've never used the code (which isn't great). Instead, I've used it as a reference to understand approaches and algorithms. Any code I've made using that book has been my own.

If you know how to code and want a leg up on understanding an algorithm, NR is a very convenient structural cheat sheet.

Issues with the license aside, there are substantial problems with the core algorithmic material of the book. The provided codes and approaches are often just plain incorrect.

Quoting from a review: "The authors of Numerical Recipes were not specialists in numerical analysis or mathematical software prior to publication of this book and its software, and this deficiency shows WHENEVER WE TAKE A CLOSE LOOK AT A TOPIC in the book [editor's emphasis]. The authors have attempted to cover a very extensive range of topics. They have basically found `some' way to approach each topic rather than finding one of the best contemporary ways. In some cases they were apparently not aware of standard theory and algorithms, and consequently devised approaches of their own. The MEDFIT code of section 14.6 is a particularly unfortunate example of this latter situation.

"One should independently check the validity of any information or codes obtained from `Numerical Recipes'...."

Source: https://www.uwyo.edu/buerkle/misc/wnotnr.html

If you're looking for a better book, I can recommend:

Matrix computations, by Golub and van Loan.

IANA numerical expert, but I have written a lot of MATLAB in the past decade.

There are NR algorithms not present in MATLAB that are very useful. When trying to find the nearest index in one very large monotonically increasing vector in another you could use interp1, but what is much faster is to do a binary search. There is no binary search function in MATLAB. It's lightning fast to use the NR-C algorithm for a binary search in MATLAB. If you want to find many nearest indexes, then interp1 likely faster.

Thought Numerical Recipes was used mostly in courses for students to learn about numerical methods rather using those pre-made routines for building large programs.

Boycott C.

In C++ or even better a language like D or Ada or whatever with built-in slices you can write assertions that explicitly encode assumptions and "what I meant to write, modulo bugs"-isms much more easily than in C. How long is that array? How big can that float be before the model is unstable? Maybe express the units in the type system etc.

The current license can be found here: http://numerical.recipes/aboutNR3license.html

> These conditions are so restrictive that of course everybody breaks them, often without knowing it. Thus, the NR license makes itself a farce.

It doesn't necessarily make the license a farce. An Oracle-like business model can be established by auditing the users of the NR library and assessing fees for license violations, including the inadvertent ones.

I used NR a bit more in grad school (late 80s/early 90s) than undergrad. I had (have) the Fortran version. I used some of the ODE integration bits in some code I worked on in grad school, though we put that project away after we switched to a different DFT MD code. FWIW, that code (with NR) still compiles/runs. Tried it a few months ago.

As an undergrad, I had a research project which used Markov methods, and I had to generate some fairly large (at the time) matrices. I wrote that code in Fortran, on a PC/AT system (1986, so ...). It also ran on an IBM 3090 VF 180 mainframe at Stony Brook.

What I recall from that time, was the professor giving me a text on Markov methods, and me learning what I needed to construct the transition probabilities (matrix elements). I remember writing extensive tests within the code, called sanity checks, which I reported on. Fortran, no assert statement at the time, and its nice to see a note "COMPUTATION IS INSANE DO NOT TRUST" error message when you ran into problems.

Honestly, that process taught me a number of things.

1st, the code you need may or may not have been written by someone else, and you really need to make sure you have sanity checks built in to establish the portions you can and cannot trust. Later when I used the NR code, I added sanity checks into their code, to make sure we weren't getting into strange spaces.

2nd, it is (often very) important to exploit mathematical features of your problems (symmetries, relation of derivatives to exact closed form equations, etc.) to simplify your calculation. Plug and crank is usually the wrong approach, and you need a solid understanding of your problem domain to have a fighting chance at not producing a garbage simulation.

I've followed these coding practices for decades at this point, through many languages and systems. NR was a great way to start with these things, but I picked up Hamming's, Hildebrant's, and many other numerical analysis books along the way, taught myself, and generated code that was hopefully sane.

I don't care so much about the license of NR, as I am not distributing code containing it. My thesis advisor warned me about that in early 90s. So that was removed from code I shared.

I guess I am a bit sad that I offered said book to offspring as she starts a PhD in MechE later this year, and I don't think she wants it. She did take my numerical analysis books though, so I'm happy about that.

I never thought that you are supposed to copy the code from a book blindly. If you do that, there must be some library that you can use that work better. So treating the code as pseudo code, I thought numerical recipes were excellent, better than text books in many aspects.

What are some good, current sources for learning about and implementing numerical algorithms like those in Numerical Recipes?

I get that for most just using Numpy, GSL or similar library would be sufficient, but what if you really want to implement them yourself.

I used NR for a gauss elimination in FORTRAN back in the 90s. Next time I needed NR, I just used GSL. Nowadays I mostly stick to numpy...

Huh glad it's not about cooking because I was ready to be ranting (ie recipes using 1 tablespoon instead of X grams)

I really can't imagine anyone finding NR now, except as a dusty relic on a lab office shelf

Don't tell me what to boycott unless it's the worst. Tell me what's good.

As somebody who came to scientific computing late in life, after spending years as a physicist then an engineer, I will step up to defend NR, if only for fun (and because I like the book for what it is).

At least one commenter here dismissed NR, saying that it is equivalent to an undergraduate numerical analysis book. Well, duh! That's one of its strengths! It is a good book for an audience like me -- somebody with scientific background but with no formal numerical analysis training who needs to write some numerical code.

The target reader is an engineer who needs to implement something specific. That's a very common scenario -- you are given the job of implementing some calculation as part of a larger project. You're not an expert in the particular mathematical discipline. You don't want to fool around trying to find an algo in some book since you don't know which of the zillions of numerical analysis books to look at. Also, most books will just give you a general algo, not ready-to-run code. Everybody knows how hard it can be to translate a wordy description of an algo into real code. If you had the budget you'd just hire somebody to do it for you. NR solves your problem -- you use the code and your problem is solved so you can move on to the next step of your project.

A different commenter here suggested using Golub and Van Loan instead. That's laughable. First off, as others pointed out, their "Matrix Compuations" book addresses numerical linear algebra exclusively. NR is a general-purpose numerical analysis book, so covers a much larger topic set. But secondly, Matrix Computations is difficult reading for non-mathematicians. Even if I needed to implement some sort of matrix decomposition, I would find it a pain to use that book as a source.

(Yes, I am sure some testosterone-soaked HN commenter will now say that they read "Matrix Computations" in 3rd grade as a way to prove how smart they are. Allow me to roll my eyes now.)

The real complaint made in the linked post involves the licensing terms of the code in NR. I agree the license stinks by today's standards. But the book was originally written back in the 1980s. My copy is from the early 1990s. Back then the GPL was basically unknown outside of a group of happy hackers at MIT. Software licensing certainly not a topic occupying mindspace of ordinary scientists or engineers. There was some DOS-based "freeware" circulating on 3.5" disks, but the whole "free (libre) software" movement was still in the future -- it needed the internet to become widespread in order to take off. Finally, I can imagine the NR authors wanted some sort of copyright or license which prevented wholesale copying of their work by folks who would try to profit off it. It's reasonable for them to have wanted to get some monetary reward for their hard work. Their license is probably an artifact of its time.

In my travels I have seen situations where code used in commercial products incorporates this or that function presented in NR. For example, I saw the FFT algo used in a test widget. (I authored none of the copied code -- I simply observe it in projects done by others.) What's funny is that NR's FFT algo is not particularly great. However, the guys who needed an FFT probably didn't know that. They probably would not have cared anyway -- they just wanted a working solution and that's what they got from NR. I am sure that little violations of the NR license happen all the time. It may not be strictly legal, but I also see people making Xerox copies of copyrighted books and journal articles all the time too.

Finally, regarding alternatives like GSL, ccmath, or whatever is available on netlib, those projects post-date NR. In any event, the internet is what made all those things widely available. I bought my copy of NR in a bookstore before that happened.

Heh, I haven't seen or read NR but added it to my wishlist after watching Terry Davis mention it [at 8:15] https://www.youtube.com/watch?v=25hLohVZdME Seemed to work well enough for him.

I don't know about a better book, but I do know how to generate a random float in [0, 1] very efficiently and accurately. A 2007 paper describing the technique is Generating Pseudo-random Floating-Point Value by Allen B. Downey. I have implemented it here in Rust: https://gist.github.com/orlp/121bb95361cc74b88aff79a7335a1a4.... This is for generating in [0, 1), but it is trivially changed to allow 1.0 by changing

    while exponent < -127 || exponent > -1 {

into

    while exponent < -127 {

It uses a single call to the RNG 255/256 = 99.61% of the time, if it has to use an extra RNG call the chance it then needs another after that is 2^-32, ad infinitum, which is negligible. More calls to the RNG could be needed if the exponent it generated was out of bounds which is also incredibly, incredibly unlikely (~2^-33 for [0, 1), ~2^-127 for [0, 1]). Thus the expected number of RNG calls is essentially just 255/256 + 2*1/256 ~= 1.00391.

A friend of mine also made a blog post that includes another implementation and discusses it: https://lokathor.github.io/prng/.

The difference between this technique and virtually all others out there on the internet is that this one is actually unbiased and can generate every single possible floating point value in [0, 1). Unbiased here means that the procedure simulates picking a real number uniformly at random from [0, 1] (with infinite precision), and then rounding it to the closest floating point value. We don't have to use infinite RNG bits to generate our infinite precision real number because we only have to generate enough bits such that it is unambiguous which floating point number we round to. The actual number of RNG bits needed is constant in expectation (see above), but unbounded, because our real number could be arbitrarily close to the boundary between two floating point values.

    from random import Random
    from math import ldexp

    class FullRandom(Random):

        def random(self):
            mantissa = 0x10_0000_0000_0000 | self.getrandbits(52)
            exponent = -53
            x = 0
            while not x:
                x = self.getrandbits(32)
                exponent += x.bit_length() - 32
            return ldexp(mantissa, exponent)

  // Returns k with probability 2^(k-1).

If you want to generate a uniform random float in [0,1), half your results have to be in the [0.5,1) range. That space is half as dense as the [0, 0.5) range. (And this applies again for [0.25,0.5) vs. [0, 0.25), etc, and so on. And this is also disregarding the whole space of denormalized floats.) So regardless of your method, you're going to have to ignore nearly all the floats in the range! (Or generate many more than 32 bits and use them non-uniformly to make a uniform output - see orlp's comment.)

There are other problems with what you say their suggestion was, e.g. INT_MAX should be RAND_MAX, and you need to take some care with rounding/conversion, and inclusive vs. exclusive upper bound requires careful attention. But with any simple approach you just can't use most of the small float space, and so the division approach with the details done correctly is pretty common!

Edit since rate-limited: orlp, I think we're only in disagreement about what qualifies as "very efficient". It's an efficient implementation of the more complex problem, but it's also much slower than the division. I wouldn't want it to be the default float generator for a general-pupose language - even most scientific modeling or statistical analysis code won't benefit from the increased output space. OP should think very hard about how they intend to use values in that range before insisting they really need it - the division approach really is sufficient and NR is right to recommend it by default!

Because of how floating-point numbers are distributed, you should expect a random number drawn uniformly from [0,1] to basically never hit the smallest floating-point number: hitting something in the range [0, 2^-n] should happen with probability 2^-n.

I agree there is much more one could want out of a uniform random float generator than just k/N where k is a integer drawn uniformly from {0, 1, …, N}, but hitting arbitrarily small floating point numbers is not on this list.

If you are intending to sample numbers in the range [0,1] such that each float occurs with equal probability then that’s fine, but it’s certainly not uniformly random on [0,1].

Floating point numbers are represented in a form mantissa*2^exponent and they are distributed non-uniformly over their dynamic range. The density of normalized 64-bit floats around zero is about 10^-320, while around one is about 10^-16. Therefore, do not expect to get the smallest positive float in sampling from a uniform distribution over [0,1]. There is a good tutorial about floating point numbers and their pitfalls written for computer scientists [0] where Figure 1 on page 7 (pdf page count 3) shows how non-uniformly floating point numbers are distributed over the number line.

[0] https://pages.cs.wisc.edu/~david/courses/cs552/S12/handouts/...

If you want a random number in [0,1] with a uniform distribution from N random bits, I don't see how you can expect an epsilon of less than 1/(2**N-1)

Choosing a random float doesn't mean very much without a distribution.

You can apply the core ideas of Hacker's Delight to floats too. Your example problem, taken from [0]:

  float frand( int *seed )
  {
      union
      {
          float fres;
          unsigned int ires;
      };
  
      seed[0] *= 16807;
      ires = ((((unsigned int)seed[0])>>9 ) | 0x3f800000);
      return fres - 1.0f;
  }

You basically assemble a float building a random mantissa and assigning the right exponent and sign bits.

It's not portable, but it's a very neat trick!

[0] https://iquilezles.org/www/articles/sfrand/sfrand.htm*

Completely agree. You shouldn’t go to NR for the code. You should go to it because, even decades after it was first published, it’s still the only reasonably comprehensive text on scientific computing that doesn’t assume you’re in a graduate program in applied mathematics. And I say that as someone who actually took those graduate level applied mathematics classes. I can pick up NR and understand the basics of some algorithm within fifteen minutes, without hurting my brain.

If there’s a better book that is accessible to the masses I’d like to know about it. But I haven’t found it yet.

> Finally, regarding alternatives like GSL, ccmath, or whatever is available on netlib, those projects post-date NR.

No, Netlib started in April 1985[1], though as a physicist or engineer you likely had access to NAG or something similar already. It must have made an impression on me as I still remember F02ABF... (I hasten to add I hit scientific computing early in life.)

1. https://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.73...

Don't roll your own numerics.

There are normally many subtle gotchas that affect numerical stability of a method, as well as great attention to detail required to handle degenerate edge cases. And that's before we get to subtle bugs that only trigger occasionally or are otherwise hard to detect.

It's not like NAG/ISML/HSL etc don't predate NR, in some cases by several decades. Libraries were always available, though back in the day you may have needed to pay for them.