OPTIONS(/e/-q/-i/-a/-d);
EXTERNAL INTEGER PROCEDURE intrea;

!    Procedure LSUM calculates the sum of the long real array LA,
     LA[1] + ... + LA[N] (almost) preserving precision.  Thus if
     LA[1] = E30, LA[2] = 1 and LA[3] = -E30 the sum will be
     correctly calculated to 1.
     EXTERNAL Procedure required: INTEGER PROCEDURE intrea;
LONG REAL PROCEDURE lsum(la,n);
LONG REAL ARRAY la;   INTEGER n;

BEGIN   INTEGER i,icell,high,low;   LONG REAL s;
    LONG REAL ARRAY cell[0:63];

    OPTIONS(/A); COMMENT START ARRAY BOUND CHECKING;
    la[n]:= la[n];

    icell:= intrea(la[1])//8R4000000000;;
    OPTIONS(/-A); COMMENT NO ARRAY BOUND CHECKING;
    IF icell < 0 THEN icell:= -icell;
    high:= low:= icell;
    cell[icell]:= la[1];
    FOR i:= 2 STEP 1 UNTIL n DO
    BEGIN   icell:= intrea(la[i])//8R4000000000;
	! Store LA[I] in CELL [ Abs(exp-part of LA[I]) ];
	! Thus E+MM and -E+MM will be accumulated in the
	! same cell;
    	IF icell < 0 THEN icell:= -icell;
	IF icell > high THEN high:= icell ELSE
	IF icell < low THEN low:= icell;
	cell[icell]:= cell[icell] + la[i];
    END loop;


    ! Sum starting with small values;
    FOR i:= low STEP 1 UNTIL high DO
    s:= s + cell[i];
    lsum:= s

END of lsum;