From owner-freebsd-standards@FreeBSD.ORG Sun Feb 21 22:44:16 2010 Return-Path: Delivered-To: freebsd-standards@FreeBSD.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:4f8:fff6::34]) by hub.freebsd.org (Postfix) with ESMTP id 1B3D61065672 for ; Sun, 21 Feb 2010 22:44:16 +0000 (UTC) (envelope-from joerg@britannica.bec.de) Received: from www.sonnenberger.org (www.sonnenberger.org [92.79.50.50]) by mx1.freebsd.org (Postfix) with ESMTP id D2C148FC1B for ; Sun, 21 Feb 2010 22:44:15 +0000 (UTC) Received: from britannica.bec.de (www.sonnenberger.org [192.168.1.10]) by www.sonnenberger.org (Postfix) with ESMTP id 77F7D6678A for ; Sun, 21 Feb 2010 23:28:42 +0100 (CET) Received: by britannica.bec.de (Postfix, from userid 1000) id EE9FA15C6C; Sun, 21 Feb 2010 23:28:11 +0100 (CET) Date: Sun, 21 Feb 2010 23:28:11 +0100 From: Joerg Sonnenberger To: freebsd-standards@FreeBSD.org Message-ID: <20100221222811.GA10638@britannica.bec.de> MIME-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline User-Agent: Mutt/1.5.20 (2009-06-14) Cc: Subject: UTF-8 and wchar_t X-BeenThere: freebsd-standards@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list List-Id: Standards compliance List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Sun, 21 Feb 2010 22:44:16 -0000 Hi all, reviewing some libarchive code I stumpled about the code that converts UTF-8 to wide strings. As done by a lot of other software, it currently blindly assumes that wchar_t ~= UCS-4. My question is whether FreeBSD intentionally makes that decision what (and therefore should define __STDC_ISO_10646__ according to ISO C99) or what correct way for reading UTF-8 it allows. Contrary to NetBSD, FreeBSD still lacks iconv(3) support in base, so the usual approach of converting to the locale charset and using mbtowc etc. is not possible. Joerg PS: Please keep me in CC. From owner-freebsd-standards@FreeBSD.ORG Sun Feb 21 23:10:45 2010 Return-Path: Delivered-To: freebsd-standards@FreeBSD.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:4f8:fff6::34]) by hub.freebsd.org (Postfix) with ESMTP id D45631065788 for ; Sun, 21 Feb 2010 23:10:45 +0000 (UTC) (envelope-from gabor@FreeBSD.org) Received: from server.mypc.hu (server.mypc.hu [87.229.73.95]) by mx1.freebsd.org (Postfix) with ESMTP id 8E35A8FC23 for ; Sun, 21 Feb 2010 23:10:45 +0000 (UTC) Received: from server.mypc.hu (localhost [127.0.0.1]) by server.mypc.hu (Postfix) with ESMTP id 981AA14DA5F2; Sun, 21 Feb 2010 23:53:22 +0100 (CET) X-Virus-Scanned: amavisd-new at server.mypc.hu Received: from server.mypc.hu ([127.0.0.1]) by server.mypc.hu (server.mypc.hu [127.0.0.1]) (amavisd-new, port 10024) with LMTP id 9U96Fbzxp+dF; Sun, 21 Feb 2010 23:53:16 +0100 (CET) Received: from [192.168.1.105] (catv-89-132-179-104.catv.broadband.hu [89.132.179.104]) (using TLSv1 with cipher DHE-RSA-CAMELLIA256-SHA (256/256 bits)) (No client certificate requested) by server.mypc.hu (Postfix) with ESMTPSA id D0CEF14DA5E9; Sun, 21 Feb 2010 23:53:16 +0100 (CET) Message-ID: <4B81B95B.9090304@FreeBSD.org> Date: Sun, 21 Feb 2010 23:53:15 +0100 From: Gabor Kovesdan User-Agent: Mozilla/5.0 (Windows; U; Windows NT 5.2; es-ES; rv:1.9.1.7) Gecko/20100111 Thunderbird/3.0.1 MIME-Version: 1.0 To: Joerg Sonnenberger References: <20100221222811.GA10638@britannica.bec.de> In-Reply-To: <20100221222811.GA10638@britannica.bec.de> Content-Type: text/plain; charset=ISO-8859-1; format=flowed Content-Transfer-Encoding: 7bit Cc: freebsd-standards@FreeBSD.org Subject: Re: UTF-8 and wchar_t X-BeenThere: freebsd-standards@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list List-Id: Standards compliance List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Sun, 21 Feb 2010 23:10:45 -0000 > > Contrary to NetBSD, FreeBSD still lacks iconv(3) > support in base, so the usual approach of converting to the locale > charset and using mbtowc etc. is not possible. > Probably, we are near to get our own iconv(3) for base. Current TODO is: - style(9) cleanups - make a working patch for HEAD; atm it is only possible to build it independently but it should go into libc - exp-run by portmgr So hopefully, we will be able to change to the usual approach soon. -- Gabor Kovesdan FreeBSD Volunteer EMAIL: gabor@FreeBSD.org .:|:. gabor@kovesdan.org WEB: http://people.FreeBSD.org/~gabor .:|:. http://kovesdan.org From owner-freebsd-standards@FreeBSD.ORG Mon Feb 22 11:07:07 2010 Return-Path: Delivered-To: freebsd-standards@FreeBSD.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:4f8:fff6::34]) by hub.freebsd.org (Postfix) with ESMTP id B129A106568B for ; Mon, 22 Feb 2010 11:07:07 +0000 (UTC) (envelope-from owner-bugmaster@FreeBSD.org) Received: from freefall.freebsd.org (freefall.freebsd.org [IPv6:2001:4f8:fff6::28]) by mx1.freebsd.org (Postfix) with ESMTP id 8651E8FC18 for ; Mon, 22 Feb 2010 11:07:07 +0000 (UTC) Received: from freefall.freebsd.org (localhost [127.0.0.1]) by freefall.freebsd.org (8.14.3/8.14.3) with ESMTP id o1MB77NC039840 for ; Mon, 22 Feb 2010 11:07:07 GMT (envelope-from owner-bugmaster@FreeBSD.org) Received: (from gnats@localhost) by freefall.freebsd.org (8.14.3/8.14.3/Submit) id o1MB76pl039838 for freebsd-standards@FreeBSD.org; Mon, 22 Feb 2010 11:07:06 GMT (envelope-from owner-bugmaster@FreeBSD.org) Date: Mon, 22 Feb 2010 11:07:06 GMT Message-Id: <201002221107.o1MB76pl039838@freefall.freebsd.org> X-Authentication-Warning: freefall.freebsd.org: gnats set sender to owner-bugmaster@FreeBSD.org using -f From: FreeBSD bugmaster To: freebsd-standards@FreeBSD.org Cc: Subject: Current problem reports assigned to freebsd-standards@FreeBSD.org X-BeenThere: freebsd-standards@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list List-Id: Standards compliance List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Mon, 22 Feb 2010 11:07:07 -0000 Note: to view an individual PR, use: http://www.freebsd.org/cgi/query-pr.cgi?pr=(number). The following is a listing of current problems submitted by FreeBSD users. These represent problem reports covering all versions including experimental development code and obsolete releases. S Tracker Resp. Description -------------------------------------------------------------------------------- o stand/143358 standards [libm] nearbyint(3) raises spurious inexact exception o stand/142803 standards j0 Bessel function inaccurate near zeros of the functi o stand/142255 standards scandir prototype in dirent.h isn't compliant with POS s stand/141705 standards [libc] [request] libc lacks cexp (and friends) o stand/135307 standards Boot Loader problem on Acer Aspire 5735 o stand/130067 standards Wrong numeric limits in system headers? o stand/129524 standards FreeBSD 7.0 isnt detecting my hardrives with raid5 o stand/129196 standards Inconsistent errno in strtol() o bin/125855 standards sh(1) allows for multiline, non-escaped control struct o stand/124860 standards flockfile(3) doesn't work when the memory has been exh o stand/123688 standards POSIX standard changes in unistd.h and grp.h o stand/121921 standards [patch] Add leap second support to at(1), atrun(8) o stand/121568 standards [patch] ln(1): wrong "ln -s" behaviour o stand/120947 standards xsm ignores system.xsm and .xsmstartup o stand/116826 standards [patch] sh support for POSIX character classes o stand/116477 standards rm(1): rm behaves unexpectedly when using -r and relat o bin/116413 standards incorrect getconf(1) handling of unsigned constants gi o stand/116081 standards make does not work with the directive sinclude p stand/107561 standards [libc] [patch] [request] Missing SUS function tcgetsid o stand/104743 standards [headers] [patch] Wrong values for _POSIX_ minimal lim o stand/100017 standards [Patch] Add fuser(1) functionality to fstat(1) o stand/96236 standards [patch] [posix] sed(1) incorrectly describes a functio o stand/96016 standards [headers] clock_getres et al should be in o stand/94729 standards [libc] fcntl() throws undocumented ENOTTY o kern/93705 standards [headers] [patch] ENODATA and EGREGIOUS (for glibc com o stand/92362 standards [headers] [patch] Missing SIGPOLL in kernel headers a stand/86484 standards [patch] mkfifo(1) uses wrong permissions o stand/83845 standards [libm] [patch] add log2() and log2f() support for libm o stand/82654 standards C99 long double math functions are missing o stand/81287 standards [patch] fingerd(8) might send a line not ending in CRL a stand/80293 standards sysconf() does not support well-defined unistd values o stand/79056 standards [feature request] [atch] regex(3) regression tests o stand/70813 standards [patch] ls(1) not Posix compliant o stand/66357 standards make POSIX conformance problem ('sh -e' & '+' command- s kern/64875 standards [libc] [patch] [request] add a system call: fdatasync( s stand/62858 standards malloc(0) not C99 compliant o stand/56476 standards cd9660 unicode support simple hack o stand/54839 standards [pcvt] pcvt deficits o stand/54833 standards [pcvt] more pcvt deficits o stand/54410 standards one-true-awk not POSIX compliant (no extended REs) o stand/46119 standards Priority problems for SCHED_OTHER using pthreads o stand/44425 standards getcwd() succeeds even if current dir has perm 000. p stand/41576 standards POSIX compliance of ln(1) o stand/39256 standards snprintf/vsnprintf aren't POSIX-conformant for strings s stand/36076 standards Implementation of POSIX fuser command o kern/27835 standards [libc] execve() doesn't conform to execve(2) spec in s a docs/26003 standards getgroups(2) lists NGROUPS_MAX but not syslimits.h s stand/24590 standards timezone function not compatible witn Single Unix Spec o bin/24390 standards ln(1) Replacing old dir-symlinks when using /bin/ln o stand/21519 standards sys/dir.h should be deprecated some more s bin/14925 standards getsubopt isn't poisonous enough 51 problems total. From owner-freebsd-standards@FreeBSD.ORG Tue Feb 23 00:21:23 2010 Return-Path: Delivered-To: freebsd-standards@FreeBSD.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:4f8:fff6::34]) by hub.freebsd.org (Postfix) with ESMTP id 43EE2106566B; Tue, 23 Feb 2010 00:21:23 +0000 (UTC) (envelope-from kargl@troutmask.apl.washington.edu) Received: from troutmask.apl.washington.edu (troutmask.apl.washington.edu [128.208.78.105]) by mx1.freebsd.org (Postfix) with ESMTP id 07C118FC15; Tue, 23 Feb 2010 00:21:23 +0000 (UTC) Received: from troutmask.apl.washington.edu (localhost.apl.washington.edu [127.0.0.1]) by troutmask.apl.washington.edu (8.14.4/8.14.4) with ESMTP id o1N0LMdT002354; Mon, 22 Feb 2010 16:21:22 -0800 (PST) (envelope-from kargl@troutmask.apl.washington.edu) Received: (from kargl@localhost) by troutmask.apl.washington.edu (8.14.4/8.14.4/Submit) id o1N0LLNb002353; Mon, 22 Feb 2010 16:21:21 -0800 (PST) (envelope-from kargl) From: "Steven G. Kargl" Message-Id: <201002230021.o1N0LLNb002353@troutmask.apl.washington.edu> In-Reply-To: <20100115110443.Y63344@delplex.bde.org> To: Bruce Evans Date: Mon, 22 Feb 2010 16:21:21 -0800 (PST) X-Mailer: ELM [version 2.4ME+ PL124d (25)] MIME-Version: 1.0 Content-Transfer-Encoding: 7bit Content-Type: text/plain; charset="US-ASCII" Cc: FreeBSD-gnats-submit@FreeBSD.org, freebsd-standards@FreeBSD.org Subject: Re: standards/142803: j0 Bessel function inaccurate near zeros of the function X-BeenThere: freebsd-standards@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list List-Id: Standards compliance List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Tue, 23 Feb 2010 00:21:23 -0000 Bruce Evans wrote: > On Thu, 14 Jan 2010, Steven G. Kargl wrote: > > > Bruce Evans wrote: > >> On Wed, 13 Jan 2010, Steven G. Kargl wrote: > >> > >>>> Description: > >>> > >>> The j0 bessel function supplied by libm is fairly inaccurate at > >>> ... > >> > >> This is a very hard and relatively unimportant problem. > > > > Yes, it is very hard, but apparently you do not use bessel > > functions in your everyday life. :) > > > > I only discover this issue because I need bessel functions > > of complex arguments and I found my routines have issues > > in the vicinity of zeros. So, I decided to look at the > > libm routines. > > It is interesting that the often-poor accuracy of almost every system's > libm matters in real life. > > Complex args are another interesting problem since even complex > multiplication is hard to do accurately (it may have large cancelation > errors). > > >> Anyway, if you can get anywhere near < 10 ulp error near all zeros using > >> only an asymptotic method, then that would be good. Then the asymptotic > >> method would also be capable of locating the zeros very accurately. But > >> I would be very surprised if this worked. I know of nothing similar for > >> reducing mod Pi for trigonometric functions, which seems a simpler problem. > >> I would expect it to at best involve thousands of binary digits in the > >> tables for the asymptotic method, and corresponding thousands of digits > >> of precision in the calculation (4000 as for mfpr enough for the 2**100th > >> zero?). > > > > The 4000-bit setting for mpfr was a hold over from testing mpfr_j0 > > against my ascending series implementation of j0 with mpfr > > primitives. As few as 128-bits is sufficient to achieve the > > following: > > > >>> x my j0f(x) libm j0f(x) MPFR j0 my err libm err > > 2.404825 5.6434398E-08 5.9634296E-08 5.6434400E-08 0.05 152824.59 > > 5.520078 2.4476657E-08 2.4153294E-08 2.4476659E-08 0.10 18878.52 > > 8.653728 1.0355303E-07 1.0359805E-07 1.0355306E-07 0.86 1694.47 > > 11.791534 -3.5291243E-09 -3.5193941E-09 -3.5301714E-09 75.93 781.53 > > Wonder why this jumps. > > > 14.930918 -6.4815082E-09 -6.3911618E-09 -6.4815052E-09 0.23 6722.88 > > 18.071064 5.1532352E-09 5.3149818E-09 5.1532318E-09 0.23 10910.50 > > 21.211637 -1.5023349E-07 -1.5002509E-07 -1.5023348E-07 2.70 56347.01 > > ... > > > > As I suspected by adding additional terms to the asymptotic > > approximation and performing all computations with double > > precision, reduces 'my err' (5th column). The value at > > x=11.7... is the best I can get. The asymptotic approximations > > contain divergent series and additional terms do not help. > > The extra precision is almost certainly necessary. Whether double > precision is nearly enough is unclear, but the error near 11.7 suggests > that it is nearly enough except there. The large error might be caused > by that zero alone (among small zeros) being very close to a representable > value. > > I forgot the mention that the error table in my previous mail is on amd64 > for comparing float precision functions with double precision ones, assuming > that the latter are correct, which they aren't, but they are hopefully > correct enough for this comparision. The errors on i386 are much larger, > due to i386 still using i387 hardware trigonometric which are extremely > inaccurate near zeros, starting at the first zero. Here are both tables: > > amd64: > %%% > j0: max_er = 0x7fffffffffdf5e07 17179869183.9960, avg_er = 5.581, #>=1:0.5 = 1593961230:1839722934 > j1: max_er = 0x7fffffffffbcd2a1 17179869183.9918, avg_er = 4.524, #>=1:0.5 = 1597678928:1856295142 > lgamma:max_er = 0x135a0b77e00000 10145883.7461, avg_er = 0.252, #>=1:0.5 = 44084256:331444835 > y0: max_er = 0x7fffffffffbcd2a0 17179869183.9918, avg_er = 2.379, #>=1:0.5 = 837057577:1437331064 > y1: max_er = 0x7fffffffffdf5e07 17179869183.9960, avg_er = 3.761, #>=1:0.5 = 865063612:1460264955 > %%% > > i386: > %%% > j0: max_er = 0x6f8686c5c9f26 3654467.3863, avg_er = 0.418, #>=1:0.5 = 671562746:1332944948 > j1: max_er = 0x449026e9c286f 2246675.4566, avg_er = 0.425, #>=1:0.5 = 674510414:1347770568 > lgamma:max_er = 0xe4cf242400000 7497618.0703, avg_er = 0.274, #>=1:0.5 = 70033452:508702222 > y0: max_er = 0x93a2340c00000 4837658.0234, avg_er = 0.380, #>=1:0.5 = 594207097:1303826516 > y1: max_er = 0x7ffa2068256f72ab 17176789825.1699, avg_er = 5.188, #>=1:0.5 = 459137173:1213136103 > %%% > > Unfortunately, most of these i386 errors, and the amd64 error for y1() > are misreported. The huge max_er of 0x7ff... (16 hex digits) is > actually a sign error misreported. Sign errors are bad enough. They > always result at a simple zero z0 (f'(z0) != 0) when the approximation > is so inaccurate that it cannot tell on which side of infinite-precision > z0 the parameter is. Methods involving a table of zeros will not have > any of these provided the table is accurate to within 1 ulp, but other > methods can easily have them, depending on the other method's ability > to locate the zeros to within 1 ulp by solving f~(z) ~= 0 where f~ is > the approximate f. > > Having no sign errors across the whole range seems too good to believe > for the amd64 functions. All of the i387 hardware trig functions have > sign errors. > Bruce, David, Since the my real-life work requires additional accuracy for the Bessel function routines in libm, I spent sometime playing with e_j0f. I've come up with the patch that follows my signature. It only deals with the first 32 zeros of j0(x), others can be added to the table. I've yet to come up with a general method to improve the accuracy around all zeros. To help you gauge the improvement, one can look at http://troutmask.apl.washington.edu/~kargl/j0f.jpg http://troutmask.apl.washington.edu/~kargl/j0f_1.jpg which compares the unpatched j0f(x) to the patched j0f(x). The patch makes use of the addition theorem for J0(x+d) where x is a zero and d is some small deviation from the zero. Using only the first 4 terms is sufficent for j0f. J0(x+d) = - J1(x)*J1(d) + J2(x)*J2(d) - J3(x)*J3(d) + J4(x)*J4(d) where J1(d), J2(d), J3(d) and J4(d) can be approximated by truncating their series representations. The functions with x, ie the known zero, are evaluated with MPFR. Note, d is determined from d = z - x = z - high(x) - low(x) where high(x) is the 1st 24 bits of the zero, low(x) is the next 24 bits of zero, and z is some value in [x-0.065,x+0.065]. Bruce (and/or David), if you think that this merits inclusion in FreeBSD's libm, I'll work up a similar patch for e_j0.c and submit both; otherwise, I keep it as a local change. -- Steve http://troutmask.apl.washington.edu/~kargl/ Index: src/e_j0f.c =================================================================== --- src/e_j0f.c (revision 204219) +++ src/e_j0f.c (working copy) @@ -22,6 +22,8 @@ static float pzerof(float), qzerof(float); static const float +invpi = 0.31830988618379067154, +pi = 3.14159265358979323846, huge = 1e30, one = 1.0, invsqrtpi= 5.6418961287e-01, /* 0x3f106ebb */ @@ -38,17 +40,157 @@ static const float zero = 0.0; +/* J1(x) for x << 1. */ +static float +_j1(float x2) +{ + return (x2 * (1 - 0.5 * x2 * x2)); +} + +/* J2(x) for x << 1. */ +static float +_j2(float x2) +{ + x2 *= x2; + return (0.5 * x2 * (1 - x2 / 3)); +} + +/* J3(x) for x << 1. */ +static float +_j3(float x2) +{ + float y; + y = x2 * x2; + return (x2 * y * (1 - 0.25 * y) / 6); +} + +/* J4(x) for x << 1. */ +static float +_j4(float x2) +{ + x2 *= x2; + return (x2 * x2 / 24); +} + +/* + * Use the addition theorem for values near a zero the Bessel function. + * J0(z) = J0(x + d) where |d| << 1. + * = - 2*J1(x)*J1(d) + 2*J2(x)*J2(d) - ... + * Note, the factor of 2 is included in the table. + */ float +j0zerof(int n, float z) +{ + static const float x[32][6] = { + {0x1.33d152p+1, 0x1.d2e368p-24, 0x1.09cdb4p+0, 0x1.ba1deep-1, + 0x1.978d44p-2, 0x1.0933ccp-3}, + {0x1.6148f6p+2, -0x1.34f46ep-24, -0x1.5c6e6p-1, -0x1.f8f72ep-3, + 0x1.00f404p-1, 0x1.9588cep-1}, + {0x1.14eb56p+3, 0x1.999bdap-22, 0x1.15f798p-1, 0x1.00f7fcp-3, + -0x1.f08b9p-2, -0x1.d8c2a8p-2}, + {0x1.79544p+3, 0x1.04e56cp-26, -0x1.dc13e6p-2, -0x1.42ff0cp-4, + 0x1.c0af7ep-2, 0x1.350edcp-2}, + {0x1.ddca14p+3, -0x1.0d8e2ep-25, 0x1.a701dp-2, 0x1.c54b94p-5, + -0x1.97d3ccp-2, -0x1.b9186p-3}, + {0x1.212314p+4, -0x1.d7982ap-26, -0x1.8077f6p-2, -0x1.5467ecp-5, + 0x1.770cdp-2, 0x1.4e26dp-3}, + {0x1.5362dep+4, -0x1.d1810ep-21, 0x1.62d93ap-2, 0x1.0ba9cep-5, + -0x1.5c8a0ap-2, -0x1.08180cp-3}, + {0x1.85a3bap+4, -0x1.9fd524p-21, -0x1.4b2a2ep-2, -0x1.b3298p-6, + 0x1.46b28cp-2, 0x1.aec26ap-4}, + {0x1.b7e54ap+4, 0x1.7f57c4p-22, 0x1.37aac8p-2, 0x1.6ac0d2p-6, + -0x1.345e5cp-2, -0x1.67dfb2p-4}, + {0x1.ea275ap+4, -0x1.c68826p-21, -0x1.27407ep-2, -0x1.34695p-6, + 0x1.24bc2ep-2, 0x1.32708ap-4}, + {0x1.0e34e2p+5, -0x1.8b3204p-20, 0x1.192f24p-2, 0x1.0a6682p-6, + -0x1.17365ap-2, -0x1.08ffd2p-4}, + {0x1.275638p+5, -0x1.5a7986p-21, -0x1.0cf3eep-2, -0x1.d242aap-7, + 0x1.0b5fc4p-2, 0x1.d0352cp-5}, + {0x1.4077a8p+5, -0x1.29d6c6p-23, 0x1.0230bap-2, 0x1.9c8084p-7, + -0x1.00e734p-2, -0x1.9af5aap-5}, + {0x1.59992cp+5, 0x1.974364p-21, -0x1.f13fbp-3, -0x1.70558ep-7, + 0x1.ef1ep-3, 0x1.6f2666p-5}, + {0x1.72bacp+5, 0x1.f02102p-20, 0x1.e018dap-3, 0x1.4b858ap-7, + -0x1.de4fp-3, -0x1.4a986ap-5}, + {0x1.8bdc62p+5, 0x1.27e0cap-20, -0x1.d09b22p-3, -0x1.2c74f6p-7, + 0x1.cf1686p-3, 0x1.2bb87cp-5}, + {0x1.a4fe0ep+5, 0x1.c8899p-20, 0x1.c28612p-3, 0x1.11f526p-7, + -0x1.c138e4p-3, -0x1.115d32p-5}, + {0x1.be1fc4p+5, 0x1.a4c606p-23, -0x1.b5a622p-3, -0x1.f645fep-8, + 0x1.b485eap-3, 0x1.f54de8p-6}, + {0x1.d7418p+5, 0x1.93c83ep-20, 0x1.a9d184p-3, 0x1.cea254p-8, + -0x1.a8d632p-3, -0x1.cdd58ap-6}, + {0x1.f06344p+5, -0x1.7b4716p-22, -0x1.9ee5eep-3, -0x1.abf28ap-8, + 0x1.9e093ap-3, 0x1.ab47cep-6}, + {0x1.04c286p+6, 0x1.0f88f2p-21, 0x1.94c6f6p-3, 0x1.8d6372p-8, + -0x1.9403e4p-3, -0x1.8cd3dp-6}, + {0x1.11536cp+6, 0x1.645ae6p-19, -0x1.8b5ccap-3, -0x1.724d02p-8, + 0x1.8aaf6p-3, 0x1.71d33p-6}, + {0x1.1de456p+6, -0x1.6bc7a4p-19, 0x1.829356p-3, 0x1.5a280ep-8, + -0x1.81f85cp-3, -0x1.59bff8p-6}, + {0x1.2a754p+6, -0x1.5f7eep-20, -0x1.7a597ep-3, -0x1.4486cp-8, + 0x1.79ce5p-3, 0x1.442d38p-6}, + {0x1.37062cp+6, -0x1.ab28bap-20, 0x1.72a09ap-3, 0x1.310f06p-8, + -0x1.72230ep-3, -0x1.30c186p-6}, + {0x1.439718p+6, 0x1.c5c6f4p-19, -0x1.6b5c04p-3, -0x1.1f766p-8, + 0x1.6aea5p-3, 0x1.1f32e8p-6}, + {0x1.502808p+6, -0x1.2cbbcep-19, 0x1.6480c4p-3, 0x1.0f7eb8p-8, + -0x1.641964p-3, -0x1.0f43acp-6}, + {0x1.5cb8f8p+6, -0x1.f0c70ap-19, -0x1.5e0544p-3, -0x1.00f3f2p-8, + 0x1.5da6f4p-3, 0x1.00c004p-6}, + {0x1.6949e8p+6, -0x1.81383ep-20, 0x1.57e12p-3, 0x1.e75414p-9, + -0x1.578accp-3, -0x1.e6f852p-7}, + {0x1.75dadap+6, -0x1.cea80cp-19, -0x1.520ceep-3, -0x1.cef722p-9, + 0x1.51bdacp-3, 0x1.cea5bap-7}, + {0x1.826bccp+6, -0x1.46c93cp-19, 0x1.4c8222p-3, 0x1.b89114p-9, + -0x1.4c392ap-3, -0x1.b8489p-7}, + {0x1.8efcbep+6, 0x1.615096p-20, -0x1.473ae6p-3, -0x1.a3eaeap-9, + 0x1.46f78ap-3, 0x1.a3aa16p-7}, + }; + + float j, d, d2; + d = z - x[n][0] - x[n][1]; + d2 = d / 2; + j = - x[n][2] * _j1(d2) + x[n][3] * _j2(d2) - + x[n][4] * _j3(d2) + x[n][5] * _j4(d2); + return (j); +} + +/* + * The zeros of J0(x) fall near x = (n + 0.75) * pi. These are offsets + * that are used to localize a value of x into a small interval about a + * zero. + */ +static const float offset[32] = { + 4.86309e-02, 2.22910e-02, 1.43477e-02, 1.05619e-02, + 8.35263e-03, 6.90623e-03, 5.88708e-03, 5.12924e-03, + 4.54305e-03, 4.07894e-03, 3.70066e-03, 3.38531e-03, + 3.11957e-03, 2.89196e-03, 2.69488e-03, 2.52450e-03, + 2.37319e-03, 2.24095e-03, 2.12016e-03, 2.01463e-03, + 1.91673e-03, 1.82645e-03, 1.75144e-03, 1.67643e-03, + 1.60904e-03, 1.54166e-03, 1.48953e-03, 1.43740e-03, + 1.38528e-03, 1.34078e-03, 1.29628e-03, 1.25179e-03 +}; + + +float __ieee754_j0f(float x) { float z, s,c,ss,cc,r,u,v; - int32_t hx,ix; + int32_t hx,ix,k; GET_FLOAT_WORD(hx,x); ix = hx&0x7fffffff; if(ix>=0x7f800000) return one/(x*x); x = fabsf(x); if(ix >= 0x40000000) { /* |x| >= 2.0 */ + + if (x < 100.) { + k = (int)(x * invpi - 0.75); + if (fabsf(x - (k + 0.75) * pi - offset[k]) < 0.065) + return (j0zerof(k, x)); + } + s = sinf(x); c = cosf(x); ss = s-c; From owner-freebsd-standards@FreeBSD.ORG Tue Feb 23 00:30:04 2010 Return-Path: Delivered-To: freebsd-standards@hub.freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:4f8:fff6::34]) by hub.freebsd.org (Postfix) with ESMTP id 870E31065672 for ; Tue, 23 Feb 2010 00:30:04 +0000 (UTC) (envelope-from gnats@FreeBSD.org) Received: from freefall.freebsd.org (freefall.freebsd.org [IPv6:2001:4f8:fff6::28]) by mx1.freebsd.org (Postfix) with ESMTP id 59D7A8FC14 for ; Tue, 23 Feb 2010 00:30:04 +0000 (UTC) Received: from freefall.freebsd.org (localhost [127.0.0.1]) by freefall.freebsd.org (8.14.3/8.14.3) with ESMTP id o1N0U4B3030836 for ; Tue, 23 Feb 2010 00:30:04 GMT (envelope-from gnats@freefall.freebsd.org) Received: (from gnats@localhost) by freefall.freebsd.org (8.14.3/8.14.3/Submit) id o1N0U4Wk030831; Tue, 23 Feb 2010 00:30:04 GMT (envelope-from gnats) Date: Tue, 23 Feb 2010 00:30:04 GMT Message-Id: <201002230030.o1N0U4Wk030831@freefall.freebsd.org> To: freebsd-standards@FreeBSD.org From: "Steven G. Kargl" Cc: Subject: Re: standards/142803: j0 Bessel function inaccurate near zeros of the function X-BeenThere: freebsd-standards@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list Reply-To: "Steven G. Kargl" List-Id: Standards compliance List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Tue, 23 Feb 2010 00:30:04 -0000 The following reply was made to PR standards/142803; it has been noted by GNATS. From: "Steven G. Kargl" To: Bruce Evans Cc: FreeBSD-gnats-submit@FreeBSD.org, freebsd-standards@FreeBSD.org Subject: Re: standards/142803: j0 Bessel function inaccurate near zeros of the function Date: Mon, 22 Feb 2010 16:21:21 -0800 (PST) Bruce Evans wrote: > On Thu, 14 Jan 2010, Steven G. Kargl wrote: > > > Bruce Evans wrote: > >> On Wed, 13 Jan 2010, Steven G. Kargl wrote: > >> > >>>> Description: > >>> > >>> The j0 bessel function supplied by libm is fairly inaccurate at > >>> ... > >> > >> This is a very hard and relatively unimportant problem. > > > > Yes, it is very hard, but apparently you do not use bessel > > functions in your everyday life. :) > > > > I only discover this issue because I need bessel functions > > of complex arguments and I found my routines have issues > > in the vicinity of zeros. So, I decided to look at the > > libm routines. > > It is interesting that the often-poor accuracy of almost every system's > libm matters in real life. > > Complex args are another interesting problem since even complex > multiplication is hard to do accurately (it may have large cancelation > errors). > > >> Anyway, if you can get anywhere near < 10 ulp error near all zeros using > >> only an asymptotic method, then that would be good. Then the asymptotic > >> method would also be capable of locating the zeros very accurately. But > >> I would be very surprised if this worked. I know of nothing similar for > >> reducing mod Pi for trigonometric functions, which seems a simpler problem. > >> I would expect it to at best involve thousands of binary digits in the > >> tables for the asymptotic method, and corresponding thousands of digits > >> of precision in the calculation (4000 as for mfpr enough for the 2**100th > >> zero?). > > > > The 4000-bit setting for mpfr was a hold over from testing mpfr_j0 > > against my ascending series implementation of j0 with mpfr > > primitives. As few as 128-bits is sufficient to achieve the > > following: > > > >>> x my j0f(x) libm j0f(x) MPFR j0 my err libm err > > 2.404825 5.6434398E-08 5.9634296E-08 5.6434400E-08 0.05 152824.59 > > 5.520078 2.4476657E-08 2.4153294E-08 2.4476659E-08 0.10 18878.52 > > 8.653728 1.0355303E-07 1.0359805E-07 1.0355306E-07 0.86 1694.47 > > 11.791534 -3.5291243E-09 -3.5193941E-09 -3.5301714E-09 75.93 781.53 > > Wonder why this jumps. > > > 14.930918 -6.4815082E-09 -6.3911618E-09 -6.4815052E-09 0.23 6722.88 > > 18.071064 5.1532352E-09 5.3149818E-09 5.1532318E-09 0.23 10910.50 > > 21.211637 -1.5023349E-07 -1.5002509E-07 -1.5023348E-07 2.70 56347.01 > > ... > > > > As I suspected by adding additional terms to the asymptotic > > approximation and performing all computations with double > > precision, reduces 'my err' (5th column). The value at > > x=11.7... is the best I can get. The asymptotic approximations > > contain divergent series and additional terms do not help. > > The extra precision is almost certainly necessary. Whether double > precision is nearly enough is unclear, but the error near 11.7 suggests > that it is nearly enough except there. The large error might be caused > by that zero alone (among small zeros) being very close to a representable > value. > > I forgot the mention that the error table in my previous mail is on amd64 > for comparing float precision functions with double precision ones, assuming > that the latter are correct, which they aren't, but they are hopefully > correct enough for this comparision. The errors on i386 are much larger, > due to i386 still using i387 hardware trigonometric which are extremely > inaccurate near zeros, starting at the first zero. Here are both tables: > > amd64: > %%% > j0: max_er = 0x7fffffffffdf5e07 17179869183.9960, avg_er = 5.581, #>=1:0.5 = 1593961230:1839722934 > j1: max_er = 0x7fffffffffbcd2a1 17179869183.9918, avg_er = 4.524, #>=1:0.5 = 1597678928:1856295142 > lgamma:max_er = 0x135a0b77e00000 10145883.7461, avg_er = 0.252, #>=1:0.5 = 44084256:331444835 > y0: max_er = 0x7fffffffffbcd2a0 17179869183.9918, avg_er = 2.379, #>=1:0.5 = 837057577:1437331064 > y1: max_er = 0x7fffffffffdf5e07 17179869183.9960, avg_er = 3.761, #>=1:0.5 = 865063612:1460264955 > %%% > > i386: > %%% > j0: max_er = 0x6f8686c5c9f26 3654467.3863, avg_er = 0.418, #>=1:0.5 = 671562746:1332944948 > j1: max_er = 0x449026e9c286f 2246675.4566, avg_er = 0.425, #>=1:0.5 = 674510414:1347770568 > lgamma:max_er = 0xe4cf242400000 7497618.0703, avg_er = 0.274, #>=1:0.5 = 70033452:508702222 > y0: max_er = 0x93a2340c00000 4837658.0234, avg_er = 0.380, #>=1:0.5 = 594207097:1303826516 > y1: max_er = 0x7ffa2068256f72ab 17176789825.1699, avg_er = 5.188, #>=1:0.5 = 459137173:1213136103 > %%% > > Unfortunately, most of these i386 errors, and the amd64 error for y1() > are misreported. The huge max_er of 0x7ff... (16 hex digits) is > actually a sign error misreported. Sign errors are bad enough. They > always result at a simple zero z0 (f'(z0) != 0) when the approximation > is so inaccurate that it cannot tell on which side of infinite-precision > z0 the parameter is. Methods involving a table of zeros will not have > any of these provided the table is accurate to within 1 ulp, but other > methods can easily have them, depending on the other method's ability > to locate the zeros to within 1 ulp by solving f~(z) ~= 0 where f~ is > the approximate f. > > Having no sign errors across the whole range seems too good to believe > for the amd64 functions. All of the i387 hardware trig functions have > sign errors. > Bruce, David, Since the my real-life work requires additional accuracy for the Bessel function routines in libm, I spent sometime playing with e_j0f. I've come up with the patch that follows my signature. It only deals with the first 32 zeros of j0(x), others can be added to the table. I've yet to come up with a general method to improve the accuracy around all zeros. To help you gauge the improvement, one can look at http://troutmask.apl.washington.edu/~kargl/j0f.jpg http://troutmask.apl.washington.edu/~kargl/j0f_1.jpg which compares the unpatched j0f(x) to the patched j0f(x). The patch makes use of the addition theorem for J0(x+d) where x is a zero and d is some small deviation from the zero. Using only the first 4 terms is sufficent for j0f. J0(x+d) = - J1(x)*J1(d) + J2(x)*J2(d) - J3(x)*J3(d) + J4(x)*J4(d) where J1(d), J2(d), J3(d) and J4(d) can be approximated by truncating their series representations. The functions with x, ie the known zero, are evaluated with MPFR. Note, d is determined from d = z - x = z - high(x) - low(x) where high(x) is the 1st 24 bits of the zero, low(x) is the next 24 bits of zero, and z is some value in [x-0.065,x+0.065]. Bruce (and/or David), if you think that this merits inclusion in FreeBSD's libm, I'll work up a similar patch for e_j0.c and submit both; otherwise, I keep it as a local change. -- Steve http://troutmask.apl.washington.edu/~kargl/ Index: src/e_j0f.c =================================================================== --- src/e_j0f.c (revision 204219) +++ src/e_j0f.c (working copy) @@ -22,6 +22,8 @@ static float pzerof(float), qzerof(float); static const float +invpi = 0.31830988618379067154, +pi = 3.14159265358979323846, huge = 1e30, one = 1.0, invsqrtpi= 5.6418961287e-01, /* 0x3f106ebb */ @@ -38,17 +40,157 @@ static const float zero = 0.0; +/* J1(x) for x << 1. */ +static float +_j1(float x2) +{ + return (x2 * (1 - 0.5 * x2 * x2)); +} + +/* J2(x) for x << 1. */ +static float +_j2(float x2) +{ + x2 *= x2; + return (0.5 * x2 * (1 - x2 / 3)); +} + +/* J3(x) for x << 1. */ +static float +_j3(float x2) +{ + float y; + y = x2 * x2; + return (x2 * y * (1 - 0.25 * y) / 6); +} + +/* J4(x) for x << 1. */ +static float +_j4(float x2) +{ + x2 *= x2; + return (x2 * x2 / 24); +} + +/* + * Use the addition theorem for values near a zero the Bessel function. + * J0(z) = J0(x + d) where |d| << 1. + * = - 2*J1(x)*J1(d) + 2*J2(x)*J2(d) - ... + * Note, the factor of 2 is included in the table. + */ float +j0zerof(int n, float z) +{ + static const float x[32][6] = { + {0x1.33d152p+1, 0x1.d2e368p-24, 0x1.09cdb4p+0, 0x1.ba1deep-1, + 0x1.978d44p-2, 0x1.0933ccp-3}, + {0x1.6148f6p+2, -0x1.34f46ep-24, -0x1.5c6e6p-1, -0x1.f8f72ep-3, + 0x1.00f404p-1, 0x1.9588cep-1}, + {0x1.14eb56p+3, 0x1.999bdap-22, 0x1.15f798p-1, 0x1.00f7fcp-3, + -0x1.f08b9p-2, -0x1.d8c2a8p-2}, + {0x1.79544p+3, 0x1.04e56cp-26, -0x1.dc13e6p-2, -0x1.42ff0cp-4, + 0x1.c0af7ep-2, 0x1.350edcp-2}, + {0x1.ddca14p+3, -0x1.0d8e2ep-25, 0x1.a701dp-2, 0x1.c54b94p-5, + -0x1.97d3ccp-2, -0x1.b9186p-3}, + {0x1.212314p+4, -0x1.d7982ap-26, -0x1.8077f6p-2, -0x1.5467ecp-5, + 0x1.770cdp-2, 0x1.4e26dp-3}, + {0x1.5362dep+4, -0x1.d1810ep-21, 0x1.62d93ap-2, 0x1.0ba9cep-5, + -0x1.5c8a0ap-2, -0x1.08180cp-3}, + {0x1.85a3bap+4, -0x1.9fd524p-21, -0x1.4b2a2ep-2, -0x1.b3298p-6, + 0x1.46b28cp-2, 0x1.aec26ap-4}, + {0x1.b7e54ap+4, 0x1.7f57c4p-22, 0x1.37aac8p-2, 0x1.6ac0d2p-6, + -0x1.345e5cp-2, -0x1.67dfb2p-4}, + {0x1.ea275ap+4, -0x1.c68826p-21, -0x1.27407ep-2, -0x1.34695p-6, + 0x1.24bc2ep-2, 0x1.32708ap-4}, + {0x1.0e34e2p+5, -0x1.8b3204p-20, 0x1.192f24p-2, 0x1.0a6682p-6, + -0x1.17365ap-2, -0x1.08ffd2p-4}, + {0x1.275638p+5, -0x1.5a7986p-21, -0x1.0cf3eep-2, -0x1.d242aap-7, + 0x1.0b5fc4p-2, 0x1.d0352cp-5}, + {0x1.4077a8p+5, -0x1.29d6c6p-23, 0x1.0230bap-2, 0x1.9c8084p-7, + -0x1.00e734p-2, -0x1.9af5aap-5}, + {0x1.59992cp+5, 0x1.974364p-21, -0x1.f13fbp-3, -0x1.70558ep-7, + 0x1.ef1ep-3, 0x1.6f2666p-5}, + {0x1.72bacp+5, 0x1.f02102p-20, 0x1.e018dap-3, 0x1.4b858ap-7, + -0x1.de4fp-3, -0x1.4a986ap-5}, + {0x1.8bdc62p+5, 0x1.27e0cap-20, -0x1.d09b22p-3, -0x1.2c74f6p-7, + 0x1.cf1686p-3, 0x1.2bb87cp-5}, + {0x1.a4fe0ep+5, 0x1.c8899p-20, 0x1.c28612p-3, 0x1.11f526p-7, + -0x1.c138e4p-3, -0x1.115d32p-5}, + {0x1.be1fc4p+5, 0x1.a4c606p-23, -0x1.b5a622p-3, -0x1.f645fep-8, + 0x1.b485eap-3, 0x1.f54de8p-6}, + {0x1.d7418p+5, 0x1.93c83ep-20, 0x1.a9d184p-3, 0x1.cea254p-8, + -0x1.a8d632p-3, -0x1.cdd58ap-6}, + {0x1.f06344p+5, -0x1.7b4716p-22, -0x1.9ee5eep-3, -0x1.abf28ap-8, + 0x1.9e093ap-3, 0x1.ab47cep-6}, + {0x1.04c286p+6, 0x1.0f88f2p-21, 0x1.94c6f6p-3, 0x1.8d6372p-8, + -0x1.9403e4p-3, -0x1.8cd3dp-6}, + {0x1.11536cp+6, 0x1.645ae6p-19, -0x1.8b5ccap-3, -0x1.724d02p-8, + 0x1.8aaf6p-3, 0x1.71d33p-6}, + {0x1.1de456p+6, -0x1.6bc7a4p-19, 0x1.829356p-3, 0x1.5a280ep-8, + -0x1.81f85cp-3, -0x1.59bff8p-6}, + {0x1.2a754p+6, -0x1.5f7eep-20, -0x1.7a597ep-3, -0x1.4486cp-8, + 0x1.79ce5p-3, 0x1.442d38p-6}, + {0x1.37062cp+6, -0x1.ab28bap-20, 0x1.72a09ap-3, 0x1.310f06p-8, + -0x1.72230ep-3, -0x1.30c186p-6}, + {0x1.439718p+6, 0x1.c5c6f4p-19, -0x1.6b5c04p-3, -0x1.1f766p-8, + 0x1.6aea5p-3, 0x1.1f32e8p-6}, + {0x1.502808p+6, -0x1.2cbbcep-19, 0x1.6480c4p-3, 0x1.0f7eb8p-8, + -0x1.641964p-3, -0x1.0f43acp-6}, + {0x1.5cb8f8p+6, -0x1.f0c70ap-19, -0x1.5e0544p-3, -0x1.00f3f2p-8, + 0x1.5da6f4p-3, 0x1.00c004p-6}, + {0x1.6949e8p+6, -0x1.81383ep-20, 0x1.57e12p-3, 0x1.e75414p-9, + -0x1.578accp-3, -0x1.e6f852p-7}, + {0x1.75dadap+6, -0x1.cea80cp-19, -0x1.520ceep-3, -0x1.cef722p-9, + 0x1.51bdacp-3, 0x1.cea5bap-7}, + {0x1.826bccp+6, -0x1.46c93cp-19, 0x1.4c8222p-3, 0x1.b89114p-9, + -0x1.4c392ap-3, -0x1.b8489p-7}, + {0x1.8efcbep+6, 0x1.615096p-20, -0x1.473ae6p-3, -0x1.a3eaeap-9, + 0x1.46f78ap-3, 0x1.a3aa16p-7}, + }; + + float j, d, d2; + d = z - x[n][0] - x[n][1]; + d2 = d / 2; + j = - x[n][2] * _j1(d2) + x[n][3] * _j2(d2) - + x[n][4] * _j3(d2) + x[n][5] * _j4(d2); + return (j); +} + +/* + * The zeros of J0(x) fall near x = (n + 0.75) * pi. These are offsets + * that are used to localize a value of x into a small interval about a + * zero. + */ +static const float offset[32] = { + 4.86309e-02, 2.22910e-02, 1.43477e-02, 1.05619e-02, + 8.35263e-03, 6.90623e-03, 5.88708e-03, 5.12924e-03, + 4.54305e-03, 4.07894e-03, 3.70066e-03, 3.38531e-03, + 3.11957e-03, 2.89196e-03, 2.69488e-03, 2.52450e-03, + 2.37319e-03, 2.24095e-03, 2.12016e-03, 2.01463e-03, + 1.91673e-03, 1.82645e-03, 1.75144e-03, 1.67643e-03, + 1.60904e-03, 1.54166e-03, 1.48953e-03, 1.43740e-03, + 1.38528e-03, 1.34078e-03, 1.29628e-03, 1.25179e-03 +}; + + +float __ieee754_j0f(float x) { float z, s,c,ss,cc,r,u,v; - int32_t hx,ix; + int32_t hx,ix,k; GET_FLOAT_WORD(hx,x); ix = hx&0x7fffffff; if(ix>=0x7f800000) return one/(x*x); x = fabsf(x); if(ix >= 0x40000000) { /* |x| >= 2.0 */ + + if (x < 100.) { + k = (int)(x * invpi - 0.75); + if (fabsf(x - (k + 0.75) * pi - offset[k]) < 0.065) + return (j0zerof(k, x)); + } + s = sinf(x); c = cosf(x); ss = s-c; From owner-freebsd-standards@FreeBSD.ORG Tue Feb 23 08:06:30 2010 Return-Path: Delivered-To: freebsd-standards@freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:4f8:fff6::34]) by hub.freebsd.org (Postfix) with ESMTP id 445DD106566B for ; Tue, 23 Feb 2010 08:06:30 +0000 (UTC) (envelope-from jdc@koitsu.dyndns.org) Received: from qmta05.emeryville.ca.mail.comcast.net (qmta05.emeryville.ca.mail.comcast.net [76.96.30.48]) by mx1.freebsd.org (Postfix) with ESMTP id DC6598FC18 for ; Tue, 23 Feb 2010 08:06:29 +0000 (UTC) Received: from omta07.emeryville.ca.mail.comcast.net ([76.96.30.59]) by qmta05.emeryville.ca.mail.comcast.net with comcast id lL191d0021GXsucA5L6WNN; Tue, 23 Feb 2010 08:06:30 +0000 Received: from koitsu.dyndns.org ([98.248.46.159]) by omta07.emeryville.ca.mail.comcast.net with comcast id lL6V1d00U3S48mS8TL6VK6; Tue, 23 Feb 2010 08:06:30 +0000 Received: by icarus.home.lan (Postfix, from userid 1000) id 5B5D41E301A; Tue, 23 Feb 2010 00:06:28 -0800 (PST) Date: Tue, 23 Feb 2010 00:06:28 -0800 From: Jeremy Chadwick To: freebsd-stable@freebsd.org Message-ID: <20100223080628.GA22397@icarus.home.lan> References: <4B6228EF.5050400@laposte.net> <201001290802.o0T829sd050043@lurza.secnetix.de> <20100129090357.GA38872@icarus.home.lan> <4B71F0A3.6020401@laposte.net> MIME-Version: 1.0 Content-Type: text/plain; charset=iso-8859-1 Content-Disposition: inline Content-Transfer-Encoding: 8bit In-Reply-To: <4B71F0A3.6020401@laposte.net> User-Agent: Mutt/1.5.20 (2009-06-14) Cc: freebsd-standards@freebsd.org Subject: Re: su password prompt to stdout instead of /dev/tty X-BeenThere: freebsd-standards@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list List-Id: Standards compliance List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Tue, 23 Feb 2010 08:06:30 -0000 On Wed, Feb 10, 2010 at 12:32:51AM +0100, Cyrille Lefevre wrote: > Jeremy Chadwick a écrit : > > > >OpenPAM is des@'s responsibility. Has anyone brought this up to him? > > still no answer from des@ ! any idea ? I don't know what's going on. He's around/active[1], but this must not be high on his priority list. [1]: http://lists.freebsd.org/pipermail/freebsd-arch/2010-February/009906.html -- | Jeremy Chadwick jdc@parodius.com | | Parodius Networking http://www.parodius.com/ | | UNIX Systems Administrator Mountain View, CA, USA | | Making life hard for others since 1977. PGP: 4BD6C0CB | From owner-freebsd-standards@FreeBSD.ORG Tue Feb 23 14:30:02 2010 Return-Path: Delivered-To: freebsd-standards@hub.freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:4f8:fff6::34]) by hub.freebsd.org (Postfix) with ESMTP id 194C5106568D for ; Tue, 23 Feb 2010 14:30:02 +0000 (UTC) (envelope-from gnats@FreeBSD.org) Received: from freefall.freebsd.org (freefall.freebsd.org [IPv6:2001:4f8:fff6::28]) by mx1.freebsd.org (Postfix) with ESMTP id D00718FC12 for ; Tue, 23 Feb 2010 14:30:01 +0000 (UTC) Received: from freefall.freebsd.org (localhost [127.0.0.1]) by freefall.freebsd.org (8.14.3/8.14.3) with ESMTP id o1NEU1Kx093188 for ; Tue, 23 Feb 2010 14:30:01 GMT (envelope-from gnats@freefall.freebsd.org) Received: (from gnats@localhost) by freefall.freebsd.org (8.14.3/8.14.3/Submit) id o1NEU1sq093182; Tue, 23 Feb 2010 14:30:01 GMT (envelope-from gnats) Resent-Date: Tue, 23 Feb 2010 14:30:01 GMT Resent-Message-Id: <201002231430.o1NEU1sq093182@freefall.freebsd.org> Resent-From: FreeBSD-gnats-submit@FreeBSD.org (GNATS Filer) Resent-To: freebsd-standards@FreeBSD.org Resent-Reply-To: FreeBSD-gnats-submit@FreeBSD.org, Axel Dörfler Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:4f8:fff6::34]) by hub.freebsd.org (Postfix) with ESMTP id 3F22B106566C for ; Tue, 23 Feb 2010 14:25:04 +0000 (UTC) (envelope-from nobody@FreeBSD.org) Received: from www.freebsd.org (www.freebsd.org [IPv6:2001:4f8:fff6::21]) by mx1.freebsd.org (Postfix) with ESMTP id 2E1498FC12 for ; Tue, 23 Feb 2010 14:25:04 +0000 (UTC) Received: from www.freebsd.org (localhost [127.0.0.1]) by www.freebsd.org (8.14.3/8.14.3) with ESMTP id o1NEP3ab077939 for ; Tue, 23 Feb 2010 14:25:03 GMT (envelope-from nobody@www.freebsd.org) Received: (from nobody@localhost) by www.freebsd.org (8.14.3/8.14.3/Submit) id o1NEP3wQ077938; Tue, 23 Feb 2010 14:25:03 GMT (envelope-from nobody) Message-Id: <201002231425.o1NEP3wQ077938@www.freebsd.org> Date: Tue, 23 Feb 2010 14:25:03 GMT From: Axel Dörfler To: freebsd-gnats-submit@FreeBSD.org X-Send-Pr-Version: www-3.1 Cc: Subject: standards/144231: bind/connect/sendto too strict about sockaddr length X-BeenThere: freebsd-standards@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list List-Id: Standards compliance List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Tue, 23 Feb 2010 14:30:02 -0000 >Number: 144231 >Category: standards >Synopsis: bind/connect/sendto too strict about sockaddr length >Confidential: no >Severity: non-critical >Priority: low >Responsible: freebsd-standards >State: open >Quarter: >Keywords: >Date-Required: >Class: sw-bug >Submitter-Id: current-users >Arrival-Date: Tue Feb 23 14:30:01 UTC 2010 >Closed-Date: >Last-Modified: >Originator: Axel Dörfler >Release: 8.0-RELEASE >Organization: >Environment: i386 >Description: bind(), connect(), and sendto() all have a socklen_t argument specifying the size of the sockaddr buffer passed in. However, with IPv4 at least, passing any other value than 16 (sizeof(sockaddr)) will result in EINVAL. There is no reason why this restriction is there, and it's quite unhandy when dealing with the sockaddr_storage structure. (I have not compiled the test program, but it should be able to reproduce the problem) >How-To-Repeat: #include #include #include int main() { sockaddr_storage buffer; ((sockaddr_in*)&buffer)->sa_family = AF_INET; ((sockaddr_in*)&buffer)->sa_addr.s_addr = INADDR_ANY; ((sockaddr_in*)&buffer)->sa_port = 0; int fd = socket(AF_INET, SOCK_DGRAM, 0); if (bind(fd, &buffer, sizeof(buffer)) != 0) perror("bind"); return 0; } >Fix: >Release-Note: >Audit-Trail: >Unformatted: From owner-freebsd-standards@FreeBSD.ORG Thu Feb 25 11:09:13 2010 Return-Path: Delivered-To: freebsd-standards@freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:4f8:fff6::34]) by hub.freebsd.org (Postfix) with ESMTP id 8DA33106566B; Thu, 25 Feb 2010 11:09:13 +0000 (UTC) (envelope-from brde@optusnet.com.au) Received: from mail10.syd.optusnet.com.au (mail10.syd.optusnet.com.au [211.29.132.191]) by mx1.freebsd.org (Postfix) with ESMTP id 22C228FC0C; Thu, 25 Feb 2010 11:09:12 +0000 (UTC) Received: from besplex.bde.org (c122-106-163-215.carlnfd1.nsw.optusnet.com.au [122.106.163.215]) by mail10.syd.optusnet.com.au (8.13.1/8.13.1) with ESMTP id o1PB990J013704 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO); Thu, 25 Feb 2010 22:09:10 +1100 Date: Thu, 25 Feb 2010 22:09:09 +1100 (EST) From: Bruce Evans X-X-Sender: bde@besplex.bde.org To: "Steven G. Kargl" In-Reply-To: <201002230021.o1N0LLNb002353@troutmask.apl.washington.edu> Message-ID: <20100225214825.E1337@besplex.bde.org> References: <201002230021.o1N0LLNb002353@troutmask.apl.washington.edu> MIME-Version: 1.0 Content-Type: TEXT/PLAIN; charset=US-ASCII; format=flowed Cc: FreeBSD-gnats-submit@freebsd.org, freebsd-standards@freebsd.org Subject: Re: standards/142803: j0 Bessel function inaccurate near zeros of the function X-BeenThere: freebsd-standards@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list List-Id: Standards compliance List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Thu, 25 Feb 2010 11:09:13 -0000 On Mon, 22 Feb 2010, Steven G. Kargl wrote: > Since the my real-life work requires additional accuracy for > the Bessel function routines in libm, I spent sometime playing > with e_j0f. I've come up with the patch that follows my signature. > It only deals with the first 32 zeros of j0(x), others can be > added to the table. I've yet to come up with a general method > to improve the accuracy around all zeros. To help you gauge > the improvement, one can look at Interesting. I will reply more later... I just browsed Abromowitz (sp?) and Stegun and saw that it gives methods for finding all the zeros and says that there are complex zeros. It seems strong on formulae and weak on practical computational methods and algorithms (even for its time; of course the best methods now are different). > +float > __ieee754_j0f(float x) > { > float z, s,c,ss,cc,r,u,v; > - int32_t hx,ix; > + int32_t hx,ix,k; > > GET_FLOAT_WORD(hx,x); > ix = hx&0x7fffffff; > if(ix>=0x7f800000) return one/(x*x); > x = fabsf(x); > if(ix >= 0x40000000) { /* |x| >= 2.0 */ > + > + if (x < 100.) { > + k = (int)(x * invpi - 0.75); > + if (fabsf(x - (k + 0.75) * pi - offset[k]) < 0.065) These and many other constants should by float constants. Otherwise you get slowness and negatively useful extra accuracy (no useful extra accuracy, due to floats elsewhere, but more complicated error analysis). This reduction is similar to the one for the first few and/or first few million zeros of trig functions in e_rem_pio2*.c. > + return (j0zerof(k, x)); > + } > + > s = sinf(x); > c = cosf(x); > ss = s-c; > Bruce From owner-freebsd-standards@FreeBSD.ORG Thu Feb 25 11:10:02 2010 Return-Path: Delivered-To: freebsd-standards@hub.freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:4f8:fff6::34]) by hub.freebsd.org (Postfix) with ESMTP id 9DC53106564A for ; Thu, 25 Feb 2010 11:10:02 +0000 (UTC) (envelope-from gnats@FreeBSD.org) Received: from freefall.freebsd.org (freefall.freebsd.org [IPv6:2001:4f8:fff6::28]) by mx1.freebsd.org (Postfix) with ESMTP id 7293A8FC12 for ; Thu, 25 Feb 2010 11:10:02 +0000 (UTC) Received: from freefall.freebsd.org (localhost [127.0.0.1]) by freefall.freebsd.org (8.14.3/8.14.3) with ESMTP id o1PBA2BA055243 for ; Thu, 25 Feb 2010 11:10:02 GMT (envelope-from gnats@freefall.freebsd.org) Received: (from gnats@localhost) by freefall.freebsd.org (8.14.3/8.14.3/Submit) id o1PBA247055242; Thu, 25 Feb 2010 11:10:02 GMT (envelope-from gnats) Date: Thu, 25 Feb 2010 11:10:02 GMT Message-Id: <201002251110.o1PBA247055242@freefall.freebsd.org> To: freebsd-standards@FreeBSD.org From: Bruce Evans Cc: Subject: Re: standards/142803: j0 Bessel function inaccurate near zeros of the function X-BeenThere: freebsd-standards@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list Reply-To: Bruce Evans List-Id: Standards compliance List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Thu, 25 Feb 2010 11:10:02 -0000 The following reply was made to PR standards/142803; it has been noted by GNATS. From: Bruce Evans To: "Steven G. Kargl" Cc: Bruce Evans , FreeBSD-gnats-submit@freebsd.org, freebsd-standards@freebsd.org Subject: Re: standards/142803: j0 Bessel function inaccurate near zeros of the function Date: Thu, 25 Feb 2010 22:09:09 +1100 (EST) On Mon, 22 Feb 2010, Steven G. Kargl wrote: > Since the my real-life work requires additional accuracy for > the Bessel function routines in libm, I spent sometime playing > with e_j0f. I've come up with the patch that follows my signature. > It only deals with the first 32 zeros of j0(x), others can be > added to the table. I've yet to come up with a general method > to improve the accuracy around all zeros. To help you gauge > the improvement, one can look at Interesting. I will reply more later... I just browsed Abromowitz (sp?) and Stegun and saw that it gives methods for finding all the zeros and says that there are complex zeros. It seems strong on formulae and weak on practical computational methods and algorithms (even for its time; of course the best methods now are different). > +float > __ieee754_j0f(float x) > { > float z, s,c,ss,cc,r,u,v; > - int32_t hx,ix; > + int32_t hx,ix,k; > > GET_FLOAT_WORD(hx,x); > ix = hx&0x7fffffff; > if(ix>=0x7f800000) return one/(x*x); > x = fabsf(x); > if(ix >= 0x40000000) { /* |x| >= 2.0 */ > + > + if (x < 100.) { > + k = (int)(x * invpi - 0.75); > + if (fabsf(x - (k + 0.75) * pi - offset[k]) < 0.065) These and many other constants should by float constants. Otherwise you get slowness and negatively useful extra accuracy (no useful extra accuracy, due to floats elsewhere, but more complicated error analysis). This reduction is similar to the one for the first few and/or first few million zeros of trig functions in e_rem_pio2*.c. > + return (j0zerof(k, x)); > + } > + > s = sinf(x); > c = cosf(x); > ss = s-c; > Bruce