Subject: CLTL compatible type contagion?

From: Raymond Toy <[email-masked]>
Date: Sun, 06 Aug 1995 15:08:30 -0400

I asked about this before, and the answer was no, clisp does not do
type contagion as specified in CLTL2.  The current behavior is to
convert everything to the smallest floating-point type instead of the
largest as most other languages typically do and as specified by CLTL.
Some sensible reasons for this are given in impnotes.txt.

However, I was doing some numerical prototyping with clisp and came
across this very annoying problem with the current type contagion:

	(+ 1d50 1.0)

does not return 1d50.  Instead, it signals a floating-point overflow.
I presume that clisp does the addition in double-float and then
coerces the result to single-float, causing the overflow.

This behavior makes numerical work more difficult than it
should be.

Does anyone know where the routines are that cause this behavior?  I'd
like to change it in my own copy of clisp.  I tried looking through
the source, but I don't read German.

Thanks for any pointers,

Ray

---

From: Marcus Daniels <[email-masked]>
Date: Sun, 6 Aug 1995 19:47:06 GMT

>>>>> "Raymond" == Raymond Toy <[email-masked]> writes:

Raymond> I presume that clisp does the addition in
Raymond> double-float and then coerces the result to single-float,
Raymond> causing the overflow.

Raymond> This behavior makes numerical work more difficult than it
Raymond> should be.

Others may disagree.  :-)

Raymond> Does anyone know where the routines are that cause this
Raymond> behavior?  I'd like to change it in my own copy of clisp.  I
Raymond> tried looking through the source, but I don't read German.

Well, you could take a look at flo_rest.d, and specifically, the GEN_F_op2
macro.  Look at the end of the "RETURN CONCAT" lines, where you see
identifiers like FF_to_SF, DF_to_SF ..  You could change these
to FF_to_SF_, DF_to_SF_, etc.  And immediately above have a
conditional like :

#ifdef CONTAGION_HACK
#define FF_to_SF_(x) x
#define DF_to_SF_(x) x
#define LF_to_SF_(x) x
#define DF_to_FF_(x) x
#define LF_to_FF_(x) x
#define LF_to_DF_(x) x
#else
#define FF_to_SF_(x) FF_to_SF(x)
#define DF_to_SF_(x) DF_to_SF(x)
#define LF_to_SF_(x) LF_to_SF(x)
#define DF_to_FF_(x) DF_to_FF(x)
#define LF_to_FF_(x) LF_to_FF(x)
#define LF_to_DF_(x) LF_to_DF(x)
#endif

This hack just addresses your particular example, though.

Personally, I like being reminded by CLISP about such careless behavior.

---

From: [email-masked] (Richard Shepherd)
Date: Mon, 7 Aug 1995 10:08:54 +1200 (NZST)

> This behavior makes numerical work more difficult than it

This is weird to me because:

(+ 1d50 1)

-> 1.0d50

and

(+ 1d50 1d0)

-> 1.0d50

but as pointed out, adding 1.0 or 1e0 causes the overflow error.  Is
this consistent?

8<--------------------------->8
Richard Shepherd. 
([email-masked])
WWW home page:- http://therat.math.waikato.ac.nz/

---

From: Marcus Daniels <[email-masked]>
Date: Sun, 6 Aug 1995 20:45:28 GMT

I'll spare another the trouble...

12.1.4.2 Rule of Float Approximation

Computations with floats are only approximate, although they are
described as if the results were mathematically accurate. Two
mathematically identical expressions may be computationally different
because of erros inherent in the floating-point approximation
process.  The precision of a float is not necessarily correlated with
the accuracy of that number.  For instance, 3.142857142857142857 is a
more precise approximation to pi than 3.14159, but the latter is more
accurate.  The precision refers to the number of bits retained in the
representation.  When a operation combines a short float with a long
float, the result will be a long float.  Common Lisp functions assume
that the accuracy of arguments to them does not exceed their
precision.  Therefore when two small floats are combined, the result
is a small float.  Common Lisp functions never convert automatically
from a larger size to a smaller size.

Yet, CLISP does, and my understanding is that it does this because:

"Common Lisp functions assume that the accuracy of arguments to them
does not exceed their precision" 

is such a wildly optimistic assumption.

Seems to me it is fairly easy to set things up to avoid the need for
contagion, and the same program will work fine on a CL that doesn't
mind playing it fast and loose with floating point precision.

---

From: [email-masked] (Bruno Haible)
Date: Mon, 7 Aug 95 19:11:28 +0200

> coerces the result to single-float, causing the overflow.

Yes. It would actually make sense to return 1d50 in that case since
1.0e0 is so much less in magnitude. But the problem remains for example
in the case

        (+ 1d40 1e37)

Returning 1.001d40 for this would be wrong (IMO) since the two numbers,
considered as intervals, add up like this:

         1d40 +/- 1e24
         1e37 +/- 1e29
     =================
     1.001d40 +/- 1e29

And if ANSI CL says that:

then of course the implementation shall do its best to guarantee that if
the arguments of the function `+' satisfy this assumption, then the
result of  (+ 1d40 1e37)  satisfies this assumption as well. If it can't,
I prefer to see an error.

My advice is: if you need numbers larger than 1e38 in your computations,
then use double floats or long floats for everything. And if you need
numbers larger than 1e600 then use long floats. Their exponent range
is virtually unlimited in CLISP.

Bruno Haible                        email: <[email-masked]>
Software Engineer                   phone: +33-1-49083585

---

From: Raymond Toy <[email-masked]>
Date: Mon, 07 Aug 1995 14:43:28 -0400

>>>>> "Bruno" == Bruno Haible <[email-masked]> writes:

    Bruno> Raymond Toy writes:
    >> 
    >> (+ 1d50 1.0)
    >> 
    >> does not return 1d50.  Instead, it signals a floating-point
    >> overflow.  I presume that clisp does the addition in

[nice example of why clisp does it this way deleted]

    Bruno> then of course the implementation shall do its best to
    Bruno> guarantee that if the arguments of the function `+' satisfy
    Bruno> this assumption, then the result of (+ 1d40 1e37) satisfies
    Bruno> this assumption as well. If it can't, I prefer to see an
    Bruno> error.

By your example, (+ 1d40 1e37) does not satisfy the assumption right?
Should clisp then signal an error?  

I don't know what the "right" answer would be.  At times I like
promoting to the bigger type, but at other times, I like instantly
seeing when I've lost precision in clisp.

However, I have discovered the "errant" part of my program.  I lose
precision because I had something like

	(sqrt 22/7)

and clisp returned a float value.  I guess the result would depend on
whether 22/7 has infinite precision/accuracy in which case the answer
should be a long-float (?) or whether 22/7 has accuracy of +/- 1/14 in
which case the answer should be short-float or less.  The assumption
that 22/7 is a single-float is a bit hazy.

Ray

---

From: Raymond Toy <[email-masked]>
Date: Mon, 07 Aug 1995 14:55:39 -0400

>>>>> "Marcus" == Marcus Daniels <[email-masked]> writes:

    Marcus> Others may disagree.  :-)

Of course!

    Raymond> Does anyone know where the routines are that cause this
    Raymond> behavior?  I'd like to change it in my own copy of clisp.
    Raymond> I tried looking through the source, but I don't read
    Raymond> German.

[instructions deleted]

    Marcus> This hack just addresses your particular example, though.

Thanks for the hint.  I haven't decided if I really want to do this or
not, though.

    Marcus> Personally, I like being reminded by CLISP about such
    Marcus> careless behavior.

I agree, it is nice, but I wish it could remind me a bit more
loudly. :-) Silently making all (following) results a single-float
makes it very difficult to find out where I made the original
mistake. :-)

Ray

---

From: Raymond Toy <[email-masked]>
Date: Mon, 07 Aug 1995 15:10:09 -0400

Bruno>         (+ 1d40 1e37)

    Bruno> Returning 1.001d40 for this would be wrong (IMO) since the two numbers,
    Bruno> considered as intervals, add up like this:

    Bruno>          1d40 +/- 1e24
    Bruno>          1e37 +/- 1e29
    Bruno>      =================
    Bruno>      1.001d40 +/- 1e29

This is a nice example of why the answer should be a single-float.

But what about this example?

	  1.000001 +/- 1e-8
	- 1.000000 +/- 1e-8
	===================
	  0.000001 +/- 2e-8

What should the type of the result be?  Something less than a
short-float?  The current scheme seems like a half-attempt at interval
arithmetic.  Either go all the way or don't go at all.

Ray

P.S.  If I sound harsh or antagonistic, I'm not.  I'm just really
curious about how clisp does it's arithmetic and the reasoning behind
it.  I have found that clisp is at least as good as gcl and cmulisp
(the only lisp implementations I used) in speed and executable size
for general purpose prototyping.  Clisp only loses when I do
heavy-duty floating point operations, and, then, only if I'm very
careful to declare everything correctly.

---

From: [email-masked] (Bruno Haible)
Date: Mon, 7 Aug 95 22:28:15 +0200

> makes it very difficult to find out where I made the original

Not so difficult: you just have to look at the number of decimal
places of your intermediate results. It would be more difficult to
find out (or even notice) if you got double-floats with the 8 trailing
decimal places filled with nearly random digits -- which is what ANSI CL
suggests.

> that 22/7 is a single-float is a bit hazy.

Agreed. This is a second flaw in CLtL, this time in the description of
the FLOAT function. 22/7, as a rational number, has infinite precision,
so the result should be a long-float of the longest available float
format. Since CLISP can deal with long-floats of more than 100000 decimal
places, this is not very practical.

For this purpose, CLISP provides a variable LISP:*DEFAULT-FLOAT-FORMAT*.
Set it to DOUBLE-FLOAT, and (sqrt 22/7) will return a float of this type.

> What should the type of the result be?

These phenomena (extinction and the growth of the interval length) are
considered to be inevitable with the usual IEEE floating point arithmetic.

> The current scheme seems like a half-attempt at interval arithmetic.

It is.

> Either go all the way or don't go at all.

One should go all the way. Follow-up to the comp.arch.arithmetic newsgroup.