Subject: iPSC Floating Point Exception

From: Cleve Moler <[email-masked]>
Date: Fri, 18 Apr 86 13:45:14 pst

4/16/86

To: Chris Grant, Paul Kautz, Ray Asbury, Peter Ross, David Scott,
    Bill Hughey, Tony Anderson, Paul Pierce, Roger Golliver, Ted Forgeron,
    Beth Ducot (MIT), Virginia Klema (MIT), Ed McDonald (Yale),
    Alan Laub (UCSB), John Bruno (UCSB), Layne Watson (GMR), 
    Paul Fredricksen (CMI), Mike Heath (Oak Ridge), Jack Dongarra (ANL)

From: Cleve Moler

Re: iPSC Floating Point Exception Handling

Here is a preliminary version of a floating point exception handler for
the iPSC.  The code does not yet do everything that I would like it
to do, but it does provide much better handling of exceptions than
results from standard Intel Fortran and the iPSC node operating system.

If you are doing floating point calculations on the iPSC with Fortran,
I invite you to start using this exception handler.  Please tell me
if it changes anything, for better or worse.  If you are doing floating
point calculations on the iPSC with some other programming language,
please switch to Fortran immediately so you can use this handler.

More importantly, I would like to invite those of you that are interested
to use this code as a starting point for a more complete handler.  We
will also continue to work on it with two goals in mind:

     1.  Exception handling consistent with draft 10 of the IEEE
         floating point standard instead of draft 8.

     2.  Provide useful debugging information for error conditions.

The software included with this memo goes part way toward achieving
the first goal, but does very little toward achieving the second.

As you know, the Intel 8087 was designed, and first shipments made, before
the IEEE floating point standard was finally adopted.  The 8087, and hence
the 80287 used in the iPSC, conforms to draft 8 of the standard, rather than
to the final, official, draft 10.  The primary difference between the two
drafts concerns the handling of denormal numbers.  For 64-bit arithmetic
(Fortran double precision), these are floating point numbers in the range
from 2.0**(-1022) to 2.0**(-1074), or roughly 2.22E-308 to 4.94E-324.
The difficulty, and controversy, occurs when denormal numbers are involved
in arithmetic operations that yield results in the normal range, greater
than 2.0**(-1022).  I call such results "promoted denormals".  For example:

            X = 1.0E-305
            Y = 1.0E+10
            Z = X/Y
            A = (X/Y)*Y
            B = 1.0E0 + X/Y

The quantity Z, whose value is 1.0E-315, is a denormal number.  It is
generated correctly by Intel, or any other, Fortran, running on any 8087
machine, with no special exception handling.  No trouble so far.  The
difficulties occur with the next two statements, both of which produce
"promoted denormals".  The quantity A should be equal to X, except that 
several bits of accuracy have been lost because the intermediate denormal 
result is represented with less relative precision.  Should the user be 
warned of this loss of accuracy?  The value stored in B should be exactly 
1.0 since X/Y is numerically negligible compared to 1.0.  Should the warning 
that might be expected for A also be generated for B?

Draft 8 provides for a "warning mode" to signal the promotion of denormals.
But it was judged to be of such esoteric utility that draft 10 makes no 
mention of warning mode.  The 8087 allows for software implementation of 
warning mode, but, as far as I know, no such software has ever been written.

The 8087 produces a NaN, a special bit pattern signifying "Not a Number",
for a promoted denormal.  If the last bit of the 8087 control word is on,
then NaN's propagate through subsequent arithmetic operations, producing
further NaN's.  In this case, both A and B would be NaN's.  It might be
argued that such a result is reasonable for A, but I think that it is hard
to justify for B.  However, Intel Fortran appears to arrange to have the
last bit of the 8087 control word turned off.  (It is not clear to me
where or why this happens, or even that it is sure to happen.)   In this
case, any attempt to store a NaN, or to do arithmetic involving a NaN,
is trapped and declared illegal.  The result is Exception 16, a dead 
program, and little helpful traceback information.

Fortunately, we can now do better than either of these cases.  The exception
handler accompanying this memo allows the generation of a result close
to X for A, an exact 1.0 for B, and continued execution of any subsequent
statements.

We're not done.  So far, we only know how to use the exception handler
in iPSC node programs; we need help getting it to work under Xenix on
the manager.  Elementary functions don't work correctly.  Taking the
square root of a denormal produces a NaN, although the result is well 
within range.  Taking the cosine of a denormal causes the node to hang,
although the result 1.0 would be accurate to over 600 decimal places.
We haven't tried single precision, although that should work.
We haven't tried C, but nobody does floating point in C.

Here is what I would like to see in a complete exception handler.

     --  Promoted denormals are quietly set to the proper normalized
         values with no warning messages or other interruptions.

     --  Underflows, less than 2.0**(-1074), are quietly set to zero
         with no warning messages or other interruptions.

     --  Other exceptional values, like Infinite or NaN, are set to the 
         values defined by the IEEE standard.  Warning messages are
         printed (via SYSLOG) the first few (say 5) times this happens,
         then the messages are turned off. 

     --  When possible, the messages should reference the relevant 
         source line, preferably in the source file, or at least in 
         some listing.

     --  System calls should provide counts of each exception and control
         the printing of warning messages.

     --  Elementary functions should return sensible results for
         exceptional arguments.  

     --  Generation of exceptional floating point values should not
         cause termination.

     --  There is no need to provide values other than the IEEE default
         values for exceptions.

Any volunteers to create such software?  

The remainder of this memo consists of:

     1.  Fortran source code for a test program.
         The program runs on one node, puts its output in SYSLOG,
         and does not require a host program.

     2.  Fortran source for dhex, which prints floating point
         numbers in hexadecimal format with reasonable byte order.

     3.  PL/M source code for the exception handler,
         (written by Paul Pierce).

     4.  A makefile to generate the test program.
         (This requires PLM286.  If you don't have PLM286 and 
         want to contribute to this effort, let me know.)

     5.  The output from the test program.

Please feel free to pass this material on to anybody who is really
interested, or who needs it, but try to keep the distribution limited for
a while.  Please return any comments, results, or improved code to me.
We'll provide a little more formal version to all iPSC users as soon 
as we've had a chance to digest this first attempt.

Thanks for your interest and participation so far.

______________________________________________________________________

      program test
      external trap_handler
      double precision x,y,z
      integer cw
      data cw /#33EH/
c
      call handler(16, trap_handler)
      call ldcw87(cw)
c
      call syslog(0,'test dhex')
      x = 1.0d0
      do 01 k = 1, 13
         x = 16.0d0*x + k
   01 continue
      call out(x)     
c
      call syslog(0,'-1/0')
      y = -1.0d0
      z = 0.0d0
      x = y/z
      call out(x)
c
      call syslog(0,'denormal')
      y = 1.0d-305
      z = y / 1.d10
      call out(z)
c
      call syslog(0,'tiny + denormal')
      x = y + z
      call out(x)
c
      call syslog(0,'promote denormal')
      x = z * 1.d10
      call out(x)
c
      call syslog(0,'1.0 + promoted denormal')
      y = 1.0d0 + x
      call out(y)
c
      call syslog(0,'sqrt(denormal)')
      x = sqrt(z)
      call out(x)
c
      call syslog(0,'overflow')
      y = 1.0d+250
      z = y*y
      call out(z)
c
      call syslog(0,'sqrt(overflow)')
      x = sqrt(z)
      call out(x)
c
      call syslog(0,'0/0')
      x = 0.0d0
      y = 0.0d0
      z = x/y
      call out(z)
c
      call syslog(0,'sqrt(0/0)')
      x = sqrt(z)
      call out(x)
      end
c
      subroutine out(x)
      double precision x
      character*16 dhex
      character*41 string
      write(string,'(1pd23.15,a18)') x,dhex(x)
      call syslog(0,string)
      return
      end

______________________________________________________________________

      character*16 function dhex(x)
      double precision x
c
c     Convert double precision value to hex string
c     Suitable for output with format A16
c
      double precision xx
      integer*4 ix(2),nohi,m
      integer sig,i,j,k,l
      character*1 d(16)
      character*16 s
      equivalence (xx,ix)
      data d/'0','1','2','3','4','5','6','7',
     >       '8','9','A','B','C','D','E','F'/
      data nohi /2147483647/
c
      xx = x
      k = 16
      do 20 i = 1, 2
         m = ix(i)
         if (m .lt. 0) then
            sig = 1
            m = m .AND. nohi
         else
            sig = 0
         endif
         do 20 j = 1, 8
            l = mod(m,16)+1
            m = m/16
            if (j.eq.8 .and. sig.eq.1) l = l+8
            s(k:k) = d(l)
            k = k-1
   20 continue
      dhex = s
      end

______________________________________________________________________

$large
$ram

HANDLER_MODULE:  DO;

DECLARE RETRY_FLAG BYTE;

DECODE: PROCEDURE (ESTATE287_PTR, ERRORS287) EXTERNAL;
	DECLARE ESTATE287_PTR POINTER, ERRORS287 WORD;
END;

ENCODE: PROCEDURE (ESTATE287_PTR, ERRORS287, CONTROL287, RETRY_FLAG) EXTERNAL;
	DECLARE ESTATE287_PTR POINTER;
	DECLARE ERRORS287 WORD;
	DECLARE CONTROL287 WORD;
	DECLARE RETRY_FLAG BYTE;
END;

NORMAL: PROCEDURE (ESTATE287_PTR, ERRORS287) BYTE EXTERNAL;
	DECLARE ESTATE287_PTR POINTER,
			ERRORS287 WORD;
END;

SIEVE: PROCEDURE (ESTATE287_PTR, ERRORS287) BYTE EXTERNAL;
	DECLARE ESTATE287_PTR POINTER,
			ERRORS287 WORD;
END;

DECLARE I$ERROR$BIT LITERALLY '0001H';

declare LIT literally 'LITERALLY';
declare DCL literally 'DECLARE';
declare PROC literally 'PROCEDURE';
declare IS literally 'LITERALLY';
declare FOREVER literally 'WHILE (1)';
declare UNTIL literally 'WHILE NOT';
declare ENDDO literally 'END';
declare ENDDOCASE literally 'END';
declare ENDMODULE literally 'END';
declare ENDPROC literally 'END';
declare ENDTHENDO literally 'END';
declare ENDELSEDO literally 'END';
declare ENDDOWHILE literally 'END';
declare ENDDOUNTIL literally 'END';
declare ENDDOFOREVER literally 'END';
declare BOOLEAN literally 'BYTE';
declare TRUE literally '0FFH';
declare FALSE literally '0';
declare ESC literally '1BH';
declare SPACE literally '20H';
declare CR literally '0DH';
declare LF literally '0AH';

$if DEBUG

declare
	syslogbuf (80) byte,
	syslognum	word initial (0);

syslog: procedure(pid, msgp, msgl) external;
	declare (pid, msgp) pointer, msgl word;
endproc syslog;

/*****************************************************************************/
/*	Write - convert message for SYSLOG
 */

Write: procedure(MsgP, MsgL);
    declare
	MsgP	pointer,
	MsgL	word;
    declare
	i	word,
	Msg	based MsgP (1) byte;

    do i = 0 to MsgL-1;
	if Msg(i) = CR then
	    /* Nothing */;
	else if Msg(i) = LF then
	    do;
		call syslog(@(0, 0), @syslogbuf, syslognum);
		syslognum = 0;
	    endthendo;
	else
	    do;
		syslogbuf(syslognum) = Msg(i);
		syslognum = syslognum + 1;
	    endelsedo;
	if syslognum = 79 then
	    do;
		call syslog(@(0, 0), @syslogbuf, syslognum);
		syslognum = 0;
	    end;
    enddo;
endproc Write;

/*****************************************************************************/

Binary_to_Hex_Ascii: proc (Buffer_Ptr,Precision);

dcl Precision                         integer;
dcl Index                             integer;
dcl Count                             integer;
dcl Buffer_Ptr                        pointer;
dcl Binary_Buffer based Buffer_Ptr (1A0H) byte;
dcl Hex_Index                         byte;
dcl Hex_Ascii_Buffer (25)             byte;

    Count = Precision - 1;
    Hex_Index = 0;
    do Index = Count to 0 by -1;
        Hex_Ascii_Buffer(Hex_Index) = shr(Binary_Buffer(Index) and 0F0H, 4) + '0';
        if Hex_Ascii_Buffer(Hex_Index) > '9' then Hex_Ascii_Buffer(Hex_Index) =
                Hex_Ascii_Buffer(Hex_Index) + 7;
        Hex_Index = Hex_Index + 1;
        Hex_Ascii_Buffer(Hex_Index) = (Binary_Buffer(Index) and 0FH) + '0';
        if Hex_Ascii_Buffer(Hex_Index) > '9' then Hex_Ascii_Buffer(Hex_Index) =
                Hex_Ascii_Buffer(Hex_Index) + 7;
        Hex_Index = Hex_Index + 1;
    endDo;
    call Write (@Hex_Ascii_Buffer, unsign((Count+1)*2));
endProc Binary_to_Hex_Ascii;

DUMP287: PROCEDURE (ESTATE287_PTR, ERRORS287);
	DECLARE ESTATE287_PTR POINTER,
			ERRORS287 WORD;
	declare estate287 based estate287_ptr structure (
		OPERATION WORD,
		ARGUMENT BYTE,
		ARG1(5) WORD, ARG1_FULL BYTE,
		ARG2(5) WORD, ARG2_FULL BYTE,
		RESULT BYTE,
		RES1(5) WORD, RES1_FULL BYTE,
		RES2(5) WORD, RES2_FULL BYTE,
 		FORMAT BYTE,
		REGISTER BYTE,
		CONTROL_WORD WORD,
	 	STATUS_WORD WORD,
		TAG_WORD WORD,
		ERROR_POINTERS(4) WORD,
		STACK_287(40) WORD);
	declare i word;

	call write(@('Op '), 3);
	call binary_to_hex_ascii(@estate287.operation, 2);
	call write(@(' Status '), 8);
	call binary_to_hex_ascii(@errors287, 2);
	call write(@(CR, LF), 2);
	call write(@('Argument '), 9);
	call binary_to_hex_ascii(@estate287.argument, 1);
	call write(@(CR, LF), 2);
	call write(@('Arg1 '), 5);
	call binary_to_hex_ascii(@estate287.arg1, 10);
	call write(@(' full '), 6);
	call binary_to_hex_ascii(@estate287.arg1_full, 1);
	call write(@(CR, LF), 2);
	call write(@('Arg2 '), 5);
	call binary_to_hex_ascii(@estate287.arg2, 10);
	call write(@(' full '), 6);
	call binary_to_hex_ascii(@estate287.arg2_full, 1);
	call write(@(CR, LF), 2);
	call write(@('Result '), 7);
	call binary_to_hex_ascii(@estate287.result, 1);
	call write(@(CR, LF), 2);
	call write(@('Res1 '), 5);
	call binary_to_hex_ascii(@estate287.res1, 10);
	call write(@(' full '), 6);
	call binary_to_hex_ascii(@estate287.res1_full, 1);
	call write(@(CR, LF), 2);
	call write(@('Res2 '), 5);
	call binary_to_hex_ascii(@estate287.res2, 10);
	call write(@(' full '), 6);
	call binary_to_hex_ascii(@estate287.res2_full, 1);
	call write(@(CR, LF), 2);
	call write(@('Format '), 7);
	call binary_to_hex_ascii(@estate287.format, 1);
	call write(@(' Reg '), 5);
	call binary_to_hex_ascii(@estate287.register, 1);
	call write(@(CR, LF), 2);
	call write(@('CW '), 3);
	call binary_to_hex_ascii(@estate287.control_word, 2);
	call write(@(' SW '), 4);
	call binary_to_hex_ascii(@estate287.status_word, 2);
	call write(@(' Tag '), 5);
	call binary_to_hex_ascii(@estate287.tag_word, 2);
	call write(@(CR, LF), 2);
	call write(@('Inst. '), 6);
	call binary_to_hex_ascii(@estate287.error_pointers(1), 2);
	call write(@(':'), 1);
	call binary_to_hex_ascii(@estate287.error_pointers(0), 2);
	call write(@('   Data '), 8);
	call binary_to_hex_ascii(@estate287.error_pointers(3), 2);
	call write(@(':'), 1);
	call binary_to_hex_ascii(@estate287.error_pointers(2), 2);
	call write(@(CR, LF), 2);
	call write(@('Stack:',CR,LF), 8);
	do i=0 to 35 by 5;
		call write(@('   '), 1);
		call binary_to_hex_ascii(@estate287.stack_287(i), 10);
		call write(@(CR, LF), 2);
	end;
	call write(@(CR, LF), 2);
	
END;

$endif

normalize: PROCEDURE(temp_p);
	DECLARE
		temp_p	POINTER;
	DECLARE
		temp	BASED temp_p (5) WORD;

	IF (temp(3) = 0) AND
	   (temp(2) = 0) AND
	   (temp(1) = 0) AND
	   (temp(0) = 0) THEN DO;
		temp(4) = 0;
	END;
	ELSE
	DO WHILE (temp(3) AND 08000H) = 0;
		temp(3) = SHL(temp(3), 1) OR SHR(temp(2), 15);
		temp(2) = SHL(temp(2), 1) OR SHR(temp(1), 15);
		temp(1) = SHL(temp(1), 1) OR SHR(temp(0), 15);
		temp(0) = SHL(temp(0), 1);
		temp(4) = temp(4) - 1;
	END;
END;

TRAP_HANDLER: PROCEDURE INTERRUPT PUBLIC;

	DECLARE ESTATE287 STRUCTURE (
		OPERATION WORD,
		ARGUMENT BYTE,
		ARG1(5) WORD, ARG1_FULL BYTE,
		ARG2(5) WORD, ARG2_FULL BYTE,
		RESULT BYTE,
		RES1(5) WORD, RES1_FULL BYTE,
		RES2(5) WORD, RES2_FULL BYTE,
 		FORMAT BYTE,
		REGISTER BYTE,
		CONTROL_WORD WORD,
	 	STATUS_WORD WORD,
		TAG_WORD WORD,
		ERROR_POINTERS(4) WORD,
		STACK_287(40) WORD);
	DECLARE ERRORS287 BYTE;
	DECLARE CONTROL287 WORD;

	ERRORS287 = get$real$error;
	CALL DECODE(@ESTATE287, ERRORS287);
	CONTROL287 = ESTATE287.CONTROL_WORD;

$if DEBUG
	CALL SYSLOG(@(0,0), @('FP Interrupt'), 12);
	CALL DUMP287(@ESTATE287, ERRORS287);
$endif

	IF RETRY_FLAG:=NORMAL(@ESTATE287, ERRORS287) THEN DO;
$if DEBUG
		CALL SYSLOG(@(0,0), @('Normalize operand'), 17);
		CALL DUMP287(@ESTATE287, ERRORS287);
$endif
	END;
	ELSE IF (ESTATE287.OPERATION = 0FH) AND
		((ERRORS287 AND 01H) = 1) THEN DO;
$if DEBUG
		CALL SYSLOG(@(0,0), @('Normalize store'), 15);
$endif
		IF (ESTATE287.ARG1(0) AND 07FFFH) <> 07FFFH THEN DO;
			CALL NORMALIZE(@ESTATE287.ARG1);
			ESTATE287.RES1(0) = ESTATE287.ARG1(0);
			ESTATE287.RES1(1) = ESTATE287.ARG1(1);
			ESTATE287.RES1(2) = ESTATE287.ARG1(2);
			ESTATE287.RES1(3) = ESTATE287.ARG1(3);
			ESTATE287.RES1(4) = ESTATE287.ARG1(4);
		END;
$if DEBUG
		CALL DUMP287(@ESTATE287, ERRORS287);
$endif
		RETRY_FLAG = TRUE;
	END;

	CALL ENCODE(@ESTATE287, ERRORS287, CONTROL287, RETRY_FLAG);
END TRAP_HANDLER;

END HANDLER_MODULE;

______________________________________________________________________

OBJS =		node.obj dhex.obj eh287n.obj

node:		$(OBJS)
		bnd286 -fl -oj node -rs '"0 to 6"' -ss '"stack(+1480)"' -- \
		$(OBJS) \
		/usr/ipsc/lib/lf286r0.lib, \
		/usr/intel/lib/f286r1.lib, \
		/usr/intel/lib/f286r2.lib, \
		/usr/intel/lib/dc287.lib, \
		/usr/intel/lib/eh287.lib, \
		/usr/intel/lib/80287.lib

node.obj:	node.f
		ftn286 -db -noty -code node.f

dhex.obj:	dhex.f
		ftn286 -db -noty dhex.f

eh287n.obj:	eh287n.plm
		/usr/prp/bin/plm286 -db -code eh287n.plm

______________________________________________________________________

   test dhex
     4.823855600872397D+15  433123456789ABCD
   -1/0
   -----------------------  FFF0000000000000
   denormal
     0.000000010000000-307  000000000C1069CD
   tiny + denormal
     1.000000000100000-305  009C16C5C53145DF
   promote denormal
     9.999999984813737-306  009C16C5C46E0000
   1.0 + promoted denormal
     1.000000000000000D+00  3FF0000000000000
   sqrt(denormal)
   .......................  FFF8000000000000
   overflow
   +++++++++++++++++++++++  7FF0000000000000
   sqrt(overflow)
   .......................  FFF8000000000000
   0/0
   .......................  FFF8000000000000
   sqrt(0/0)
   .......................  FFF8000000000000