Re: [MIT-Scheme-devel] floating-point environment

[Matt Birkholz]

> From: Taylor R Campbell <address@hidden>
> Date: Thu, 14 Oct 2010 02:12:03 +0000
> 
> [...] But if you mean readers and writers for the state which one
> can use in DYNAMIC-WIND, well, that's what I proposed --
> FLO:MASKED-EXCEPTIONS, FLO:SET-MASKED-EXCEPTIONS!, &c., with a
> thread-local floating-point environment.

The latter, yes, thank you.  Dyno-wind the proverbial save/restore
context switch around a thunk, but stay to THIS side of the native /
alien divide -- the two stack&register save/restore swaps that get us
in and out of (shudder) libc.

> [...]
>    Do you expect the condition system to be too costly, thus deferred
>    exceptions (which can be checked and possibly cleared later) are
>    preferable?
> 
> The machines and C libraries give us two choices in exceptional
> situations: trap; or give a well-defined, if undesirable, result (0,
> NaN, +/-inf) on which further arithmetic can happen, and set the
> exception flags.  The condition system can be used to handle traps;
> using it to substitute results deep inside a computation is not as
> straightforward.

I peeked at the SIGFPE handler, but it will take more than a peek.

Does the SIGFPE handler not do the usual "request interrupt" thing?
It is a "trap", not "interrupt" handler, so its continuation is... ?

My simple-minded theory (all I can fit in my head :-) is that I might
get my flo:distance procedure

    (define (flo:distance x1 y1 x2 y2)
      (flo:sqrt (flo:+ (flo:square (flo:- x1 x2))
                       (flo:square (flo:- y1 y2)))))

to compile down to something like

        (entry (arity 4 0 #f) "flo:distance-1")
    flo:distance-1:
        (interrupt-test procedure)
        (load (fr  6) (arg-loc 0))     ;x1
        (load (fr  7) (arg-loc 1))     ;y1
        (load (fr  8) (arg-loc 2))     ;x2
        (load (fr  9) (arg-loc 3))     ;y2
        (- (fr 10) (fr  6) (fr  8))
        (- (fr 11) (fr  7) (fr  9))    ;(flo:- x1 x2)
        (* (fr 12) (fr 10) (fr 10))    ;(flo:square (flo:- x1...))
    Oy!:
        (* (fr 13) (fr 11) (fr 11))    ;(flo:square (flo:- y1...))
        (+ (fr 14) (fr 12) (fr 13))    ;(flo:+ (flo:square...)...)
        (cons-float (wr 5) (fr 14))
        (store (wr 5) (arg-loc 0)
        (pop 3)
        (tail-call flo:sqrt)

i.e. enjoying all manner of float instructions and temporary
registers, consing just one float, etc. -- an ideal case.  In this
simple-minded ideal, any floating-point exceptions (SIGFPEs) will be
noted (e.g. at "Oy!" above), but the Scheme machine will not invoke an
interrupt (trap?)  handler until the next interrupt check point, thus
the continuation of the trap will be pseudo-random (e.g. the
application of flo:sqrt).

Assuming my simple mind has grasped any of The Truth, there is no
especially "Scheme friendly" mode of operation -- just a default
expectation that any floating-point exception will be signaled sooner
rather than later.  And the continuation of such a condition will be
pseudo-random -- e.g. the pop-return of a general-purpose procedure
like flo:distance.  Yes?

> As for performance costs, I don't know what they will be.

I suppose sticking explicit tests and jumps (over
internal-continuations) in the code would trash any processor
pipelines and put the costs in the "Ouch!" magnitude.  I am sorry I
even brought it up. :-}

> [...]  But if we instead had primitives to read floating-point
> values from and write them to octet vectors, which could be
> relatively easily open-coded, I imagine the difference between
> frobbing the machine's floating-point exception mask versus
> involving the condition system would be pretty substantial.

We have "floating-point vector primitives"

        flo:vector-cons
        flo:vector-ref
        flo:vector-set!
        flo:vector-length

and they appear to already be open-coded!  I am looking forward to
using them.  Bonus question: how do I write the following so that it
uses fancy SIMD instructions? :-)

    (define (transform-3d-points transform points)
      (define-integrable (%+ x y z) (flo:+ x (flo:+ y z)))
      (map
       (lambda (point)
         (let ((new (flo:vector-cons 3)))
           (flo:vector-set! new 0 (%+
                                   (flo:* (flo:vector-ref point 0)
                                          (flo:vector-ref transform 0))
                                   (flo:* (flo:vector-ref point 1)
                                          (flo:vector-ref transform 1))
                                   (flo:* (flo:vector-ref point 2)
                                          (flo:vector-ref transform 2))))
           (flo:vector-set! new 0 (%+
                                   (flo:* (flo:vector-ref point 0)
                                          (flo:vector-ref transform 3))
                                   (flo:* (flo:vector-ref point 1)
                                          (flo:vector-ref transform 4))
                                   (flo:* (flo:vector-ref point 2)
                                          (flo:vector-ref transform 5))))
           (flo:vector-set! new 0 (%+
                                   (flo:* (flo:vector-ref point 0)
                                          (flo:vector-ref transform 6))
                                   (flo:* (flo:vector-ref point 1)
                                          (flo:vector-ref transform 7))
                                   (flo:* (flo:vector-ref point 2)
                                          (flo:vector-ref transform 8))))
           new))
       points))