From owner-freebsd-bugs@FreeBSD.ORG Wed Nov 10 12:52:59 2010 Return-Path: Delivered-To: freebsd-bugs@FreeBSD.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:4f8:fff6::34]) by hub.freebsd.org (Postfix) with ESMTP id 6B783106566C for ; Wed, 10 Nov 2010 12:52:59 +0000 (UTC) (envelope-from brde@optusnet.com.au) Received: from mail08.syd.optusnet.com.au (mail08.syd.optusnet.com.au [211.29.132.189]) by mx1.freebsd.org (Postfix) with ESMTP id E20CB8FC0A for ; Wed, 10 Nov 2010 12:52:58 +0000 (UTC) Received: from c122-107-121-73.carlnfd1.nsw.optusnet.com.au (c122-107-121-73.carlnfd1.nsw.optusnet.com.au [122.107.121.73]) by mail08.syd.optusnet.com.au (8.13.1/8.13.1) with ESMTP id oAACqrZu021196 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO); Wed, 10 Nov 2010 23:52:56 +1100 Date: Wed, 10 Nov 2010 23:52:53 +1100 (EST) From: Bruce Evans X-X-Sender: bde@besplex.bde.org To: Ulrich =?utf-8?B?U3DDtnJsZWlu?= In-Reply-To: <201011092010.oA9KABNt076837@freefall.freebsd.org> Message-ID: <20101110225059.M1461@besplex.bde.org> References: <201011092010.oA9KABNt076837@freefall.freebsd.org> MIME-Version: 1.0 Content-Type: MULTIPART/MIXED; BOUNDARY="0-1469788072-1289393573=:1461" Cc: freebsd-bugs@FreeBSD.org Subject: Re: bin/144306: [libm] [patch] Nasty bug in jn(3) X-BeenThere: freebsd-bugs@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list List-Id: Bug reports List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Wed, 10 Nov 2010 12:52:59 -0000 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--