MODULE VARIABLE_PRECISION_ARITHMETIC
! (c) copyright John Lawrence Schonfelder (1999)
! version 1.1 produced Jan-1999
!-----------------------------------------------------------------------------!
! Provides facilities for representing and manipulating real values of        !
! arbitrary precision and extreme range.                                      !
! This version is produced specifically for a PC                              !                                            
! Default INTEGER 32bit, HUGE=2147483647, RANGE=9                             !
!-----------------------------------------------------------------------------!

PRIVATE

INTEGER,PARAMETER ::  maxdecdig=9  ! max number of decimal digits that 
                                   ! can be represented in a default integer
                                   ! maxdecdig=RANGE(0)
INTEGER,PARAMETER :: radd=4  ! number of decimal digits representable in the 
                             ! number radix, radd=maxdecdig/2
INTEGER,PARAMETER :: rad=10000 ! number system radix chosen so that rad**2
                               ! is representable in a default integer
                               ! rad=10**radd
INTEGER,PARAMETER :: maxd=HUGE(0)/rad ! used in normalisation

INTEGER :: ndig=(30+radd-1)/radd ! controls the current accuracy
                                 ! initially set to provide at least 30D

TYPE NUMBER
  INTEGER         :: exp=0     ! holds the base rad exponent
  INTEGER,POINTER :: sig(:)=>NULL()  ! holds the significand
ENDTYPE NUMBER
!-----------------------------------------------------------------------------!
! the value of a NUMBER is given by                                           !
!       rad(**exp)*SUM(sig(i)*rad**(-i))  i=0:ndig                            !
! where ABS(sig(i))<rad, sig(0) /= 0 and either sig(i) = 0 exactly            !
! for i > ndig, or the significand has been truncated to retain only          !
! this number of digits                                                       !
! exact zero is represented with sig(0)=0                                     !
!-----------------------------------------------------------------------------!
! N.B. operations that can cause overflow like * and / flag an error if an    !
!       exponent greater than rad is produced. Numbers larger than            !
!       rad**(rad), 10**(radd*rad), are considered too large to be of         !
!       use even though very much larger numbers are representable            ! 
!-----------------------------------------------------------------------------! 

INTEGER,PARAMETER :: ichzero = ICHAR("0")
INTEGER,PARAMETER :: nbits = BIT_SIZE(0)-1
INTEGER,PARAMETER :: bigint0=HUGE(0)/(rad*rad)      ! digits of a NUMBER
INTEGER,PARAMETER :: remdigs=HUGE(0)-bigint0*rad**2 ! with value HUGE(0)
INTEGER,PARAMETER :: bigint1=remdigs/rad            ! for use in tests for
INTEGER,PARAMETER :: bigint2=remdigs-bigint1*rad    ! integer overflow
  ! INT(HUGE(0)) = NUMBER(2,(/bigint0,bigint1,bigint2/))

!- PRECISION set -------------------------------------------------------------!
INTERFACE PRECISION
 MODULE PROCEDURE set_num_prec ! set or inquire on current working precision
ENDINTERFACE  PRECISION

!- Conversion procedure interfaces -------------------------------------------!
INTERFACE NUM
 MODULE PROCEDURE int_to_num, &    ! default integer to NUMBER, exact
                 char_to_num       ! value represented as literal to NUMBER
                                   ! conversion exact
ENDINTERFACE  NUM

INTERFACE INT
 MODULE PROCEDURE num_to_int ! truncate NUMBER to default integer
ENDINTERFACE  INT

INTERFACE CHAR
 MODULE PROCEDURE num_to_char ! NUMBER to character literal, exact
ENDINTERFACE  CHAR

!- Assignment operator interface ---------------------------------------------!
INTERFACE ASSIGNMENT(=)
  MODULE PROCEDURE num_ass_num ! NUMBER to NUMBER assignment
ENDINTERFACE  ASSIGNMENT(=)

!- Formatted conversion procedures -------------------------------------------!
INTERFACE EFCHAR
 MODULE PROCEDURE num_to_ef  ! NUMBER chopped to E format character literal
ENDINTERFACE  EFCHAR

!- arithmetic operator interfaces --------------------------------------------!
INTERFACE OPERATOR(+)
 MODULE PROCEDURE num_plus_num,plus_num
ENDINTERFACE  OPERATOR(+)

INTERFACE OPERATOR(-)
 MODULE PROCEDURE num_minus_num,minus_num
ENDINTERFACE  OPERATOR(-)

INTERFACE OPERATOR(*)
 MODULE PROCEDURE num_mult_num
ENDINTERFACE  OPERATOR(*)

INTERFACE OPERATOR(/)
 MODULE PROCEDURE num_div_num
ENDINTERFACE  OPERATOR(/)

INTERFACE OPERATOR(**) ! raises a NUMBER to an INTEGER Power
 MODULE PROCEDURE num_power_int
ENDINTERFACE OPERATOR(**)

!- logical comparison operator interfaces ------------------------------------!
INTERFACE OPERATOR(<)
 MODULE PROCEDURE num_lt_num
ENDINTERFACE  OPERATOR(<)

INTERFACE OPERATOR(<=)
 MODULE PROCEDURE num_le_num
ENDINTERFACE  OPERATOR(<=)

INTERFACE OPERATOR(>=)
 MODULE PROCEDURE num_ge_num
ENDINTERFACE  OPERATOR(>=)

INTERFACE OPERATOR(>)
 MODULE PROCEDURE num_gt_num
ENDINTERFACE  OPERATOR(>)

INTERFACE OPERATOR(==)
 MODULE PROCEDURE num_eq_num
ENDINTERFACE  OPERATOR(==)

INTERFACE OPERATOR(/=)
 MODULE PROCEDURE num_ne_num
ENDINTERFACE  OPERATOR(/=)

!--Rounding procedures--------------------------------------------------------!
INTERFACE TRUNC
 MODULE PROCEDURE num_trunc ! truncates a NUMBER towards zero to a secified 
                            ! decimal precision
ENDINTERFACE  TRUNC

!-Overload of intrinsic procedures--------------------------------------------!
INTERFACE ABS
 MODULE PROCEDURE num_abs
ENDINTERFACE ABS

INTERFACE SIGN
 MODULE PROCEDURE num_sign
ENDINTERFACE SIGN

!-----------------------------------------------------------------------------!
PUBLIC :: maxdecdig,rad,radd,ndig, &
          NUMBER,PRECISION,NUM,INT,CHAR,ASSIGNMENT(=),EFCHAR,&
          OPERATOR(+),OPERATOR(-),OPERATOR(*),OPERATOR(/), &
          OPERATOR(<),OPERATOR(<=),OPERATOR(>=),OPERATOR(>), &
          OPERATOR(==),OPERATOR(/=),OPERATOR(**), &
          TRUNC,ABS,SIGN
!-----------------------------------------------------------------------------!

CONTAINS

!- private auxilliary procedures ---------------------------------------------!
SUBROUTINE fault(message)
 CHARACTER(LEN=*) :: message ! intended to identify the procedure and
                             ! nature of the error
 WRITE(*,*) "**** ERROR in EXTENDED_PRECISION MODULE *****"
 WRITE(*,*) "**** ",message
 STOP
ENDSUBROUTINE fault

FUNCTION IDIG(cval)
  INTEGER :: IDIG
  CHARACTER(*)::cval
  ! converts digits contained in cval to equivalent decimal integer
  ! greater than zero, returns a negative value if a nonvalid char included
  INTEGER :: dwk
  INTEGER,PARAMETER :: chzero = ICHAR("0")
  IDIG=0
  DO i = 1,LEN(cval)
    dwk = ICHAR(cval(i:i))-chzero
    IF(dwk<0 .OR. dwk>9)THEN
      IDIG=-1; RETURN
    ELSE
      IDIG=10*IDIG + dwk
    ENDIF
  ENDDO
ENDFUNCTION IDIG

FUNCTION DIGS(intval,ndgc)
 INTEGER, PARAMETER :: iczero=ICHAR("0")
 INTEGER :: intval,ndgc
 CHARACTER(LEN=ndgc) :: DIGS
 ! returns the digits of intval as characters right justified
 ! with least significant digit in DIGS(ndgc:ndgc)
 INTEGER :: wka,wkb
 wka=ABS(intval)
 DO i=ndgc,1,-1
   wkb = wka
   wka = wkb/10
   wkb = wkb - wka*10
   DIGS(i:i)=CHAR(wkb+iczero)
 ENDDO
ENDFUNCTION DIGS

!-PRECISION set -------------------------------------------------------------!
FUNCTION set_num_prec(dec_acc) ! generic PRECISION
 INTEGER,INTENT(IN) :: dec_acc
 INTEGER :: set_num_prec
 ! sets or inquires on current working precision
 ! dec_acc specifies the required decimal accuracy to be used
 ! returns the precision actually set
 ! if dec_acc <= 0 no change is made in the current accuracy setting
 IF(dec_acc > 0)THEN
   ndig =CEILING(REAL(dec_acc)/radd) 
 ENDIF
 set_num_prec = ndig*radd
ENDFUNCTION set_num_prec

!- Conversion procedures -----------------------------------------------------!
FUNCTION int_to_num(val) ! generic NUM, default INTEGER to NUMBER
 INTEGER,INTENT(IN) :: val
 type(NUMBER) :: int_to_num
 ! converts from a default integer value to a number representation
 ! of exactly the same value
 INTEGER :: iwk
 IF(ABS(val) < rad)THEN
   ALLOCATE(int_to_num%sig(0:0))
   int_to_num%exp=0
   int_to_num%sig(0)=val
 ELSEIF(ABS(val) < rad*rad)THEN
   ALLOCATE(int_to_num%sig(0:1))
   int_to_num%exp=1
   int_to_num%sig(0)=val/rad
   int_to_num%sig(1)=val - int_to_num%sig(0)*rad
 ELSE
   ALLOCATE(int_to_num%sig(0:2))
   int_to_num%exp=2
   int_to_num%sig(0)=val/(rad*rad)
   iwk=val - int_to_num%sig(0)*rad*rad
   int_to_num%sig(1)=iwk/rad
   int_to_num%sig(2)=iwk - int_to_num%sig(1)*rad
 ENDIF
ENDFUNCTION int_to_num

FUNCTION char_to_num(val) ! generic NUM, value represented as literal to NUMBER
 CHARACTER(LEN=*),INTENT(IN) :: val
 type(NUMBER) :: char_to_num
 ! a value represented as a character literal, possibly with leading or trailing 
 ! blanks, in the form
 ! sdddd.ddddEsdddd
 ! where s is the sign + - or nothing
 !       dddd is a string of one or more decimal digits
 ! is converted to the equivalent NUMBER
 INTEGER :: lv,&      ! length of input string less leading and trailing blanks
            sgn,&     ! sign of number
            lsi,  &   ! start and end of integer string
            lef,&     ! start and end of fraction string
            sex,&     ! sign of exponent
            lse,    & ! start and end of exponent string
            we,  &     ! working result exponent and shift to change base.
            iwk,jwk,dwk, & ! working integers
            inz,lnz,idp ! first and last nonzero digits, decimal point 
 CHARACTER(LEN=LEN(val)) :: wch
 wch=ADJUSTL(val); lv=LEN_TRIM(wch)
 IF(lv<5)call fault("invalid form of value, too few characters, in NUM("//val//")")
 ! sort out decimal exponent, scan backwards to find E 
 lse=SCAN(wch(1:lv),set="Ee",back=.TRUE.)
 IF(lse<4 .OR. lv-lse>maxdecdig+2) &
      call fault("invalid exponent in NUM("//val//")")
 ! fix the sign of the exponent and start of exponent string
 lef=lse-1 ! end of fraction
 lse=lse+1 ! sign or start of exponent
 IF(wch(lse:lse) == "-")THEN
   sex = -1; lse = lse+1
 ELSEIF(wch(lse:lse) == "+")THEN
   sex = +1; lse = lse+1
 ELSEIF(wch(lse:lse) >= "0" .AND. wch(lse:lse) <= "9")THEN
   sex = +1
 ELSE
 call fault("invalid sign of exponent in NUM("//val//")")
 ENDIF 
 we = IDIG(wch(lse:lv)) 
 IF(we < 0)call fault("invalid character in exponent value in NUM("//val//")")
 we=we*sex ! we will now contain the decimal exponent of val
 ! fix the sign of the number and start of integer string
 IF(wch(1:1) == "-")THEN
   sgn = -1; lsi = 2
 ELSEIF(wch(1:1) == "+")THEN
   sgn = +1; lsi = 2
 ELSEIF(wch(1:1) >= "0" .AND. wch(1:1) <= "9")THEN
   sgn = +1; lsi = 1
 ELSE
   call fault("invalid sign of value in NUM("//val//")")
 ENDIF
! scan forward to find decimal point
 idp=SCAN(wch(lsi:lef),set=".")
 IF(idp == 0)call fault("invalid form of value, no decimal point, in NUM("//val//")")
 idp=idp+lsi-1
 ! scan forward to find first non-zero digit 
 inz=SCAN(wch(lsi:lef),set="123456789")
 IF(inz == 0)THEN ! either format invalid or number is all zero
   iwk=VERIFY(wch(lsi:idp-1)//wch(idp+1:lef),set="0")
   IF(iwk /= 0) call fault("invalid form of value in NUM("//val//")")
   ! result is zero
   char_to_num%exp=0
   ALLOCATE(char_to_num%sig(0:0))
   char_to_num%sig(0)=0
   RETURN
 ENDIF
 ! value is nonzero
 inz=inz+lsi-1
 ! scan backwards to find last nonzero digit
 lnz=SCAN(wch(lsi:lef),set="123456789",back=.TRUE.)+lsi-1
 IF(inz<idp)THEN ! nonzero integer part
   we=we+idp-inz  !adjust exponent so first digit of fraction nonzero
   IF(lnz>idp)THEN ! nonzero fraction copy significant digits into low end of wch
     wch(1:idp-inz)=wch(inz:idp-1)
     lv=lnz-inz
     wch(idp-inz+1:lv)=wch(idp+1:lnz)
   ELSE ! zero fraction copy significant digits into low end of wch
     lv=lnz-inz+1
     wch(1:lv)=wch(inz:lnz)
   ENDIF
 ELSE !only fraction nonzero
   we=we+idp-inz+1
   lv=lnz-inz+1
   wch(1:lv)=wch(inz:lnz)
 ENDIF 
 wch(lv+1:)="000000" !must be at least radd-1 zeros
 ! we is the normalised decimal exponent and wch(1:lv) contains
 ! the significant digits as a fraction
 iwk=we/radd
 jwk=we - iwk*radd
 iwk=iwk-1
 IF(jwk > 0)THEN
   iwk=iwk+1; jwk=jwk-radd
 ENDIF
 IF(ABS(iwk)>rad) call fault("exponent too large in NUM("//val//")")
 we=CEILING(REAL(lv-radd-jwk)/radd)
 ALLOCATE(char_to_num%sig(0:we))
 char_to_num%exp=iwk
 dwk = IDIG(wch(1:radd+jwk))
 IF(dwk < 0)call fault("invalid character in mantissa value in NUM("//val//")")
 char_to_num%sig(0)=sgn*dwk
 iwk=1
 DO i=radd+jwk+1,lv,radd
   dwk=IDIG(wch(i:i+radd-1)) 
   IF(dwk < 0)call fault("invalid character in mantissa value in NUM("//val//")")
   char_to_num%sig(iwk)=dwk*sgn
   iwk=iwk+1
 ENDDO
ENDFUNCTION char_to_num

FUNCTION num_to_int(val) ! generic INT
  type(NUMBER),INTENT(IN) :: val
  INTEGER :: num_to_int
  ! truncate a NUMBER value to a default integer
  IF(val%exp>2) call fault("value too large for conversion in INT")
  IF(val%exp<0)THEN
    num_to_int=0; RETURN
  ELSEIF(val%exp==0)THEN
    num_to_int=val%sig(0); RETURN
  ELSEIF(val%exp==1)THEN
    num_to_int=val%sig(0)*rad + val%sig(1); RETURN
  ELSEIF(val%exp==2 .AND. ABS(val%sig(0))<bigint0)THEN
    num_to_int=(val%sig(0)*rad + val%sig(1))*rad + val%sig(2); RETURN
  ELSEIF(ABS(val%sig(0))==bigint0 .AND. ABS(val%sig(1))<bigint1)THEN
    num_to_int=(val%sig(0)*rad + val%sig(1))*rad + val%sig(2); RETURN
  ELSEIF(ABS(val%sig(0))==bigint0 .AND. ABS(val%sig(1))==bigint1 &
                        .AND.ABS(val%sig(2))<=bigint2)THEN
    num_to_int=(val%sig(0)*rad + val%sig(1))*rad + val%sig(2); RETURN
  ELSE
    call fault("value too large for conversion in INT")
  ENDIF
ENDFUNCTION num_to_int

FUNCTION num_to_char(val) ! generic CHAR
  type(NUMBER) :: val
  CHARACTER(LEN=radd*(SIZE(val%sig)+2)+7) :: num_to_char
  ! converts a NUMBER value to a character string representing the 
  ! value exactly in the format ^^^sd.dddddEsdddd 
  ! where ^^^ is zero or more blancs
  !       s is + or - 
  !       ddd are the digits of the value
  ! The value is normalised so that there is exactly one digit before
  ! the decimal point
  CHARACTER(LEN=radd) :: fsig 
  CHARACTER,PARAMETER :: plus="+",minus="-"
  INTEGER :: wka,wkb,dex,nres,ndi
  IF(val%sig(0)==0)THEN  ! value zero
    num_to_char = " +0.0E+0"
    num_to_char = ADJUSTR(num_to_char)
    RETURN
  ELSE ! value non-zero
    num_to_char=" " !fill result with blancs
    nres = LEN(num_to_char)
    ! work out number of digits before point in mantissa
    wka=ABS(val%sig(0))
    DO i=1,radd
       wka=wka/10
       IF(wka==0)EXIT
    ENDDO
    ndi=i
    ! work out the decimal equivalent of the normalised exponent
    dex=val%exp*radd+ndi-1 
    num_to_char(nres-2*radd+1:nres)=DIGS(dex,2*radd)
    wkb=SCAN(num_to_char(nres-2*radd+1:nres),set="123456789")
    IF(wkb==0)THEN !exponent zero
      i=nres 
    ELSE ! exponent nonzero
      i=nres-2*radd+wkb
    ENDIF
    IF(dex >= 0)THEN
      num_to_char(i-2:i-1)="E+"
    ELSE
      num_to_char(i-2:i-1)="E-"
    ENDIF
    nres=i-3
    wka=UBOUND(val%sig,1)
    DO i=nres,1,-radd ! add digits to fraction from the right
      IF(wka==0)EXIT
      num_to_char(i-radd+1:i)=DIGS(val%sig(wka),radd)
      wka=wka-1
    ENDDO
    fsig=DIGS(val%sig(0),radd)
    num_to_char(i-ndi+2:i)=fsig(radd-ndi+2:radd)
    i=i-ndi+1
    num_to_char(i:i)="."; i=i-1
    num_to_char(i:i)=fsig(radd-ndi+1:radd-ndi+1); i=i-1
    IF(val%sig(0)>0)THEN
      num_to_char(i:i)="+"
    ELSE
      num_to_char(i:i)="-"
    ENDIF
  ENDIF
ENDFUNCTION num_to_char

FUNCTION num_to_ef(val,w,d) ! formatted conversion, generic EFCHAR
  type(NUMBER),INTENT(IN) :: val
  INTEGER,INTENT(IN) :: w,d
  CHARACTER(LEN=w) :: num_to_ef
  ! converts to value character string with format 1PEw.d
  CHARACTER(LEN=(SIZE(val%sig)+2)*radd+7) :: chval
  INTEGER :: pp,ep,nc, & ! point, E position and nonblank representation length
             ev, &       ! 
             i,j        ! working variables
  chval = ADJUSTR(CHAR(val)) ! places the character representation of the value
  nc = LEN(chval)            ! in chval(1:nc) rt justified
  ! scan forward to locate point
  pp = SCAN(chval,".") 
  ! scan backwards to locate E
  ep = SCAN(chval,"E",back=.TRUE.)
    IF(nc-ep+d+7 >= w )THEN ! number overflows format
      DO i=1,w
        num_to_ef(i:i) = "*"
      ENDDO
      RETURN
    ENDIF
  i=w-nc+ep
  num_to_ef(i:w)=chval(ep:nc);i=i-1
  DO j=1,d
    IF(pp+j < ep)THEN
      num_to_ef(i-d+j:i-d+j)=chval(pp+j:pp+j)
    ELSE
      num_to_ef(i-d+j:i-d+j)="0"
    ENDIF
  ENDDO 
  i=i-d
  num_to_ef(i-2:i)=chval(pp-2:pp)
  num_to_ef(1:i-3)=" "
ENDFUNCTION num_to_ef

!- Assignment procedure ------------------------------------------------------!
SUBROUTINE num_ass_num(var,expr)
  type(NUMBER),INTENT(IN) :: expr
  type(NUMBER),INTENT(INOUT) :: var
  ! assigns a number value expression to a variable making any change in 
  ! accuracy required by the current precision
  INTEGER :: wex,n ! working registers
  INTEGER :: wdig(0:SIZE(expr%sig)-1)
  wex = expr%exp
  wdig = expr%sig
  n = MIN(ndig,UBOUND(wdig,1))
  IF(ASSOCIATED(var%sig)) DEALLOCATE(var%sig)
  ALLOCATE(var%sig(0:n))
  var%exp = wex
  var%sig = wdig(0:n)
ENDSUBROUTINE num_ass_num

!- arithmetic operator procedures --------------------------------------------!
!- Operator (+) --------------------------------------------------------------!
FUNCTION num_plus_num(l,r)  
  type(NUMBER),INTENT(IN) :: l,r
  type(NUMBER) :: num_plus_num
  ! performs an algebraic addition of the two operands and returns a 
  ! result correctly truncated to the current precision
  INTEGER::wkr(-1:ndig+1), & ! working mantissa register
           wke, &            ! working exponent register
           it,is,ir,&        ! working indices
           sgn,srad          ! sign and signed radix of result
  wkr=0 ! clear working mantissa
  IF(l%exp >= r%exp)THEN ! l has larger abs value
    wke=l%exp; it=wke-r%exp ! shift to align digits
    ir=MIN(ndig,UBOUND(l%sig,1))
    IF(it>ndig)THEN !smaller value too small to effect result
      num_plus_num%exp=wke
      ALLOCATE(num_plus_num%sig(0:ir))
      num_plus_num%sig=l%sig(0:ir)
      RETURN
    ENDIF
    ! copy l%sig to register
    wkr(0:ir)=l%sig(0:ir)
    ir=MIN(ndig+1,it+UBOUND(r%sig,1))
    ! add relevant digits shifted as necessary to register
    wkr(it:ir)=wkr(it:ir)+r%sig(0:ir-it)
  ELSE ! r has larger abs value
    wke=r%exp; it=wke-l%exp ! shift to align digits
    ir=MIN(ndig,UBOUND(r%sig,1))
    IF(it>ndig)THEN !smaller value too small to effect result
      num_plus_num%exp=wke
      ALLOCATE(num_plus_num%sig(0:ir))
      num_plus_num%sig=r%sig(0:ir)
      RETURN
    ENDIF
    ! copy r to register
    wkr(0:ir)=r%sig(0:ir)
    ir=MIN(ndig+1,it+UBOUND(l%sig,1))
    wkr(it:ir)=wkr(it:ir)+l%sig(0:ir-it) 
  ENDIF
  ! wke and wkr(0:ndig+1) contain unnormalised result
  ! sign of result is same as sign of first non-zero digit
  ir=ndig+1
  DO i=0,ir
    IF(wkr(i)>0)THEN
      sgn=+1;srad=rad;EXIT
    ELSEIF(wkr(i)<0)THEN
      sgn=-1;srad=-rad;EXIT
    ENDIF
    ! digit is zero cycle loop
  ENDDO
  is=i ! set index of most significant non-zero digit
  IF(is>ir)THEN ! result is zero
    num_plus_num%exp=0
    ALLOCATE(num_plus_num%sig(0:0))
    num_plus_num%sig=0
    RETURN
  ENDIF
  ! scan backwards to normalise digits
  DO i=ir,is,-1
    IF(wkr(i)*sgn<0)THEN     ! sign of digit wrong
      wkr(i)=wkr(i)+srad     ! carry radix from higher digit
      wkr(i-1)=wkr(i-1)-sgn  ! to switch sign
    ENDIF
    IF(ABS(wkr(i))>=rad)THEN  ! digit larger than radix
      wkr(i)=wkr(i)-srad     ! reduce by radix and carry
      wkr(i-1)=wkr(i-1)+sgn  ! to higher digit
    ENDIF
    IF(wkr(i)==0.AND.i==ir) ir=i-1 ! least significant digit zero
                                   ! reduce length of number
  ENDDO
  IF(wkr(is-1)/=0)is=is-1
 IF(ABS(wke-is)>rad) call fault("exponent too large in + ")
  num_plus_num%exp=wke-is
  ALLOCATE(num_plus_num%sig(0:ir-is))
  num_plus_num%sig=wkr(is:ir)
ENDFUNCTION num_plus_num

FUNCTION plus_num(r)
  type(NUMBER),INTENT(IN) :: r
  type(NUMBER) :: plus_num
  ALLOCATE(plus_num%sig(0:UBOUND(r%sig,1)))
  plus_num%exp = r%exp
  plus_num%sig = r%sig
ENDFUNCTION plus_num

!- Operator (-) --------------------------------------------------------------!
FUNCTION num_minus_num(l,r)  
  type(NUMBER),INTENT(IN) :: l,r
  type(NUMBER) :: num_minus_num
  ! performs an algebraic subtraction of the two operands and returns a 
  ! result correctly truncated to the current precision
  INTEGER::wkr(-1:ndig+1), & ! working mantissa register
           wke, &            ! working exponent register
           it,is,ir,&        ! working indices
           sgn,srad          ! sign and signed radix of result
  wkr=0 ! clear working mantissa
  IF(l%exp >= r%exp)THEN ! l has larger abs value
    wke=l%exp; it=wke-r%exp ! shift to align digits
    ir=MIN(ndig,UBOUND(l%sig,1))
    IF(it>ndig+1)THEN !smaller value too small to effect result
      num_minus_num%exp=wke
      ALLOCATE(num_minus_num%sig(0:ir))
      num_minus_num%sig=l%sig(0:ir)
      RETURN
    ENDIF
    ! copy l to register
    wkr(0:ir)=l%sig(0:ir)
    ir=MIN(ndig+1,it+UBOUND(r%sig,1))
    wkr(it:ir)=wkr(it:ir)-r%sig(0:ir-it)
  ELSE ! r has larger abs value
    wke=r%exp; it=wke-l%exp ! shift to align digits
    ir=MIN(ndig,UBOUND(r%sig,1))
    IF(it>ndig)THEN !smaller value too small to effect result
      num_minus_num%exp=wke
      ALLOCATE(num_minus_num%sig(0:ir))
      num_minus_num%sig=-r%sig(0:ir)
      RETURN
    ENDIF
    ! copy r to register
    wkr(0:ir)=-r%sig(0:ir)
    ir=MIN(ndig+1,it+UBOUND(l%sig,1))
    wkr(it:ir)=wkr(it:ir)+l%sig(0:ir-it) 
  ENDIF
  ! wke and wkr(0:ir) contain unnormalised result
  ! sign of result is same as sign of first non-zero digit
  ir=ndig+1
  DO i=0,ir
    IF(wkr(i)>0)THEN
      sgn=+1;srad=rad;EXIT
    ELSEIF(wkr(i)<0)THEN
      sgn=-1;srad=-rad;EXIT
    ENDIF
    ! digit is zero cycle loop
  ENDDO
  is=i ! set index of most significant non-zero digit
  IF(is>ir)THEN ! result is zero
    num_minus_num%exp=0
    ALLOCATE(num_minus_num%sig(0:0))
    num_minus_num%sig=0
    RETURN
  ENDIF
  ! scan backwards to normalise digits
  DO i=ir,is,-1
    IF(wkr(i)*sgn<0)THEN     ! sign of digit wrong
      wkr(i)=wkr(i)+srad     ! carry radix from higher digit
      wkr(i-1)=wkr(i-1)-sgn  ! to switch sign
    ENDIF
    IF(ABS(wkr(i))>=rad)THEN  ! digit larger than radix
      wkr(i)=wkr(i)-srad     ! reduce by radix and carry
      wkr(i-1)=wkr(i-1)+sgn  ! to higher digit
    ENDIF
    IF(wkr(i)==0.AND.i==ir) ir=i-1 ! least significant digit zero
                                   ! reduce length of number
  ENDDO
  IF(wkr(is-1)/=0)is=is-1
 IF(ABS(wke-is)>rad) call fault("exponent too large in - ")
  num_minus_num%exp=wke-is
  ALLOCATE(num_minus_num%sig(0:ir-is))
  num_minus_num%sig=wkr(is:ir)
ENDFUNCTION num_minus_num

FUNCTION minus_num(r)
  type(NUMBER),INTENT(IN) :: r
  type(NUMBER) :: minus_num
  ALLOCATE(minus_num%sig(0:UBOUND(r%sig,1)))
  minus_num%exp = r%exp
  minus_num%sig = -r%sig
ENDFUNCTION minus_num

!- Operator(*) ---------------------------------------------------------------!
FUNCTION num_mult_num(l,r)
  type(NUMBER),INTENT(IN) :: l,r
  type(NUMBER) :: num_mult_num
  ! performs an algebraic multiplicaltion of the two operands and returns a 
  ! result correctly truncated to the current precision
  INTEGER,DIMENSION(-1:SIZE(l%sig)+SIZE(r%sig)-2) :: pd
  INTEGER :: nl,nr, & ! upper bounds of operand arrays
             i,j,wd   ! working integers
  IF( ABS(l%exp+r%exp)>rad )THEN
    CALL fault("danger of exponent overflow/underflow in multiply")
  ENDIF 
  pd = 0 ! clear result working array
  nl = UBOUND(l%sig,1); nr = UBOUND(r%sig,1)
  DO i=nl,0,-1
    DO j=nr,0,-1
      wd = l%sig(i)*r%sig(j) + pd(i+j)
      pd(i+j) = MOD(wd,rad)  ! normalises (i+j) result digit
      pd(i+j-1) = pd(i+j-1) + wd/rad  ! adds any carry 
    ENDDO
  ENDDO
  ! pd now contains the correct result digits
  ! locate the last nonzero digit
  DO wd=UBOUND(pd,1),0,-1
     IF(pd(wd)/=0)EXIT
  ENDDO
  IF( pd(-1)==0 )THEN 
    num_mult_num%exp = l%exp + r%exp
    wd = MAX(0,MIN(ndig,wd))
    ALLOCATE(num_mult_num%sig(0:wd))
    num_mult_num%sig = pd(0:wd)
  ELSE
    num_mult_num%exp = l%exp + r%exp + 1
    wd = MIN(ndig,wd+1)
    ALLOCATE(num_mult_num%sig(0:wd))
    num_mult_num%sig = pd(-1:wd-1)
  ENDIF
ENDFUNCTION num_mult_num

!- Operator(/) ---------------------------------------------------------------!
FUNCTION num_div_num(l,r)
  type(NUMBER),INTENT(IN) :: l,r
  type(NUMBER) :: num_div_num
  ! performs an algebraic division of the two operands and returns a 
  ! result correctly truncated to the current precision
  type(NUMBER) :: xn,tm  ! working number variables
  INTEGER :: cacc, wacc   ! current accuracy & working accuracy
  REAL :: rx ! initial real approx to reciprocal of r
  IF( ABS(l%exp-r%exp-1)>rad )THEN
    CALL fault("danger of exponent overflow/underflow in divide")
  ENDIF 
  IF(r%sig(0)==0)THEN
    CALL fault("zero denominator in divide")
  ENDIF
  cacc = PRECISION(-1) ! record current working precision
  wacc = PRECISION(radd) ! set initial internal working precision
  ! calculate initial real value for 1/r
  IF(SIZE(r%sig)==1)THEN
    rx = REAL(rad)/REAL(r%sig(0))
  ELSE
    rx = REAL(rad*rad)/REAL(rad*r%sig(0)+r%sig(1)) 
  ENDIF
  xn%exp = -r%exp-1
  ALLOCATE(xn%sig(0:1))
  xn%sig(0) = INT(rx)
  xn%sig(1) = INT((rx - xn%sig(0))*rad)
  ! xn now contains a first approximation (1/r) 
  DO ! perform Newton iterations to improve the accuracy until this greater than cacc
    xn = xn*(NUM(2)-r*xn)
    IF( wacc>cacc ) EXIT
    wacc = PRECISION(2*wacc)
  ENDDO
  xn=l*xn ! result = l*(1/r)
  wacc = PRECISION(cacc)
  num_div_num = xn ! truncate to required accuracy
ENDFUNCTION num_div_num

!--OPERATOR(**)---------------------------------------------------------------!
FUNCTION num_power_int(l,r) 
  type(NUMBER),INTENT(IN) :: l
  INTEGER,INTENT(IN) :: r
  type(NUMBER) :: num_power_int
  ! calculates l**r by repeated squaring of l and multiplying into the result
  ! when the corresponding bit position in r is set
  type(NUMBER) :: wkr,msq
  INTEGER :: bi
  IF(r==0)THEN
    num_power_int=NUM(1)
    RETURN
  ELSEIF(ABS(l%exp*r)>rad)THEN
    call fault("exponent overflow in ** operation")
  ELSE
    wkr=NUM(1); bi=ABS(r); msq=l
    DO i=0,nbits
      IF(BTEST(bi,i))THEN
        wkr=wkr*msq
        bi=IBCLR(bi,i)
      ENDIF
      IF(bi==0)EXIT
      msq=msq*msq
    ENDDO
    IF(r<0)THEN
      num_power_int=NUM(1)/wkr
    ELSE
      num_power_int=wkr
    ENDIF
  ENDIF
ENDFUNCTION num_power_int

!- logical relational operators ----------------------------------------------!
FUNCTION num_lt_num(l,r) ! OPERATOR(<)
  type(NUMBER),INTENT(IN) :: l,r
  LOGICAL :: num_lt_num
  ! compares the values, if l<r algebraically the result is true
  INTEGER :: i ! working variables
  IF( l%sig(0)*r%sig(0)<=0 )THEN ! signs are different or one at least zero
    IF( l%sig(0)<0 )THEN
      num_lt_num = .TRUE.; RETURN
    ELSEIF( l%sig(0)==0 .AND. r%sig(0)>0 )THEN
      num_lt_num = .TRUE.; RETURN
    ELSE 
      num_lt_num = .FALSE.; RETURN
    ENDIF
  ELSE ! signs same and neither is zero
    IF(l%exp/=r%exp)THEN
      num_lt_num = (l%sig(0)<0.AND.l%exp>r%exp).OR. &
                   (l%sig(0)>0.AND.l%exp<r%exp); RETURN
    ENDIF
    ! signs and exponents are equal 
    DO i=0,MIN(UBOUND(l%sig,1),UBOUND(r%sig,1))
      IF(l%sig(i)/=r%sig(i))THEN
         num_lt_num = l%sig(i)<r%sig(i); RETURN
      ENDIF
    ENDDO
    IF(UBOUND(l%sig,1)>UBOUND(r%sig,1))THEN ! l has larger abs value
      num_lt_num = l%sig(0)<0; RETURN
    ELSEIF(UBOUND(l%sig,1)==UBOUND(r%sig,1))THEN ! values are equal
      num_lt_num = .FALSE.; RETURN
    ELSE ! r has larger abs value
      num_lt_num = r%sig(0)>0; RETURN
    ENDIF
  ENDIF
ENDFUNCTION num_lt_num

FUNCTION num_le_num(l,r) ! OPERATOR(<=)
  type(NUMBER),INTENT(IN) :: l,r
  LOGICAL :: num_le_num
  ! compares the values, if l<=r algebraically the result is true
  INTEGER :: i ! working variables
  IF( l%sig(0)*r%sig(0)<=0 )THEN ! signs are different or one at least zero
    IF( l%sig(0)<0 )THEN
      num_le_num = .TRUE.; RETURN
    ELSEIF( l%sig(0)==0 .AND. r%sig(0)>=0 )THEN
      num_le_num = .TRUE.; RETURN
    ELSE 
      num_le_num = .FALSE.; RETURN
    ENDIF
  ELSE ! signs same and neither is zero
    IF(l%exp/=r%exp)THEN
      num_le_num = (l%sig(0)<0.AND.l%exp>r%exp).OR. &
                   (l%sig(0)>0.AND.l%exp<r%exp); RETURN
    ENDIF
    ! signs and exponents are equal 
    DO i=0,MIN(UBOUND(l%sig,1),UBOUND(r%sig,1))
      IF(l%sig(i)/=r%sig(i))THEN
         num_le_num = l%sig(i)<r%sig(i); RETURN
      ENDIF
    ENDDO
    IF(UBOUND(l%sig,1)>UBOUND(r%sig,1))THEN ! l has larger abs value
      num_le_num = l%sig(0)<0; RETURN
    ELSEIF(UBOUND(l%sig,1)==UBOUND(r%sig,1))THEN ! values are equal
      num_le_num = .TRUE.; RETURN
    ELSE ! r has larger abs value
      num_le_num = r%sig(0)>0; RETURN
    ENDIF
  ENDIF
ENDFUNCTION num_le_num

FUNCTION num_ge_num(l,r) ! OPERATOR(>=)
  type(NUMBER),INTENT(IN) :: l,r
  LOGICAL :: num_ge_num
  ! compares the values, if l>=r algebraically the result is true
  INTEGER :: i ! working variables
  IF( l%sig(0)*r%sig(0)<=0 )THEN ! signs are different or one at least zero
    IF( l%sig(0)>0 )THEN
      num_ge_num = .TRUE.; RETURN
    ELSEIF( l%sig(0)==0 .AND. r%sig(0)<=0 )THEN
      num_ge_num = .TRUE.; RETURN
    ELSE 
      num_ge_num = .FALSE.; RETURN
    ENDIF
  ELSE ! signs same and neither is zero
    IF(l%exp/=r%exp)THEN
      num_ge_num = (l%sig(0)>0.AND.l%exp>r%exp).OR. &
                   (l%sig(0)<0.AND.l%exp<r%exp); RETURN
    ENDIF
    ! signs and exponents are equal 
    DO i=0,MIN(UBOUND(l%sig,1),UBOUND(r%sig,1))
      IF(l%sig(i)/=r%sig(i))THEN
         num_ge_num = l%sig(i)>r%sig(i); RETURN
      ENDIF
    ENDDO
    IF(UBOUND(l%sig,1)>UBOUND(r%sig,1))THEN ! l has larger abs value
      num_ge_num = l%sig(0)>0; RETURN
    ELSEIF(UBOUND(l%sig,1)==UBOUND(r%sig,1))THEN ! values are equal
      num_ge_num = .TRUE.; RETURN
    ELSE ! r has larger abs value
      num_ge_num = r%sig(0)<0; RETURN
    ENDIF
  ENDIF
ENDFUNCTION num_ge_num

FUNCTION num_gt_num(l,r) ! OPERATOR(>)
  type(NUMBER),INTENT(IN) :: l,r
  LOGICAL :: num_gt_num
  ! compares the values, if l>r algebraically the result is true
  INTEGER :: i ! working variables
  IF( l%sig(0)*r%sig(0)<=0 )THEN ! signs are different or one at least zero
    IF( l%sig(0)>0 )THEN
      num_gt_num = .TRUE.; RETURN
    ELSEIF( l%sig(0)==0 .AND. r%sig(0)<0 )THEN
      num_gt_num = .TRUE.; RETURN
    ELSE 
      num_gt_num = .FALSE.; RETURN
    ENDIF
  ELSE ! signs same and neither is zero
    IF(l%exp/=r%exp)THEN
      num_gt_num = (l%sig(0)>0.AND.l%exp>r%exp).OR. &
                   (l%sig(0)<0.AND.l%exp<r%exp); RETURN
    ENDIF
    ! signs and exponents are equal 
    DO i=0,MIN(UBOUND(l%sig,1),UBOUND(r%sig,1))
      IF(l%sig(i)/=r%sig(i))THEN
         num_gt_num = l%sig(i)>r%sig(i); RETURN
      ENDIF
    ENDDO
    IF(UBOUND(l%sig,1)>UBOUND(r%sig,1))THEN ! l has larger abs value
      num_gt_num = l%sig(0)>0; RETURN
    ELSEIF(UBOUND(l%sig,1)==UBOUND(r%sig,1))THEN ! values are equal
      num_gt_num = .FALSE.; RETURN
    ELSE ! r has larger abs value
      num_gt_num = r%sig(0)<0; RETURN
    ENDIF
  ENDIF
ENDFUNCTION num_gt_num

FUNCTION num_ne_num(l,r) ! OPERATOR(/=)
  type(NUMBER),INTENT(IN) :: l,r
  LOGICAL :: num_ne_num
  ! compares the values, if l/=r the result is true
  INTEGER :: i ! working variables
  IF( l%sig(0)*r%sig(0)<0 )THEN ! signs are different and neither zero
    num_ne_num = .TRUE.; RETURN
  ELSEIF( l%sig(0)==0 .AND. r%sig(0)==0 )THEN ! both zero
    num_ne_num = .FALSE.; RETURN
  ELSEIF(l%exp/=r%exp)THEN ! exponents differ
    num_ne_num = .TRUE.; RETURN
  ELSEIF(UBOUND(l%sig,1)/=UBOUND(r%sig,1))THEN !number of digits differ
    num_ne_num = .TRUE.; RETURN
  ELSE ! signs, exponents and number of digits are equal 
    DO i=0,UBOUND(l%sig,1)
      IF(l%sig(i)/=r%sig(i))THEN
         num_ne_num = .TRUE.; RETURN
      ENDIF
    ENDDO
    num_ne_num = .FALSE.; RETURN
  ENDIF
ENDFUNCTION num_ne_num

FUNCTION num_eq_num(l,r) ! OPERATOR(==)
  type(NUMBER),INTENT(IN) :: l,r
  LOGICAL :: num_eq_num
  ! compares the values, if l==r the result is true
  INTEGER :: i ! working variables
  IF( l%sig(0)*r%sig(0)<0 )THEN ! signs are different and neither zero
    num_eq_num = .FALSE.; RETURN
  ELSEIF( l%sig(0)==0 .AND. r%sig(0)==0 )THEN ! both zero
    num_eq_num = .TRUE.; RETURN
  ELSEIF(l%exp/=r%exp)THEN ! exponents differ
    num_eq_num = .FALSE.; RETURN
  ELSEIF(UBOUND(l%sig,1)/=UBOUND(r%sig,1))THEN !number of digits differ
    num_eq_num = .FALSE.; RETURN
  ELSE ! signs, exponents and number of digits are equal 
    DO i=0,UBOUND(l%sig,1)
      IF(l%sig(i)/=r%sig(i))THEN
         num_eq_num = .FALSE.; RETURN
      ENDIF
    ENDDO
    num_eq_num = .TRUE.; RETURN
  ENDIF
ENDFUNCTION num_eq_num

!--precision rounding procedures----------------------------------------------!
FUNCTION num_trunc(val,dp) ! generic TRUNC
  type(NUMBER),INTENT(IN) :: val
  INTEGER,INTENT(IN),OPTIONAL :: dp
  type(NUMBER) :: num_trunc
  ! truncates val towards zero to the value less than or equal to val
  ! with dp decimal digits after the decimal point
  ! if dp is absent zero is assumed and the result is the
  ! the largest integer value less than or equal to ABS(val) in absolute value
  INTEGER :: nd,nsigd,nrp,ntop,nbot
  nbot=-val%exp
  ntop=nbot+UBOUND(val%sig,1)
  IF(PRESENT(dp))THEN
    nd=dp
  ELSE
    nd=0
  ENDIF
  nsigd=nd/radd
  nrp=nd-nsigd*radd
  IF(nrp>0)THEN
    nsigd=nsigd+1
  ELSE
    nrp=radd
  ENDIF 
  IF(nbot>nsigd)THEN ! result is zero correct to nd places
   num_trunc%exp=0
   ALLOCATE(num_trunc%sig(0:0))
   num_trunc%sig(0)=0
  ELSEIF(ntop<nsigd)THEN ! result is already truncated 
   num_trunc%exp=val%exp
   ALLOCATE(num_trunc%sig(0:UBOUND(val%sig,1)))
   num_trunc%sig=val%sig
  ELSE ! result requires truncation
   num_trunc%exp=val%exp
   ALLOCATE(num_trunc%sig(0:-nbot+nsigd))
   num_trunc%sig(0:-nbot+nsigd)=val%sig(0:-nbot+nsigd)
   nrp=10**(radd-nrp)
   num_trunc%sig(-nbot+nsigd)=(num_trunc%sig(-nbot+nsigd)/nrp)*nrp
  ENDIF
ENDFUNCTION num_trunc

!--Overloads of intrinsic procedures------------------------------------------!
FUNCTION num_abs(a)
  type(NUMBER),INTENT(IN) :: a
  type(NUMBER) :: num_abs
  ! returns the absolute value of its argument
  IF(ASSOCIATED(a%sig))THEN
    ALLOCATE(num_abs%sig(0:UBOUND(a%sig,1)))
    num_abs%sig=ABS(a%sig)
    num_abs%exp=a%exp
  ELSE
    NULLIFY(num_abs%sig)
  ENDIF
ENDFUNCTION num_abs

FUNCTION num_sign(a,b)
  type(NUMBER),INTENT(IN) :: a,b
  type(NUMBER) :: num_sign
  ! returns the absolute value of a times the sign of b
  ALLOCATE(num_sign%sig(0:UBOUND(a%sig,1)))
  num_sign%exp=a%exp
  IF(b%sig(0)>=0)THEN
    num_sign%sig=ABS(a%sig)
   ELSE
    num_sign%sig=-ABS(a%sig)
  ENDIF
ENDFUNCTION num_sign

ENDMODULE VARIABLE_PRECISION_ARITHMETIC