gamma_isf nan about owl HOT 3 OPEN

I pulled the relevant C code out of owl, compiled it independently, and ran the following:

#include<math.h>
#include "owl_maths.h"
#include "igami.c"
#include <stdio.h>


int main() {
  double t = -1;
  double x = 0.9;
  double y = 0.292000000000000037;
  t = igami(x,y);
  printf("%.10f\n",t);
  return(0);
}

This returns the correct value of 1.1000096747.
This would seem to disqualify the C library itself being buggy, I think. I tried it at the -O3 optimization level and it still worked correctly. I can post the standalone C code if anyone wants to try it.

I'm at a loss, this is really bizarre.

from owl.

I edited the igami function in-place so it would run in the same context (being called from OCaml) and

• isnormal(x) returns false.
• x == x returns true
• isnan(x) returns false
• fpclassify(x) returns FP_NORMAL
• Returning to OCaml, Float.classify_float says nan
• the hexadecimal value is 0xFFF8000000000000

from owl.

I have identified the source of the bug.

As documented in the INSTALL.md, by default configure.ml sets the default C compiler flags to

OWL_CFLAGS="-g -O3 -Ofast -funroll-loops -ffast-math -DSFMT_MEXP=19937 -msse2 -fno-strict-aliasing"`

Removing -Ofast and -ffast-math from the configuration causes the bug to disappear. The problem originates from -ffinite-math-only, one of the flags set by -ffast-math. The igami algorithm contains multiple checks for isnan and to check whether a value is infinite, see here, here and here.

When the -ffinite-math-only flag is passed to the compiler, these checks are ignored (they are optimized away under the assumption that your program contains no nan's). Checking for whether a number is infinite is likely a core part of the algorithm which cannot be ignored.

Because the code relies on numerous calls to owl_isnan, and references owl_posinf and other related values, it seems imprudent to use the -ffast-math flag which eliminates all such function calls and references. x<OWL_POSINF will be optimized to true. It is at least an inconsistency: if one does not need the isnan checks, why write them in the first place?

I propose a git pr of changing the default OWL_CFLAGS to omit -ffast-math and -Ofast.

I believe that beyond the immediate bug, this change would be more in keeping with the spirit of the OCaml community, which values the correctness guarantees offered by static typing and would also probably appreciate compliance with IEEE 754 standards about floating point numbers. -ffast-math includes -funsafe-math-optimizations. The gcc manual says that this option “allows optimizations for floating-point arithmetic that (a) assume that arguments and results are valid and (b) may violate IEEE or ANSI standards.” This is somewhat embarrassing given the goal stated in the README to "provide both researchers and industry programmers a powerful framework to write concise, fast and safe analytical code" (emphasis mine).

@jzstark @mseri any thoughts?

from owl.