80-bit extended precision floating point in OCaml

fish001

80-bit extended precision floating point in OCaml [Copy link]

Is there an OCaml library that can take advantage of the 80-bit extended precision floating point types on the IA-32 and x86-64 architectures?
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 language's ffi support.
The library must be split into two parts: a native OCaml source part, and a C runtime part.
The OCaml source must contain data type declarations, as well as declarations of all imported functions. For example, adding an 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 representation to a native C value, use binary operators on them, and return a new OCaml value. The ml_float80_copy function performs the assignment of that runtime representation.

Similarly, the C implementations of the sub, mul, and div functions should also be defined. You can notice the similarity in the signatures and implementations of these functions, which are abstracted by using 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 modules should follow.

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

type t = string
On the C side, we define functions that convert OCaml values back and forth to C values:

#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, this implementation does not support the polymorphic comparison functionality built into OCaml Pervasive.compare, as well as some other features. Using that function on the float80 type above will mislead the comparison function, because the values are strings, and their contents are compared lexicographically.

Supporting these special features is straightforward though. We redefine the OCaml type as an abstract class and change the C code to create and handle a custom structure for 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 recommend building the whole thing using ocamlbuild and a small bash script.