test

Sun Studio 12: Fortran プログラミングガイド

第 6 章 浮動小数点演算

この章では、浮動小数点演算について説明し、数値計算の誤差を回避および検出するための方針を示します。

SPARC および x86 プロセッサ上の浮動小数点計算についての詳細は、『数値計算ガイド』を参照してください。

6.1 はじめに

SPARC プロセッサ上での Fortran 95 の浮動小数点環境は、『IEEE Standard 754 for Binary Floating Point Arithmetic』で指定された演算モデルを実装しています。この環境では、安定した、高性能の、移植可能な数値計算アプリケーションを開発できます。また、数値計算プログラムによる異常な動作を調査するためのツールも提供します。

数値計算プログラムでは、計算誤差を引き起こす次のような潜在的な原因がたくさんあります。

間違った数値計算の誤差の原因を見つけることはたいへん困難です。市販の検査済みライブラリパッケージを使用することで、コーディングの誤りの可能性を減らすことはできます。アルゴリズムを選択することも重要な問題です。また、コンピュータの特性を考慮した優れた算術演算法を使用することも重要です。

この章では、数値誤差の解析については説明していません。ここでは、Fortran 95 によって実装された IEEE 浮動小数点モデルを紹介します。

6.2 IEEE 浮動小数点演算

IEEE 演算は、無効演算、ゼロ除算、オーバーフロー、アンダーフロー、結果不正確などの問題を発生させる算術演算を取り扱う、比較的新しい方式です。丸め、ゼロに近い数の扱い方、その機種で取り扱える最大数に近い数の扱い方に違いがあります。

IEEE 規格は、例外、丸め、精度のユーザー処理をサポートします。その結果、IEEE 規格は、区間演算と変則性の診断をサポートします。 IEEE 規格 754 は、exp や cos などの基本関数の標準化、高精度の演算、数値計算と記号代数計算の結合を実現します。

IEEE 演算では、ほかの浮動小数点演算と比べると、ユーザーが計算を大きく制御できます。IEEE 規格は、数値計算的に洗練された移植性のあるプログラムを作成する作業を簡単にします。浮動小数点演算についての多くの質問は、基本的な数の演算に関連します。たとえば、次のようにします。

もう 1 つのクラスの質問は、浮動小数点例外や例外処理に関連するものです。たとえば、次のとおりです。

旧式の演算モデルでは、最初のクラスの質問は期待された答えにならず、2 番目のク ラスの例外ケースはすべて同じ結果になります。つまり、プログラムは該当する箇所で異常終了するか、演算は続行されるが結果は意味のないものとなります。

IEEE 規格は、演算が数学的に期待される結果を、期待される特性で出すことを保証します。また、ユーザーが特にほかを選択しないかぎり、例外ケースが指定された結果を出すことも保証します。

たとえば、例外値 +Inf、-Inf、NaN は次のように直感的に理解できます。

big*big = +Inf 正の無限大

big*(-big) = -Inf 負の無限大

num/0.0 = +Inf num > 0.0 のとき

num/0.0 = -Inf num < 0.0 のとき

0.0/0.0 = NaN 非数

また、次のような 5 つの種類の浮動小数点例外も発生します。

IEEE 規格の実装については、『数値計算ガイド』を参照してください。

6.2.1 -ftrap=mode コンパイラオプション

-ftrap=mode オプションは、浮動小数点例外用のトラップを有効にします。ieee_handler() 呼び出しによってシグナルハンドラが設定されていない場合、例外はプログラムを終了し、メモリーダンプコアファイルを作成します。このコンパイラオプションに関する詳細は、『Fortran ユーザーズガイド』を参照してください。たとえば、オーバーフロー、ゼロ除算、無効演算のトラップを有効にするには、-ftrap=common を付けてコンパイルします。これが f95 のデフォルトです。

注 –

トラップを有効にするには、アプリケーションの主プログラムを -ftrap= を付けてコンパイルする必要があります。

6.2.2 浮動小数点演算の例外

f95 プログラムは例外を自動的に報告しません。プログラムの終了時に、発生した浮動小数点演算の例外リストを表示するには、ieee_retrospective(3M) を明示的に呼び出す必要があります。一般的には、無効演算、ゼロ除算、オーバーフローのいずれかの例外が発生すると、メッセージが表示されます。実際のプログラムでは結果不正確の例外は多く発生するため、結果不正確の例外のメッセージは表示されません。

6.2.2.1 発生した例外の通知

ieee_retrospective 関数は、浮動小数点のステータスレジスタを調べて、どの例外が発生したのかを突き止め、標準エラーにメッセージを表示し、どの例外が発生してクリアーされていないのかをプログラマに知らせます。 メッセージは通常、次のようになります (リリースによって若干異なります)。


注意: 次の IEEE 浮動小数点例外フラグが立っています:
ゼロ除算;
次の IEEE 浮動小数点例外トラップが有効になっています:
    inexact;  underflow;  overflow;  invalid operation;
『数値計算ガイド』、ieee_flags(3M) のマニュアルページを参照してください
    ieee_handler(3M)

Fortran 95 プログラムでは、ieee_retrospective を明示的に呼び出し、-xlang=f77 を使用してコンパイルして、f77 互換性ライブラリを使用してリンクする必要があります。

-f77 互換性フラグを使用してコンパイルすると、プログラムの終了時に自動的に ieee_retrospective を呼び出す FORTRAN 77 の規則が有効になります。

ieee_retrospective を呼び出す前に例外ステータスフラグをクリアーすると、ieee_flags() を使用して、メッセージの一部また全部を表示させないようにすることができます。

6.2.3 例外処理

IEEE 規格準拠の例外処理は、SPARC および x86 プロセッサ上ではデフォルトで行われます。ただし、浮動小数点例外の検出と、浮動小数点例外に対するシグナル (SIGFPE) の生成には、違いがあります。

浮動小数点演算中にトラップしていない例外が発生すると、IEEE 規格に従って、次の 2 つの処理が行われます。

6.2.4 浮動小数点演算の例外のトラップ

浮動小数点演算の例外を処理する方法は、f95 と以前の f77 では大きく異なります。

f95 のデフォルトでは、ゼロ除算、オーバーフロー、無効演算の場合に自動的にトラップします。f77 のデフォルトでは、浮動小数点演算例外において実行プログラム を中断するためのシグナルを自動的に生成しません。これは、トラップはパフォーマンスを低下させるということと、期待された値が戻ってくれば、ほとんどの例外は重要でないという前提に基づいていました。

f95 のコマンド行オプション -ftrap を使用すると、このデフォルトを変更できます。f95 のデフォルトは -ftrap=common です。以前の f77 のデフォルトに従うには、-ftrap=%none を使用して主プログラムをコンパイルします。

6.2.5 非標準の算術演算

段階的なアンダーフローと呼ばれる、標準の IEEE 算術演算の 1 つは手作業で無効にできます。無効にしたとき、プログラムは非標準の算術演算で実行していると考えられます。

算術演算に関する IEEE 規格では、アンダーフローとなった結果は、有効桁の小数部の基点を動的に調整することにより、段階的に扱うように指定しています。 IEEE の浮動小数点の形式では、有効桁の前に小数点が現れ、暗黙的な 1 の先行ビットがあります。段階的なアンダーフローでは、浮動小数点の演算結果がアンダーフローとなったときに、段階的なアンダーフローでは、暗黙的な先行ビットをゼロにクリアーし、小数点を有効桁方向にシフトさせるようになっています。これは、SPARC プロセッサではハードウェアではなく、ソフトウェアで行われます。このため、プログラムにアンダーフローが多数発生すると (アルゴリズムに問題があることを示す可能性がある)、パフォーマンスが低下することになります。

段階的なアンダーフローを無効にするには、-fns オプションを付けてコンパイルします。または、プログラムの中からライブラリルーチン nonstandard_arithmetic() を呼び出して、段階的なアンダーフローを無効にします。standard_arithmetic() を呼び出して、段階的なアンダーフローを有効に戻します。

注 –

有効にするためには、アプリケーションの主プログラムを -fns を使用してコンパイルする必要があります。『Fortran ユーザーズガイド』を参照してください。

古いアプリケーションの場合、次のことに注意してください。

6.3 IEEE ルーチン

次のインタフェースは、IEEE 演算を使用するユーザーの役に立ちます。これらのほとんどは、数学ライブラリ libsunmath と、いくつかの .h ファイルに置かれています。

6.3.1 フラグと ieee_flags()

ieee_flags 関数は、例外ステータスフラグの問い合わせやクリアーに使用します。この関数は、Sun コンパイラとともに出荷される libsunmath ライブラリに組み込まれていて、次の作業を行うことができます。

ieee_flags の一般的な呼び出し方法は次のとおりです。


      
flags = ieee_flags( action, mode, in, out )

4 つの引数はすべて文字列です。入力は、action、mode、および in です。出力は、out と flags です。ieee_flags は、整数値の関数です。flags は 1 ビットフラグの集合で、ここには有用な情報が返されます。詳細は、ieee_flags(3m) のマニュアルページを参照してください。

パラメータの取り得る値は次のとおりです。

表 6–1 ieee_flags( action, mode, in, out ) の引数の値

引数の 

使用可能な値 

action

get、set、clear、clearall のいずれか

mode

direction、exception

in、out

nearest、tozero、negative、positive、extended、double single、inexact、division、underflow、overflow、invalid all、common

これらはリテラル文字列であること、出力パラメータ out は最低でも CHARACTER*9 である必要があることに注目してください。in と out の取り得る値の意味は、使用時の動作とモードによって異なります。これらの関係を次の表に要約します。

表 6–2 ieee_flags の引数 in と out の意味

in と out の値

意味 

nearest、tozero、negative、positive

丸めの方向 

extended、double、single

丸めの精度 

inexact、division、underflow、overflow、invalid

例外 

all

5 種類すべての例外 

common

一般的な例外: invalid (無効演算)、division (ゼロ除算)、overflow (オーバーフロー) 

たとえば、フラグが立てられている例外のうちで何がもっとも優先順位が高いかを判別するときは、引数の in に NULL 文字列を渡します。


      CHARACTER *9, out
      ieeer = ieee_flags( 'get', 'exception', '', out )
      PRINT *, out, ' flag raised'

また、overflow 例外フラグが立てられているかどうかを判別するときは、引数の in に overflow を設定します。復帰時に、out が overflow になっていれば、overflow 例外のフラグが立てられているということです。


      ieeer = ieee_flags( 'get', 'exception', 'overflow', out )
      IF ( out.eq. 'overflow') PRINT *, 'overflow flag raised'

例: invalid 例外をクリアーします。


      ieeer = ieee_flags( 'clear', 'exception', 'invalid', out )

例: すべての例外をクリアーします。


      ieeer = ieee_flags( 'clear', 'exception', 'all', out )

例: ゼロに近づけるように丸めます。


      ieeer = ieee_flags( 'set', 'direction', 'tozero', out )

例: 丸めの精度を double に設定します。


      ieeer = ieee_flags( 'set', 'precision', 'double', out )

6.3.1.1 ieee_flags を使用して警告メッセージを抑制する

次の例のように、動作 clear で ieee_flags を呼び出すと、クリアーされてない例外すべてがリセットされます。 プログラムが終了する前にこの呼び出しを置くと、プログラム終了時の浮動小数点例外に関するシステム警告メッセージを抑制します。

例: ieee_flags() を使用して、発生していたすべての例外をクリアーします。


      i = ieee_flags('clear', 'exception', 'all', out )

6.3.1.2 ieee_flags を使用して例外を検出する

次の例は、前の計算によってどの浮動小数点例外フラグが立ったかを決定する方法を示しています。システム include ファイル floatingpoint.h に定義されたビットマスクは、ieee_flags により返された値に適用されます。

この例 DetExcFlg.F では、インクルードファイルは #include 前処理部指令を使用して導入されます。 この指令には、.F 接尾辞のソースファイルを指定しなければいけません。アンダーフローの原因は、最小の倍精度数を 2 で除算しているためです。

例: ieee_flags を使用して例外を検出し、復号化します。


#include "floatingpoint.h"
         CHARACTER*16 out
         DOUBLE PRECISION d_max_subnormal, x
         INTEGER div, flgs, inv, inx, over, under

       x = d_max_subnormal() / 2.0                ! アンダーフローを発生

       flgs=ieee_flags('get','exception','',out)  ! どの例外が発生したのか ?

       inx   = and(rshift(flgs, fp_inexact)  , 1)  ! ieee_flags
       div   = and(rshift(flgs, fp_division) , 1)   ! によって
       under = and(rshift(flgs, fp_underflow), 1)    ! 返される
       over  = and(rshift(flgs, fp_overflow) , 1)     ! 値を
       inv   = and(rshift(flgs, fp_invalid)  , 1)      ! 復号化

       PRINT *, "もっとも優先順位が高い例外:", out
       PRINT *, ' invalid  divide  overflo underflo inexact'
       PRINT '(5i8)', inv, div, over, under, inx
       PRINT *, '(1 = 例外フラグが立っている; 0 = 立っていない)'
       i = ieee_flags('clear', 'exception', 'all', out)    ! すべてクリアー
       END

例: 前出の例 (DetExcFlg.F) をコンパイルして実行します。


demo% f95 DetExcFlg.F
demo% a.out
 もっとも優先順位が高い例外: underflow
  invalid  divide  overflo underflo inexact
       0       0       0       1       1
 (1 = 例外フラグが立っている; 0 = 立っていない)
demo%

6.3.2 IEEE 極値関数

コンパイラは、特別な IEEE 極値を返すために呼び出すことができる関数のセットを提供します。無限大や最小の正規化数などの値は、アプリケーションプログラムの中で直接使用できます。

例: 収束のテストはハードウェアによってサポートされる最小数に基づき、次のようになります。


      IF ( delta .LE. r_min_normal() ) RETURN

利用できる値を次の表にリストします。

表 6–3 IEEE の値を使用する関数

IEEE の値 

倍精度 

単精度 

infinity (無限大)

d_infinity()

r_infinity()

quiet NaN (シグナルを発しない非数)

d_quiet_nan()

r_quiet_nan()

signaling NaN (シグナルを発する非数)

d_signaling_nan()

r_signaling_nan()

min normal (最小の正規化数)

d_min_normal()

r_min_normal()

min subnormal (最小の非正規化数)

d_min_subnormal()

r_min_subnormal()

max subnormal (最大の非正規化数)

d_max_subnormal()

r_max_subnormal()

max normal (最大の正規化数)

d_max_normal()

r_max_normal()

2 つの NaN 値 (quiet と signaling) は順序がないものなので、IF(X.ne.r_quiet_nan())THEN... のように比較の中で使用してはいけません。値が NaN であるかどうかを判別するときは、ir_isnan(r) か id_isnan(d) 関数を使用します。

これらの関数の Fortran 名は、次のマニュアルページにリストされています。

また、次も参照してください。

6.3.3 例外ハンドラと ieee_handler()

通常、IEEE 例外について、次のことを知っておく必要があります。

ユーザールーチンへの例外トラップは、システムの浮動小数点例外に対するシグナルの生成から始まります。「浮動小数点例外のシグナル」を表す UNIX の公式名は SIGFPE です。 SPARC プラットフォームは、デフォルトでは、例外が発生しても SIGFPE を生成しません。システムが SIGFPE を生成するためには、まず、例外トラップを有効にしなければいけません。これは通常は、ieee_handler() への呼び出しによって行います。

6.3.3.1 例外ハンドラ関数を設定する

例外ハンドラとして関数を設定するときは、監視する例外や対応動作といっしょに、関数の名前を ieee_handler() に渡します。 ハンドラを一度設定すると、特定の浮動小数点の例外が発生するたびに、SIGFPE シグナルが生成され、指定した関数が呼び出されます。

ieee_handler() を起動する形式を次の表に示します。

表 6–4 ieee_handler( action , exception , handler) の引数

引数 

種類 

可能な値 

action

文字列

get、set、clear のいずれか

exception

文字列

invalid、division、overflow、underflow、inexact のいずれか

handler

関数名 

ユーザーハンドラ関数の名前、または、SIGFPE_DEFAULT、SIGFPE_IGNORE、SIGFPE_ABORT のいずれか

戻り値 

整数

0 =OK

f95 でコンパイルする Fortran 95 ルーチンで ieee_handler() を呼び出す場合は、次の宣言も必要です。

#include 'floatingpoint.h'

特別な引数 SIGFPE_DEFAULT、SIGFPE_IGNORE、および SIGFPE_ABORT は、これらのインクルードファイルで定義され、特定の例外に対するプログラムの動作を変更するのに使用できます。

SIGFPE_DEFAULT または SIGFPE_IGNORE

指定した例外が発生しても何も動作しない。 

SIGFPE_ABORT

例外発生時にはプログラムが異常終了し、おそらくダンプファイルを生成する。 

6.3.3.2 ユーザー例外ハンドラ関数を作成する

ユーザーの例外ハンドラが行う動作はユーザーが自由に設定できます。しかし、ルーチンは、次に示す 3 つの引数を取る、整数型関数でなければいけません。

handler_name( sig, sip, uap )

例: 例外ハンドラ関数


      INTEGER FUNCTION hand( sig, sip, uap )
      INTEGER sig, location
      STRUCTURE /fault/
           INTEGER address
           INTEGER trapno
      END STRUCTURE
      STRUCTURE /siginfo/
           INTEGER si_signo
           INTEGER si_code
           INTEGER si_errno
           RECORD /fault/ fault
      END STRUCTURE
      RECORD /siginfo/ sip
      location = sip.fault.address
      ... ユーザーの行う処理 ...
      END

64 ビット SPARC アーキテクチャーで実行するには、STRUCTURE 内のすべての INTEGER 宣言を INTEGER*8 で置き換えて、この例を変更する必要があります。

ieee_handler() によって有効にされるハンドラルーチンが例で示すように Fortran で書かれている場合、ハンドラルーチンは 1 番目の引数 (sig) に対しては、一切の参照を行なってはいけません。1 番目の引数は値でルーチンに渡され、loc(sig) としてのみ参照できます。この値はシグナル番号です。この値はシグナル番号です。

ハンドラを使用して例外を検出する

次の例では、浮動小数点の例外を検出するためのハンドラルーチンを作成する方法を示します。

例: 例外を検出して、異常終了します。 


demo% cat DetExcHan.f
      EXTERNAL myhandler
      REAL :: r = 14.2 , s = 0.0
      i = ieee_handler ('set', 'division', myhandler )
      t = r/s
      END

      INTEGER FUNCTION myhandler(sig,code,context)
      INTEGER sig, code, context(5)
      CALL abort()
      END
demo% f95 DetExcHan.f
demo% a.out
異常終了
demo%

SIGFPE は、浮動小数点演算の例外が発生すると必ず生成されます。SIGFPE が検出されると、制御が myhandler 関数に渡され、即座に異常終了します。-g を付けてコンパイルし、dbx を使用すると、例外箇所を突き止めることができます。

ハンドラを使用して例外箇所を突き止める

例: 例外箇所を突き止め (アドレスを出力して)、異常終了します。


demo% cat LocExcHan.F
#include "floatingpoint.h"
      EXTERNAL Exhandler
      INTEGER Exhandler, i, ieee_handler
      REAL:: r = 14.2 ,  s = 0.0 , t
C  ゼロ除算の検出
      i = ieee_handler( 'set', 'division', Exhandler )
      t = r/s
      END

      INTEGER FUNCTION Exhandler( sig, sip, uap)
      INTEGER sig
      STRUCTURE /fault/
                  INTEGER address
      END STRUCTURE
      STRUCTURE /siginfo/
                  INTEGER si_signo
                  INTEGER si_code
                  INTEGER si_errno
                  RECORD /fault/ fault
      END STRUCTURE
      RECORD /siginfo/ sip
      WRITE (*,10)  sip.si_signo, sip.si_code, sip.fault.address
10        FORMAT(' シグナル ',i4,' コード ',i4,' 16 進アドレス ', Z8 )
      Exhandler=1
      CALL abort()
      END
demo%  f95 -g LocExcHan.F
demo%  a.out
シグナル 8 コード 3 16 進アドレス 11230
異常終了
demo%

64 ビット SPARC 環境では、各 STRUCTURE 内の INTEGER 宣言を INTEGER*8 で置き換えて、書式中の i4 を i8 で置き換えます。この例では、f95 コンパイラへの拡張により、VAX Fortran の STRUCTURE 文を受け付けられるようにしています。

ほとんどの場合、例外の実際のアドレスを知るということは、dbx だけに意味があります。


demo% dbx a.out
(dbx) stopi at 0x11230     ブレークポイントをアドレスに設定する
(2) stopi at &MAIN+0x68
(dbx) run               プログラムを実行する
実行中: a.out
(プロセス id 18803)
MAIN で 0x11230 で停止しました
MAIN+0x68:      fdivs   %f3, %f2, %f2
(dbx) where           例外が発生した行番号を表示する
=>[1] MAIN()、"LocExcHan.F" の 7 行目
(dbx) list 7          ソースコード行を表示する
    7         t = r/s
(dbx) cont            ブレークポイント後、継続してハンドラルーチンに入る
シグナル 8 コード 3 16 進アドレス 11230
異常終了:
_kill 0xef6e18a4 でシグナル ABRT (異常終了) が呼び出されました
_kill+0x8:      bgeu    _kill+0x30
現関数: exhandler
   24         CALL abort()
(dbx) quit
demo%

もちろん、エラーの原因となるソース行を決定するためのより簡単な方法があります。しかし、この例は、例外処理の基本を示すのが目的です。

6.4 IEEE の例外のデバッグ

どこで 例外が発生したかを突き止めるためには、トラップを有効にした例外が必要です。そのためには、 -ftrap=common オプション (f95 でコンパイルする場合のデフォルト) を付けてコンパイルするか、ieee_handler() により例外ハンドラルーチンを呼び出します。例外トラップを有効にして、dbx からプログラムを実行し、dbx catch FPE コマンドを使用すれば、どこでエラーが発生するかを見つけることができます。

-ftrap=common を付けてコンパイルすることの利点は、例外をトラップするようにソースコードを変更する必要がないことです。しかし、ieee_handler() を呼び出すことによって、探すべき例外をより限定できます。

例: コンパイルし直して、dbx を使用します。


demo% f95 -g myprogram.f
demo% dbx a.out
a.out の読み込み中
...
(dbx) catch FPE
(dbx) run
実行中: a.out
(プロセス id 19739)
シグナル FPE (ゼロによる浮動小数点除算) 関数 MAIN 行番号 212 ファイル "myprogram.f"
  212               Z = X/Y
(dbx)  print Y
y = 0.0
(dbx)

プログラムが終了したときオーバーフローとほかの例外が発生していた場合、ieee_handler() を呼び出してオーバーフローだけをトラップすることによって、最初のオーバーフローを突き止めることができます。このためには、次の例に示すように、主プログラムのソースコードを必要最小限だけ修正する必要があります。

例: ほかの例外も発生しているときにオーバーフローを突き止めます。


demo% 


cat myprog.F
#include "floatingpoint.h"
        program myprogram
...
      ier = ieee_handler('set','overflow',SIGFPE_ABORT)
...
demo% f95 -g myprog.F
demo% dbx a.out
a.out の読み込み中
...
(dbx) catch FPE
(dbx) run
実行中: a.out
(プロセス id 19793)
シグナル FPE (浮動小数点オーバーフロー) 関数 MAIN 行番号 55 ファイル "myprog.F"
   55               w = rmax * 200.                     ! オーバーフローの原因
(dbx) cont                                   ! 実行を継続して終了する
実行完了。終了コードは、0 です
(dbx)

例外を選択するため、この例では、#include を導入しています。 このため、ソースファイルの名前を .F 接尾辞に変更し ieee_handler() を呼び出す必要があります。さらに一歩進んで、オーバーフロー例外時に呼び出されるユーザー独自のハンドラ関数を作成し、アプリケーション特有な解析を行い、異常終了する前に中間結果またはデバッグ結果を出力することもできます。

6.5 数値に関連したその他の問題

この節では、無効演算、ゼロ除算、オーバーフロー、アンダーフロー、結果不正確の例外を予想外に生成する算術演算に関連する、より実際的な問題について説明します。

たとえば、IEEE 規格より前には、2 つの極めて小さな数をコンピュータで乗算すると、結果はゼロになっていました。メインフレームやミニコンピュータのほとんどが、この方式でした。これに対し、IEEE の算術演算では、段階的なアンダーフローが計算の範囲を動的に拡張します。

たとえば、マシンのイプシロン、つまり表現可能な最小値が 1.0E-38 である 32 ビットプロセッサを想定し、2 つの小さな数値を乗算してみます。


      a = 1.0E-30
      b = 1.0E-15
      x = a * b

旧式の算術演算では、0.0 になりますが、IEEE の算術演算では (同じワード長で) 結果は 1.40130E-45 となります。アンダーフローは、マシンが本来表せる値よりも小さい結果になったということを示します。この結果は、いくつかのビットが仮数部から「盗まれ」、指数部へシフトされたことによってできました。 結果 (非正規化数)は、ある場合には精度が低くなりますが、逆に高くなる場合もあります。これに関する詳細な説明はこのマニュアルの範囲外です。詳細に興味のある方は、1980 年 1 月に発行された『Computer』誌 (第 13 巻 No.1)、特に J. Coonen 氏の論文「Underflow and the Denormalized Numbers」を参照してください。

科学技術プログラムのほとんどは、方程式を解いたり、行列の因数分解など、丸めに敏感に反応するコードセクションがあります。段階的なアンダーフローがなければ、プログラマは不正確なしきい値への接近を検出する独自の方法を実装します。実装できなければ、独自のアルゴリズムの安定した実装はあきらめるかしかありません。

これらの話題に関する詳細は、『数値計算ガイド』を参照してください。

6.5.1 単純なアンダーフローを防ぐ

アプリケーションの中には、実際に、非常にゼロに近いところで多くの処理をしているものがあります。これは、誤差や差分補正のアルゴリズムでよく発生します。数値的に安全で最大のパフォーマンスを得るには、重要な演算は拡張精度で行うべきです。アプリケーションが単精度の場合は、重要な演算を倍精度にすればかまいません。

例: 単精度の、簡単なドット積の演算


      sum = 0
      DO i = 1, n
         sum = sum + a(i) * b(i)
      END DO

a(i) と b(i) が極めて小さい数であれば、多数のアンダーフローが発生することになります。演算を倍精度に移行すると、ドット積をより正確に計算でき、アンダーフローもなくなります。


      DOUBLE PRECISION sum
      DO i = 1, n
         sum = sum + dble(a(i)) * dble(b(i))
      END DO
      result = sum

アンダーフローに関して、SPARC プロセッサを旧式のシステムのように動作 (Store Zero) させることができます。このためには、ライブラリルーチン nonstandard_arithmetic() への呼び出しを追加するか、アプリケーションの主プログラムを -fns オプションを付けてコンパイルします。

6.5.2 間違った答えのまま継続する

結果が明らかに間違っている場合でも処理が継続されることを、疑問に思われるかもしれません。IEEE の算術演算では、NaN や Inf など、どの種別の間違った答えを無視できるかをユーザーが指定できます。そして、そのような区別に基づいて決定が行われます。

たとえば、回路シミュレーションを想定してみましょう。特有の 50 行の演算の中で、引数の目的となる重要な変数は電圧量だけであり、この変数のとり得る値は +5v、0、-5v だけであると仮定します。

途中の演算結果が正当な範囲にくるように、演算の各部分を詳細に調整できます。

さらに、Inf は許可されない値であるので、大きな数を乗算しないような特別なロジックが必要です。

IEEE 算術演算ではロジックはかなり単純にできます。演算を明確な形式で記述し、最終結果を正しい値に導くだけでかまいません。 なぜなら、Inf の発生する可能性があり、それは簡単にテストできるからです。

さらに、0/0 の特殊なケースも検出でき、ユーザーの望む形で処理できます。不必要な比較を行わなくてもすむため、結果は読みやすく、実行も高速です。

6.5.3 アンダーフローの頻発

2 つの極めて小さな数を乗算すると、結果はアンダーフローとなります。

あらかじめ、乗算 (または減算) の被演算子が小さくなり、アンダーフローしそうであることがわかっていれば、計算を倍精度で行い、そのあとで結果を単精度に変換します。

たとえば、ドット積ループは次のようになります。


  real sum, a(maxn), b(maxn)
  ...
  do i =1, n
      sum = sum + a(i)*b(i)
  enddo

a(*) と b(*) は、小さな要素が入ることがわかっているので、倍精度で実行し、数値の正確性を保持します。


  real a(maxn), b(maxn)
  double sum
  ...
  do i =1, n
      sum = sum + a(i)*dble(b(i))
  enddo

このようにすることによって、オリジナルのループが原因であるアンダーフロー頻発をソフトウェア的に解決し、パフォーマンスを上げることができます。しかし、これに関しては、確実で手間のかからない方法はありません。計算が多いコードを使用するユーザー自身の経験のほうが、より適切な解決策を決定することでしょう。

6.6 区間演算

注: 現在、区間演算は SPARC プラットフォームのみ利用できます。

Fortran 95 コンパイラの f95 では、組み込みデータ型として intervals をサポートしています。区間は、[a, b] ={z | a≤ z≤ b} で表され、1 組の数字は a ≤ b となります。

区間を組み込みデータ型として Fortran 95 に導入したので、開発者は Fortran 95 の適用可能な構文や記号をすべてすぐに使用できるようになります。INTERVAL データ型のほかに、f95 により次の区間拡張機能が Fortran 95 に取り込まれます。

f95 コマンド行オプション -xinterval により、コンパイラの区間演算機能が使用できるようになります。『Fortran ユーザーズガイド』を参照してください。

Fortran 95 の区間演算の詳細については、『Fortran 95 Interval Arithmetic Programming Reference』を参照してください。