log(3) と 1ulp

log(3.0) は 1.0986122886681096913952452369225 くらいの値になるはずです (Windows 付録の電卓調べ)。これがC言語の log(3) で log(3) した場合、double に丸められるので、実際の値は 1.0986122886681098 になります。

言い換えると、だいたい 0x1.193ea7aad030a97p+0 あたりの値を 0x1.193ea7aad030bp+0 に丸めるはずなのです。しかし、netlibFreeBSDNetBSD (i386 ではアセンブラを使うので影響しない) では log(3) は 1.0986122886681096、つまり、0x1.193ea7aad030ap+0 になります。
原因は

return dk*ln2_hi-((s*(f-R)-dk*ln2_lo)-f);

の計算において、(s*(f-R)-dk*ln2_lo) の値が 0.037682 (0x1.34b10865f9eacp-5) で、f が -0.250000 (-0x1p-2) になります。で、((s*(f-R)-dk*ln2_lo)-f) の結果は 0.287682 (0x1.2696210cbf3d58p-2) になり、丸めて、0.287682 0x1.2696210cbf3d6p-2 になります。
dk*ln2_hi は 1.386294 (0x1.62e42feep+0) なので、結果である、dk*ln2_hi-((s*(f-R)-dk*ln2_lo)-f) は 1.098612 (0x1.193ea7aad030a8p+0) になり、丸めて (銀行型丸めなのでこっちでは切り捨てる)、0x1.193ea7aad030ap+0 になります。これは結局、2回の丸めが両方とも結果を少なくする方向に働くために、1ulp 分少なくなってしまう。

log(3.0) に限れば、これは以下の通り書きかえると、f が -0.25 なこともあって丸めが1回だけになり精度が上がるのだが、副作用がある気もして怖い。

return (dk*ln2_hi+f)-(s*(f-R)-dk*ln2_lo);