/************************************************************************************
*   The subroutines NUM, ADD, MULT, DIV enable basic arithmetic operations of numbers
*   with exponents up to ±1015. 
*   NUM separates usual floating point numbers (X) into a mantissa (X) and exponent (XEX),
*   and ADD, MULT, DIV provide the basic arithmetic operations of these representations,
*   where XM,XEX,YM,YEX are the two input numbers and XYM, XYEX the result.
*
*   These routines have worked well and correct for me in the past (using IBM PL/1 compilers),  
*   but I can't guarantee this for all compilers, especially if the code is translated
*   into a different language.
*
*   If there are any problems, please let me know and I mention the issue here.
*   (for contact details see  http://www.plasmaphysics.org.uk/feedback.htm ).
*
*                                               Thomas Smid, May 2007

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

 NUM: PROCEDURE(X,XEX) ;						
 DEFAULT RANGE(A:Z) DEC FLOAT VALUE(DEC FLOAT(16)) ;
 DCL X DEC FLOAT(16) ;
 DCL XEX DEC FIXED(15,0) ;
 LAX=LOG10(ABS(X)) ;
 XEX=FLOOR(LAX) ;
 X=10**(LAX-XEX)*SIGN(X) ;
 END NUM ;



 ADD: PROCEDURE(XM,XEX,YM,YEX,XYM,XYEX) ;
 DEFAULT RANGE(A:Z) DEC FLOAT VALUE(DEC FLOAT(16)) ;
 DCL (XM,YM,XYM) DEC FLOAT(16) ;
 DCL (XEX,YEX,XYEX) DEC FIXED(15,0) ;
 IF XM=0&YM=0 THEN DO;
 XYM=0; XYEX=0; 
 RETURN;
 END ;
 IF XM=0 THEN DO;
 XYM=YM; XYEX=YEX;
 RETURN;
 END ;
 IF YM=0 THEN DO;
 XYM=XM; XYEX=XEX;
 RETURN;
 END ;

 IF XEX>=YEX THEN DO ;
 XYEX=XEX ;
 XYM=YEX-XEX; XYME=XYM; 
 IF XYM<-60 THEN XYM=0; ELSE XYM=10**XYME ;
 XYM=XM+YM*XYM ;
 END ;
 ELSE DO ;
 XYEX=YEX;
 XYM=XEX-YEX; XYME=XYM; 
 IF XYM<-60 THEN XYM=0; ELSE XYM=10**XYME ;
 XYM=XM*XYM+YM ;
 END ;
 IF ABS(XYM)>=10 THEN DO ;
 XYM=XYM/10; XYEX=XYEX+1 ;
 END ;
 IF ABS(XYM)<1 THEN DO; 
 XYM=XYM*10; XYEX=XYEX-1;
 END ;

 END ADD ;



 MULT: PROCEDURE(XM,XEX,YM,YEX,XYM,XYEX) ;
 DEFAULT RANGE(A:Z) DEC FLOAT VALUE(DEC FLOAT(16)) ;
 DCL (XM,YM,XYM) DEC FLOAT(16) ;
 DCL (XEX,YEX,XYEX) DEC FIXED(15,0) ;
 IF XM=0|YM=0 THEN DO ; 
 XYM=0; XYEX=0; 
 RETURN;
 END ;
 XYM=XM*YM ;
 XYEX=XEX+YEX ; 
 IF ABS(XYM)>=10 THEN DO ;
 XYM=XYM/10; XYEX=XYEX+1 ;
 END ;
 END MULT ;


 DIV: PROCEDURE(XM,XEX,YM,YEX,XYM,XYEX) ;
 DEFAULT RANGE(A:Z) DEC FLOAT VALUE(DEC FLOAT(16)) ;
 DCL (XM,YM,XYM) DEC FLOAT(16) ;
 DCL (XEX,YEX,XYEX) DEC FIXED(15,0) ;
 IF XM=0 THEN DO;
 XYM=0; XYEX=0; 
 RETURN;
 END ;
 XYM=XM/YM ;
 XYEX=XEX-YEX ; 
 IF ABS(XYM)<1 THEN DO ;
 XYM=XYM*10; XYEX=XYEX-1 ;
 END ;
 END DIV ;