しゅどうです。

正木さんは、x86 プロセッサが載ったマシンをお使いだということが判ります。
また、お使いの OS は FreeBSD 以外、かと。

正木さん wrote:

> 1.0/n は n=2731,4623,5462,5607,5963,6147,9246,・・・
> の場合に正しい(Float の中で最も近い)値を返しません。

僕も n = 2731 の場合に限って調査をしました。
1.0 / 2731 は、x86 の特殊な仕様にちょうど引っかかる演算でした。
x86 が他のプロセッサ (SPARC, PowerPC, MIPS, Alpha など) と
どう異なるのかは、次のウェブページに貼ってあるスライドを御覧ください:

  厳密な浮動小数点演算セマンティクスの Java 実行時コンパイラへの実装
  http://www.shudo.net/publications/swopp0107/


> これらの場合に、1/n の仮数部を64桁まで(と 1.0/n の仮数部を52
> 桁まで)調べてみると
> 1/2731
> 1.0111111111110100000000000101111111111101000000000001011111111111
> 1.0/2731
> 1.0111111111110100000000000101111111111101000000000010

x86 以外 (SPARC 等) では、これとは異なる結果が得られます↓。

  1.0/2731
  1.0111111111110100000000000101111111111101000000000001
                                                      ^^ 異なる箇所
また、x86 でも、
丸め時の精度を拡張精度 (80 bit) ではなくて倍精度 (64 bit) に設定して
演算すると、他のプロセッサと同じ結果になります。
FreeBSD 4.X では、はじめから倍精度 (64 bit) になっています。


x86 で妙な結果が得られる原因は、正木さんの推測通りです。
丸めが 2度起きてしまっています。

演算結果が FPU のレジスタに格納される時点で、
仮数部は 64 bit の精度をもってしまっています。
レジスタには、次の仮数部が入ります (1度目の丸め):

  1. 01111111 11110100 00000000 01011111 11111101 00000000 00011000 0000000
                                                               ^^^^ ^^^^^^^(1)
次に、この仮数部がメモリにストアされる時点で、
53 bit 精度 (IEEE 754 倍精度) に丸められます (2度目の丸め):

  1. 01111111 11110100 00000000 01011111 11111101 00000000 0010

> これから想像すると FPU の設計ミスかなにかで(本当は切捨てにすべき)
> 64 bit 目を 0 捨 1 入して
> 100000000000

その通りでして、上記 (1) のようになっています。

> この bug のせいで Float の誤差は少し大きくなりますので、
> 次のように再訂正します。

この x86 の挙動は、バグではなくて、
仕様であって、しかも IEEE 754 に準拠したものです。

(これが本当にプロセッサのバグだったら、
 世の中大変な騒ぎ (祭り?) になります。
 Pentium の除算用テーブルにバグが発見されたときは、
 盛り上がり (失礼) ましたねー。)


ではどうすればよいか、ですが、
丸め時の精度を、拡張精度のままにしておかず、明示的に倍精度に設定することで
今回のような 2度丸めは防止できます。
すると、x86 でも他の種類のプロセッサと同じ結果が得られるようになります。

(実は、x86 では、
 丸め時の精度を倍精度に設定しても、2度丸めは完全に防げません。
 他のプロセッサと異なる結果が得られてしまうという演算があります。)


Kazuyuki Shudo/首藤一幸   私をたばねないで あらせいとうの花のように
  shudo / computer.org   http://www.shudo.net/