From moler%isc.intel.com@CSNET-RELAY.ARPA Sat Apr 19 08:18:21 1986 Received: from CSNET-RELAY.ARPA (csnet-relay.arpa.ARPA) by anl-mcs.ARPA (4.12/4.9) id AA16445; Sat, 19 Apr 86 08:18:01 cst Message-Id: <8604191418.AA16445@anl-mcs.ARPA> Received: from isc.intel.com by csnet-relay.csnet id af07655; 19 Apr 86 9:14 EST Received: by isc.intel.com (4.12/2.0.iSC) id AA16284; Fri, 18 Apr 86 13:45:14 pst Date: Fri, 18 Apr 86 13:45:14 pst From: Cleve Moler Posted-Date: Fri, 18 Apr 86 13:45:14 pst To: dongarra@anl-mcs.arpa, na.laub@su-score.arpa, mth%ornl-msr.uucp@CSNET-RELAY.ARPA Subject: iPSC Floating Point Exception Status: RO 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