80-bit extended precision floating-point in OCaml

Is there an OCaml library to take advantage of the 80-bit extended precision floating-point type on the IA-32 and x86-64 architectures?

I know about the MPFR bindings but my ideal library would be more lightweight. Taking advantage of the historical floating-point instructions would be ideal.

The implementation of such a library is possible outside the compiler, thanks to the ffi support of the language.

The library must be split in two parts: the native ocaml source part, and the C runtime part. the OCaml source must contain the datatype declaration, as well as the declaration of all the imported functions. For instance, the add operation would be:

(** basic binary operations on long doubles *)
external add : t -> t -> t = "ml_float80_add"
external sub : t -> t -> t = "ml_float80_sub"
external mul : t -> t -> t = "ml_float80_mul"
external div : t -> t -> t = "ml_float80_div"

in the C code, the ml_float80_add function should be defined, as described in the OCaml manual:

CAMLprim value ml_float80_add(value l, value r){
   float80 rlf = Float80_val(l);
   float80 rrf = Float80_val(r);
   float80 llf = rlf + rrf;
   value res = ml_float80_copy(llf);
   return res;
}

Here we convert the OCaml value runtime representations to native C values, use the binary operator on them, and return a new OCaml value. the ml_float80_copy function does the allocation of that runtime representation.

Likewise,the C implementations of sub, mul and div functions should be defined there too. One can notice the similarity in signature and implementation of these functions, and abstract away through the use of C macros:

#define FLOAT80_BIN_OP(OPNAME,OP)                   \
  CAMLprim value ml_float80_##OPNAME(value l, value r){     \
    float80 rlf = Float80_val(l);                           \
    float80 rrf = Float80_val(r);                           \
    float80 llf = rlf OP rrf;                               \
    value res = ml_float80_copy(llf);           \
    return res;                     \
  }


FLOAT80_BIN_OP(add,+);
FLOAT80_BIN_OP(sub,-);
FLOAT80_BIN_OP(mul,*);
FLOAT80_BIN_OP(div,/);

The rest of the OCaml and C module should follow.

There are many possibilities as to how encode the float80 C type into an OCaml value. The simplest choice is to use a string, and store in it the raw long double.

type t = string

On the C side, we define the functions to convert an OCaml value back and forth to a C value:

#include <caml/mlvalues.h>
#include <caml/alloc.h>
#include <caml/misc.h>
#include <caml/memory.h>


#define FLOAT80_SIZE 10  /* 10 bytes */

typedef long double float80;

#define Float80_val(x) *((float80 *)String_val(x))

void float80_copy_str(char *r, const char *l){
   int i;
   for (i=0;i<FLOAT80_SIZE;i++)
      r[i] = l[i];
}

void store_float80_val(value v,float80 f){
   float80_copy_str(String_val(v), (const char *)&f);
}

CAMLprim value ml_float80_copy(value r, value l){
   float80_copy_str(String_val(r),String_val(l));
   return Val_unit;
}

However, that implementation doesn't bring support for the polymorphic comparison functions built into OCaml Pervasive.compare, and a few other features. Using that function on the above float80 type will mislead the comparison function into beleiving that the values are strings, and do a lexicographical comparison on their content.

Supporting these special features is simple enough though. We redefine the OCaml type as abstract, and change the C code to create and handle custom structs for our float80:

#include <caml/mlvalues.h>
#include <caml/alloc.h>
#include <caml/misc.h>
#include <caml/memory.h>
#include <caml/custom.h>
#include <caml/intext.h>

typedef struct {
   struct custom_operations *ops;
   float80 v;  
} float80_s;

#define Float80_val(x) *((float80 *)Data_custom_val(x))

inline int comp(const float80 l, const float80 r){
   return l == r ? 0: (l < r ? -1: 1); 
}

static int float80_compare(value l, value r){
   const float80 rlf = Float80_val(l);
   const float80 rrf = Float80_val(r);
   const int llf = comp(rlf,rrf);
   return llf;
}

/* other features implementation here */

CAMLexport struct custom_operations float80_ops = {
  "float80", custom_finalize_default, float80_compare, float80_hash,
  float80_serialize, float80_deserialize, custom_compare_ext_default
};

CAMLprim value ml_float80_copy(long double ld){
  value res = caml_alloc_custom(&float80_ops, FLOAT80_SIZE, 0, 1);  
  Float80_val(res) = ld;
  return res;
}

We then propose to build the whole thing using ocamlbuild and a small bash script.