bin/144306: [libm] [patch] Nasty bug in jn(3)
Bruce Evans
brde at optusnet.com.au
Wed Nov 10 12:52:59 UTC 2010
On Tue, 9 Nov 2010, Ulrich [utf-8] Spörlein wrote:
> Hello Steve, I saw your updated patch for jnf(3) on the mailing list and
> was wondering what the correct values for your test case would look like?
> I'm asking because on Ubuntu I get slightly different values. Especially
> disturbing is, that for FreeBSD the jn/jnf values are identical in this
> case.
The printfs (in the test program for doubles) don't have enough precision
for doubles.
> Ubuntu 10.04:
>
> # cc -o testjn -lm testjn.c && ./testjn
> 2 4.317548e-01
> 3 1.850071e-01
> 4 6.121505e-02
> 5 1.568141e-02
> 6 3.251921e-03
> 7 5.737757e-04
> 8 8.808225e-05
> 9 1.195869e-05
Ubuntu is unlikely to be more correct than FreeBSD :-), and is in fact
off by about 7% (~2**49 ulps ~= 500 tulps (tera-ulps)) for most except
the first according to hopefully-correct values (but excessively rounded
to be directly comparable with the above) produced by pari:
2 4.317548e-1
3 1.989999e-1
4 6.474667e-2
5 1.638924e-2
6 3.404818e-3
7 6.006884e-4
8 9.216579e-5
9 1.251727e-5
> # cc -o testjnf -lm testjnf.c && ./testjnf
> 2 4.317548e-01
> 3 2.126431e-01
> 4 6.348598e-02
> 5 1.606404e-02
> 6 3.341151e-03
> 7 5.894695e-04
> 8 9.044447e-05
> 9 1.228349e-05
Single precision is more accurate than double precision -- now the errors
are a measly couple of hundred tulps.
FreeBSD's maximum errors for the j* family are more like a few gulps
giga-ulps). That is poor, but is hard to fix, and the average error
is small. The above seems to show an enormous average error.
> FreeBSD -CURRENT:
>
> freebsd64# cc -o testjn -lm testjn.c && ./testjn
> 2 4.317548e-01
> 3 1.989999e-01
> 4 6.474667e-02
> 5 1.638924e-02
> 6 3.404818e-03
> 7 6.006884e-04
> 8 9.216579e-05
> 9 1.251727e-05
FreeBSD agrees perfectly with pari. This is not surprising, since only
7 significant decimal digits are printed and almost 15 should be correct
-- any error here would be an error of about 100 million decimal ulps
or a couple of binary gulps.
> freebsd64# cc -o testjnf -lm testjnf.c && ./testjnf
> 2 4.317548e-01
> 3 1.989999e-01
> 4 6.474666e-02
> 5 1.638924e-02
> 6 3.404818e-03
> 7 6.006881e-04
> 8 9.216577e-05
> 9 1.251727e-05
This is harder to explain. Now we should only expect 6 digits of accuracy,
but we are getting 7. Anyway, these are all perfectly rounded according to
pari.
> (both systems are amd64, I tried gcc and clang on FreeBSD, giving the
> same results.)
pari program:
---
\\ \r ftoe.gp
ftoe(x, n) =
{
local(fmt, k, len, ofmt, s, v);
ofmt = default(format);
fmt = concat("e0.", Str(n));
default(format, fmt);
s = Str(x + 0.0);
default(format, ofmt);
v = Vec(s);
len = matsize(v)[2];
s = "";
for (k = 1, len,
if (v[k] == "E",
v[k] = "e";
);
if (v[k] != " ",
s = concat(s, v[k]);
);
);
return (s);
}
\\ \r roundn.gp
roundn(x, n) =
{
local(er, ex, y);
if (x == 0, return (x););
ex = floor(log(abs(x)) / log(2));
y0 = x * 2^(n - ex - 1);
y = round(y0, &er);
if (y % 2 == 1,
if (y > 0 && y == y0 + 0.5, y = y - 1;);
if (y < 0 && y == y0 - 0.5, y = y + 1;);
);
return (y / 2^(n - ex - 1) + 0.0);
}
FLT_MANT_DIG = 24;
FLT_DIG = 6;
DBL_MANT_DIG = 53;
DBL_DIG = 15;
z = 2.4048255576957729;
print("Double precision with bad decimal field width:");
for (n = 2, 9, print(n," ",ftoe(roundn(besselj(n,z),DBL_MANT_DIG),7)));
print("Double precision:");
for (n = 2, 9, print(n," ",ftoe(roundn(besselj(n,z),DBL_MANT_DIG),DBL_DIG)));
print("Single precision with bad decimal field width:");
for (n = 2, 9, print(n," ",ftoe(roundn(besselj(n,z),DBL_MANT_DIG),7)));
print("Single precision:");
for (n = 2, 9, print(n," ",ftoe(roundn(besselj(n,z),FLT_MANT_DIG),FLT_DIG)));
---
Description:
- pari besselj(n,z) is C jn(n,z) (but arbitary precision, and complex z works)
- roundn(x,n) is a bde library function that rounds real x to n binary digits.
Only rounding to nearest is supported.
- ftoe(x,n) is a bde library function that formats real x in C %e format with
n decimal digits. This involves rounding to a decimal number;
only rounding to nearest is supported.
Output:
---
Double precision with bad decimal field width:
2 4.317548e-1
3 1.989999e-1
4 6.474667e-2
5 1.638924e-2
6 3.404818e-3
7 6.006884e-4
8 9.216579e-5
9 1.251727e-5
Double precision:
2 4.31754807019680e-1
3 1.98999905357691e-1
4 6.47466661641780e-2
5 1.63892432048059e-2
6 3.40481847202783e-3
7 6.00688365732954e-4
8 9.21657867053449e-5
9 1.25172709779615e-5
Single precision with bad decimal field width:
2 4.317548e-1
3 1.989999e-1
4 6.474667e-2
5 1.638924e-2
6 3.404818e-3
7 6.006884e-4
8 9.216579e-5
9 1.251727e-5
Single precision:
2 4.31755e-1
3 1.99000e-1
4 6.47467e-2
5 1.63892e-2
6 3.40482e-3
7 6.00688e-4
8 9.21658e-5
9 1.25173e-5
---
To examine for perfect or good enough (< 1 ulps of error and preferably
< 0.51 ulps of error), the binary values must be examined. Here they are:
pari program:
---
\\ \r ftoa.gp
ftoa(x, n) =
{
local(er, ex, hex, k, s, y);
if (x == 0, return (" 0"););
ex = floor(log(abs(x)) / log(2));
y = round(abs(x) * 2^(n - ex - 1), &er);
if (y == 2^n,
y = y / 2;
ex = ex + 1;
);
if (sign(x) == -1, s = "-", s = " ");
s = concat(s, "0x");
hex = ["0", "1", "2", "3", "4", "5", "6", "7",
"8", "9", "a", "b", "c", "d", "e", "f"];
forstep (k = 4 * truncate((n - 1) / 4), 0, -4,
s = concat(s, hex[1 + (y >> k) % 16]);
);
s = concat(s, ".0p");
return (concat(s, 1 + ex - n));
}
FLT_MANT_DIG = 24;
DBL_MANT_DIG = 53;
z = 2.4048255576957729;
print("Double precision:");
for (n = 2, 9, print(n," ",ftoa(besselj(n,z),DBL_MANT_DIG+4)));
print("Single precision:");
for (n = 2, 9, print(n," ",ftoa(besselj(n,z),FLT_MANT_DIG+4)));
---
Description:
- ftoa(x,n) is a bde library function that formats real x in a modified C %a
format with n binary digits, after first rounding x to n binary
digits. Only rounding to nearest is supported.
- Note that 4 extra binary digits are printed, so that you can see the rounding
error to within 1/16 ulps instead of to within 1 ulp after completing the
rounding manually and looking at the lost digits.
- The modified %a format involves lining up the nybbles so that the values
look like raw binary ones. I think C allows this variant and the usual
variant is a bad one. In both, there is not anough control over the
location of the (binary or decimal) point, so the value in each nybble
tends to shift too often. I get enough control for purposes like the
above by printing a number of extra digits that is the multiple of 4.
Output:
---
Double precision:
2 0x1ba1deea029493f.0p-58
3 0x1978d432b58d635.0p-59
4 0x10933ccdb33e9d1.0p-60
5 0x10c8577e4888049.0p-62
6 0x1be46bff8ec1777.0p-65
7 0x13aef0716585c9f.0p-67
8 0x18292428a9a0003.0p-70
9 0x1a40289fa016035.0p-73
Single precision:
2 0xdd0ef75.0p-29
3 0xcbc6a19.0p-30
4 0x8499e67.0p-31
5 0x8642bbf.0p-33
6 0xdf23600.0p-36
7 0x9d77839.0p-38
8 0xc149214.0p-41
9 0xd201450.0p-44
---
for (n = 2, 9, print(n," ",ftoa(besselj(n,z),FLT_MANT_DIG+4)));
Bruce
More information about the freebsd-bugs
mailing list