! 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
! 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
Return to contents page