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".

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

Given this representation, a very simple interval

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.

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

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

Another issue that the code above has, is that the two

where

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

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)

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(:constructorint(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(defunint-sum(i1 i2) (int(+ (int-li1) (int-li2)) (+ (int-hi1) (int-hi2))))

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 (int0.0 1.0)) (defvar i2 (int42.0 most-positive-double-float) (int-sumi1 i2) (int-sumi2 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,(defunint-sum(i1 i2) (int(+↓ (int-li1) (int-li2)) (+↑ (int-hi1) (int-hi2))))

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)

Note that infinities and such are pretty portable with the float-features library. https://shinmera.github.io/portability/#float-features

ReplyDeleteBTW. Would you like to join the fun?

DeleteThanks. Writing a layer for infinities is straightforward. Even NaNs are relatively easy to hanlde. It is their use that requires a more coordinated approach and one of the points I am probably not really making is that implementations should provide much of the ground work as many things implied by the LIA specs are not easily packaged in a library.

ReplyDeleteI would add that, in signal and image processing domains, the use of denormals becomes problematic. We aren't so concerned there with retaining bits of precision (which are fading anyway), but with the consequent destruction of processing speed which falls by orders of magnitude with denormal numbers. We need real-time floating point processing and denormals in that domain can safely be truncated toward zero. We should at least have the option to specify that DSP-behavior.

ReplyDeleteHi David

Deletesorry but I saw this only now.

I understand your concerns. I agree that an API should be provided for that.

I am not sure how though, as there isn't anything in the LIA spec, AFAIU, that would allow you that level of control.