Date: Wed, 10 Nov 2010 23:52:53 +1100 (EST) From: Bruce Evans <brde@optusnet.com.au> To: Ulrich =?utf-8?B?U3DDtnJsZWlu?= <uqs@spoerlein.net> Cc: freebsd-bugs@FreeBSD.org Subject: Re: bin/144306: [libm] [patch] Nasty bug in jn(3) Message-ID: <20101110225059.M1461@besplex.bde.org> In-Reply-To: <201011092010.oA9KABNt076837@freefall.freebsd.org> References: <201011092010.oA9KABNt076837@freefall.freebsd.org>
next in thread | previous in thread | raw e-mail | index | archive | help
This message is in MIME format. The first part should be readable text,
while the remaining parts are likely unreadable without MIME-aware tools.
--0-1469788072-1289393573=:1461
Content-Type: TEXT/PLAIN; charset=X-UNKNOWN; format=flowed
Content-Transfer-Encoding: QUOTED-PRINTABLE
On Tue, 9 Nov 2010, Ulrich [utf-8] Sp=C3=B6rlein 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 ~=3D 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) =3D
{
=09local(fmt, k, len, ofmt, s, v);
=09ofmt =3D default(format);
=09fmt =3D concat("e0.", Str(n));
=09default(format, fmt);
=09s =3D Str(x + 0.0);
=09default(format, ofmt);
=09v =3D Vec(s);
=09len =3D matsize(v)[2];
=09s =3D "";
=09for (k =3D 1, len,
=09=09if (v[k] =3D=3D "E",
=09=09=09v[k] =3D "e";
=09=09);
=09=09if (v[k] !=3D " ",
=09=09=09s =3D concat(s, v[k]);
=09=09);
=09);
=09return (s);
}
\\ \r roundn.gp
roundn(x, n) =3D
{
=09local(er, ex, y);
=09if (x =3D=3D 0, return (x););
=09ex =3D floor(log(abs(x)) / log(2));
=09y0 =3D x * 2^(n - ex - 1);
=09y =3D round(y0, &er);
=09if (y % 2 =3D=3D 1,
=09=09if (y > 0 && y =3D=3D y0 + 0.5, y =3D y - 1;);
=09=09if (y < 0 && y =3D=3D y0 - 0.5, y =3D y + 1;);
=09);
=09return (y / 2^(n - ex - 1) + 0.0);
}
FLT_MANT_DIG=09=3D=0924;
FLT_DIG=09=09=3D=096;
DBL_MANT_DIG=09=3D=0953;
DBL_DIG=09=09=3D=0915;
z =3D 2.4048255576957729;
print("Double precision with bad decimal field width:");
for (n =3D 2, 9, print(n," ",ftoe(roundn(besselj(n,z),DBL_MANT_DIG),7)));
print("Double precision:");
for (n =3D 2, 9, print(n," ",ftoe(roundn(besselj(n,z),DBL_MANT_DIG),DBL_DIG=
)));
print("Single precision with bad decimal field width:");
for (n =3D 2, 9, print(n," ",ftoe(roundn(besselj(n,z),DBL_MANT_DIG),7)));
print("Single precision:");
for (n =3D 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 wor=
ks)
- roundn(x,n) is a bde library function that rounds real x to n binary digi=
ts.
Only rounding to nearest is supported.
- ftoe(x,n) is a bde library function that formats real x in C %e format wi=
th
n decimal digits. This involves rounding to a decimal number;
=09 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) =3D
{
=09local(er, ex, hex, k, s, y);
=09if (x =3D=3D 0, return (" 0"););
=09ex =3D floor(log(abs(x)) / log(2));
=09y =3D round(abs(x) * 2^(n - ex - 1), &er);
=09if (y =3D=3D 2^n,
=09=09y =3D y / 2;
=09=09ex =3D ex + 1;
=09);
=09if (sign(x) =3D=3D -1, s =3D "-", s =3D " ");
=09s =3D concat(s, "0x");
=09hex =3D ["0", "1", "2", "3", "4", "5", "6", "7",
=09 "8", "9", "a", "b", "c", "d", "e", "f"];
=09forstep (k =3D 4 * truncate((n - 1) / 4), 0, -4,
=09=09s =3D concat(s, hex[1 + (y >> k) % 16]);
=09);
=09s =3D concat(s, ".0p");
=09return (concat(s, 1 + ex - n));
}
FLT_MANT_DIG=09=3D=0924;
DBL_MANT_DIG=09=3D=0953;
z =3D 2.4048255576957729;
print("Double precision:");
for (n =3D 2, 9, print(n," ",ftoa(besselj(n,z),DBL_MANT_DIG+4)));
print("Single precision:");
for (n =3D 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 binar=
y
=09 digits. Only rounding to nearest is supported.
- Note that 4 extra binary digits are printed, so that you can see the roun=
ding
error to within 1/16 ulps instead of to within 1 ulp after completing th=
e
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 =3D 2, 9, print(n," ",ftoa(besselj(n,z),FLT_MANT_DIG+4)));
Bruce
--0-1469788072-1289393573=:1461--
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20101110225059.M1461>