gamma_isf nan

[patrick-nicodemus]

I am trying to see what needs to be done to port Owl to OCaml 5.2 and I found a problem in the gamma_isf function, which depends upon the igami function from the cephes mathematics library.
https://github.com/owlbarn/owl/blob/main/src/owl/maths/cephes/igami.c

I am in OCaml 5.2. I am using gcc (Ubuntu 13.2.0-23ubuntu4) 13.2.0 on Ubuntu 24.04.
I reproduced the error on OCaml 4.14.2.
It works with the main branch of owl as well as the version of owl on opam.

To reproduce the error, open utop, #require "owl" (disabling Float16 warnings if necessary), and

#use "./owl/src/owl/stats/owl_stats_dist.ml";;

Then run

gamma_isf 0.292000000000000037 ~shape:0.9 ~scale:1.;;

(where this number is 1.0 - 0.708)

This returns nan and it really shouldn't, the correct value is not extreme. It's somewhere in the vicinity of 1.10000967467090516

I am using Ubuntu 24.04, which is very recent so may have bugs.

I was unable to reproduce the error on Ubuntu 22.04 and gcc 11. It may be a problem with my configuration.

I was able to reproduce the error on a second Ubuntu 24.04 + gcc (Ubuntu 13.2.0-23ubuntu4) 13.2.0 machine.

[patrick-nicodemus]

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.

[patrick-nicodemus]

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

[patrick-nicodemus]

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?