80-bit extended precision floating-point in OCaml

Question

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.

@didierc A simple C library with bindings is exactly what I am looking for, but it cannot be very simple: there is 32- and 64-bit OCaml to handle and various conversion functions (to and from float, string, big_int). It is all perfectly doable and the performance would be fine but it's enough work to be worth avoiding duplication. — Pascal Cuoq, Nov 21 '12 at 07:31
@didierc Good point. I did not think about that. Indeed, if such a library was to be written today, custom types would be an option to consider, although that would mean one more word of overhead. Now, for the particular application I intended, I would always need four `long double`s and four booleans together. It seems it makes sense for me to pack these tightly together. I will be writing my library in C then. — Pascal Cuoq, Nov 25 '12 at 00:49
@didierc To answer your question more directly, custom types are nice-ish. Serialization is a plus in general but I couldn't use it if it was there (hashconsing forces me to use a more sophisticated system http://blog.frama-c.com/public/unmarshal.pdf ). Still, serialization is useful for most programmers, as opposed to generic comparison and hashing, which it is bad style to rely on because they break for structured types anyway (hashtables, maps, sets, …). — Pascal Cuoq, Nov 25 '12 at 00:55
@subodh Not one of the words you put between backticks in your edit was an OCaml expression or keyword. I reverted. — Pascal Cuoq, Feb 21 '13 at 12:55
@didierc Github is good for this sort of thing, but I would start toying with it even if you put it on pastebin. — Pascal Cuoq, Mar 02 '13 at 07:49

didierc · Accepted Answer · 2013-03-19T18:06:57.827

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.

finally managed to get that repository online, and here's my answer. It's not ideal, since it's a "userland" support, and compiler native support would of course be much much better. But I guess it can work as a substitute until someone provides a better alternative. — didierc, Mar 18 '13 at 16:04

80-bit extended precision floating-point in OCaml

1 Answers1