!   EXPERIMENTAL MODULES FOR ARITHMETIC WITH ERROR BOUNDS  V0.1
    !   
    !   Author:       Abraham Agay  (agay@vms.huji.ac.il)
    !   Date:         May 1999
    !   
    !   Disclaimer:  The following software may have bugs!

    !   Limitations:
    !
    !     1)  Only the most important overloadings were implemented,
    !         e.g. operations involving integers or reals other than
    !         the PRECISE type were not, also some intrinsic functions.
    !
    !     2)  The routines in this package doesn't check input vars,
    !         so bad input, e.g. negative error bounds will produce
    !         ridiculous results.
    !
    !     3)  No attempt was made to code for efficiency, and speed
    !         was willingly sacrificed for clarity/elegance.

    !   Usage: 
    !
    !     0)  Call "num_check()" to verify you have PRECISE REALs.
    !     1)  Add "use err_arithmetic" in every compilation unit
    !         that does computations with REALs.
    !     2)  Replace all REAL declarations by "type(err_real) :: ".
    !     3)  Convert integers to "real(kind = PRECISE)" in all mixed 
    !         REAL/INTEGER operations, e.g  N --> real(N, PRECISE).
    !     4)  Add "use err_io" in compilation units doing REAL I/O,
    !         replace all I/O statements involving REALs with 
    !         "err_write_line" and "...".
    !     5)  Provide input with error bounds.

    !   Remarks:
    !
    !   A strange feature of DEC Fortran 90: PARAMETERs can't be SAVEd!
    !   According to the Fortran 90 Handbook this is a non-standard feature.
    !   Does it matter when using module sub-programs?
    !
    !   Another problem of DEC Fortran 90: too many optional procedure 
    !   arguments breaks the code.

    !   Personal Remark:  Benjamin Netanyahu was a good prime-minister.


MODULE NUM_CONSTS 
  implicit none

    !   The following module selects the kind of real used according
    !   to specified precision and range parameters:
    !   
    !      N_SIG_DIGITS       Number of significant digits required
    !      EXP_RANGE          Exponent range required
    !
    !   Various useful constants are then defined in the selected
    !   real type, and two useful procedures.

    integer, parameter  ::                                      &
      N_SIG_DIGITS =      10,                                   &
      EXP_RANGE =         100,                                  &
      PRECISE =  SELECTED_REAL_KIND(N_SIG_DIGITS, EXP_RANGE)

    real(kind=PRECISE),parameter  ::                            &
      ZERO =              0.0_PRECISE,                          &
      ONE =               1.0_PRECISE,                          &
      TWO =               2.0_PRECISE,                          &
      THREE =             3.0_PRECISE,                          &
      HALF =              ONE / TWO                               ! 0.5_PRECISE

    real(kind=PRECISE),parameter  ::                            &
      PI =                3.141592653589793238462643_PRECISE,   &
      HALFPI =            PI / TWO,                             &
      DEG2RAD =           PI / 180.0_PRECISE,                   &
      RAD2DEG =           180.0_PRECISE / PI,                   &
      RAD2MIN =           10800.0_PRECISE / PI,                 &
      RAD2SEC =           648000.0_PRECISE / PI


  contains

    SUBROUTINE NUM_CHECK()
      select case (PRECISE)
      case (-1)
        stop 'num_check: Insufficient precision '
      case (-2)
        stop 'num_check: Insufficient exponent range '
      case (-3)
        stop 'num_check: Insufficient precision and exponent range '
      case default
        write(*,*) 'num_check: suitable type of REAL is available '
!       call num_status()
      end select
      return
    END SUBROUTINE NUM_CHECK

    SUBROUTINE NUM_STATUS()
      real(kind=PRECISE) :: x
      write (unit=*, fmt='(1x,a,i12)')                                        &
            'num_status: The base of the number model is:  ', radix(x)
      write (unit=*, fmt='(1x,a,i12)')                                        &
            'num_status: Number of digits in this base:    ', digits(x)
      write (unit=*, fmt='(1x,a,i12)')                                        &
            'num_status: Minimum exponent in this base:    ', minexponent(x)
      write (unit=*, fmt='(1x,a,i12)')                                        &
            'num_status: Maximum exponent in this base:    ', maxexponent(x)
      write (unit=*, fmt='(1x,a,i12)')                                        &
            'num_status: The PRECISE real type is kind:    ', kind(x)
      write (unit=*, fmt='(1x,a,i12)')                                        & 
            'num_status: Decimal exponent range is:        ', range(x)
      write (unit=*, fmt='(1x,a,es12.3)')                                     &
            'num_status: Smallest positive number is:      ', tiny(x)
      write (unit=*, fmt='(1x,a,es12.3)')                                     &
            'num_status: Largest positive number is:       ', huge(x)
      write (unit=*, fmt='(1x,a,i12,a)')                                      & 
            'num_status: Decimal precision is:  ', precision(x), ' digits '
      write (unit=*, fmt='(1x,a,es12.3)')                                     &
            'num_status: Epsilon is:                       ', epsilon(x)
      write (unit=*, fmt='(1x)')
    END SUBROUTINE NUM_STATUS

END MODULE NUM_CONSTS 



MODULE ERR_ARITHMETIC
  use NUM_CONSTS
  implicit none

    !   We define 3 new compound types of real numbers: intervals,
    !   absolute error-ranges, and relative error-ranges.
    !
    !     IVL_REAL      (a, b)
    !
    !     ERR_REAL      (A +|- a)  =  (A - a, A + a)
    !
    !     REL_REAL      (A +|- a/A)

    !   Err-arithmetic is performed by converting to intervals, 
    !   implicitly or explicitly, performing the required operation, 
    !   rounding outwards, and converting back to err-arithmetic.

    !   Arithmetical operations near a singularity, for example,
    !   when taking tan() near the point  PI / TWO, give rise to
    !   multiply-connected error ranges.  We can define, for each  
    !   real type two sub-types:
    !
    !     I.   Normal, contiguous, one-part interval
    !          IVL_REAL:   LOWER  <= UPPER
    !          ERR_REAL:   ABSERR >= ZERO
    !
    !     II.  A two-part interval which is the set-theoretic
    !          complement of a type I interval.
    !          IVL_REAL:   LOWER  > UPPER
    !          ERR_REAL:   ABSERR < ZERO
    !
    !   However, this doesn't close interval operations, as more 
    !   complicated error ranges may be produced.  We will use 
    !   only TYPE I intervals in this package, and flag an error 
    !   when a singularity is encountered.  The following two 
    !   parameters control error handling:

    integer, parameter     ::  MAXWARN = 20, MAXERR = 1

    TYPE IVL_REAL
      real(kind = PRECISE) ::  LOWER, UPPER
    END TYPE IVL_REAL

    TYPE ERR_REAL
      real(kind = PRECISE) ::  NUMBER, ABSERR
    END TYPE ERR_REAL

!   TYPE REL_REAL
!     real(kind = PRECISE) ::  NUMBER, RELERR
!   END TYPE REL_REAL

    !   Extending assignment to convert between the 2 new types 
    !   of reals and between the new types and the old ones.
    !   This technique will be used extensively, it's probably 
    !   bad programming (see Cooper Redwine) but it's tempting. 
    !
    !   Conversion to real is not supported on purpose.

    INTERFACE ASSIGNMENT(=) 
      MODULE PROCEDURE ERR_2_IVL
      MODULE PROCEDURE IVL_2_ERR
      MODULE PROCEDURE REAl_2_ERR
      MODULE PROCEDURE REAL_2_IVL
    END INTERFACE

    !   N e w   o p e r a t o r:  i n t e r v a l   i n c l u s i o n

    INTERFACE OPERATOR(.in.)
      MODULE PROCEDURE ERR_IN_IVL_IVL
      MODULE PROCEDURE ERR_IN_REAL_IVL
    END INTERFACE

    !   R e l a t i o n a l   o p e r a t o r s 

    INTERFACE OPERATOR(.gt.)
      MODULE PROCEDURE ERR_GT_ERR_ERR
      MODULE PROCEDURE ERR_GT_REAL_ERR
      MODULE PROCEDURE ERR_GT_ERR_REAL
    END INTERFACE

    INTERFACE OPERATOR(.ge.)
      MODULE PROCEDURE ERR_GE_ERR_ERR
      MODULE PROCEDURE ERR_GE_REAL_ERR
      MODULE PROCEDURE ERR_GE_ERR_REAL
    END INTERFACE

    INTERFACE OPERATOR(.lt.)
      MODULE PROCEDURE ERR_LT_ERR_ERR
      MODULE PROCEDURE ERR_LT_REAL_ERR
      MODULE PROCEDURE ERR_LT_ERR_REAL
    END INTERFACE

    INTERFACE OPERATOR(.le.)
      MODULE PROCEDURE ERR_LE_ERR_ERR
      MODULE PROCEDURE ERR_LE_REAL_ERR
      MODULE PROCEDURE ERR_LE_ERR_REAL
    END INTERFACE

!   INTERFACE OPERATOR(.eq.)             !  Not very useful, problematic
!     MODULE PROCEDURE ERR_EQ_ERR_ERR
!     MODULE PROCEDURE ERR_EQ_REAL_ERR
!     MODULE PROCEDURE ERR_EQ_ERR_REAL
!   END INTERFACE

!   INTERFACE OPERATOR(.ne.)             !  Not very useful, problematic
!     MODULE PROCEDURE ERR_NE_ERR_ERR
!     MODULE PROCEDURE ERR_NE_REAL_ERR
!     MODULE PROCEDURE ERR_NE_ERR_REAL
!   END INTERFACE

    !   I n t e r f a c e s   t o   b a s i c   o p e r a t i o n s

    INTERFACE OPERATOR(+)
      MODULE PROCEDURE ERR_ADD_ERR_ERR
      MODULE PROCEDURE ERR_ADD_REAL_ERR
      MODULE PROCEDURE ERR_ADD_ERR_REAL
    END INTERFACE

    INTERFACE OPERATOR(-)
      MODULE PROCEDURE ERR_SUB_ERR_ERR
      MODULE PROCEDURE ERR_SUB_REAL_ERR
      MODULE PROCEDURE ERR_SUB_ERR_REAL
      MODULE PROCEDURE ERR_SUB_UNARY_ERR
    END INTERFACE

    INTERFACE OPERATOR(*)
      MODULE PROCEDURE ERR_MULT_ERR_ERR
      MODULE PROCEDURE ERR_MULT_REAL_ERR
      MODULE PROCEDURE ERR_MULT_ERR_REAL
    END INTERFACE

    INTERFACE OPERATOR(/)
      MODULE PROCEDURE ERR_DIV_ERR_ERR
      MODULE PROCEDURE ERR_DIV_ERR_REAL
      MODULE PROCEDURE ERR_DIV_REAL_ERR
    END INTERFACE

    !   I n t e r f a c e s   t o   e l e m e n t a r y   f u n c t i o n s

    INTERFACE OPERATOR(**)
      MODULE PROCEDURE ERR_POWER_ERR_INT
      MODULE PROCEDURE ERR_POWER_ERR_REAL
    END INTERFACE

    INTERFACE ABS
      MODULE PROCEDURE ERR_ABS
    END INTERFACE

    INTERFACE SQRT
      MODULE PROCEDURE ERR_SQRT
    END INTERFACE

    INTERFACE EXP
      MODULE PROCEDURE ERR_EXP
    END INTERFACE

    INTERFACE LOG
      MODULE PROCEDURE ERR_LOG
    END INTERFACE

    INTERFACE SIN
      MODULE PROCEDURE ERR_SIN
    END INTERFACE

    INTERFACE COS
      MODULE PROCEDURE ERR_COS
    END INTERFACE

    INTERFACE TAN
      MODULE PROCEDURE ERR_TAN
    END INTERFACE

!   INTERFACE ALT_TAN
!     MODULE PROCEDURE ERR_ALT_TAN
!   END INTERFACE

!   INTERFACE COT
!     MODULE PROCEDURE ERR_COT
!   END INTERFACE

    INTERFACE ATAN
      MODULE PROCEDURE ERR_ATAN
    END INTERFACE

!   INTERFACE ACOT
!     MODULE PROCEDURE ERR_ACOT
!   END INTERFACE

    INTERFACE ASIN
      MODULE PROCEDURE ERR_ASIN
    END INTERFACE

    INTERFACE ACOS
      MODULE PROCEDURE ERR_ACOS
    END INTERFACE

!   INTERFACE SINH
!     MODULE PROCEDURE ERR_SINH
!   END INTERFACE


  contains

    !   These important checking routines are not used yet.

    SUBROUTINE ERR_CHECK_ERR(A)
      type(err_real)               :: a
      if (a%abserr < ZERO)   &
         call err_error(2, 'err_check_err: bad error range ')
      return
    END SUBROUTINE ERR_CHECK_ERR

    SUBROUTINE ERR_CHECK_IVL(A)
      type(ivl_real)               :: a
      if (a%lower > a%upper)   &
         call err_error(2, 'err_check_ivl: bad error range ')
      return
    END SUBROUTINE ERR_CHECK_IVL

    !   A useful routine, rounds an interval outwards.
    !
    !   Using the "spacing" function instead of "nearest" is 
    !   probably incorrect. 
    !
    !      err_round%lower = a%lower - spacing(a%lower)
    !      err_round%upper = a%upper + spacing(a%upper)
    !
    !   One example that comes to mind: with DEC floating-point 
    !   numbers, the spacing is actually asymmetrical (?).

    FUNCTION ERR_ROUND(A)
      type(ivl_real)               :: err_round, a
      err_round%lower = nearest(a%lower, -1.0)
      err_round%upper = nearest(a%upper, +1.0)
      return
    END FUNCTION ERR_ROUND

    !   E x t e n d i n g   a s s i g n m e n t   (see remark above) 

    SUBROUTINE ERR_2_IVL(A, B)
      type(ivl_real)               :: a
      type(err_real), intent(in)   :: b
      a%lower = b%number - b%abserr
      a%upper = b%number + b%abserr
      return
    END SUBROUTINE ERR_2_IVL

    SUBROUTINE IVL_2_ERR(A, B)
      type(err_real)               :: a
      type(ivl_real), intent(in)   :: b
      a%number = (b%upper + b%lower) / TWO
      a%abserr = (b%upper - b%lower) / TWO
      return
    END SUBROUTINE IVL_2_ERR

    SUBROUTINE REAL_2_ERR(A, B)
      type(err_real)                     :: a
      real(kind = PRECISE), intent(in)   :: b
      a%number = b 
      a%abserr = ZERO 
      return
    END SUBROUTINE REAL_2_ERR

    SUBROUTINE REAL_2_IVL(A, B)
      type(ivl_real)                     :: a
      real(kind = PRECISE), intent(in)   :: b
      a%lower = b 
      a%upper = b 
      return
    END SUBROUTINE REAL_2_IVL

    !   D e f i n i n g   i n t e r v a l   i n c l u s i o n  (new operator)

    FUNCTION ERR_IN_IVL_IVL(A, B)
      logical                       :: err_in_ivl_ivl
      type(ivl_real), intent(in)    :: a, b
      if ((a%lower >= b%lower) .and. (a%upper <= b%upper)) then
        err_in_ivl_ivl = .true.
      else
        err_in_ivl_ivl = .false.
      endif
      return
    END FUNCTION ERR_IN_IVL_IVL

    FUNCTION ERR_IN_REAL_IVL(A, B)
      logical                             :: err_in_real_ivl
      real(kind = PRECISE), intent(in)    :: a
      type(ivl_real), intent(in)          :: b
      type(ivl_real)                      :: tmp
      tmp = a
      err_in_real_ivl = tmp .in. b
      return
    END FUNCTION ERR_IN_REAL_IVL

    !   E x t e n d i n g   r e l a t i o n a l   o p e r a t o r s 

    !   .gt.

    FUNCTION ERR_GT_ERR_ERR(A, B)
      logical                      :: err_gt_err_err
      type(err_real), intent(in)   :: a, b
      type(ivl_real)               :: x, y
      x = a
      y = b
      if (x%lower > y%upper) then
        err_gt_err_err = .true.
      else
        err_gt_err_err = .false.
      endif
      return
    END FUNCTION ERR_GT_ERR_ERR

    FUNCTION ERR_GT_REAL_ERR(A, B)
      logical                            :: err_gt_real_err
      real(kind = PRECISE), intent(in)   :: a
      type(err_real),       intent(in)   :: b
      type(err_real)                     :: tmp
      tmp = a
      err_gt_real_err = (tmp .gt. b)
      return
    END FUNCTION ERR_GT_REAL_ERR

    FUNCTION ERR_GT_ERR_REAL(A, B)
      logical                            :: err_gt_err_real
      type(err_real),       intent(in)   :: a
      type(err_real)                     :: tmp
      real(kind = PRECISE), intent(in)   :: b
      tmp = b
      err_gt_err_real = (a .gt. tmp)
      return
    END FUNCTION ERR_GT_ERR_REAL

    !   .ge.

    FUNCTION ERR_GE_ERR_ERR(A, B)
      logical                      :: err_ge_err_err
      type(err_real), intent(in)   :: a, b
      type(ivl_real)               :: x, y
      x = a
      y = b
      if (x%lower >= y%upper) then
        err_ge_err_err = .true.
      else
        err_ge_err_err = .false.
      endif
      return
    END FUNCTION ERR_GE_ERR_ERR

    FUNCTION ERR_GE_REAL_ERR(A, B)
      logical                            :: err_ge_real_err
      real(kind = PRECISE), intent(in)   :: a
      type(err_real),       intent(in)   :: b
      type(err_real)                     :: tmp
      tmp = a
      err_ge_real_err = (tmp .ge. b)
      return
    END FUNCTION ERR_GE_REAL_ERR

    FUNCTION ERR_GE_ERR_REAL(A, B)
      logical                            :: err_ge_err_real
      type(err_real),       intent(in)   :: a
      type(err_real)                     :: tmp
      real(kind = PRECISE), intent(in)   :: b
      tmp = b
      err_ge_err_real = (a .ge. tmp)
      return
    END FUNCTION ERR_GE_ERR_REAL

    !   .lt.

    FUNCTION ERR_LT_ERR_ERR(A, B)
      logical                      :: err_lt_err_err
      type(err_real), intent(in)   :: a, b
      type(ivl_real)               :: x, y
      x = a
      y = b
      if (x%lower < y%upper) then
        err_lt_err_err = .true.
      else
        err_lt_err_err = .false.
      endif
      return
    END FUNCTION ERR_LT_ERR_ERR

    FUNCTION ERR_LT_REAL_ERR(A, B)
      logical                            :: err_lt_real_err
      real(kind = PRECISE), intent(in)   :: a
      type(err_real),       intent(in)   :: b
      type(err_real)                     :: tmp
      tmp = a
      err_lt_real_err = (tmp .lt. b)
      return
    END FUNCTION ERR_LT_REAL_ERR

    FUNCTION ERR_LT_ERR_REAL(A, B)
      logical                            :: err_lt_err_real
      type(err_real),       intent(in)   :: a
      type(err_real)                     :: tmp
      real(kind = PRECISE), intent(in)   :: b
      tmp = b
      err_lt_err_real = (a .lt. tmp)
      return
    END FUNCTION ERR_LT_ERR_REAL

    !   .le.

    FUNCTION ERR_LE_ERR_ERR(A, B)
      logical                      :: err_le_err_err
      type(err_real), intent(in)   :: a, b
      type(ivl_real)               :: x, y
      x = a
      y = b
      if (x%lower <= y%upper) then
        err_le_err_err = .true.
      else
        err_le_err_err = .false.
      endif
      return
    END FUNCTION ERR_LE_ERR_ERR

    FUNCTION ERR_LE_REAL_ERR(A, B)
      logical                            :: err_le_real_err
      real(kind = PRECISE), intent(in)   :: a
      type(err_real),       intent(in)   :: b
      type(err_real)                     :: tmp
      tmp = a
      err_le_real_err = (tmp .le. b)
      return
    END FUNCTION ERR_LE_REAL_ERR

    FUNCTION ERR_LE_ERR_REAL(A, B)
      logical                            :: err_le_err_real
      type(err_real),       intent(in)   :: a
      type(err_real)                     :: tmp
      real(kind = PRECISE), intent(in)   :: b
      tmp = b
      err_le_err_real = (a .le. tmp)
      return
    END FUNCTION ERR_LE_ERR_REAL

    !   E x t e n d i n g   a d d i t i o n

    !   Master addition definition, if one operand is real we convert.
    !
    !   Why not use just:
    !
    !      err_add_err_err%number = a%number + b%number
    !      err_add_err_err%abserr = a%abserr + b%abserr
    !
    !   The idea is to have the same procedure of rounding as in the
    !   multiplication and division procedures.  In both cases we go
    !   via interval arithmetic, and perform there the rounding.
    !
    !   Another alternative is: 
    !
    !      tmp%lower = a%number + b%number - a%abserr - b%abserr
    !      tmp%upper = a%number + b%number + a%abserr + b%abserr
    !      err_add_err_err = err_round(tmp) 
    !
    !   This is o.k., but we will go for a more elegant method.

    FUNCTION ERR_ADD_ERR_ERR(A,B)
      type(err_real)               :: err_add_err_err
      type(err_real), intent(in)   :: a, b
      type(ivl_real)               :: x, y, tmp
      x = a
      y = b
      tmp%lower = x%lower + y%lower
      tmp%upper = x%upper + y%upper
      err_add_err_err = err_round(tmp) 
      return
    END FUNCTION ERR_ADD_ERR_ERR

    FUNCTION ERR_ADD_REAL_ERR(A,B)   
      type(err_real)                      :: err_add_real_err, tmp
      real(kind = PRECISE), intent(in)    :: a
      type(err_real), intent(in)          :: b
      tmp = a
      err_add_real_err = tmp + b
      return
    END FUNCTION ERR_ADD_REAL_ERR

    FUNCTION ERR_ADD_ERR_REAL(A,B)
      type(err_real)                      :: err_add_err_real, tmp
      type(err_real), intent(in)          :: a
      real(kind = PRECISE), intent(in)    :: b
      tmp = b
      err_add_err_real = a + tmp 
      return
    END FUNCTION ERR_ADD_ERR_REAL

    !   E x t e n d i n g   s u b s t r a c t i o n

    !   Master substraction definition, if one operand is real we convert.
    !
    !      err_sub_err_err%number = a%number - b%number
    !      err_sub_err_err%abserr = a%abserr + b%abserr
    !
    !   Again, we don't use these formulae, as we want to get a standard 
    !   method of rounding.
    !
    !   Another method:
    !
    !      tmp%lower = a%number - b%number - a%abserr - b%abserr
    !      tmp%upper = a%number - b%number + a%abserr + b%abserr
    !      err_sub_err_err = err_round(tmp) 
    !
    !   This is o.k., but we will go for a more elegant method,
    !   via the operation of unary minus and addition.

    FUNCTION ERR_SUB_UNARY_ERR(A)
      type(err_real)                      :: err_sub_unary_err
      type(err_real), intent(in)          :: a
      err_sub_unary_err = err_real(- a%number, a%abserr) 
      return
    END FUNCTION ERR_SUB_UNARY_ERR

    FUNCTION ERR_SUB_ERR_ERR(A,B)
      type(err_real)               :: err_sub_err_err
      type(err_real), intent(in)   :: a, b
      err_sub_err_err = a + (- b) 
      return
    END FUNCTION ERR_SUB_ERR_ERR

    FUNCTION ERR_SUB_REAL_ERR(A,B)
      type(err_real)                      :: err_sub_real_err, tmp
      real(kind = PRECISE), intent(in)    :: a
      type(err_real), intent(in)          :: b
      tmp = a
      err_sub_real_err = tmp - b 
      return
    END FUNCTION ERR_SUB_REAL_ERR

    FUNCTION ERR_SUB_ERR_REAL(A,B)
      type(err_real)                      :: err_sub_err_real, tmp
      type(err_real), intent(in)          :: a
      real(kind = PRECISE), intent(in)    :: b
      tmp = b
      err_sub_err_real = a - tmp 
      return
    END FUNCTION ERR_SUB_ERR_REAL

    !   Not used, a different algorithm.  Is it o.k.? 

    FUNCTION ERR_ALT_MULT(A,B)
      type(err_real)             :: err_alt_mult
      type(err_real), intent(in) :: a, b
      err_alt_mult%number = a%number * b%number + a%abserr * b%abserr
      err_alt_mult%abserr = abs(a%number) * b%abserr + a%abserr * abs(b%number)
      return
    END FUNCTION ERR_ALT_MULT

    !   E x t e n d i n g   m u l t i p l i c a t i o n

    FUNCTION ERR_MULT_ERR_ERR(A,B)
      type(err_real)               :: err_mult_err_err
      type(err_real), intent(in)   :: a, b
      type(ivl_real)               :: x, y, tmp
      x = a   ;   y = b
      tmp%lower = min(x%lower * y%lower, x%lower * y%upper,      &  
                      x%upper * y%lower, x%upper * y%upper)
      tmp%upper = max(x%lower * y%lower, x%lower * y%upper,      &  
                      x%upper * y%lower, x%upper * y%upper)
      err_mult_err_err = err_round(tmp) 
      return
    END FUNCTION ERR_MULT_ERR_ERR

    FUNCTION ERR_MULT_REAL_ERR(A,B)
      type(err_real)               :: err_mult_real_err, tmp
      real(kind = PRECISE), intent(in)   :: a
      type(err_real), intent(in)         :: b
      tmp = a
      err_mult_real_err = tmp + b
      return
    END FUNCTION ERR_MULT_REAL_ERR

    FUNCTION ERR_MULT_ERR_REAL(A,B)
      type(err_real)               :: err_mult_err_real, tmp
      type(err_real), intent(in)         :: a
      real(kind = PRECISE), intent(in)   :: b
      tmp = b
      err_mult_err_real = a + tmp 
      return
    END FUNCTION ERR_MULT_ERR_REAL

    !   E x t e n d i n g   d i v i s i o n

    FUNCTION ERR_DIV_ERR_ERR(A,B)
      type(err_real)               :: err_div_err_err
      type(err_real), intent(in)   :: a, b
      type(ivl_real)               :: x, y, tmp
      x = a
      y = b
      if (ZERO .in. y)  then
        err_div_err_err = err_real(ZERO, huge(ZERO))      !  Stupid handling?
        call err_error(1, 'err_div_err_err:  infinite error range ')
      else
        tmp%lower = min(x%lower / y%lower, x%lower / y%upper,      &  
                        x%upper / y%lower, x%upper / y%upper)
        tmp%upper = max(x%lower / y%lower, x%lower / y%upper,      &  
                        x%upper / y%lower, x%upper / y%upper)
        err_div_err_err = err_round(tmp) 
      endif
      return
    END FUNCTION ERR_DIV_ERR_ERR

    FUNCTION ERR_DIV_REAL_ERR(A,B)
      type(err_real)                      :: err_div_real_err, tmp
      real(kind = PRECISE), intent(in)    :: a
      type(err_real),       intent(in)    :: b
      tmp = a
      err_div_real_err = tmp / b
      return
    END FUNCTION ERR_DIV_REAL_ERR

    FUNCTION ERR_DIV_ERR_REAL(A,B)
      type(err_real)                      :: err_div_err_real, tmp
      type(err_real),       intent(in)    :: a
      real(kind = PRECISE), intent(in)    :: b
      tmp = b
      err_div_err_real = a / tmp 
      return
    END FUNCTION ERR_DIV_ERR_REAL

    !   E x t e n d i n g   e l e m e n t a l   f u n c t i o n s

    !   The ABS intrinsic doesn't need an outward rounding. 

    FUNCTION ERR_ABS(A)
      type(err_real)                      :: err_abs
      type(err_real), intent(in)          :: a
      type(ivl_real)                      :: x, tmp
      x = a
      tmp%upper = max(abs(x%lower), abs(x%upper))
      tmp%lower = min(abs(x%lower), abs(x%upper))     ! Should be in "else"
      if (ZERO .in. x) then
        err_abs = ivl_real(ZERO, tmp%upper)
      else
        err_abs = tmp
      endif
      return
    END FUNCTION ERR_ABS

    !   V a r i o u s   p o w e r   o p e r a t i o n s

    !   How to treat a sqrt() of a range which is partly negative? 

    FUNCTION ERR_SQRT(A)
      type(err_real)                      :: err_sqrt 
      type(err_real), intent(in)          :: a
      type(ivl_real)                      :: x, tmp
      x = a
      tmp%lower = sqrt(max(x%lower, ZERO))      !  o.k. ???
      tmp%upper = sqrt(x%upper)
      err_sqrt = err_round(tmp) 
      return
    END FUNCTION ERR_SQRT

    FUNCTION ERR_EXP(A)
      type(err_real)                      :: err_exp
      type(err_real), intent(in)          :: a
      type(ivl_real)                      :: x, tmp
      x = a
      tmp%lower = exp(x%lower) 
      tmp%upper = exp(x%upper)  
      err_exp = err_round(tmp) 
      return
    END FUNCTION ERR_EXP

    FUNCTION ERR_LOG(A)
      type(err_real)                      :: err_log
      type(err_real), intent(in)          :: a
      type(ivl_real)                      :: x, tmp
      x = a
      tmp%lower = log(x%lower) 
      tmp%upper = log(x%upper)  
      err_log = err_round(tmp) 
      return
    END FUNCTION ERR_LOG

    !   Defining the exponentiation operator.
    !   The integer case is not treated as a private case of real, 
    !   to allow the compiler to do its accuracy optimizations

    FUNCTION ERR_POWER_ERR_INT(A,B)
      type(err_real)                      :: err_power_err_int
      type(err_real), intent(in)          :: a
      integer,        intent(in)          :: b
      type(ivl_real)                      :: x, tmp
      x = a
      tmp%lower = min(x%lower ** b, x%upper ** b) 
      tmp%upper = max(x%lower ** b, x%upper ** b) 
      err_power_err_int = err_round(tmp) 
      return
    END FUNCTION ERR_POWER_ERR_INT

    FUNCTION ERR_POWER_ERR_REAL(A,B)
      type(err_real)                      :: err_power_err_real
      type(err_real),       intent(in)    :: a
      real(kind = PRECISE), intent(in)    :: b
      type(ivl_real)                      :: x, tmp
      x = a
      tmp%lower = min(x%lower ** b, x%upper ** b) 
      tmp%upper = max(x%lower ** b, x%upper ** b) 
      err_power_err_real = err_round(tmp) 
      return
    END FUNCTION ERR_POWER_ERR_REAL

    !   E x t e n d i n g   t r i g o n o m e t r i c   f u n c t i o n s

    FUNCTION ERR_SIN(A)
      type(err_real)              :: err_sin
      type(err_real), intent(in)  :: a
      type(ivl_real)              :: x, tmp
      real(kind = PRECISE)        :: offset
      if (a%abserr >= TWO * PI) then
        err_sin = err_real(ZERO, ONE)
      else 
        x = a
        offset = x%lower - mod(x%lower, TWO * PI)
        if ((offset + HALFPI) .in. x) then
          tmp%upper = ONE
          if ((offset + THREE * HALFPI) .in. x) then
            tmp%lower = - ONE
          else
            tmp%lower = min(sin(x%lower), sin(x%upper))
          endif
        else
          tmp%upper = max(sin(x%lower), sin(x%upper))
          if ((offset + THREE * HALFPI) .in. x) then
            tmp%lower = - ONE
          else
            tmp%lower = min(sin(x%lower), sin(x%upper))
          endif
        endif
        err_sin = err_round(tmp)  
      endif
      return
    END FUNCTION ERR_SIN 

    FUNCTION ERR_COS(A)
      type(err_real)              :: err_cos
      type(err_real), intent(in)  :: a
      err_cos = sin(err_real(PI / TWO - a%number, a%abserr))
      return
    END FUNCTION ERR_COS

    FUNCTION ERR_TAN(A)
      type(err_real)              :: err_tan
      type(err_real), intent(in)  :: a
      err_tan = sin(a) / cos(a) 
      return
    END FUNCTION ERR_TAN

    !   This direct algorithm seems to give identical results.

    FUNCTION ERR_ALT_TAN(A)
      type(err_real)              :: err_alt_tan
      type(err_real), intent(in)  :: a
      type(ivl_real)              :: x, tmp
      real(kind = PRECISE)        :: offset 
      if (a%abserr >= TWO * PI) then
        err_alt_tan = err_real(ZERO, huge(ZERO))
        call err_error(1, 'err_alt_tan:  a singularity ')
      else
        x = a
        offset = x%lower - mod(x%lower, PI)
        if ((offset + HALFPI) .in. x) then
!         err_alt_tan = ivl_real(tan(x%lower), tan(x%upper))   ! Reversed range
          err_alt_tan = err_real(ZERO, huge(ZERO))
          call err_error(1, 'err_alt_tan:  a singularity ')
        else
          tmp%lower = min(tan(x%lower), tan(x%upper))
          tmp%upper = max(tan(x%lower), tan(x%upper))
          err_alt_tan = err_round(tmp) 
        endif
      endif
      return
    END FUNCTION ERR_ALT_TAN

    !   ASIN is monotonic increasing, and hence easy to implement. 

    FUNCTION ERR_ASIN(A)
      type(err_real)              :: err_asin
      type(err_real), intent(in)  :: a
      type(ivl_real)              :: x, tmp
      x = a
      tmp%lower = asin(x%lower)
      tmp%upper = asin(x%upper)
      err_asin = err_round(tmp) 
    END FUNCTION ERR_ASIN

    !   We use the identity:   asin(z) + acos(z) = PI / TWO

    FUNCTION ERR_ACOS(A)
      type(err_real)              :: err_acos
      type(err_real), intent(in)  :: a
      err_acos = PI / TWO - asin(a)
    END FUNCTION ERR_ACOS

    !   ATAN is monotonic increasing on the whole real line.

    FUNCTION ERR_ATAN(A)
      type(err_real)              :: err_atan
      type(err_real), intent(in)  :: a
      type(ivl_real)              :: x, tmp
      x = a
      tmp%lower = atan(x%lower)
      tmp%upper = atan(x%upper)
      err_atan = err_round(tmp) 
    END FUNCTION ERR_ATAN

    SUBROUTINE ERR_ERROR(N,STRING)
      integer,           intent(in)   :: N
      character(len =*), intent(in)   :: string
      integer, save                   :: nwarn = 0, nerr = 0
      select case (N)
      case (0)
      case (1)
        write(*,*) string
        nwarn = nwarn + 1
      case (2)
        write(*,*) string
        nerr = nerr + 1
      case default
        stop 'err_error: bad error number '
      end select
      if (nwarn >= MAXWARN) stop 'err_error: too many warnings, aborting '
      if (nerr  >= MAXERR)  stop 'err_error: too many errors, aborting '
      return
    END SUBROUTINE ERR_ERROR

END MODULE ERR_ARITHMETIC


MODULE ERR_IO
  use NUM_CONSTS,     only  : PRECISE
  use ERR_ARITHMETIC, only  : ivl_real, err_real
  implicit none

  contains

    !   Reading one err-real without advancing

    SUBROUTINE ERR_READ(LUN, A)
      integer              :: lun, length
      type(err_real)       :: a
      character(len = 80)  :: buf
      read(unit = lun, fmt = '(a)', advance = 'NO')  buf 
    END SUBROUTINE ERR_READ

    !   Advancing output

    SUBROUTINE ERR_EOR(LUN)
      integer              :: lun
      write (unit = lun, fmt = '(1x)')
      return 
    END SUBROUTINE ERR_EOR

    !   Writing one err-real without advancing

    SUBROUTINE ERR_WRITE(LUN, TEXT, PRECISION, A, MARK, B)
      integer                :: lun, precision, length
      character(len = *)     :: text
      real(kind = PRECISE)   :: a, b
      character(len = 3)     :: mark
      character(len = 80)    :: buf
      buf = '(1x,a,a1,sp,es' 
      length = 3 + precision + 2 + 3
      write (unit = buf(15:17), fmt = '(i3)') length
      buf(18:18) = '.'
      write (unit = buf(19:21),  fmt = '(i3)') precision 
      buf(22:23) = 'e3'
      buf(24:) = ',a3,es' // buf(15:23) // ',a1)'
!     write (*,*) '|', buf, '|'
      write (unit = lun, fmt = buf, advance = 'NO')    &
            text, '(', a, mark, b, ')'
      return
    END SUBROUTINE ERR_WRITE

    !   Only 2 variable sets are supported.

    SUBROUTINE ERR_WRITE_ERR_ERR(LUN, TEXT1, PRECISION1, A,   & 
                                   TEXT2, PRECISION2, B)
      integer                        :: lun
      integer,            optional, volatile   :: precision1, precision2
      character(len = *), optional, volatile   :: text1, text2
      type(err_real),     optional, volatile   :: a, b
      if (present(a))     &
         call ERR_WRITE(lun, text1, precision1, a%number, ' # ', a%abserr)
      if (present(b))     &
         call ERR_WRITE(lun, text2, precision2, b%number, ' # ', b%abserr)
      call ERR_EOR(LUN)
      return
    END SUBROUTINE ERR_WRITE_ERR_ERR

    SUBROUTINE ERR_WRITE_ERR_IVL(LUN, TEXT1, PRECISION1, A,   & 
                                   TEXT2, PRECISION2, B)
      integer                        :: lun
      integer,            optional, volatile   :: precision1, precision2
      character(len = *), optional, volatile   :: text1, text2
      type(err_real),     optional, volatile   :: a
      type(ivl_real),     optional, volatile   :: b
      if (present(a))     &
         call ERR_WRITE(lun, text1, precision1, a%number, ' # ', a%abserr)
      if (present(b))     &
         call ERR_WRITE(lun, text2, precision2, b%lower, ' , ', b%upper)
      call ERR_EOR(LUN)
      return
    END SUBROUTINE ERR_WRITE_ERR_IVL

    SUBROUTINE ERR_WRITE_IVL_ERR(LUN, TEXT1, PRECISION1, A,   & 
                                   TEXT2, PRECISION2, B)
      integer                        :: lun
      integer,            optional, volatile   :: precision1, precision2
      character(len = *), optional, volatile   :: text1, text2
      type(ivl_real),     optional, volatile   :: a
      type(err_real),     optional, volatile   :: b
      if (present(a))     &
         call ERR_WRITE(lun, text1, precision1, a%lower, ' , ', a%upper)
      if (present(b))     &
         call ERR_WRITE(lun, text2, precision2, b%number, ' # ', b%abserr)
      call ERR_EOR(LUN)
      return
    END SUBROUTINE ERR_WRITE_IVL_ERR

    SUBROUTINE ERR_WRITE_IVL_IVL(LUN, TEXT1, PRECISION1, A,   & 
                                   TEXT2, PRECISION2, B)
      integer                        :: lun
      integer,            optional, volatile   :: precision1, precision2
      character(len = *), optional, volatile   :: text1, text2
      type(ivl_real),     optional, volatile   :: a, b
      if (present(a))     &
         call ERR_WRITE(lun, text1, precision1, a%lower, ' , ', a%upper)
      if (present(b))     &
         call ERR_WRITE(lun, text2, precision2, b%lower, ' , ', b%upper)
      call ERR_EOR(LUN)
      return
    END SUBROUTINE ERR_WRITE_IVL_IVL

END MODULE ERR_IO