Ziele Untersuchung
mit Columbo Integrität von
Datenbanken Interaktion und
Portierbarkeit Ergonomie der
Schnittstellen

Angebot Produkte Projekt Beratung

Mittel Analytik Modellierung Sprachen Algebra Logik Hardware Denken Kreativität

Zusammenhänge Gesellschaft Wirtschaft Branche Firma


products/sources/formale Sprachen/Fortran/ Datei vom 4.1.2008 mit Größe 38 kB

Quelle Arithmetic.ftn Sprache: Fortran

    !   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
    !         "write_err" and "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  :: read_err,   write_err,    &
               write_txt,  write_nl,    &
               write_txt_nl,   write_one

  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 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 READ_ERR

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

    SUBROUTINE 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 WRITE_ERR

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

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

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

    !   New line

    SUBROUTINE WRITE_NL(LUN)
      integer              :: lun

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

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

      call write_txt(lun, text)
      call write_nl(lun)
      return
    END SUBROUTINE WRITE_TXT_NL

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

      call write_txt(lun, text)
      call write_err(lun, precision, a)
      call write_nl(lun)
      return
    END SUBROUTINE WRITE_ONE

    !   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

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 write_txt_nl(6, "*** z = x / x")     !  Write a line of text
   x = err_real(1.0, 0.1)                    !  Assign:  x = <1.0,0.1>
   call write_one(6, 5, "x = ", x)         !  Write X
   y = x                                     !  Assign:  y = x
   call write_one(6, 5, "y = ", y)         !  Write Y
   z = x / y                                 !  Error analytic division
   call write_one(6, 5, "x/y = ", z)         !  Write result
   call write_nl(6)                          !  Write a newline

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

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

   call write_txt_nl(6, "*** Interactive")   !  Write a line of text
   call read_err(5, x, "Enter x: ")          !  Input X
   call read_err(5, y, "Enter y: ")          !  Input Y
   z = x / y                                 !  Error analytic division
   call write_one(6, 5, "x/y = ", z)         !  Write result
   call write_nl(6)                          !  Write a newline

end program test

quality99%

¤ Dauer der Verarbeitung: 0.9 Sekunden (vorverarbeitet) ¤

Normalansicht

Suchen

Beweissystem der NASA

Beweissystem Isabelle

NIST Cobol Testsuite

Cephes Mathematical Library

Wiener Entwicklungsmethode

Haftungshinweis

Die Informationen auf dieser Webseite wurden nach bestem Wissen sorgfältig zusammengestellt. Es wird jedoch weder Vollständigkeit, noch Richtigkeit, noch Qualität der bereit gestellten Informationen zugesichert.