IEEE 規格では、アンダーフローおよびオーバーフローがトラップされた場合、指数がラップされた結果 (つまり、指数がその通常の範囲の終わりでラップされていることを除けば、オーバーフローまたはアンダーフローした演算を丸めた結果と一致する値) をトラップハンドラが置換する方法をシステムが提供して、結果を 2 のべき乗によって位取りすることを推奨しています。以降の計算でアンダーフローやオーバーフローが発生しないように、指数範囲の中央のできるだけ近くにアンダーフローまたはオーバーフローした結果を割り当てるような位取りが選択されます。発生したアンダーフローやオーバーフローの回数を追跡することによって、プログラムは最終的な結果を位取りし、ラップされた指数を補正できます。このアンダーフロー/オーバーフローの「カウントモード」は、有効な浮動小数点形式の範囲を超えてしまうような計算において、正確な結果を生成するために使用できます。詳細は、『Floating-Point Computation』、P. Sterbenz 著を参照してください。
SPARC ベースのシステムでは、浮動小数点命令がトラップされた例外の原因である場合、システムは宛先レジスタを変更しません。このため、指数がラップされた結果を置換するには、アンダーフロー/オーバーフローのハンドラが命令をデコードし、オペランドレジスタを検査して、位取りされた結果自体を生成する必要があります。次の例では、この手順を実行するハンドラを示します。
使用例 4-4 SPARC ベースのシステムでの、IEEE トラップされたアンダーフロー/オーバーフローハンドラの結果の置換#include <stdio.h>
#include <ieeefp.h>
#include <math.h>
#include <sunmath.h>
#include <siginfo.h>
#include <ucontext.h>
#ifdef V8PLUS
/* The upper 32 floating point registers are stored in an area
pointed to by uap->uc_mcontext.xrs.xrs_prt. Note that this
pointer is valid ONLY when uap->uc_mcontext.xrs.xrs_id ==
XRS_ID (defined in sys/procfs.h). */
#include <assert.h>
#include <sys/procfs.h>
#define FPxreg(x) ((prxregset_t*)uap->uc_mcontext.xrs.xrs_ptr)
->pr_un.pr_v8p.pr_xfr.pr_regs[(x)]
#endif
#define FPreg(x) uap->uc_mcontext.fpregs.fpu_fr.fpu_regs[(x)]
/*
* Supply the IEEE 754 default result for trapped under/overflow
*/
void
ieee_trapped_default(int sig, siginfo_t *sip, ucontext_t *uap)
{
unsigned instr, opf, rs1, rs2, rd;
long double qs1, qs2, qd, qscl;
double ds1, ds2, dd, dscl;
float fs1, fs2, fd, fscl;
/* get the instruction that caused the exception */
instr = uap->uc_mcontext.fpregs.fpu_q->FQu.fpq.fpq_instr;
/* extract the opcode and source and destination register
numbers */
opf = (instr >> 5) & 0x1ff;
rs1 = (instr >> 14) & 0x1f;
rs2 = instr & 0x1f;
rd = (instr >> 25) & 0x1f;
/* get the operands */
switch (opf & 3) {
case 1: /* single precision */
fs1 = *(float*)&FPreg(rs1);
fs2 = *(float*)&FPreg(rs2);
break;
case 2: /* double precision */
#ifdef V8PLUS
if (rs1 & 1)
{
assert(uap->uc_mcontext.xrs.xrs_id == XRS_ID);
ds1 = *(double*)&FPxreg(rs1 & 0x1e);
}
else
ds1 = *(double*)&FPreg(rs1);
if (rs2 & 1)
{
assert(uap->uc_mcontext.xrs.xrs_id == XRS_ID);
ds2 = *(double*)&FPxreg(rs2 & 0x1e);
}
else
ds2 = *(double*)&FPreg(rs2);
#else
ds1 = *(double*)&FPreg(rs1);
ds2 = *(double*)&FPreg(rs2);
#endif
break;
case 3: /* quad precision */
#ifdef V8PLUS
if (rs1 & 1)
{
assert(uap->uc_mcontext.xrs.xrs_id == XRS_ID);
qs1 = *(long double*)&FPxreg(rs1 & 0x1e);
}
else
qs1 = *(long double*)&FPreg(rs1);
if (rs2 & 1)
{
assert(uap->uc_mcontext.xrs.xrs_id == XRS_ID);
qs2 = *(long double*)&FPxreg(rs2 & 0x1e);
}
else
qs2 = *(long double*)&FPreg(rs2);
#else
qs1 = *(long double*)&FPreg(rs1);
qs2 = *(long double*)&FPreg(rs2);
#endif
break;
}
/* set up scale factors */
if (sip->si_code == FPE_FLTOVF) {
fscl = scalbnf(1.0f, -96);
dscl = scalbn(1.0, -768);
qscl = scalbnl(1.0, -12288);
} else {
fscl = scalbnf(1.0f, 96);
dscl = scalbn(1.0, 768);
qscl = scalbnl(1.0, 12288);
}
/* disable traps and generate the scaled result */
fpsetmask(0);
switch (opf) {
case 0x41: /* add single */
fd = fscl * (fscl * fs1 + fscl * fs2);
break;
case 0x42: /* add double */
dd = dscl * (dscl * ds1 + dscl * ds2);
break;
case 0x43: /* add quad */
qd = qscl * (qscl * qs1 + qscl * qs2);
break;
case 0x45: /* subtract single */
fd = fscl * (fscl * fs1 - fscl * fs2);
break;
case 0x46: /* subtract double */
dd = dscl * (dscl * ds1 - dscl * ds2);
break;
case 0x47: /* subtract quad */
qd = qscl * (qscl * qs1 - qscl * qs2);
break;
case 0x49: /* multiply single */
fd = (fscl * fs1) * (fscl * fs2);
break;
case 0x4a: /* multiply double */
dd = (dscl * ds1) * (dscl * ds2);
break;
case 0x4b: /* multiply quad */
qd = (qscl * qs1) * (qscl * qs2);
break;
case 0x4d: /* divide single */
fd = (fscl * fs1) / (fs2 / fscl);
break;
case 0x4e: /* divide double */
dd = (dscl * ds1) / (ds2 / dscl);
break;
case 0x4f: /* divide quad */
qd = (qscl * qs1) / (qs2 / dscl);
break;
case 0xc6: /* convert double to single */
fd = (float) (fscl * (fscl * ds1));
break;
case 0xc7: /* convert quad to single */
fd = (float) (fscl * (fscl * qs1));
break;
case 0xcb: /* convert quad to double */
dd = (double) (dscl * (dscl * qs1));
break;
}
/* store the result in the destination */
if (opf & 0x80) {
/* conversion operation */
if (opf == 0xcb) {
/* convert quad to double */
#ifdef V8PLUS
if (rd & 1)
{
assert(uap->uc_mcontext.xrs.xrs_id == XRS_ID);
*(double*)&FPxreg(rd & 0x1e) = dd;
}
else
*(double*)&FPreg(rd) = dd;
#else
*(double*)&FPreg(rd) = dd;
#endif
} else
/* convert quad/double to single */
*(float*)&FPreg(rd) = fd;
} else {
/* arithmetic operation */
switch (opf & 3) {
case 1: /* single precision */
*(float*)&FPreg(rd) = fd;
break;
case 2: /* double precision */
#ifdef V8PLUS
if (rd & 1)
{
assert(uap->uc_mcontext.xrs.xrs_id == XRS_ID);
*(double*)&FPxreg(rd & 0x1e) = dd;
}
else
*(double*)&FPreg(rd) = dd;
#else
*(double*)&FPreg(rd) = dd;
#endif
break;
case 3: /* quad precision */
#ifdef V8PLUS
if (rd & 1)
{
assert(uap->uc_mcontext.xrs.xrs_id == XRS_ID);
*(long double*)&FPxreg(rd & 0x1e) = qd;
}
else
*(long double*)&FPreg(rd & 0x1e) = qd;
#else
*(long double*)&FPreg(rd & 0x1e) = qd;
#endif
break;
}
}
}
int
main()
{
volatile float a, b;
volatile double x, y;
ieee_handler("set", "underflow", ieee_trapped_default);
ieee_handler("set", "overflow", ieee_trapped_default);
a = b = 1.0e30f;
a *= b; /* overflow; will be wrapped to a moderate number */
printf( "%g\n", a );
a /= b;
printf( "%g\n", a );
a /= b; /* underflow; will wrap back */
printf( "%g\n", a );
x = y = 1.0e300;
x *= y; /* overflow; will be wrapped to a moderate number */
printf( "%g\n", x );
x /= y;
printf( "%g\n", x );
x /= y; /* underflow; will wrap back */
printf( "%g\n", x );
ieee_retrospective(stdout);
return 0;
}
この例で変数 a、b、x、および y が volatile と宣言されているのは、コンパイラがコンパイル時に a * b などを評価することを防ぐためにすぎません。通常の使用では、volatile 宣言は必要ありません。
上記のプログラムの出力は、次のとおりです。
159.309
1.59309e-28
1
4.14884e+137
4.14884e-163
1
Note: IEEE floating-point exception traps enabled:
underflow; overflow;
See the Numerical Computation Guide, ieee_handler(3M)
x86 ベースのシステムでは、浮動小数点命令によってアンダーフローまたはオーバーフローがトラップされ、その宛先がレジスタである場合、浮動小数点ハードウェアによって指数がラップされた結果が提供されます。ただし、トラップされたアンダーフローまたはオーバーフローが浮動小数点のストア命令で発生した場合、ハードウェアはストアを完了させずにトラップを行い、そのストア命令がストアおよびポップである場合は、スタックのポップも行いません。このため、カウントモードを実装するには、アンダーフロー/オーバーフローのハンドラが位取りされた結果を生成し、ストア命令でトラップが発生した場合は、スタックを修正する必要があります。Example 4–5 はそのようなハンドラを示しています。
使用例 4-5 x86 ベースのシステムでの IEEE トラップされたアンダーフロー/オーバーフローハンドラの結果の置換#include <stdio.h>
#include <ieeefp.h>
#include <math.h>
#include <sunmath.h>
#include <siginfo.h>
#include <ucontext.h>
/* offsets into the saved fp environment */
#define CW 0 /* control word */
#define SW 1 /* status word */
#define TW 2 /* tag word */
#define OP 4 /* opcode */
#define EA 5 /* operand address */
#define FPenv(x) uap->uc_mcontext.fpregs.fp_reg_set.
fpchip_state.state[(x)]
#define FPreg(x) *(long double *)(10*(x)+(char*)&uap->
uc_mcontext.fpregs.fp_reg_set.fpchip_state.state[7])/*
* Supply the IEEE 754 default result for trapped under/overflow
*/
void
ieee_trapped_default(int sig, siginfo_t *sip, ucontext_t *uap)
{
double dscl;
float fscl;
unsigned sw, op, top;
int mask, e;
/* preserve flags for untrapped exceptions */
sw = uap->uc_mcontext.fpregs.fp_reg_set.fpchip_state.status;
FPenv(SW) |= (sw & (FPenv(CW) & 0x3f));
/* if the excepting instruction is a store, scale the stack
top, store it, and pop the stack if need be */
fpsetmask(0);
op = FPenv(OP) >> 16;
switch (op & 0x7f8) {
case 0x110:
case 0x118:
case 0x150:
case 0x158:
case 0x190:
case 0x198:
fscl = scalbnf(1.0f, (sip->si_code == FPE_FLTOVF)?
-96 : 96);
*(float *)FPenv(EA) = (FPreg(0) * fscl) * fscl;
if (op & 8) {
/* pop the stack */
FPreg(0) = FPreg(1);
FPreg(1) = FPreg(2);
FPreg(2) = FPreg(3);
FPreg(3) = FPreg(4);
FPreg(4) = FPreg(5);
FPreg(5) = FPreg(6);
FPreg(6) = FPreg(7);
top = (FPenv(SW) >> 10) & 0xe;
FPenv(TW) |= (3 << top);
top = (top + 2) & 0xe;
FPenv(SW) = (FPenv(SW) & ~0x3800) | (top << 10);
}
break;
case 0x510:
case 0x518:
case 0x550:
case 0x558:
case 0x590:
case 0x598:
dscl = scalbn(1.0, (sip->si_code == FPE_FLTOVF)?
-768 : 768);
*(double *)FPenv(EA) = (FPreg(0) * dscl) * dscl;
if (op & 8) {
/* pop the stack */
FPreg(0) = FPreg(1);
FPreg(1) = FPreg(2);
FPreg(2) = FPreg(3);
FPreg(3) = FPreg(4);
FPreg(4) = FPreg(5);
FPreg(5) = FPreg(6);
FPreg(6) = FPreg(7);
top = (FPenv(SW) >> 10) & 0xe;
FPenv(TW) |= (3 << top);
top = (top + 2) & 0xe;
FPenv(SW) = (FPenv(SW) & ~0x3800) | (top << 10);
}
break;
}
}
int main()
{
volatile float a, b;
volatile double x, y;
ieee_handler("set", "underflow", ieee_trapped_default);
ieee_handler("set", "overflow", ieee_trapped_default);
a = b = 1.0e30f;
a *= b;
printf( "%g\n", a );
a /= b;
printf( "%g\n", a );
a /= b;
printf( "%g\n", a );
x = y = 1.0e300;
x *= y;
printf( "%g\n", x );
x /= y;
printf( "%g\n", x );
x /= y;
printf( "%g\n", x );
ieee_retrospective(stdout);
return 0;
}
SPARC ベースのシステムの場合および –xarch=386 を指定してコンパイルされている場合と同様に、x86 での上記のプログラムの出力は次のようになります。
159.309
1.59309e-28
1
4.14884e+137
4.14884e-163
1
Note: IEEE floating-point exception traps enabled:
underflow; overflow;
See the Numerical Computation Guide, ieee_handler(3M)
C/C++ プログラムでは、libm にある fex_set_handling 関数を使用すると、アンダーフローおよびオーバーフローに対する FEX_CUSTOM ハンドラをインストールできます。SPARC ベースのシステムでは、このようなハンドラに渡される情報には、例外の原因である演算およびオペランドが常に含まれており、上記に示したように、ハンドラは、この情報を使用して IEEE の指数がラップされた結果を計算できます。x86 ベースのシステムでは、例外の原因である特定の演算、および超越命令の 1 つが例外を発生させる時点 (たとえば、info->op パラメータが fex_other に設定された) を、提供される情報が常に示しているとは限りません(定義については、fenv.h ファイルを参照してください)。また、x86 のハードウェアは指数がラップされた結果を自動的に提供するため、例外を発生させている命令の宛先が浮動小数点レジスタである場合は、オペランドの 1 つが上書きされる場合があります。
幸いなことに、fex_set_handling 機能は、FEX_CUSTOM モードでインストールされているハンドラがアンダーフローまたはオーバーフローした演算を IEEE の指数がラップされた結果に置換する簡単な方法を提供しています。これらの例外のいずれかがトラップされたときに、ハンドラは次のように設定できます。
info->res.type = fex_nodata;
これにより、指数がラップされた結果を提供することが示されます。次に、そのようなハンドラの例を示します。
#include <stdio.h>
#include <fenv.h>
void handler(int ex, fex_info_t *info) {
info->res.type = fex_nodata;
}
int main()
{
volatile float a, b;
volatile double x, y;
fex_set_log(stderr);
fex_set_handling(FEX_UNDERFLOW | FEX_OVERFLOW, FEX_CUSTOM,
handler);
a = b = 1.0e30f;
a *= b; /* overflow; will be wrapped to a moderate number */
printf("%g\n", a);
a /= b;
printf("%g\n", a);
a /= b; /* underflow; will wrap back */
printf("%g\n", a);
x = y = 1.0e300;
x *= y; /* overflow; will be wrapped to a moderate number */
printf("%g\n", x);
x /= y;
printf("%g\n", x);
x /= y; /* underflow; will wrap back */
printf("%g\n", x);
return 0;
}
上記のプログラムの出力は、次のようになります。
Floating point overflow at 0x00010924 main, handler: handler 0x00010928 main 159.309 1.59309e-28 Floating point underflow at 0x00010994 main, handler: handler 0x00010998 main 1 Floating point overflow at 0x000109e4 main, handler: handler 0x000109e8 main 4.14884e+137 4.14884e-163 Floating point underflow at 0x00010a4c main, handler: handler 0x00010a50 main 1
上記のプログラムの例は、SPARC、–xarch=sse2 を指定した x86、および –xarch=386 を指定した x86 上で動作します。