通常、累積例外ビットの検査またはクリアはユーザープログラムで行います。次の例は、累積例外フラグを 検査する C プログラムです。
使用例 15 累積例外フラグの検査#include <sunmath.h>
#include <sys/ieeefp.h>
int main()
{
int code, inexact, division, underflow, overflow, invalid;
double x;
char *out;
/* cause an underflow exception */
x = max_subnormal() / 2.0;
/* this statement insures that the previous */
/* statement is not optimized away */
printf("x = %g\n",x);
/* find out which exceptions are raised */
code = ieee_flags("get", "exception", "", &out);
/* decode the return value */
inexact = (code >> fp_inexact) & 0x1;
underflow = (code >> fp_underflow) & 0x1;
division = (code >> fp_division) & 0x1;
overflow = (code >> fp_overflow) & 0x1;
invalid = (code >> fp_invalid) & 0x1;
/* "out" is the raised exception with the highest priority */
printf(" Highest priority exception is: %s\n", out);
/* The value 1 means the exception is raised, */
/* 0 means it isn't. */
printf("%d %d %d %d %d\n", invalid, overflow, division,
underflow, inexact);
ieee_retrospective_();
return 0;
}
上記のプログラムを実行した場合の出力は次のようになります。
demo% a.out
x = 1.11254e-308
Highest priority exception is: underflow
0 0 0 1 1
Note:IEEE floating-point exception flags raised:
Inexact; Underflow;
See the Numerical Computation Guide, ieee_flags(3M)
使用例 16 累積例外フラグの検査 – Fortran/*
A Fortran example that:
* causes an underflow exception
* uses ieee_flags to determine which exceptions are raised
* decodes the integer value returned by ieee_flags
* clears all outstanding exceptions
Remember to save this program in a file with the suffix .F, so that
the c preprocessor is invoked to bring in the header file
floatingpoint.h.
*/
#include <floatingpoint.h>
program decode_accrued_exceptions
double precision x
integer accrued, inx, div, under, over, inv
character*16 out
double precision d_max_subnormal
c Cause an underflow exception
x = d_max_subnormal() / 2.0
c Find out which exceptions are raised
accrued = ieee_flags('get', 'exception', '', out)
c Decode value returned by ieee_flags using bit-shift intrinsics
inx = and(rshift(accrued, fp_inexact) , 1)
under = and(rshift(accrued, fp_underflow), 1)
div = and(rshift(accrued, fp_division) , 1)
over = and(rshift(accrued, fp_overflow) , 1)
inv = and(rshift(accrued, fp_invalid) , 1)
c The exception with the highest priority is returned in "out"
print *, "Highest priority exception is ", out
c The value 1 means the exception is raised; 0 means it is not
print *, inv, over, div, under, inx
c Clear all outstanding exceptions
i = ieee_flags('clear', 'exception', 'all', out)
end
出力は次のようになります。
Highest priority exception is underflow 0 0 0 1 1
通常、ユーザープログラムでは例外フラグを設定 しませんが、設定は可能です。次の C の例でこれを 示します。
#include <sunmath.h>
int main()
{
int code;
char *out;
if (ieee_flags("clear", "exception", "all", &out) != 0)
printf("could not clear exceptions\n");
if (ieee_flags("set", "exception", "division", &out) != 0)
printf("could not set exception\n");
code = ieee_flags("get", "exception", "", &out);
printf("out is: %s , fp exception code is: %X \n",
out, code);
return 0;
}
SPARC では、上記のプログラムの出力は次のようになります。
out is: division , fp exception code is: 2
x86 では、出力は次のようになります。
out is: division , fp exception code is: 4
次の例は、例外を特定するためのシグナルハンドラをインストールする Fortran プログラムです (SPARC ベースシステムの場合のみ)。
使用例 17 アンダーフローでのトラップ (SPARC) program demo
c declare signal handler function
external fp_exc_hdl
double precision d_min_normal
double precision x
c set up signal handler
i = ieee_handler('set', 'common', fp_exc_hdl)
if (i.ne.0) print *, 'ieee trapping not supported here'
c cause an underflow exception (it will not be trapped)
x = d_min_normal() / 13.0
print *, 'd_min_normal() / 13.0 = ', x
c cause an overflow exception
c the value printed out is unrelated to the result
x = 1.0d300
x = x * x
print *, '1.0d300*1.0d300 = ', x
end
c
c the floating-point exception handling function
c
integer function fp_exc_hdl(sig, sip, uap)
integer sig, code, addr
character label*16
c
c The structure /siginfo/ is a translation of siginfo_t
c from <sys/siginfo.h>
c
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
c See <sys/machsig.h> for list of FPE codes
c Figure out the name of the SIGFPE
code = sip.si_code
if (code.eq.3) label = 'division'
if (code.eq.4) label = 'overflow'
if (code.eq.5) label = 'underflow'
if (code.eq.6) label = 'inexact'
if (code.eq.7) label = 'invalid'
addr = sip.fault.address
c Print information about the signal that happened
write (*,77) code, label, addr
77 format ('floating-point exception code ', i2, ',',
* a17, ',', ' at address ', z8 )
end
前のコードを –f77 を指定してコンパイルすると、出力は次のようになります。
d_min_normal() / 13.0 = 1.7115952757748-309
floating-point exception code 4, overflow , at address 1131C
1.0d300*1.0d300 = 1.0000000000000+300
Note: IEEE floating-point exception flags raised:
Inexact; Underflow;
IEEE floating-point exception traps enabled:
overflow; division by zero; invalid operation;
See the Numerical Computation Guide, ieee_flags(3M),
ieee_handler(3M)
次の例は、SPARC ベースシステムでの より複雑な C の例です。
使用例 18 無効、ゼロ除算、オーバーフロー、アンダーフロー、および不正確でのトラップ (SPARC)/*
* Generate the 5 IEEE exceptions: invalid, division,
* overflow, underflow and inexact.
*
* Trap on any floating point exception, print a message,
* and continue.*
* Note that you could also inquire about raised exceptions by
* i = ieee("get","exception","",&out);* where out contains the name of the highest exception
* raised, and i can be decoded to find out about all the
* exceptions raised.
*/
#include <sunmath.h>
#include <signal.h>
#include <siginfo.h>
#include <ucontext.h>
extern void trap_all_fp_exc(int sig, siginfo_t *sip,
ucontext_t *uap);
int main()
{
double x, y, z;
char *out;
/*
* Use ieee_handler to establish "trap_all_fp_exc"
* as the signal handler to use whenever any floating
* point exception occurs.
*/
if (ieee_handler("set", "all", trap_all_fp_exc) != 0)
printf(" IEEE trapping not supported here.\n");
/* disable trapping (uninteresting) inexact exceptions */
if (ieee_handler("set", "inexact", SIGFPE_IGNORE) != 0)
printf("Trap handler for inexact not cleared.\n");
/* raise invalid */
if (ieee_flags("clear", "exception", "all", &out) != 0)
printf(" could not clear exceptions\n");
printf("1. Invalid: signaling_nan(0) * 2.5\n");
x = signaling_nan(0);
y = 2.5;
z = x * y;
/* raise division */
if (ieee_flags("clear", "exception", "all", &out) != 0)
printf(" could not clear exceptions\n");
printf("2. Div0: 1.0 / 0.0\n");
x = 1.0;
y = 0.0;
z = x / y;
/* raise overflow */
if (ieee_flags("clear", "exception", "all", &out) != 0)
printf(" could not clear exceptions\n");
printf("3. Overflow: -max_normal() - 1.0e294\n");
x = -max_normal();
y = -1.0e294;
z = x + y;
/* raise underflow */
if (ieee_flags("clear", "exception", "all", &out) != 0)
printf(" could not clear exceptions\n");
printf("4. Underflow: min_normal() * min_normal()\n");
x = min_normal();
y = x;
z = x * y;
/* enable trapping on inexact exception */
if (ieee_handler("set", "inexact", trap_all_fp_exc) != 0)
printf("Could not set trap handler for inexact.\n");
/* raise inexact */
if (ieee_flags("clear", "exception", "all", &out) != 0)
printf(" could not clear exceptions\n");
printf("5. Inexact: 2.0 / 3.0\n");
x = 2.0;
y = 3.0;
z = x / y;
/* don't trap on inexact */
if (ieee_handler("set", "inexact", SIGFPE_IGNORE) != 0)
printf(" could not reset inexact trap\n");
/* check that we're not trapping on inexact anymore */
if (ieee_flags("clear", "exception", "all", &out) != 0)
printf(" could not clear exceptions\n");
printf("6. Inexact trapping disabled; 2.0 / 3.0\n");
x = 2.0;
y = 3.0;
z = x / y;
/* find out if there are any outstanding exceptions */
ieee_retrospective_();
/* exit gracefully */
return 0;
}
void trap_all_fp_exc(int sig, siginfo_t *sip, ucontext_t *uap) {
char *label = "undefined";
/* see /usr/include/sys/machsig.h for SIGFPE codes */
switch (sip->si_code) {
case FPE_FLTRES:
label = "inexact";
break;
case FPE_FLTDIV:
label = "division";
break;
case FPE_FLTUND:
label = "underflow";
break;
case FPE_FLTINV:
label = "invalid";
break;
case FPE_FLTOVF:
label = "overflow";
break;
}
printf(
" signal %d, sigfpe code %d: %s exception at address %x\n",
sig, sip->si_code, label, sip->__data.__fault.__addr);
}
出力は、次のようになります。
1. Invalid: signaling_nan(0) * 2.5 signal 8, sigfpe code 7: invalid exception at address 10da8 2. Div0: 1.0 / 0.0 signal 8, sigfpe code 3: division exception at address 10e44 3. Overflow: -max_normal() - 1.0e294 signal 8, sigfpe code 4: overflow exception at address 10ee8 4. Underflow: min_normal() * min_normal() signal 8, sigfpe code 5: underflow exception at address 10f80 5. Inexact: 2.0 / 3.0 signal 8, sigfpe code 6: inexact exception at address 1106c 6. Inexact trapping disabled; 2.0 / 3.0 Note: IEEE floating-point exception traps enabled: underflow; overflow; division by zero; invalid operation; See the Numerical Computation Guide, ieee_handler(3M)
次のコードは、SPARC で ieee_handler およびインクルードファイルを使用して、特定の例外状況のデフォルトの結果を変更する方法を示しています。
使用例 19 例外状況のデフォルトの結果の変更/*
* Cause a division by zero exception and use the
* signal handler to substitute MAXDOUBLE (or MAXFLOAT)
* as the result.
*
* compile with the flag -Xa
*/
#include <values.h>
#include <siginfo.h>
#include <ucontext.h>
void division_handler(int sig, siginfo_t *sip, ucontext_t *uap);
int main() {
double x, y, z;
float r, s, t;
char *out;
/*
* Use ieee_handler to establish division_handler as the
* signal handler to use for the IEEE exception division.
*/
if (ieee_handler("set","division",division_handler)!=0) {
printf(" IEEE trapping not supported here.\n");
}
/* Cause a division-by-zero exception */
x = 1.0;
y = 0.0;
z = x / y;
/*
* Check to see that the user-supplied value, MAXDOUBLE,
* is indeed substituted in place of the IEEE default
* value, infinity.
*/
printf("double precision division: %g/%g = %g \n",x,y,z);
/* Cause a division-by-zero exception */
r = 1.0;
s = 0.0;
t = r / s;
/*
* Check to see that the user-supplied value, MAXFLOAT,
* is indeed substituted in place of the IEEE default
* value, infinity.
*/
printf("single precision division: %g/%g = %g \n",r,s,t);
ieee_retrospective_();
return 0;
}
void division_handler(int sig, siginfo_t *sip, ucontext_t *uap) {
int inst;
unsigned rd, mask, single_prec=0;
float f_val = MAXFLOAT;
double d_val = MAXDOUBLE;
long *f_val_p = (long *) &f_val;
/* Get instruction that caused exception. */
inst = uap->uc_mcontext.fpregs.fpu_q->FQu.fpq.fpq_instr;
/*
* Decode the destination register. Bits 29:25 encode the
* destination register for any SPARC floating point
* instruction.
*/
mask = 0x1f;
rd = (mask & (inst >> 25));
/*
* Is this a single precision or double precision
* instruction? Bits 5:6 encode the precision of the
* opcode; if bit 5 is 1, it's sp, else, dp.
*/
mask = 0x1;
single_prec = (mask & (inst >> 5));
/* put user-defined value into destination register */
if (single_prec) {
uap->uc_mcontext.fpregs.fpu_fr.fpu_regs[rd] =
f_val_p[0];
} else {
uap->uc_mcontext.fpregs.fpu_fr.fpu_dregs[rd/2] = d_val;
}
}
次の出力が予測されます。
double precision division: 1/0 = 1.79769e+308 single precision division: 1/0 = 3.40282e+38 Note: IEEE floating-point exception traps enabled: division by zero; See the Numerical Computation Guide, ieee_handler(3M)
特定の浮動小数点例外の場合は、 ieee_handler を使用して強制的にプログラムを中止させることができます。
#include <floatingpoint.h>
program abort
c
ieeer = ieee_handler('set', 'division', SIGFPE_ABORT)
if (ieeer .ne. 0) print *, ' ieee trapping not supported'
r = 14.2
s = 0.0
r = r/s
c
print *, 'you should not see this; system should abort'
c
end
次の例は、libm によって提供される例外処理機能のいくつかを使用する方法を示しています。最初の例は、数値 x と係数 a0、a1、...、aN、および b0、b1、...、bN-1 の場合に、関数 f(x) とその一次導関数 f'(x) を評価するタスクに基づいています。ここで f() は次の連分数です。
f(x) =a0 + b0 / (x + a1 + b1 / (x + ... / (x + aN-1 + bN-1 / (x + aN))...))
IEEE 演算では、f() の計算は単純です。中間除算の 1 つがオーバーフローしたり、ゼロ除算を行う場合でも、標準規格によって指定されているデフォルトの値 (正しい符号が付けられた無限大) から正しい結果が算出されます。一方、f'() の計算は、これを評価するもっとも簡潔な形式が、削除可能な特異点を持つことができるため、f の計算よりも難しくなります。計算中にこれらの特異点の 1 つが出現すると、不定形 0/0、0*無限大、または無限大/無限大のいずれか 1 つの評価が試みられます (これらはすべて無効な演算例外を発生します)。W. Kahan は、前置換という機能によってこれらの例外を処理する方法を提唱しています。
前置換は、例外に対する IEEE のデフォルトの応答を拡張したものであり、例外演算の結果を置換する値をユーザーが前もって指定できます。libm の例外処理機能を使用すると、FEX_CUSTOM 例外処理モードでハンドラをインストールすることによって、プログラムで簡単に前置換を実装できます。このモードでは、ハンドラに渡された info パラメータが差しているデータ構造体に値を格納するだけで、ハンドラは例外演算の結果に任意の値を設定できます。次の例は、FEX_CUSTOM ハンドラによって実装した前置換を使用して連分数とその導関数を計算するプログラムの例を示しています。
使用例 20 FEX_CUSTOM ハンドラを使用した連分数とその導関数の計算#include <stdio.h>
#include <sunmath.h>
#include <fenv.h>
volatile double p;
void handler(int ex, fex_info_t *info)
{
info->res.type = fex_double;
if (ex == FEX_INV_ZMI)
info->res.val.d = p;
else
info->res.val.d = infinity();
}
/*
* Evaluate the continued fraction given by coefficients a[j] and
* b[j] at the point x; return the function value in *pf and the
* derivative in *pf1
*/
void continued_fraction(int N, double *a, double *b,
double x, double *pf, double *pf1)
{
fex_handler_t oldhdl; /* for saving/restoring handlers */
volatile double t;
double f, f1, d, d1, q;
int j;
fex_getexcepthandler(&oldhdl, FEX_DIVBYZERO | FEX_INVALID);
fex_set_handling(FEX_DIVBYZERO, FEX_NONSTOP, NULL);
fex_set_handling(FEX_INV_ZDZ | FEX_INV_IDI | FEX_INV_ZMI,
FEX_CUSTOM, handler);
f1 = 0.0;
f = a[N];
for (j = N - 1; j >= 0; j--) {
d = x + f;
d1 = 1.0 + f1;
q = b[j] / d;
/* the following assignment to the volatile variable t
is needed to maintain the correct sequencing between
assignments to p and evaluation of f1 */
t = f1 = (-d1 / d) * q;
p = b[j-1] * d1 / b[j];
f = a[j] + q;
}
fex_setexcepthandler(&oldhdl, FEX_DIVBYZERO | FEX_INVALID);
*pf = f;
*pf1 = f1;
}
/* For the following coefficients, x = -3, 1, 4, and 5 will all
encounter intermediate exceptions */
double a[] = { -1.0, 2.0, -3.0, 4.0, -5.0 };
double b[] = { 2.0, 4.0, 6.0, 8.0 };
int main()
{
double x, f, f1;
int i;
feraiseexcept(FE_INEXACT); /* prevent logging of inexact */
fex_set_log(stdout);
fex_set_handling(FEX_COMMON, FEX_ABORT, NULL);
for (i = -5; i <= 5; i++) {
x = i;
continued_fraction(4, a, b, x, &f, &f1);
printf("f(% g) = %12g, f'(% g) = %12g\n", x, f, x, f1);
}
return 0;
}
このプログラムについては、順を追って解説します。入口で、関数 continued_fraction は、ゼロ除算およびすべての無効な演算例外に対する現在の例外処理モードを保存します。続いて、ゼロ除算に対して無停止の例外処理、および 3 つの不定形に対して FEX_CUSTOM ハンドラを設定します。このハンドラは 0/0 と無限大/無限大の両方を無限大で置換しますが、0*無限大は大域変数 p の値で置換します。後続の 0*無限大の無効な演算の置換には、正しい値を設定できるように p は関数を評価するループで毎回再計算する必要があります。また、p はループ内で明示的に記述されていないため、コンパイラによって削除されないように volatile と宣言する必要があります。最後に、コンパイラが p に対する代入を、例外を発生させる可能性のある計算 (これに対して p が前置換の値を提供します) の上または下に移動させないように、その計算の結果も volatile 変数 (このプログラムでは t) に代入しています。fex_setexcepthandler の最後の呼び出しによって、ゼロ除算および無効な演算に対する元の処理モードが復元されます。
メインプログラムは、fex_set_log 関数を呼び出して、遡及診断のロギングを有効にしています。この呼び出しの前に、メインプログラムは不正確フラグを発生させ、これによって不正確演算がロギングされなくなります。遡及診断セクションで説明しているように、FEX_NONSTOP モードでは例外フラグが発生しても例外がログに記録されません。また、メインプログラムは、一般的な例外に対して FEX_ABORT モードを設定し、continued_fraction によって明示的に処理されない特異な例外によってプログラムが終了しないようにしています。最後に、プログラムは複数のポイントで特定の連分数を評価しています。次の出力例が示しているように、計算で実際に中間例外が発生しています。
f(-5) = -1.59649, f'(-5) = -0.1818 f(-4) = -1.87302, f'(-4) = -0.428193 Floating point division by zero at 0x08048dbe continued_fraction, nonstop mode 0x08048dc1 continued_fraction 0x08048eda main Floating point invalid operation (inf/inf) at 0x08048dcf continued_fraction, handler: handler 0x08048dd2 continued_fraction 0x08048eda main Floating point invalid operation (0*inf) at 0x08048dd2 continued_fraction, handler: handler 0x08048dd8 continued_fraction 0x08048eda main f(-3) = -3, f'(-3) = -3.16667 f(-2) = -4.44089e-16, f'(-2) = -3.41667 f(-1) = -1.22222, f'(-1) = -0.444444 f( 0) = -1.33333, f'( 0) = 0.203704 f( 1) = -1, f'( 1) = 0.333333 f( 2) = -0.777778, f'( 2) = 0.12037 f( 3) = -0.714286, f'( 3) = 0.0272109 f( 4) = -0.666667, f'( 4) = 0.203704 f( 5) = -0.777778, f'( 5) = 0.0185185
x = 1、4、および 5 の場合の f'(x) の計算で発生する例外は、プログラム内で x = -3 のときに起きる例外と同じサイトで発生するため、遡及診断メッセージに出力されません。
上記のプログラムは、連分数とその導関数の評価で発生する例外を処理する場合の、もっとも効率がよい方法ではない可能性があります。その理由の 1 つは、必要であるかどうかにかかわらず、ループが繰り返されるごとに前置換の値を再計算する必要があることです。この場合、前置換の値の計算で浮動小数点除算が行われますが、最新の SPARC および x86 プロセッサでは浮動小数点除算は比較的遅い演算です。そのうえ、ループ自体にすでに 2 つの除算が含まれていますが、ほとんどの SPARC および x86 プロセッサは 2 つの除算演算の実行をオーバーラップできないため、除算がループ内のボトルネックとなりやすくなります。別の除算を追加するとボトルネックはさらに悪化します。
1 つの除算のみを必要とするようにループを記述しなおすことができます。実際、前置換値の計算では除算を必要としません。このようにループを記述しなおすには、b 配列内の係数の隣接要素比を前もって計算する必要があります。これにより、複数の除算の演算でボトルネックは排除されますが、前置換値の計算に関連するすべての算術演算が除外されるわけではありません。さらに、前置換値と演算結果の両方が volatile 変数に前置換されるように割り当てる必要があるため、プログラムの速度を低下させるメモリー演算が増えます。これらの代入は、コンパイラが特定のキー操作を並べ替えないようにするために必要ですが、コンパイラがほかの無関係な演算を並べ替えることも防止することになります。このため、この例のように前置換を使用して例外を処理すると、メモリー演算が増え、通常では可能な最適化が行われなくなります。これらの例外を、さらに効率よく処理することは可能でしょうか。
高速な前置換のための特別なハードウェアがサポートされていない場合、この例の例外をもっとも効率よく処理する方法は、次のバージョンに示すようにフラグを使用することです。
使用例 21 例外処理へのフラグの使用 (続き)#include <stdio.h>
#include <math.h>
#include <fenv.h>
/*
* Evaluate the continued fraction given by coefficients a[j] and
* b[j] at the point x; return the function value in *pf and the
* derivative in *pf1
*/
void continued_fraction(int N, double *a, double *b,
double x, double *pf, double *pf1)
{
fex_handler_t oldhdl;
fexcept_t oldinvflag;
double f, f1, d, d1, pd1, q;
int j;
fex_getexcepthandler(&oldhdl, FEX_DIVBYZERO | FEX_INVALID);
fegetexceptflag(&oldinvflag, FE_INVALID);
fex_set_handling(FEX_DIVBYZERO | FEX_INV_ZDZ | FEX_INV_IDI |
FEX_INV_ZMI, FEX_NONSTOP, NULL);
feclearexcept(FE_INVALID);
f1 = 0.0;
f = a[N];
for (j = N - 1; j >= 0; j--) {
d = x + f;
d1 = 1.0 + f1;
q = b[j] / d;
f1 = (-d1 / d) * q;
f = a[j] + q;
}
if (fetestexcept(FE_INVALID)) {
/* recompute and test for NaN */
f1 = pd1 = 0.0;
f = a[N];
for (j = N - 1; j >= 0; j--) {
d = x + f;
d1 = 1.0 + f1;
q = b[j] / d;
f1 = (-d1 / d) * q;
if (isnan(f1))
f1 = b[j] * pd1 / b[j+1];
pd1 = d1;
f = a[j] + q;
}
}
fesetexceptflag(&oldinvflag, FE_INVALID);
fex_setexcepthandler(&oldhdl, FEX_DIVBYZERO | FEX_INVALID);
*pf = f;
*pf1 = f1;
}
このバージョンでは、最初のループはデフォルトの無停止モードで f(x) および f'(x) の計算を試みています。無効フラグが発生すると、2 番目のループは NaN の形態をテストして f(x) および f'(x) を明示的に再計算します。通常、無効な演算例外は発生しないため、プログラムは最初のループのみを実行します。このループには、volatile 変数に対する参照も余分な算術演算も含まれていないため、コンパイラが許すかぎりの速度で実行されます。この効率を得るためには、例外が発生した場合の処理について、2 番目のループを最初のループとほとんど同じように記述する必要があります。このトレードオフは、浮動小数点例外処理において生じることがある典型的なジレンマです。
libm の例外処理機能は主に C/C++ プログラムからの使用を想定していますが、Sun の Fortran 言語の相互運用性機能を使用すると、Fortran プログラムからも一部の libm 関数を呼び出すことができます。
次の例では、前置換を使用して連分数とその導関数を評価するための Fortran バージョンのプログラムを示します (SPARC のみ)。
使用例 22 前置換を使用した連分数とその導関数の評価 (SPARC)c
c Presubstitution handler
c
subroutine handler(ex, info)
structure /fex_numeric_t/
integer type
union
map
integer i
end map
map
integer*8 l
end map
map
real f
end map
map
real*8 d
end map
map
real*16 q
end map
end union
end structure
structure /fex_info_t/
integer op, flags
record /fex_numeric_t/ op1, op2, res
end structure
integer ex
record /fex_info_t/ info
common /presub/ p
double precision p, d_infinity
volatile p
c 4 = fex_double; see <fenv.h> for this and other constants
info.res.type = 4
c x'80' = FEX_INV_ZMI
if (loc(ex) .eq. x'80') then
info.res.d = p
else
info.res.d = d_infinity()
endif
return
end
c
c Evaluate the continued fraction given by coefficients a(j) and
c b(j) at the point x; return the function value in f and the
c derivative in f1
c
subroutine continued_fraction(n, a, b, x, f, f1)
integer n
double precision a(*), b(*), x, f, f1
common /presub/ p
integer j, oldhdl
dimension oldhdl(24)
double precision d, d1, q, p, t
volatile p, t
data ixff2/x'ff2'/
data ix2/x'2'/
data ixb0/x'b0'/
external fex_getexcepthandler, fex_setexcepthandler
external fex_set_handling, handler
c$pragma c(fex_getexcepthandler, fex_setexcepthandler)
c$pragma c(fex_set_handling)
c x'ff2' = FEX_DIVBYZERO | FEX_INVALID
call fex_getexcepthandler(oldhdl, %val(ixff2))
c x'2' = FEX_DIVBYZERO, 0 = FEX_NONSTOP
call fex_set_handling(%val(ix2), %val(0), %val(0))
c x'b0' = FEX_INV_ZDZ | FEX_INV_IDI | FEX_INV_ZMI, 3 = FEX_CUSTOM
call fex_set_handling(%val(ixb0), %val(3), handler)
f1 = 0.0d0
f = a(n+1)
do j = n, 1, -1
d = x + f
d1 = 1.0d0 + f1
q = b(j) / d
f1 = (-d1 / d) * q
c
c the following assignment to the volatile variable t
c is needed to maintain the correct sequencing between
c assignments to p and evaluation of f1
t = f1
p = b(j-1) * d1 / b(j)
f = a(j) + q
end do
call fex_setexcepthandler(oldhdl, %val(ixff2))
return
end
c Main program
c
program cf
integer i
double precision a, b, x, f, f1
dimension a(5), b(4)
data a /-1.0d0, 2.0d0, -3.0d0, 4.0d0, -5.0d0/
data b /2.0d0, 4.0d0, 6.0d0, 8.0d0/
data ixffa/x'ffa'/
external fex_set_handling
c$pragma c(fex_set_handling)
c x'ffa' = FEX_COMMON, 1 = FEX_ABORT
call fex_set_handling(%val(ixffa), %val(1), %val(0))
do i = -5, 5
x = dble(i)
call continued_fraction(4, a, b, x, f, f1)
write (*, 1) i, f, i, f1
end do
1 format('f(', I2, ') = ', G12.6, ', f''(', I2, ') = ', G12.6)
end
–f77 フラグを指定してコンパイルされたこのプログラムの出力は次のようになります。
f(-5) = -1.59649 , f'(-5) = -.181800
f(-4) = -1.87302 , f'(-4) = -.428193
f(-3) = -3.00000 , f'(-3) = -3.16667
f(-2) = -.444089E-15, f'(-2) = -3.41667
f(-1) = -1.22222 , f'(-1) = -.444444
f( 0) = -1.33333 , f'( 0) = 0.203704
f( 1) = -1.00000 , f'( 1) = 0.333333
f( 2) = -.777778 , f'( 2) = 0.120370
f( 3) = -.714286 , f'( 3) = 0.272109E-01
f( 4) = -.666667 , f'( 4) = 0.203704
f( 5) = -.777778 , f'( 5) = 0.185185E-01
Note: IEEE floating-point exception flags raised:
Inexact; Division by Zero; Underflow; Invalid Operation;
IEEE floating-point exception traps enabled:
overflow; division by zero; invalid operation;
See the Numerical Computation Guide, ieee_flags(3M),
ieee_handler(3M)