double exp(double) unexpectedly returns NaN (no float overflow)

Question

I've been working on implementing black-body radiation according to Planck's law with the following:

double BlackBody(double T, double wavelength) {
  wavelength /= 1e9; // pre-scale wavelength to meters

  static const double h = 6.62606957e-34; // Planck constant
  static const double c = 299792458.0; // speed of light in vacuum
  static const double k = 1.3806488e-23; // Boltzmann constant

  double exparg = h*c / (k*wavelength*T);
  double exppart = std::exp(exparg) - 1.0;
  double constpart = (2.0*h*c*c);
  double powpart = pow(wavelength, -5.0);
  double v = constpart * powpart / exppart;
  return v;
}

I have a float[max-min+1] array, where static const int max=780, static const int min = 380. I simply iterate over the array, and put in what the BlackBody gives for the wavelength (wavelength = array-index + min). The IntensitySpectrum::BlackBody performs this iteration, while both min and max are static member vars, and the array is also inside IntensitySpectrum.

IntensitySpectrum spectrum;

Vec3 rgb = spectrum.ToRGB();
rgb /= std::max(rgb.x, std::max(rgb.y, rgb.z));
for (int xc = 0; xc < grapher.GetWidth(); xc++) {
  if (xc % 10 == 0) {
    spectrum.BlackBody(200.f + xc * 200.f);
    spectrum.Scale(1.0f / 1e+14f);
    rgb = spectrum.ToRGB();
    rgb /= std::max(rgb.x, std::max(rgb.y, rgb.z));
  }
  for (int yc = 20; yc < 40; yc++) {
    grapher(xc, yc) = grapher.FloatToUint(rgb.x, rgb.y, rgb.z);
  }
}

The problem is that, the line spectrum.BlackBody() sets the 0th element of the array to NaN, and only the 0th. Also it does not happen for the very first iteration, but all the following ones where xc>=10.

The text from the VS debugger: spectrum = {intensity=0x009bec50 {-1.#IND0000, 520718784., 537559104., 554832896., 572547904., 590712128., 609333504., ...} }

I tracked the error down, and exppart in the ::BlackBody() function becomes NaN, basically exp() returns NaN, even though it's argument is near 2.0, so definetely not overflow. But only for array index 0. It magically starts working for the rest 400 indices.

I know memory overruns might cause things like that. That's why I double checked my memory handling. I'm linking Vec3 from another self-made library, which is much bigger, and might contain errors, but what I use from Vec3 has nothing to do with memory.

After many hours I'm completely clueless. What else can cause this? Is the optimizer or WINAPI fooling me...? (Uhm, yes, the program creates a window, with WINAPI, and uses a nearly empty WndProc that calls my code on WM_PAINT.)

Thanks for you help in advance.

Sorry for making it unclear. This is the layout:

// member
class IntensitySpectrum {
public:
    void BlackBody(float temperature) {
        // ...
        this->intensity[i] = ::BlackBody(temperature, wavelength(i));
        // ...
    }
private:
    static const int min = 380;
    static const int max = 780;
    float intensity[max-min+1];
}

// global
double BlackBody(double T, double wavelength);

Answer 1

If you happen to be using MSVC 2013, one possible explanation is that you have some code somewhere that is trying to convert a float infinity to int. A bug in MSVC 2013 causes an unbalanced push on the x87 FPU stack when this happens. Trigger that bug 8 times and your FPU stack is totally full, and any subsequent attempt to push a value (such as calling 'exp()') will result in an 'invalid operation' and return an indefinite (like 1.#IND). Note that even if you are compiling with SSE2 floating point instructions, this bug can still bite because the calling convention dictates that floating point return values are returned on the top of the FPU stack.

To check if this is your issue, have a look at your FPU registers just prior to the bad call to 'exp()'. If your TAGS register is all zero, then your FPU stack is full.

• http://connect.microsoft.com/VisualStudio/feedback/details/806362/vc12-pollutes-the-floating-point-stack-when-casting-infinity-nan-to-unsigned-long

MS claims this will be fixed in update 2 for MSVC 2013.

Answer 2

The following function call only has 1 parameter:

spectrum.BlackBody(200.f + xc * 200.f);

So it cannot be calling the function you defined as

double BlackBody(double T, double wavelength)

If you look at the ::BlackBody implementation, I'm betting you have a divide by 0 error somewhere.