!   EXPERIMENTAL MODULES FOR ARITHMETIC WITH ERROR BOUNDS - V0.2
    !
    !   Contents:
    !
    !     NUM_CONSTS        Module defining the PRECISE real kind,
    !                       some useful constants and routines
    !     ERR_ARITHMETIC    Module implementing error analysis 
    !     test              An example program
    !     
    !  -----------------------------------------------------------------
    !     
    !   Author:       Abraham Agay  (agay@vms.huji.ac.il)
    !   Date:         May 1999, revised on June 2001
    !
    !   Disclaimer:  The following software may contain 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)  Add "use num_consts" in the main program,
    !         and 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)  replace all I/O statements involving REALs with
    !         "err_write_err" and "err_read_err", etc.
    !     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 break the code.
    !


MODULE NUM_CONSTS
  implicit none

    !   The following module selects the kind of real used according
    !   to specified precision and range parameters:
    !
    !      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  ::                                      &
      SIG_DIGITS =         10,                                  &
      EXP_RANGE =         100,                                  &
      PRECISE =  SELECTED_REAL_KIND(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,   &
      TWOPI   =           PI * TWO,                             &
      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: Base of 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
  PRIVATE

    !   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)

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

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

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

    !   Err-arithmetic is performed by: 
    !
    !     1.  Converting to intervals, implicitly or explicitly, 
    !     2.  Performing the required operation
    !     3.  Rounding outwards
    !     4.  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
    integer, parameter     ::  INFO    = 1,               &
                               WARN    = 2,               &
                               ERROR   = 3,               &
                               FATAL   = 4

    !   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 too tempting.
    !
    !   Conversion to real is not supported on purpose.

    public  :: assignment(=)

    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

    public  :: operator(.gt.), operator(.ge.),    &
               operator(.lt.), operator(.le.)

    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

    public  :: operator(+), operator(-),     &
               operator(*), operator(/)

    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

    public  :: operator(**)

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

    public  :: abs,  sqrt, exp, log,      &
               sin,  cos,  tan,           &
               atan, asin, acos

    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

    public  :: err_read_err,   err_write_err,    &
               err_write_txt,  err_write_eor,    &
               err_write_tl,   err_write_ter

  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(ERROR, '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(ERROR, '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(WARN, '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

    !   SIN is our primitive.  We will reduce COS and TAN to SIN.  
    !   There is an alternative primitive implementation of TAN,
    !   but it is not used.

    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 >= TWOPI) then
        err_sin = err_real(ZERO, ONE)
      else
        x = a
        offset = x%lower - mod(x%lower, TWOPI)
        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(HALFPI - 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 >= TWOPI) then
        err_alt_tan = err_real(ZERO, huge(ZERO))
        call err_error(WARN, '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(WARN, '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) = HALFPI

    FUNCTION ERR_ACOS(A)
      type(err_real)              :: err_acos
      type(err_real), intent(in)  :: a

      err_acos = HALFPI - 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

    !   We are trying to introduce here a new convention:
    !
    !     Name           Type          Format
    !     =========      ========      ======
    !     Intervals      IVL_REAL      (x,y)
    !     Abs error      ERR_REAL      <x,y>
    !     Rel error      REL_REAL      [x,y]

    !   Reading one err-real without advancing

    SUBROUTINE ERR_READ_ERR(LUN, A, PROMPT)
      integer                               :: lun, l, l1, l2, i, j
      type(err_real)                        :: a
      character(len=*), optional, volatile  :: prompt
      character(len=80)                     :: buf1, buf2
      character(len=1)                      :: c, begin, end

      if (present(prompt))   &
         write(unit=*, fmt='(1x,a)', advance='NO') prompt
      buf1 = ' '
      read(unit=lun, fmt='(a)')  buf1
      buf2 = ' '
      l = len_trim(buf1)
      j = 1
      do i = 1, l
        c = buf1(i:i)
        if (c .ne. ' ') then
          buf2(j:j) = c
          j = j + 1
        end if
      end do
      begin = buf2(1:1)
      if (begin .ne. '<')      &
         call err_error(ERROR, 'err_read:  number should begin with "<"')
      l1 = index(buf2, ',')
      if (l1 .eq. 0)           &
         call err_error(ERROR, 'err_read:  no comma separator found')
      l2 = len_trim(buf2)
      end = buf2(l2:l2)
      if (end .ne. '>')        &
         call err_error(ERROR, 'err_read:  number should end with ">"')
      read(unit=buf2(2:l1 - 1),      fmt=*) a%number
      read(unit=buf2(l1 + 1:l2 - 1), fmt=*) a%abserr
    END SUBROUTINE ERR_READ_ERR

    !   Writing one err-real without advancing
    !   Prepends a space, may be eaten by carriage control.  

    SUBROUTINE ERR_WRITE_ERR(LUN, PRECISION, A)
      integer                :: lun, precision, length
      type(err_real)         :: a
      character(len = 80)    :: rtf
      character(len=3)       :: l, p

      length = 2 + precision + 2 + 3
      write (unit=l, fmt='(i3)') length
      write (unit=p, fmt='(i3)') precision
      rtf = '(1x,a1,es' // l // '.' // p // 'e3' //        &
               ',a1,es' // l // '.' // p // 'e3' // ',a1)'
      write (unit=lun, fmt=rtf, advance='NO')   &
            '<', a%number, ',', a%abserr, '>'
      return
    END SUBROUTINE ERR_WRITE_ERR

    !   Prepends a space, may be eaten by carriage control.  

    SUBROUTINE ERR_WRITE_TXT(LUN, TEXT)
      integer              :: lun
      character(len=*)     :: text

      write (unit=lun, fmt='(1x,a)', advance='NO') text
      return
    END SUBROUTINE ERR_WRITE_TXT

    !   New line 

    SUBROUTINE ERR_WRITE_EOR(LUN)
      integer              :: lun

      write (unit=lun, fmt='(1x)')
      return
    END SUBROUTINE ERR_WRITE_EOR

    SUBROUTINE ERR_WRITE_TL(LUN, TEXT)
      integer              :: lun
      character(len=*)     :: text

      call err_write_txt(lun, text)
      call err_write_eor(lun)
      return
    END SUBROUTINE ERR_WRITE_TL

    SUBROUTINE ERR_WRITE_TER(LUN, PRECISION, TEXT, A)
      integer              :: lun, precision
      character(len=*)     :: text
      type(err_real)       :: a

      call err_write_txt(lun, text)
      call err_write_err(lun, precision, a)
      call err_write_eor(lun)
      return
    END SUBROUTINE ERR_WRITE_TER


    !   E r r o r   h a n d l i n g

    SUBROUTINE ERR_ERROR(N, STRING)
      integer,           intent(in)   :: N
      character(len =*), intent(in)   :: string
      integer, save                   :: nwarn = 0, nerr = 0

      select case (N)
      case (INFO)
        write(*,*) string
      case (WARN)
        write(*,*) string
        nwarn = nwarn + 1
      case (ERROR)
        write(*,*) string
        nerr  = nerr  + 1
      case (FATAL)
        write(*,*) string
        stop 'Fatal error'
      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


    !  The following simple test program can be used to demonstrate 
    !  a major weakness of our simple-minded approach.  
    !
    !  Suppose we want to compute the expression:  z = x / x, 
    !  given the value of x.
    !
    !  Since the numerator and denominator are the same variable,
    !  whatever value one of them has, the other has the same one,
    !  and the answer should be something like this:  <1.0,0.0>
    !  
    !  Our approach implicitly assumes that all operation arguments
    !  are independent, which is obviously not true if variables 
    !  occur more than once in any expression.
    !
    !  This example is very artificial, but  z = (x + 1) / (x + 2)
    !  is not, and it suffers from the same problem.
    !

 program test
   use num_consts                            !  Numerics module
   use err_arithmetic                        !  Error analytic module
   type(err_real) :: x, y, z                 !  Declare our vars

   call num_check()                          !  Verify we have PRECISE reals

   call err_write_tl(6, "*** z = x / x")     !  Write a line of text
   x = err_real(1.0, 0.1)                    !  Assign:  x = <1.0,0.1>
   call err_write_ter(6, 5, "x   = ", x)     !  Write X
   y = x                                     !  Assign:  y = x
   call err_write_ter(6, 5, "y   = ", y)     !  Write Y
   z = x / y                                 !  Error analytic division
   call err_write_ter(6, 5, "x/y = ", z)     !  Write result
   call err_write_eor(6)                     !  Write a newline

   call err_write_tl(6, "*** singularity")   !  Write a line of text
   call err_write_ter(6, 5, "x   = ", x)     !  Write X
   y = err_real(0.0, 0.1)                    !  Assign:  y = <0.0,0.1>
   call err_write_ter(6, 5, "y   = ", y)     !  Write Y
   z = x / y                                 !  Error analytic division
   call err_write_ter(6, 5, "x/y = ", z)     !  Write result
   call err_write_eor(6)                     !  Write a newline

   call err_write_tl(6, "*** trigo... ")     !  Write a line of text
   x = err_real(0.0, 0.1)                    !  Assign:  y = <0.0,0.1>
   call err_write_ter(6, 5, "x   = ", x)     !  Write X
   z = sin(x)                                !  Error analytic sin(x)
   call err_write_ter(6, 5, "sin(x) = ", z)  !  Write sin(x)
   z = cos(x)                                !  Error analytic cos(x)
   call err_write_ter(6, 5, "cos(x) = ", z)  !  Write cos(x)
   z = tan(x)                                !  Error analytic tan(x)
   call err_write_ter(6, 5, "tan(x) = ", z)  !  Write tan(x)
   call err_write_eor(6)                     !  Write a newline

   call err_write_tl(6, "*** Interactive")   !  Write a line of text
   call err_read_err(5, x, "Enter x: ")      !  Input X
   call err_read_err(5, y, "Enter y: ")      !  Input Y
   z = x / y                                 !  Error analytic division
   call err_write_ter(6, 5, "x/y = ", z)     !  Write result
   call err_write_eor(6)                     !  Write a newline

 end program test