Go to Qiita Advent Calendar 2024 Top
LoginSignup
7
2

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 5 years have passed since last update.

FortranAdvent Calendar 2018Day 16
@implicit_none

Fortranで区間演算

Last updated at Posted at 2018-12-15

概要

先日の記事において，Fortranで浮動小数点の丸めモードを変更する方法を述べました．浮動小数点の丸めモードが変更できると，精度保証付き数値計算の基礎となる区間演算が行えるようになります1．

この記事では，Fortranのオブジェクト指向プログラミングの機能を用いて，区間演算を簡単に実行できる型を作ります．区間演算については，最近出版された書籍2を参考にしてください．

使用環境

コンパイラ バージョン
intel Parallel Studio XE Composer Edition for Fortran 17.0.4.210
gfortran 7.3.0 (for Ubuntu on WSL)

PGI Fortranでは，ユーザ定義派生型入出力という，本記事にとって本質的でないところでコンパイルエラーが出て，簡単には修正できなかったので使用しませんでした．

更新履歴

  • 2018年12月16日
    • ソースのスペーシングを変更
    • ソースにコメントを追加
    • 除数の区間が0を含んでいたら計算を止めるように変更

区間演算

区間演算2は，区間とよばれる上端と下端の集合として表現し，その区間に対して演算を行う方法です．浮動小数点のように一つの値で実数を近似的に表現する方法とは異なります．浮動小数点を用いた演算は，丸め誤差の蓄積や情報落ちなどで間違った答えを返す場合があります．一方で，区間演算は真の値を含んだ区間として答えが得られ，少なくとも間違った答えを返すことはありません．

区間$\boldsymbol{a}$は，下端$\underline{a}$と上端$\overline{a}$の集合$[\underline{a},\overline{a}]$として定義されます．下端・上端に浮動小数点を用いる場合には機械区間とよばれます2．

機械区間を用いた機械区間演算が，次のように定義されています2．

\begin{eqnarray}
    \boldsymbol{a}+\boldsymbol{b} &\subset& [\rm{fl_\nabla}(\underline{a}+\underline{b}), \rm{fl_\Delta}(\overline{a}+\overline{b})] \quad\quad(1)\\

    \boldsymbol{a}-\boldsymbol{b} &\subset& [\rm{fl_\nabla}(\underline{a}-\overline{b}), \rm{fl_\Delta}(\overline{a}-\underline{b})] \quad\quad(2)\\

    \boldsymbol{a}\times\boldsymbol{b} &\subset& [\rm{fl_\nabla}(\min(\underline{a}\underline{b},\underline{a}\overline{b},\overline{a}\underline{b},\overline{a}\overline{b})), \rm{fl_\Delta}(\max(\underline{a}\underline{b},\underline{a}\overline{b},\overline{a}\underline{b},\overline{a}\overline{b}))] \quad\quad(3)\\
    \frac{\boldsymbol{a}}{\boldsymbol{b}} &\subset& \left[\rm{fl_\nabla}\left(\min\left(\frac{\underline{a}}{\underline{b}},\frac{\underline{a}}{\overline{b}},\frac{\overline{a}}{\underline{b}},\frac{\overline{a}}{\overline{b}}\right)\right), \rm{fl_\Delta}\left(\max\left(\frac{\underline{a}}{\underline{b}},\frac{\underline{a}}{\overline{b}},\frac{\overline{a}}{\underline{b}},\frac{\overline{a}}{\overline{b}}\right)\right)\right] \quad (0 \notin \boldsymbol{b}) \quad\quad(4)\\
\end{eqnarray}

$\rm{fl_\nabla()}$は，カッコ内の演算を切り下げのモードで実行することを，$\rm{fl_\Delta()}$は，切り上げのモードで実行することを表しています．

話だけを聞くとややこしそうに思いますが，データ構造と演算がわかっていればプログラムを作るのはさほど難しくありません．

Fortranによる区間演算の実装

派生型

データとして下端と上端をセットで取り扱うので，派生型を使うとプログラムが簡単になりそうです．

モジュールを作成し，その中で倍精度の浮動小数点を成分に持つ派生型intervalFP64と，intervalFP64に対して演算を行う手続（関数とサブルーチン）を定義します．

成分は下端と上端のみでよいでしょう．演算以外で値が変わっては困るので，private属性を付与します．

    type, public :: intervalFP64
        real(real64),private :: inf ! 下端
        real(real64),private :: sup ! 上端
    end type intervalFP64

演算

式(1)~(4)に沿って演算を行う関数を実装します．難しいことは何もありません．

加算

    function addIntervalFP64(a,b) result(add)
        implicit none
        type(intervalFP64),intent(in) :: a
        type(intervalFP64),intent(in) :: b
         
        type(intervalFP64) :: add ! 戻り値
        
        call ieee_set_rounding_mode(ieee_down   ); add%inf = a%inf + b%inf ! 下端
        call ieee_set_rounding_mode(ieee_up     ); add%sup = a%sup + b%sup ! 上端
        call ieee_set_rounding_mode(ieee_nearest); return
    end function addIntervalFP64

減算

    function subtractIntervalFP64(a,b) result(sub)
        implicit none
        type(intervalFP64),intent(in) :: a
        type(intervalFP64),intent(in) :: b
         
        type(intervalFP64) :: sub ! 戻り値
        
        call ieee_set_rounding_mode(ieee_down   ); sub%inf = a%inf - b%sup ! 下端
        call ieee_set_rounding_mode(ieee_up     ); sub%sup = a%sup - b%inf ! 上端
        call ieee_set_rounding_mode(ieee_nearest); return
    end function subtractIntervalFP64

乗算

    function multiplyByIntervalFP64(a,b) result(mul)
        implicit none
        type(intervalFP64),intent(in) :: a
        type(intervalFP64),intent(in) :: b
        
        type(intervalFP64) :: mul ! 戻り値
        
        !下端
        call ieee_set_rounding_mode(ieee_down)
        mul%inf = min(a%inf*b%inf, a%inf*b%sup, a%sup*b%inf, a%sup*b%sup)

        !上端
        call ieee_set_rounding_mode(ieee_up)
        mul%sup = max(a%inf*b%inf, a%inf*b%sup, a%sup*b%inf, a%sup*b%sup)
        
        call ieee_set_rounding_mode(ieee_nearest)
    end function multiplyByIntervalFP64

除算

    function divideByIntervalFP64(a,b) result(div)
        implicit none
        type(intervalFP64),intent(in) :: a
        type(intervalFP64),intent(in) :: b
        
        type(intervalFP64) :: div ! 戻り値

        if (    (b%inf <= 0d0 .and. b%sup >= 0d0)&
            .or.(b%inf >= 0d0 .and. b%sup <= 0d0)) stop "error : divide by zero in interval operation"

        !下端
        call ieee_set_rounding_mode(ieee_down)
        div%inf = min(a%inf/b%inf, a%inf/b%sup, a%sup/b%inf, a%sup/b%sup)
            
        !上端
        call ieee_set_rounding_mode(ieee_up)
        div%sup = max(a%inf/b%inf, a%inf/b%sup, a%sup/b%inf, a%sup/b%sup)

        call ieee_set_rounding_mode(ieee_nearest)
    end function divideByIntervalFP64

平方根

平方根は，浮動小数点の丸めモードを変更すると，そのモードに従って丸められることが保証されている2ようです．そのため，区間intervalFP64の平方根は，丸めモードを変更して，下端と上端の平方根を計算します．

    function dIntervalSqrt(a) result(sqrt)
        implicit none
        
        type(intervalFP64) :: a
        type(intervalFP64) :: sqrt ! 戻り値

        if(a%inf < 0d0) stop "error: negative value in interval sqrt"

        call ieee_set_rounding_mode(ieee_down   ); sqrt%inf = dsqrt(a%inf)
        call ieee_set_rounding_mode(ieee_up     ); sqrt%sup = dsqrt(a%sup)
        call ieee_set_rounding_mode(ieee_nearest); return
    end function dIntervalSqrt

Fortranのこれまでの関数の命名規則に従って，倍精度を表すdを関数名の先頭につけています．ただし，これでは実数型に対する平方根sqrt()との互換性に欠けるので，総称名を利用して，dIntervalSqrtをsqrtとして呼べるようにします．

    interface sqrt
        procedure :: dIntervalSqrt
    end interface

オブジェクト指向実装

区間を表す型や，それらに対する演算を行う関数は上記のように簡単に作ることができます．しかし，浮動小数点のように演算子を用いて演算を記述できないと，見通しが悪くなり，かつプログラムが煩雑になります．

そこで，オブジェクト指向プログラミングの機能を利用して，ユーザ定義派生型に型束縛手続きを追加し，かつ演算子をオーバーロードします．代入演算子はオーバーロードする必要はありません．

演算は，以下の3パターンを実装しました．

  • invervalFP64 op invervalFP64
  • invervalFP64 op 倍精度浮動小数
  • 倍精度浮動小数 op invervalFP64

加えて，初期化のパターンを考えて4個のコンストラクタを定義し，総称名で統一的に呼び出せるようにします．4個のコンストラクタは，それぞれ

  • 引数なしで呼び出すと0で初期化
  • 倍精度浮動小数一つを引数にすると，その値で下端と上端を初期化
  • 倍精度浮動小数二つを引数にすると，下端と上端を個別に初期化
  • intervalFP64型を引数にすると，その下端と上端の値を引き継いで初期化

を実装しています．

最後に，print/write文でprivateな成分を表示できるように，ユーザ定義派生型の出力用サブルーチンを作成しています．

class_intervalFP64.f90
module class_interval
    use, intrinsic :: iso_fortran_env
    use, intrinsic :: ieee_arithmetic
    implicit none
    private
    public interval
    public sqrt
    public write (formatted)
    
    ! コンストラクタ
    interface interval
        procedure :: initializeByZero,&
                     initializeByValue,&
                     initializeByInterval,&
                     initializeByIntervalFP64
    end interface

    interface write (formatted)
        procedure printIntervalFP64
    end interface

    ! 平方根
    interface sqrt
        procedure :: dIntervalSqrt
    end interface

    type, public :: intervalFP64
        real(real64),private :: inf ! 下端
        real(real64),private :: sup ! 上端
    contains
        procedure,public,pass       :: setInf                      ! 下端のsetter
        procedure,public,pass       :: setSup                      ! 上端のsetter
        procedure,public,pass       :: addIntervalFP64             ! intervalFP64 + intervalFP64 
        procedure,public,pass       :: addFP64                     ! intervalFP64 + FP64
        procedure,public,pass(this) :: addIntervalFP64toFP64       ! FP64 + intervalFP64
        procedure,public,pass       :: subtractIntervalFP64        ! intervalFP64 - intervalFP64 
        procedure,public,pass       :: subtractFP64                ! intervalFP64 - FP64
        procedure,public,pass(this) :: subtractIntervalFP64toFP64  ! FP64 - intervalFP64
        procedure,public,pass       :: multiplyByIntervalFP64      ! intervalFP64 * intervalFP64 
        procedure,public,pass       :: multiplyByFP64              ! intervalFP64 * FP64
        procedure,public,pass(this) :: multiplyFP64ByIntervalFP64  ! FP64 * intervalFP64
        procedure,public,pass       :: divideByIntervalFP64        ! intervalFP64 / intervalFP64 
        procedure,public,pass       :: divideByFP64                ! intervalFP64 / FP64
        procedure,public,pass(this) :: divideFP64ByIntervalFP64    ! FP64 / intervalFP64

        ! 演算子のオーバーロード
        generic :: operator(+) => addIntervalFP64,&
                                  addFP64,&
                                  addIntervalFP64toFP64
        generic :: operator(-) => subtractIntervalFP64,&
                                  subtractFP64,&
                                  subtractIntervalFP64toFP64
        generic :: operator(*) => multiplyByIntervalFP64,&
                                  multiplyByFP64,&
                                  multiplyFP64ByIntervalFP64
        generic :: operator(/) => divideByIntervalFP64,&
                                  divideByFP64,&
                                  divideFP64ByIntervalFP64
    end type intervalFP64

    contains

    !-----------------------------------------------------------------!
    !constructor
    function initializeByZero() result(retValue) ! 0で初期化
        implicit none
       type(intervalFP64) :: retValue
       
       retValue%inf = 0d0
       retValue%sup = 0d0
    end function initializeByZero
    
    function initializeByValue(value) result(retValue) ! 実数で初期化
        implicit none
        real(real64),intent(in) :: value
        type(intervalFP64) :: retValue
       
       retValue%inf = value
       retValue%sup = value
    end function initializeByValue

    function initializeByInterval(inf, sup) result(retValue) ! 下端と上端を個別に初期化
        implicit none
        real(real64),intent(in) :: inf
        real(real64),intent(in) :: sup
        type(intervalFP64) :: retValue
       
       retValue%inf = inf
       retValue%sup = sup
    end function initializeByInterval

    function initializeByIntervalFP64(intvlFP64) result(retValue) ! 引数の区間を引き継いで初期化
        implicit none
        type(intervalFP64),intent(in) :: intvlFP64
        type(intervalFP64) :: retValue
       
       retValue%inf = intvlFP64%inf
       retValue%sup = intvlFP64%sup
    end function initializeByIntervalFP64
    !-----------------------------------------------------------------!
    
    !-----------------------------------------------------------------!
    !setter
    subroutine setInf(this,inf)
        implicit none
        class(intervalFP64),intent(inout) :: this
        real(real64),intent(in) :: inf
        
        this%inf = inf
    end subroutine setInf

    subroutine setSup(this,sup)
        implicit none
        class(intervalFP64),intent(inout) :: this
        real(real64),intent(in) :: sup
        
        this%sup = sup
    end subroutine setSup
    !-----------------------------------------------------------------!
    
    !-----------------------------------------------------------------!
    !add
    function addIntervalFP64(this,intvlFP64) result(add)
        implicit none
        class(intervalFP64),intent(in) :: this
         type(intervalFP64),intent(in) :: intvlFP64
         
         type(intervalFP64) :: add
        
        call ieee_set_rounding_mode(ieee_down   ); add%inf = this%inf + intvlFP64%inf
        call ieee_set_rounding_mode(ieee_up     ); add%sup = this%sup + intvlFP64%sup
        call ieee_set_rounding_mode(ieee_nearest); return
    end function addIntervalFP64
    
    function addFP64(this,FP64) result(add)
        implicit none
        class(intervalFP64),intent(in) :: this
        real(real64),intent(in) :: FP64
         
         type(intervalFP64) :: add
        
        call ieee_set_rounding_mode(ieee_down   ); add%inf = this%inf + FP64
        call ieee_set_rounding_mode(ieee_up     ); add%sup = this%sup + FP64
        call ieee_set_rounding_mode(ieee_nearest); return
    end function addFP64

    function addIntervalFP64toFP64(FP64,this) result(add)
        implicit none
        real(real64),intent(in) :: FP64
        class(intervalFP64),intent(in) :: this
         
         type(intervalFP64) :: add

         add = this + FP64
    end function addIntervalFP64toFP64
    !-----------------------------------------------------------------!

    !-----------------------------------------------------------------!
    !subtract
    function subtractIntervalFP64(this,intvlFP64) result(sub)
        implicit none
        class(intervalFP64),intent(in) :: this
         type(intervalFP64),intent(in) :: intvlFP64
         
         type(intervalFP64) :: sub
        
        call ieee_set_rounding_mode(ieee_down   ); sub%inf = this%inf - intvlFP64%sup
        call ieee_set_rounding_mode(ieee_up     ); sub%sup = this%sup - intvlFP64%inf
        call ieee_set_rounding_mode(ieee_nearest); return
    end function subtractIntervalFP64

    function subtractFP64(this,FP64) result(sub)
        implicit none
        class(intervalFP64),intent(in) :: this
        real(real64),intent(in) :: FP64
         
         type(intervalFP64) :: sub
        
        call ieee_set_rounding_mode(ieee_down   ); sub%inf = this%inf - FP64
        call ieee_set_rounding_mode(ieee_up     ); sub%sup = this%sup - FP64
        call ieee_set_rounding_mode(ieee_nearest); return
    end function subtractFP64

    function subtractIntervalFP64toFP64(FP64,this) result(sub)
        implicit none
        real(real64),intent(in) :: FP64
        class(intervalFP64),intent(in) :: this
         
         type(intervalFP64) :: sub
        
        call ieee_set_rounding_mode(ieee_down   ); sub%inf = min(FP64 - this%inf, FP64 - this%sup)
        call ieee_set_rounding_mode(ieee_up     ); sub%sup = max(FP64 - this%sup, FP64 - this%inf)
        call ieee_set_rounding_mode(ieee_nearest); return
    end function subtractIntervalFP64toFP64
    !-----------------------------------------------------------------!

    !-----------------------------------------------------------------!
    !multiply
    function multiplyByIntervalFP64(this,FP64) result(mul)
        implicit none
        class(intervalFP64),intent(in) :: this
         type(intervalFP64),intent(in) :: FP64
        
        type(intervalFP64) :: mul
        
        call ieee_set_rounding_mode(ieee_down)
        mul%inf = min(this%inf*FP64%inf, this%inf*FP64%sup, this%sup*FP64%inf, this%sup*FP64%sup)
        call ieee_set_rounding_mode(ieee_up)
        mul%sup = max(this%inf*FP64%inf, this%inf*FP64%sup, this%sup*FP64%inf, this%sup*FP64%sup)
        call ieee_set_rounding_mode(ieee_nearest)
    end function multiplyByIntervalFP64
    
    function multiplyByFP64(this,FP64) result(mul)
        implicit none
        class(intervalFP64),intent(in) :: this
        real(real64),intent(in) :: FP64
        
        type(intervalFP64) :: mul

        call ieee_set_rounding_mode(ieee_down)
        mul%inf = min(this%inf*FP64, this%sup*FP64)
        call ieee_set_rounding_mode(ieee_up)
        mul%sup = max(this%inf*FP64, this%sup*FP64)
        call ieee_set_rounding_mode(ieee_nearest)
    end function multiplyByFP64

    function multiplyFP64ByIntervalFP64(FP64,this) result(mul)
        implicit none
        real(real64),intent(in) :: FP64
        class(intervalFP64),intent(in) :: this
         
         type(intervalFP64) :: mul
        
        mul = this * FP64
    end function multiplyFP64ByIntervalFP64
    !-----------------------------------------------------------------!

    !-----------------------------------------------------------------!
    !division
    function divideByIntervalFP64(this, FP64) result(div)
        implicit none
        class(intervalFP64),intent(in) :: this
         type(intervalFP64),intent(in) :: FP64
        
        type(intervalFP64) :: div

        if (    (FP64%inf <= 0d0 .and. FP64%sup >= 0d0)&
            .or.(FP64%inf >= 0d0 .and. FP64%sup <= 0d0)) stop "error : divide by zero in interval operation"

        call ieee_set_rounding_mode(ieee_down)
        div%inf = min(this%inf/FP64%inf, this%inf/FP64%sup, this%sup/FP64%inf, this%sup/FP64%sup)
        call ieee_set_rounding_mode(ieee_up)
        div%sup = max(this%inf/FP64%inf, this%inf/FP64%sup, this%sup/FP64%inf, this%sup/FP64%sup)
        call ieee_set_rounding_mode(ieee_nearest)
        
    end function divideByIntervalFP64
        
    function divideByFP64(this, FP64) result(div)
        implicit none
        class(intervalFP64),intent(in) :: this
        real(real64),intent(in) :: FP64
        
        type(intervalFP64) :: div

        if (FP64 /= 0d0) then
            call ieee_set_rounding_mode(ieee_down)
            div%inf = min(this%inf/FP64, this%sup/FP64)
            call ieee_set_rounding_mode(ieee_up)
            div%sup = max(this%inf/FP64, this%sup/FP64)
            call ieee_set_rounding_mode(ieee_nearest)
            return
        end if
        stop "error : divide by zero in interval operation"
    end function divideByFP64

    function divideFP64ByIntervalFP64(FP64,this) result(div)
        implicit none
        class(intervalFP64),intent(in) :: this
        real(real64),intent(in) :: FP64
        
        type(intervalFP64) :: div

        if (    (this%inf <= 0d0 .and. this%sup >= 0d0)&
            .or.(this%inf >= 0d0 .and. this%sup <= 0d0)) stop "error : divide by zero in interval operation"

        call ieee_set_rounding_mode(ieee_down)
        div%inf = min(FP64/this%inf, FP64/this%sup)
        call ieee_set_rounding_mode(ieee_up)
        div%sup = max(FP64/this%inf, FP64/this%sup)
        call ieee_set_rounding_mode(ieee_nearest)
    end function divideFP64ByIntervalFP64
    !-----------------------------------------------------------------!
    
    !-----------------------------------------------------------------!
    !square root
    function dIntervalSqrt(intvlFP64) result(sqrt)
        implicit none
        
        type(intervalFP64) :: intvlFP64
        type(intervalFP64) :: sqrt

        if(intvlFP64%inf < 0d0) stop "error: negative value in interval sqrt"

        call ieee_set_rounding_mode(ieee_down   ); sqrt%inf = dsqrt(intvlFP64%inf)
        call ieee_set_rounding_mode(ieee_up     ); sqrt%sup = dsqrt(intvlFP64%sup)
        call ieee_set_rounding_mode(ieee_nearest); return
    end function dIntervalSqrt
    !-----------------------------------------------------------------!

    !-----------------------------------------------------------------!
    subroutine printIntervalFP64(this, Unit, IOType, argList, IOStatus, IOMessage)
        implicit none
        class(intervalFP64),intent(in   ) :: this
        integer            ,intent(in   ) :: Unit !writeで指定した装置番号，*の場合は-1
        character(*)       ,intent(in   ) :: IOType !*の時はLISTDIRECTED
        integer            ,intent(in   ) :: argList(:)
        integer            ,intent(  out) :: IOStatus
        character(*)       ,intent(inout) :: IOMessage

        if(IOType == "LISTDIRECTED" .or. size(argList)<2)then !書式指定なしか，書式指定の数字が不足している場合
            write(unit=Unit, fmt = '(A,E27.19,A,E27.19,A)', iostat = IOStatus, iomsg = IOMessage) '[',this%inf,',',this%sup,']'
            return
        else
            if(IOType(3:) /= "intervalFP64")then
                print *,"error : type mismatch"
                return
            end if
            block
                character(2) :: width_tol, width_dec
                character(6) :: RealSpec
                character(:),allocatable :: fmt
                character(4),parameter :: comma = "', '"
            
                write(width_tol,'(I2)') argList(1)
                write(width_dec,'(I2)') argList(2)
                RealSpec = 'E'//width_tol//'.'//width_dec
                fmt = "(' ['"//RealSpec//","//comma//","//RealSpec//"']')"
            
                write(unit=Unit, fmt = fmt, iostat = IOStatus, iomsg = IOMessage) this%inf,this%sup
            end block
            return
        end if
        
    end subroutine printIntervalFP64
    !-----------------------------------------------------------------!
end module class_interval

初期化

初期値をどう与えるか非常に重要なはずですが，論文や書籍ではそういった実用的な項目は解説されません．どうするのが正しいのかわからなかったため，実数型の変数を用意して，丸めモードを変えて演算することで下端と上端の値を設定しました．

    type(intervalFP64) :: a
    real(real64) :: x,y
    
    x = 1d0
    y = 10d0

    ! 下端
    call ieee_set_rounding_mode(ieee_down)
    call a%setInf(x/y)
    ! 上端
    call ieee_set_rounding_mode(ieee_up)
    call a%setSup(x/y)
    call ieee_set_rounding_mode(ieee_nearest)

しかし，初期値が単純な演算で表せない場合にどうすればよいか，さっぱりわかりません．
最初は区間を考えずに下端と上端を同じとすればよいのかもしれませんが，そういった記述は見つかりませんでした．

メインルーチン

愚痴を言っても実装は生えてこないので，適当に与えた初期値で区間演算を実行してみました．

main.f90
program main
    use,intrinsic :: iso_fortran_env
    use,intrinsic :: ieee_arithmetic
    use class_interval
    implicit none
    
    type(intervalFP64) :: a,b
    real(real64) :: x,y
    
    x = 1d0
    y = 10d0

    ! 下端
    call ieee_set_rounding_mode(ieee_down)
    call a%setInf(x/y)
    ! 上端
    call ieee_set_rounding_mode(ieee_up)
    call a%setSup(x/y)
    call ieee_set_rounding_mode(ieee_nearest)

    b = a ! 全成分のコピー

    ! 書式なし表示
    print *,"a = ",a
    print *,"b = ",b

    ! 区間 op 区間
    print *,"a+b = ",a+b
    print *,"a-b = ",a-b
    print *,"a*b = ",a*b
    print *,"a/b = ",a/b

    print *,"sqrt(a) = ", sqrt(a)

    ! 書式付き表示
    print '(dt"intervalFP64"(28,20))',a+b

    ! 区間 op 実数
    print *,"a+x = ",a+x
    print *,"a-x = ",a-x
    print *,"a*x = ",a*x
    print *,"a/x = ",a/x

    ! 実数 op 区間
    print *,"y+a = ",y+a
    print *,"y-a = ",y-a
    print *,"y*a = ",y*a
    print *,"y/a = ",y/a
end program main

これが正しいかどうか自信はありませんが，すべての演算において，真の値を内包する区間が得られています．

書式なし表示
 a = [  0.9999999999999999167E-01,  0.1000000000000000056E+00]
 b = [  0.9999999999999999167E-01,  0.1000000000000000056E+00]

区間 op 区間
 a+b = [  0.1999999999999999833E+00,  0.2000000000000000111E+00]
 a-b = [ -0.1387778780781445676E-16,  0.1387778780781445676E-16]
 a*b = [  0.9999999999999996739E-02,  0.1000000000000000194E-01]
 a/b = [  0.9999999999999997780E+00,  0.1000000000000000222E+01]
 sqrt(a) = [  0.3162277660168378857E+00,  0.3162277660168379967E+00]

書式付き表示
 [  0.19999999999999998335E+00,   0.20000000000000001110E+00]
 
区間 op 実数
 a+x = [  0.1099999999999999867E+01,  0.1100000000000000089E+01]
 a-x = [ -0.9000000000000000222E+00, -0.8999999999999999112E+00]
 a*x = [  0.9999999999999999167E-01,  0.1000000000000000056E+00]
 a/x = [  0.9999999999999999167E-01,  0.1000000000000000056E+00]

実数 op 区間
 y+a = [  0.1009999999999999964E+02,  0.1010000000000000142E+02]
 y-a = [  0.9899999999999998579E+01,  0.9900000000000000355E+01]
 y*a = [  0.9999999999999998890E+00,  0.1000000000000000222E+01]
 y/a = [  0.9999999999999998579E+02,  0.1000000000000000142E+03]

まとめ

実数の丸めモードの変更とオブジェクト指向プログラミングを利用して，区間演算を実装しました．型とそれに対する演算だけを実装するのは非常に簡単なので，オブジェクト指向プログラミングの練習にもよいと思います．

本来は区間演算を行うのが目的ではなく，これを使って精度保証を行うのですが，それはFortranとは別の話になってくるので本記事では取り上げません．

  1. 先日の記事を書いたのは，この記事を書くための布石でした． ↩

  2. 大石進一 編著，精度保証付き数値計算の基礎，コロナ社，2018 ↩ ↩2 ↩3 ↩4 ↩5

7
2
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
7
2

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?