Why You Cannot (Yet) (Portably) Write an "Interval Arithmetic" Library in Common Lisp

Hi

it has been a while since I posted anything here.  Common Lisp is a nice "comfort zone" for me and I essentially play around with it, in ways that eventually make me wander in the wide ocean of knowledge out there.

Here is the takeaway: I spent some (too much?) time thinking about parts of Common Lisp related to numerical computations and the result is this arXiv paper with the same title of this blog post.

Everything started when I was looking for a cute exercise/project to give students. I had settled on "write an Interval Arithmetic library".  Interval Arithmetic (IA) is an interesting topic per se, with many not-so-trivial mathematical underpinning.  Much work have been done to provide implementation and standards on current processor architectures and, especially, on current IEEE-754 floating standard libraries.

A very simple interval representation in Common Lisp is the following.

(defstruct (int (:constructor int (l h))
  (l 0.0 :type real)
  (h 0.0 :type real))

Given this representation, a very simple interval sum operation could be written as

(defun int-sum (i1 i2)
  (int (+ (int-l i1) (int-l i2))
       (+ (int-h i1) (int-h i2))))

The other operations pretty much follow the same pattern until you get to division, which opens up a can of worms. But let's forget about division and see what other issues pop up.

IEEE-754 infinities etc...

First of all you have the issue of what happens in "extreme" case. What should the following sums return?

(defvar i1 (int 0.0 1.0))

(defvar i2 (int 42.0 most-positive-double-float)

(int-sum i1 i2)

(int-sum i2 i2)

Common Lisp does not quite specify how to handle these cases, therefore implementations may differ about what they return or how they handle overflows. Note also that having the feature :ieee-floating-point set does not mean that your implementation provides all the bells and whistles of IEEE-754. In particular you do not have access to a portable representation of IEEE-754 infinities (and NaNs).

Control on Rounding Direction

Another issue that the code above has, is that the two + yielding the lower and upper bound of the result interval should not behave in the same way; i.e., as we are talking floating point numbers, there are issues of rounding to be considered. Most IA specifications assume IEEE-754, or, more recently, the Language Independent Arithmetic ISO/IEC-10967 (LIA) standard, which define rounding modes and operations. Common Lisp implementations usually round to nearest, while the IA specifications require a more fine grained control on rounding directions; the int-sum function above should actually be written as follows,

(defun int-sum (i1 i2)
  (int (+↓ (int-l i1) (int-l i2))
       (+↑ (int-h i1) (int-h i2))))

where +↓ and +↑ compute results rounding downward and upward.

Again, Common Lisp implementations vary quite a lot in providing access to IEEE rounding modes, if at all, leaving he programmer in a bind.

So, What's Next?

This has been a learning experience for me. Once I found out about LIA, I took some time to try to grok it. The result is the arXiv paper I mentioned, where I advocate a community effort to come up with a specification for a LIA-compliant arithmetic library/API for Common Lisp. Note that I dare to claim that LIA is a superset of any facility provided in Common Lisp and this is the reason why I would like to suggest a different approach to "language evolution" in this case: specify first, implement later.

Maybe some of the Common Lisp regulars out there will like to join the fun. Will you?

---

(cheers)