Boycott Numerical Recipes (2007) | Hacker News Reader

Boycott Numerical Recipes (2007)

mingus.as.arizona.edu

155 pointsarthurmorgan107 comments

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.

This 100%.

20+ years ago I took my NR class as part of my physics PhD program, it was a required course for all physics and astronomy students.

I don't think anybody in the class was planning to use the actual recipes from the book as a drop-in placement for their work, but it was to explain the methods behind the software we all used (Matlab in my group, IDL and IRAF for the astronomers, etc), so users can make the best choice for which curve fitting or regression or optimisation routines to call, and avoid potential pitfalls.

To me it seems like an OS class that uses minix to explain behind-the-scenes operation system functionality. Hardly anyone in the class will go on to use Minix as their primary OS in their future career. (I've not done any formal CS study, so forgive the analogy).

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.

Oh, yes. The arrays in the C code were indexed from 1, following the FORTRAN code, with functions to convert the subscripts. Pointers to arrays were offset by one element size, so they didn't point to the array at all. They pointed to an address outside of the array.

I rewrote parts of NR in C++, but ancient, 1990s C++, with pre-STL Microsoft collection classes.

> the code, not so much

Exactly. And what's nice about that approach is that it sidesteps that ludicrous license.

As many people point out, you shouldn't roll-your-own numeric code. Unless you are solely responsible for knowing it's what you want and are cognizant of all the tricky edge cases and can make informed engineering decisions of course.

As the article says "Numerical Recipes is a good book for learning about numerical methods." I agree. If I need a quick algorithm in my language of choice, I might refer to NR in Fortran or C as a guide, but the actual implementation must be cleaned up and integrated correctly. Probably extended and specialized and combined with whatever tricks are required for the purpose. The books are (were) useful as books, but the code was never more than executable pseudo-code. I find it helpful to understand the algorithm, and once I do, I immediately begin seeing where I need to adjust it too.

And you can't copyright algorithms. So as long as I don't make a NR.py lib (which no-one needs) that contains the entire Structure, sequence and organization (SSO) of the book, the publishers can stick their crappy license where they like. I don't have to boycott the book.

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.

Open source code was a very new idea, perhaps going back no further than the first SHARE meeting 31 years earlier in 01955, but imposing legal restrictions on source code distribution was an even newer idea.

(But what about open-source numerical methods? Well, EISPACK had been published in 01974, 12 years before NR came out.)

Typing in code from a book was perfectly reasonable.

The issue today, though, is not what restrictions the book's copyright holders imposed, or purported to impose, in 01986. The issue is the restrictions they impose today.

> 01974

People from the year 100000 might be offended by this notation. But in seriousness, what is the extra 0 for? Even without it, 1974 and 10000 are perfectly distinguishable, and for computers good methods to store long dates use binary storage formats, so amount of ASCII decimal digits doesn't matter there

It is for Y10k compliance, https://longnow.org/ideas/02013/12/31/long-now-years-five-di...

To support the sibling: If we're talking history, the assumption in the early 80s was that software was shared, at least in the technical computing circles I lived in and to the extent it was portable enough. We were suspicious of people who kept code under wraps -- with good reason in some cases that dubiously reproduced data post hoc, which was an early lesson. Free Software per se was a reaction to that culture being eroded around rms earlier than I saw it, and it's basically what I've done since I started.

I'm not convinced about "black art" given, for instance, the Harwell library and later NAG, although they weren't free software. (I knew about Ian Grant's work as an undergraduate, though he didn't teach me.) Later I knew how to get code from Netlib by email before we had general trans-Atlantic network access.

I found Nick Higham has a soft spot for the undergraduate textbook we had, from rather earlier, Acton's Numerical Methods that Work.

> Source code distributions were awkward;

Tapes by mail was a thing ...

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 just like this phrase from the Usenet comments in the first link:

> Their code has the implicit assumption that all elliptic problems have homogeneous Dirichlet boundary conditions!

I think the basic problem is that the algorithms were for exposition, not production quality (they would not guard against under/overflow, edge cases, etc.)

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.

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 is not the same as choosing a random number.

What applications call for a "random float in [0,1]" with distribution defined in the space of fp representation and not the usual kind?

For example you when you want to do fuzzing of some algorithm that uses floats, a random float might come in handy.

In [0,1]?

(I consider this answer nearly disingenuous. It's obviously not what OP was talking about.)

Sure, but who said N has to be sizeof(int)? And besides that, the distance between x and successor(x) is not constant throughout the range of [0, 1).

N is sizeof(float), which is conveniently the same as sizeof(int).

> And besides that, the distance between x and successor(x) is

That is approximately balanced by using random() int / MAXINT

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!

You can make this approach very efficient using bitwise operations, and the count trailing zeros instruction common on many CPUs. See my other comment for a reference to the paper describing the technique and a link to a low level implementation by me.

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*

It's a neat trick, but can not generate all floating point values. See my other comment for that.

If you want max entropy and normally distributed values then, I think, there's no other way to go. Nice post btw

What do you mean by "normally distributed values"? Your reply also doesn't follow a normal distribution by the way. Did you mean uniform? For uniform values there is a way to go and my other reply exactly describes how.

Of course, it's uniform. I always swap the two names randomly hehe. I understand your code. The difference with the snippet I posted (which, btw, is not mine) is that the difference between successive floats is always a constant, i.e. there is no rounding. It's fast, has 23 bits of entropy and it's suited for constrained platforms/small demos [0], since it can be coded with very few asm instructions.

I'm not saying that your solution hasn't its merits, well done.

[0] The author of the code is a known member of the demoscene

Oh, I know very well who iq is :) And for a visual demo it might make perfect sense to sacrifice completeness of a solution for speed / code size. But to me, "pick a random number from {k/2^23 : 0 <= k < 2^23, k ∈ ℕ}" and "pick a random number uniformly at random from [0, 1)" are not the same.

All operations (add/sub, mul/div, sin, exp, log, etc...), with IEEE754 floats are defined using real arithmetic, and then assumed to round to the nearest representable floating point value. I don't see why random number generation, at least as a default before micro-optimizing, should be any different.

You are right, I was uncomfortable with the uneven density of the floats approaching zero, but it's just an inherent property of their representation. Have a nice day :)

The source of disagreement is: his code generates floating point approximations of uniformly-distributed real values in [0,1]; not random floats greater than 0 and smaller than 1.

This distinction has little to do with the subtlety of floats. You can e.g. generate numbers up to 100 with

- 5d20 (5 dice of 20 sides),

- 20d5 or even

- 16d6+1d4.

[Edit: the point of this exercise: these are different representations of the numbers up to 100.]

But even assuming each die with D sides has equiprobable results with prob 1/D, these choices don't have the same probability distribution. To convince yourself of this, compare 1d12 (the probability of getting 1 is 1/12) with 2d6 (the probability of getting 1 is 0).

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.

That can't be right:

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

Did you mean 2^(-k-1)?

Yes, woops, fixed.

The Python docs¹ include an implementation of Allen Downey's algorithm:

    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)

¹ https://docs.python.org/3/library/random.html#recipes

That's almost correct, but not quite (depending on your definition of correct). It fails to account for the problem that Allen Downey's paper so carefully took care of. I think it's best shown with an image (exaggerated to have a 2-bit mantissa and 2-bit exponent):

https://i.imgur.com/5TerZkm.png

Top is Python, bottom is correctly rounding to nearest float. Python's implementation does not round correctly to the nearest floating point value, it floors. This means it oversamples those floating point numbers with mantissa 0, and can not sample an inclusive range at all.

FWIW that algorithm for uniform floating-point sampling in [0,1] is actually originally due to: Walker, Alastair J. “Fast generation of uniformly distributed pseudorandom numbers with floating-point representation.” Electronics Letters 10 (1974): 533-534.

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

I think the OP knows that, but complains about the fact that randint()/INT_MAX (which my mind autocorrects to rand()/RAND_MAX) takes on at most RAND_MAX+1 different values (that _may_ be as low as 32,768), and returns each one with equal probability, so that it returns zero 1 in (RAND_MAX+1) times, on average, and won’t ever return most possible floats in [0,1]

Even if you decide that returning only RAND_MAX+1 different values, I think you should divide the [0,1] interval into RAND_MAX+1 equal sized intervals and return the floating point value at the center of one of such intervals. The given solution has both zero and one stand for a value in an interval of size 1/(2RAND_MAX), but picks them with probability 1/RAND_MAX, so it picked both of them way too often.

Aside from that there’s:

- rand(), on many systems, is a very poor random number generator.

- the API of rand()* is poor. Because there’s only one global state, you can’t run multiple independent generators, and it doesn’t help you much as soon as you want to generate anything other than a discrete distribution with RAND_MAX+1 items (a number that even is implementation dependent, so that portable code needs special care even if it may not be needed on most platforms). For example, because RAND_MAX may be even, you can’t even simulate a fair coin toss with a single call of rand().

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

Steelmanning their argument, an ideal output would consist of every possible float in [0,1] occurring with an appropriate probability -- some of those vanishingly unlikely. Kind of like how randint() should be capable of returning every possible integer rather than, say, every multiple of 64.

It's fine to counter that most of the time the difference wouldn't matter or that it might be problematic to compute for some reason or another, but I don't think it's reasonable to critique the idea on the grounds of the result not being uniform when they didn't actually ask for each float to have equal probability.

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

Randomness is not a property of numbers but of distributions.

Actually I believe it is better to say randomness is a property of a sequence of numbers. So a "random sequence" has members chosen independently of each other, but all from the same probability distribution. And random number generators create random sequences which follow a given distribution.

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.

HTTP response headers say:

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

I always had the impression that NR was about getting a broad overview of the field. I imagine each chapter covers at least a whole field of numerical research. If you’re working on applications of numerics, it is indeed important to know what’s available out there. The NR code itself I’ve never used, really.

It's more of an intro level book, aimed at scientists rather than numerics geeks or or people into theoretical math. So it does a good job of presenting things simply and in ways that don't require navigating too many fine points. That approach mostly works, though it leaves openings for things to go wrong. The topics (root finding, integration, etc.) are similar to what you'd see in an undergrad numerics course. It specifically doesn't have anything about PDE's, since it defines that as the boundary between introductory and advanced topics.

FWIW, there is a whole chapter (17 in the first edition, 20 in the third/current edition, don’t have a second edition in front of me but would be surprised if it didn’t have one) on Partial Differential Equations. It’s about 50+ pages in the first edition and 70+ in the third.

But I take your point, the chapter does admit to being intended as “the briefest possible useful introduction” and notes that the topic would require another volume entirely dedicated to it.

Ah thanks, I must have mis-remembered. I think it did say something about PDE's being the dividing line though.

I don't get what the big deal is. I remember this book. It was useful back in the day, but really, has anyone really ever been pinched as an NR license violator? I expect some orgs paid for licenses, it's not an exorbitant amount. It probably amounted to little more than pocket change for the professor/author over the years. Is the guy lawsuit-crazy, I didn't see evidence of that.

If some grad student, postdoc or independent uses it and doesn't quite follow the license terms, I expect nothing will happen. Though they would probably do well to update their code base, eventually, for technical as well as license reasons if they "scale up" beyond using the code for a paper or two. But the tone of the OP's call to "Boycott" is a bit over the top, RMS-style, and not in a good way.

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.

But its algorithms are often incorrect, inefficient and numerically unstable. Just use a library written by an expert in whatever your area of interest is.

Isn't this book essentially a library written by experts[1] in numerical processing?

You can't just use a library by experts. You have to know or learn how to evaluate said library for correctness, efficiency, and stability (numerical and operational). If you know how to do that, you may be able to use a flawed library as a tool to get the shape of your processing quickly and fill in the rest of the details yourself.

[1] They've published a major book through three editions and the book cover calls them leading scientists.

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.

Oh, the bible. It is a great and thorough book (maybe not ideal for beginners). What it covers is quite different from NR, though.

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.

Note that the article is not about denigrating the algorithms, and is instead entirely about the license restrictions.

I know but a lot of HN comments about the algorithm quality. It seems like a moot point because even if the algos were obsolete, the license issues remain.

I don't think the algos are obsolete though.

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.

Most people using NR aren't building large programs.

Think about a student who learned from NR, then in grad school needs a method, and in the spirit of "I'm smart, I can do this myself", copies the method out of the book.

In larger projects, someone is more likely to be aware of alternatives, and the issues related to using source covered under the NR copyright.

... hesitantly raises his hand

Back in ... 1991? I did use the C version of the Simplex Method almost straight from the book. I needed a concise linear programming optimization module that would be called in the middle of a not-really-concise production programming system handling (automotive) engines and gearboxes production for a couple of FIAT plants.

Basically a short-period scheduler to decide what to produce (slave to an immense system running on IBM mainframes; COBOL, again).

So I quickly put together the C version (with minimal modifications) and it chugged along happily on a VAX for at least a decade.

Thanks for the counter-example. :)

Did your company also license the code?

To explain my frame of mine, I was thinking how it was used in 2007 or later.

In the early 1990s, there wouldn't have been much off-the-shelf/free&open-source options. The GNU Scientific Library didn't start until 1996, for example (according to Wikipedia).

No. I was aware of the issue at the time, and this is part of the reason for me to alter the source a bit.

In my defense, the people I was working for at the time did not even understood the topic when I talked about that.

I think also that 30+ years ago software copyright was farther in the background.

The "Don't Copy that Floppy" propaganda, for example, didn't start until 1992.

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

Seems nr.com is no more. Wonder why they decided (or should I say how much they're offered) to switch away from such a premium domain.

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

The Eigen library (https://eigen.tuxfamily.org) has great attention to numerical robustness and accuracy in it's implementation. The documentation is good, but for the numerical discussions, you'd have to read the source code.

The main problem is each subsection of a chapter in NR is really a subfield in and of itself, and to adequately describe a decent modern implementation is at least one if not multiple books.

First make sure you're familiar with background issues in numerical analysis/computation such that you're generally familiar with IEEE floating point and its foibles, numerical stability, basic error analysis etc.

Figure out which algorithm you're interested, find a couple of open source implementations (academics generally generate a vast array of these, though many may become abandonware), look at any papers they are associated with/mention and read those. Now you're ready to start writing.

Once you have a basic implementation, most fields have a set of standard problems they benchmark on. Find those and see how you compare on numerical accuracy and performance against the libraries you looked at before.

> First make sure you're familiar with background issues in numerical analysis/computation such that you're generally familiar with IEEE floating point and its foibles, numerical stability, basic error analysis etc.

So, what's the book for learning that?

You might find the book Matters Computational: Ideas, Algorithms, Source Code by Jorg Arndt (the book is freely available) and the accompanying FXT library useful - https://www.jjj.de/fxt/fxtpage.html

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)

Funny you mention that. I was ready to rant the opposite (ie, that recipes should not include exact measurements).

I really, strongly believe that our obsession with measurement and repeatability is a net negative that has taught a whole generation to treat cooking like a standardised test rather than a creative form of human expression.

You look at the older generations from cultures that value cooking and they simply don't give a fuck about what the recipe says to do.

> that recipes should not include exact measurements

Wouldn't that make cooking much less accessible?

It’s tricky - people have different tastes, and ingredients may change from region to region.

Furthermore, making people measure things with a scale rather than a cup adds significant effort and time in the scheme of things.

It takes less time to get the scale out of the drawer than figuring out which bloody cup I have to use.

>figuring out which bloody cup I have to use

Unless you're baking, you could literally use any cup you have in your house (yes, even a drinking glass!). Unless you're baking, it really doesn't matter, and that's the beauty of it.

I don't think so. You'll fuck up a few times catastrophically in the beginning, but very quickly get a feel for what works and what doesn't (for example, you can go crazy on the garlic without ruining a dish, but you need to be a bit more careful with salt).

If anything, the belief that you must boil the potatoes for exactly 15 minutes is more likely to make you hate cooking than teach you how to enjoy it/excel at it.

Most people never have been and never will be interested in cooking as a creative outlet, and trying to convince them to do so is a highly arbitrary goal only fueled by your own subjective interests.

The vast majority of people just want to make some shit to eat everyday. And precise measurements are great for that.

Precise measurements make the problem of cooking something tasty feel much harder than it actually is.

Nothing hard in putting stuff in a bowl over a scale until the scale says 750g. On the contrary, it’s easier than trying to guess what the author had in mind when they wrote “a pinch” or “a handful”.

>On the contrary, it’s easier than trying to guess what the author had in mind when they wrote “a pinch” or “a handful”.

Again, I think you're missing the fact that whatever the author had in mind isn't "right". Do whatever feels right to you and as long as your senses are calibrated (there'll be a short period of adjustment while you learn how ingredients work) it'll come out better than if you'd measured and you'll enjoy it more too.

> Do whatever feels right to you and as long as your senses are calibrated (there'll be a short period of adjustment while you learn how ingredients work) it'll come out better than if you'd measured and you'll enjoy it more too.

You’re missing that this is a lot of work, and far from everybody is even remotely interested in doing that. Most people would rather be doing useful or entertaining things, not cooking the same thing a thousand times to discover how much salt “to taste” means. Do not assume that everybody’s hobby ought to be cooking.

>You’re missing that this is a lot of work, and far from everybody is even remotely interested in doing that

It's not. It's a couple of dozen disappointments in your teens and then you know how to cook everything you want for the rest of your life.

>Most people would rather be doing useful

I actually crunched the numbers on this once, and worked out that having the skills to cook and adapt to whatever ingredients are cheap/in season is worth somewhere in the order of $250,000 over the lifetime of an adult.

>to discover how much salt “to taste” means

"to taste" is however damn much you want it to be. That's the point I'm trying to make. It isn't that cooking needs to be everybody's hobby; it's that:

- It doesn't need to be precise

- The obsession people have with precision makes things harder and more stressful; not easier.

> It's not. It's a couple of dozen disappointments in your teens and then you know how to cook everything you want for the rest of your life.

But not everybody wants to do that. Why should everybody be forced to get good at cooking? It’s like those “Oh, everybody should learn to play a musical instrument” people. No we should not.

> I actually crunched the numbers on this once

You’re still arguing that everybody should be forced to learn to cook.

> "to taste" is however damn much you want it to be.

But if I don’t know what a typical, reasonable value is, I’m forced to guess, probably ruining the food, and now I have no food today. Repeat for multiple days. Why do you want to force people to do that? It actually feels like gatekeeping, with using weird in-group terms and assuming knowledge which can’t even be taught and must be learned personally over a long period of time.

> Again, I think you're missing the fact that whatever the author had in mind isn't "right". Do whatever feels right to you and as long as your senses are calibrated

That requires experience. I have some and have no problem cooking, but a lot of people don’t and have.

For many people precise instructions imply that you need to be precise when executing them, i.e. if you deviate you ruined it. Most recipes however are not like that and the tolerances are actually huge.

I would say it really depends. Some times exact quantities are necessary, some times they aren't.

I hate that almost no recipe specifies salt in grams. How much salt fits in half a teaspoon is pretty arbitrary since it strongly depends on the type of salt. For many dishes it doesn't matter, since you just salt according to taste. But if you are making meatloaf, getting the amount of salt right is essential, and you can't really taste it unless you like tasting raw meat, so an amount in gram would be appreciated.

I was also confused, but about an apparent gripe against cooking in sequential order. To be fair, tbsp is universally acceptable, but the apparent apparent ignorance of metric measures among US-centric authors is ludicrous.

To continue this off-topic rant, recipes have been receiving a lot of attention lately, and the medium is well overdue an overhaul to cut down on the fluff. A comprehensive wiki, with cooking markup (eg. Cooklang) would be amazing and would fill a long-overdue gap that hasn’t been filled since Larousse Gastronomic was released. The idea that recipes are ‘owned’ by someone who happened to publish it first is crazy, but could nonetheless be addressed with an attribution license (eg. the new MIT-0?)

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.

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.

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.

> Don't roll your own numerics.

This is empty advice without some context. Many of us work precisely in rolling new numerics.

Like "don't roll your own crypto", the implicit caveat is clear.

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

Correction: That date is when the Netlib email service for the existing collection started.

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.

Reply on news.ycombinator.com