[Maxima-discuss] Careful arithmetic ... was Re: +-Inf and NaN

[Richard F. <fa...@gm...>]

If we know how to extract and re-pack stuff into a NaN,  we can
elaborate on something like this.  (in SBCL Maxima)

  ... where I have defined cexpt for "careful expt" ...

(%i2) cexpt(2.0,1000000);

#<FLOATING-POINT-OVERFLOW{10038EB283}>   /* a printed message */

(%o2) 2.0^1000000

 /* see the result, but internally it is...  ((MEXPT SIMP
#<FLOATING-POINT-OVERFLOW{10038EB283}>) 2.0 1000000)  */

Here is the lisp definition of cexpt

(defun $cexpt (x y); careful expt
  (catch 'float-catch  (handler-bind
            ((error
               #'(lambda(c)
                    (print c)  ;; just for debugging ...
                    (throw 'float-catch  `((mexpt simp ,c) ,x  ,y)))))
              (expt x y))))

HOWEVER, running this example in SBCL sets some flags so that it
doesn't work a second time, returning instead
#.SB-EXT:DOUBLE-FLOAT-POSITIVE-INFINITY
Suggestions from SBCL fans?
(to be specific, there are 2 areas )
1. how to pack a pointer to, say,  the result 2.0^1000000 into a NaN.
    so cexpt will then return that NaN as its result...
2. fix sbcl flags so that the handler still works after the first
floating-point error..

The idea then is to make "careful" versions of all floating-point
arithmetic,
include sin/cos/tan.

Instead of using NaNs, in Maxima
we can also intercept the float errors and put in our
own symbolic Maxima names like inf, minf, und,

or perhaps we can  return intervals.
(I'm not entirely in agreement with this abuse of 'intervals', but
Mathematica says Sin[Infinity] is Interval[{-1,1}] .
Of course then you need to be able to compute Sin[Interval[...]]   )

RJF


On Wed, Feb 28, 2024 at 11:14 AM Richard Fateman <fa...@gm...> wrote:

> Yes for your 2nd and 3rd points.  I was thinking something like...
> m:ident(4)$
> m[2,3]:eps$
> m[1,2]: x$
> m^^2;
>
> where the vast majority of operations could be done with floats, full
> speed.
> Only 4 entries are non-numbers...
>
>
> [1, 2*x, eps*x, 0],
> [0, 1, 2*eps, 0],
> [0, 0, 1, 0],
> [0, 0, 0, 1]
>
> This kind of approach could be used to establish how some small parameter
> like eps
> might enter into the algorithm that has been in use for purely numeric
> calcs.
> RJF
>
>
> On Wed, Feb 28, 2024 at 10:58 AM Stavros Macrakis <mac...@gm...>
> wrote:
>
>> I'm not sure I see the point of encoding symbols into NaN payloads.
>> When would you want to do this?
>>
>>    - If it is mostly a symbolic calculation, this seems like a very
>>    roundabout way of invoking simplification (doesn't trap handling have a
>>    much higher overhead than a procedure call?).
>>    - If it is mostly a numerical calculation, do you expect that only a
>>    small proportion of numbers will become symbolic if (say) just one of
>>    the inputs is made symbolic? Or will most of them sooner or later become
>>    symbolic by contagion?
>>    - I suppose that you could use this approach to replace float64 with
>>    bfloat. The overhead would be high, but it might be useful for seeing the
>>    effect of higher precision on the calculation.
>>
>> How many numerical algorithms lend themselves to this approach? Anything
>> with iterative convergence presumably will fail, for example.
>>
>>
>>
>> On Wed, Feb 28, 2024 at 1:39 PM Richard Fateman <fa...@gm...>
>> wrote:
>>
>>> Years (decades?) ago I gave up on clever use of NaNs combined with
>>> speed.  If trapping
>>> is enabled, but an exception sends the processor into some painful
>>> handler, that's
>>> still a prospect for diagnostics. Presumably the handler is given some
>>> return address
>>> and maybe some other bits, and this info could be stuffed into some
>>> payload of the
>>> NaN's mantissa.  I have not seen it done.
>>> One prospect that I think I've mentioned is that one could encode a
>>> computer
>>> algebra system into a NaN-based framework.  Floating-point numbers are
>>> whatever
>>> numbers they represent. A NaN is a pointer into a table that says.. this
>>> NaN is the
>>> symbol X,   or perhaps a Lisp S-expression.  The trap handling would
>>> call a simplifier
>>> program when the operation is ADD of two NaNs representing X and Y, to
>>> make a
>>> new NaN pointing to X+Y.   etc.
>>> It seems to me this would be a fun project, but I have not studied
>>> current facilities, and
>>> probably higher-level language support for this is lacking. But maybe
>>> this is an opportunity
>>> for Lisp.
>>> .
>>> rjf
>>>
>>> On Wed, Feb 28, 2024 at 9:51 AM Henry Baker <hb...@pi...>
>>> wrote:
>>>
>>>> Re: use the mantissa of a NaN for storing info
>>>>
>>>>
>>>>
>>>> Great idea, but doesn't work in practice.
>>>>
>>>>
>>>>
>>>> The reason?
>>>>
>>>>
>>>>
>>>> The stuff you might want to store in the mantissa isn't normally
>>>> available to the FP HW,
>>>>
>>>> so you would normally consider trapping.  But as I found to my
>>>> consternation, trapping
>>>>
>>>> NaN's doesn't work/is useless, so you end up with an if-then-else
>>>> expression instead of
>>>>
>>>> a simple HW FP op.  You've now lost most -- if not all -- of your HW FP
>>>> performance.
>>>>
>>>>
>>>>
>>>> Part of the problem is that someone needed to do at least one PhD
>>>> thesis on what sorts
>>>>
>>>> of things one might want to store in said mantissa, but to my
>>>> knowledge, no such PhD
>>>>
>>>> thesis was ever written.  There's not much point in optimizing HW for
>>>> some task that's
>>>>
>>>> never been researched before in SW.
>>>>
>>>>
>>>>
>>>> Re: uninitialized/missing data:
>>>>
>>>>
>>>>
>>>> Flagging uninitialized data might be about the *only* useful use for
>>>> trapping NaN's.
>>>>
>>>>
>>>>
>>>> Missing data flags probably want to propagate transitively w/o trapping
>>>> -- e.g., what
>>>>
>>>> happens in spreadsheets when some entry is blank.
>>>>
>>>>
>>>>
>>>> -----Original Message-----
>>>> From: Richard Fateman <fa...@gm...>
>>>> Sent: Feb 28, 2024 8:11 AM
>>>> To: Camm Maguire <ca...@ma...>
>>>> Cc: Stavros Macrakis <mac...@gm...>, <
>>>> max...@li...>, Robert Dodier <
>>>> rob...@gm...>, <gcl...@gn...>
>>>> Subject: Re: [Maxima-discuss] +-Inf and NaN
>>>>
>>>>
>>>> Not sure about* most* users.  You can read my old article about IEEE
>>>> 754 and programming languages.
>>>> At least the intent for NaNs was that they would propagate through a
>>>> computation and
>>>> emerge only if relevant. At least under some conditions.
>>>> I personally have no realistic use cases for using traps for numerics,
>>>> or for propagating NaNs,
>>>> but I'm not a professional error analyst; few if any Lisp programmers
>>>> are... (Ray?)
>>>> If all that is being questioned is the default enabling of trapping, it
>>>> doesn't seem much
>>>> of an issue.
>>>>   It would be nice to see if we can (say) fill an "uninitialized" float
>>>> array by
>>>> initializing it with NaNs, or use the mantissa of a NaN for storing
>>>> info.
>>>> RJF
>>>>
>>>> On Wed, Feb 28, 2024 at 5:44 AM Camm Maguire <ca...@ma...>
>>>> wrote:
>>>>
>>>>> Greetings!
>>>>>
>>>>> Raymond Toy <toy...@gm...> writes:
>>>>>
>>>>> > On Mon, Feb 26, 2024 at 5:28 AM Camm Maguire <ca...@ma...>
>>>>> wrote:
>>>>> >
>>>>> >  Greetings, and thanks to all for the very helpful feedback!
>>>>> >
>>>>> >  I'd also like to note in passing that the excellent SBCL out of the
>>>>> box
>>>>> >  triggers an error, which we separately refer to as an 'exception
>>>>> trap',
>>>>> >  when NaN is passed to the comparison functions, e.g. '=.  And low
>>>>> and
>>>>> >  behold, gcl does the same when 'trapping' is turned on, but gcl
>>>>> turns
>>>>> >  trapping off 'out of the box'.  NaN does not trigger an exception
>>>>> nor
>>>>> >
>>>>> > I hope you will consider changing the default to trap on invalid
>>>>> > operations so that instead of returning NaN, an error is signaled.  I
>>>>> > certainly hate this behaviour, because after you do a long set of
>>>>> > numerical computations, all the results are NaN.  I've just wasted a
>>>>> > whole bunch of time.  This is super-annying in Javascript where you
>>>>> > don't have an option of disabling traps.  Now you get to go through
>>>>> > all the code and either put prints or checks for NaN to see what
>>>>> > caused it.  What a colossal waste of time when the computer could
>>>>> have
>>>>> > told me.
>>>>> >
>>>>>
>>>>> This is good to hear, as I was of the opposite impression that most
>>>>> users wanted NaNs to propagate freely.
>>>>>
>>>>>
>>>>> >  error when passed to '+,'*,'cos etc. but returns NaN.  This is the
>>>>> >  current state of play.  What follows is a discussion of what it
>>>>> means
>>>>> >  and the best way to fit NaN into the common-lisp type lattice so it
>>>>> does
>>>>> >  not get in the way of the compiler.
>>>>> >
>>>>> >  In summary, NaN seems indistinguishable to me from an ignorable
>>>>> error
>>>>> >  indicator at the hardware level.  It is not conceptually a 'valid
>>>>> input'
>>>>> >  to any of these functions, but is a clever design to propagate the
>>>>> error
>>>>> >  properly when ignored.  Conceptually, it seems no different from
>>>>> >  'ignore-errors in common-lisp.  So really the spec is not all that
>>>>> far
>>>>> >  away from IEEE.
>>>>> >
>>>>> > I suppose that that is true, but I don't think any Lisp programmer
>>>>> > would wrap their entire programs in a ignore-errors form.  That's
>>>>> > usually one done on a very limited scope.  Having traps off to get
>>>>> > NaN, is basically the same idea.
>>>>> >
>>>>>
>>>>> Wrapping code in ignore-errors is totally unnecessary.  I can insert a
>>>>> specific handler within GCL's type-error mechanism to skip the error
>>>>> when testing NaN depending on a global variable setting.  One question
>>>>> is whether such a setting governing skipping a lisp error should be the
>>>>> same as or distinct from the setting enabling hardware traps.
>>>>>
>>>>> >  As for common-lisp types, the obvious path is what GCL and others
>>>>> >  currently implement, double precision NaN is of type 'long-float and
>>>>> >  hence 'number.  This actually means that lisp will *stray* from
>>>>> IEEE and
>>>>> >
>>>>> > I think one thing that confuses me, is that for gcl, single-float =
>>>>> > double-float = double precision.  short-float is single precision.  I
>>>>> > think long-float is double-precision.  I think all other lisps have
>>>>> > single-float = single precision, double-float = double precision.
>>>>> > Short-float is usually the same as single-float.  Long-float is
>>>>> > usually double-float, if they don't have a separate long-float (like
>>>>> > clisp).
>>>>> >
>>>>> > So whenever you say long-float, I sometimes think of actual
>>>>> > long-floats like clisp, but you're actually talking about
>>>>> > double-precision floats.
>>>>>
>>>>> Thanks, this needs to be cleaned up.  I think single-float should be
>>>>> short-float (e.g. 4 bytes).
>>>>>
>>>>> I suppose in principle each of these could be distinct, long-float
>>>>> could
>>>>> be long double, supported on most machines, and short float could be a
>>>>> new 16bit float :-).  It looks like the next GCL version will be able
>>>>> to
>>>>> pass 4 unboxed arguments of the same size on 64bit machines, 'object'
>>>>> (e.g. void *), unsigned long, double, and float complex.  I have space
>>>>> in the argument descriptor for 4 such formats, but on 32bit machines,
>>>>> we
>>>>> only have three useful ones, given the absence of 'short float
>>>>> complex'.
>>>>>
>>>>> Take care,
>>>>> --
>>>>> Camm Maguire
>>>>> ca...@ma...
>>>>>
>>>>> ==========================================================================
>>>>> "The earth is but one country, and mankind its citizens."  --
>>>>> Baha'u'llah
>>>>
>>>>
>>>>
>>>