[Caml-list] Problem with Mpfr rounding modes (FFI interface)

[Caml-list] Problem with Mpfr rounding modes (FFI interface)

[Andrej Bauer]

I seem to be hitting some sort of a bug in the interface between OCaml
and the MPFR library for multiple-precision floating point.

I am on macOS 10.13.1, Ocaml 4.04.1 (but also 4.02.1 behaves the same
way), OPAM package mlgmpidl 1.2.4. MPFR was installed via Homebrew and
it is version 3.1. The following complete program shows that rounding
modes do not work correctly. Let me describe the problem.

-----------test_rounding.ml----------------
(* Testing rounding modes in Mpfr. *)
(* Requires OPAM package mlgmpidl *)
(* Compile with:
   ocamlbuild -use-ocamlfind -pkg gmp test_rounding.native *)

let test_modes prec x =
   (* Create MFPR variables with the given precision, one per rounding mode. *)
   let a_Near : Mpfr.t = Mpfr.init2 prec
   and a_Zero : Mpfr.t = Mpfr.init2 prec
   and a_Up : Mpfr.t = Mpfr.init2 prec
   and a_Down : Mpfr.t = Mpfr.init2 prec
   and a_Away : Mpfr.t = Mpfr.init2 prec
   and a_Faith : Mpfr.t = Mpfr.init2 prec
   and a_NearAway : Mpfr.t = Mpfr.init2 prec
   in
   (* Set the variables with various rounding modes *)
   ignore (Mpfr.set_d a_Near x  Mpfr.Near) ;
   ignore (Mpfr.set_d a_Zero x  Mpfr.Zero) ;
   ignore (Mpfr.set_d a_Up x  Mpfr.Up) ;
   ignore (Mpfr.set_d a_Down x  Mpfr.Down) ;
   ignore (Mpfr.set_d a_Away x  Mpfr.Away) ;
   ignore (Mpfr.set_d a_Faith x  Mpfr.Faith) ;
   ignore (Mpfr.set_d a_NearAway x  Mpfr.NearAway) ;
   (* Print the results *)
   Format.printf "x = %F@." x ;
   Format.printf "a_Near = %t@." (fun ppf -> Mpfr.print ppf a_Near) ;
   Format.printf "a_Zero = %t@." (fun ppf -> Mpfr.print ppf a_Zero) ;
   Format.printf "a_Up = %t@." (fun ppf -> Mpfr.print ppf a_Up) ;
   Format.printf "a_Down = %t@." (fun ppf -> Mpfr.print ppf a_Down) ;
   Format.printf "a_Away = %t@." (fun ppf -> Mpfr.print ppf a_Away) ;
   Format.printf "a_Faith = %t@." (fun ppf -> Mpfr.print ppf a_Faith) ;
   Format.printf "a_NearAway = %t@." (fun ppf -> Mpfr.print ppf a_NearAway)
;;

test_modes 11 1024.3333333333333 ;;
-------------------------------------------------

The output is:

x = 1024.33333333
a_Near = 0.10240E4
a_Zero = 0.10240E4
a_Up = 0.10240E4
a_Down = 0.10240E4
a_Away = 0.10250E4
a_Faith = 0.10250E4
a_NearAway = 0.10250E4

The problem is that a_Up is smaller than x, and that it is equal to
a_Down. In fact, the rounding modes behave in strange ways (try
plugging in other values of x). The problem persists if I try to
convert from string, or copy one Mpfr float into another, etc.

I am not familiar with the FFI interface, but I can tell you where the
rounding modes are defined.
In the mlgmpidl library there is in mpfr.ml the type round, defined as

type round =
  | Near
  | Zero
  | Up
  | Down
  | Away
  | Faith
  | NearAway

The corresponding enum type in the C library is defined in
/usr/local/include/mpfr.h, defined as

typedef enum {
  MPFR_RNDN=0,  /* round to nearest, with ties to even */
  MPFR_RNDZ,    /* round toward zero */
  MPFR_RNDU,    /* round toward +Inf */
  MPFR_RNDD,    /* round toward -Inf */
  MPFR_RNDA,    /* round away from zero */
  MPFR_RNDF,    /* faithful rounding (not implemented yet) */
  MPFR_RNDNA=-1 /* round to nearest, with ties away from zero (mpfr_round) */
} mpfr_rnd_t;

Conversion from the OCaml round to the C mpfr_rnd_t is in mlgmpidl
library in gmp_caml.h and it looks like this:

static inline
void camlidl_mpfr_rnd_t_ml2c(value val, mpfr_rnd_t* rnd)
{
  int arg = Int_val(val);
#if MPFR_VERSION_MAJOR >= 3
  arg = (arg==6) ? (int)MPFR_RNDNA : arg;
#endif
  *rnd = (mpfr_rnd_t)arg;
}

So, it's using Int_val to convert the OCalml tag to the enum value. Is
this the way to do it? Has there been a change in recent OCaml
versions? Or am I just doing something else wrong? Any help is
appreciated.

With kind regards,

Andrej