From owner-freebsd-standards@FreeBSD.ORG  Mon Jan 11 11:07:09 2010
Return-Path: <owner-freebsd-standards@FreeBSD.ORG>
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 D41FE1065698
	for <freebsd-standards@FreeBSD.org>;
	Mon, 11 Jan 2010 11:07:09 +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 C36C48FC30
	for <freebsd-standards@FreeBSD.org>;
	Mon, 11 Jan 2010 11:07:09 +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 o0BB79Jj034798
	for <freebsd-standards@FreeBSD.org>; Mon, 11 Jan 2010 11:07:09 GMT
	(envelope-from owner-bugmaster@FreeBSD.org)
Received: (from gnats@localhost)
	by freefall.freebsd.org (8.14.3/8.14.3/Submit) id o0BB79Uf034796
	for freebsd-standards@FreeBSD.org; Mon, 11 Jan 2010 11:07:09 GMT
	(envelope-from owner-bugmaster@FreeBSD.org)
Date: Mon, 11 Jan 2010 11:07:09 GMT
Message-Id: <201001111107.o0BB79Uf034796@freefall.freebsd.org>
X-Authentication-Warning: freefall.freebsd.org: gnats set sender to
	owner-bugmaster@FreeBSD.org using -f
From: FreeBSD bugmaster <bugmaster@FreeBSD.org>
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 <freebsd-standards.freebsd.org>
List-Unsubscribe: <http://lists.freebsd.org/mailman/listinfo/freebsd-standards>, 
	<mailto:freebsd-standards-request@freebsd.org?subject=unsubscribe>
List-Archive: <http://lists.freebsd.org/pipermail/freebsd-standards>
List-Post: <mailto:freebsd-standards@freebsd.org>
List-Help: <mailto:freebsd-standards-request@freebsd.org?subject=help>
List-Subscribe: <http://lists.freebsd.org/mailman/listinfo/freebsd-standards>, 
	<mailto:freebsd-standards-request@freebsd.org?subject=subscribe>
X-List-Received-Date: Mon, 11 Jan 2010 11:07:10 -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/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 <time.h>
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

49 problems total.


From owner-freebsd-standards@FreeBSD.ORG  Thu Jan 14 01:50:02 2010
Return-Path: <owner-freebsd-standards@FreeBSD.ORG>
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 0C7F6106566B
	for <freebsd-standards@hub.freebsd.org>;
	Thu, 14 Jan 2010 01:50: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 BF8168FC16
	for <freebsd-standards@hub.freebsd.org>;
	Thu, 14 Jan 2010 01:50: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 o0E1o1qR080601
	for <freebsd-standards@freefall.freebsd.org>;
	Thu, 14 Jan 2010 01:50:01 GMT
	(envelope-from gnats@freefall.freebsd.org)
Received: (from gnats@localhost)
	by freefall.freebsd.org (8.14.3/8.14.3/Submit) id o0E1o1Z6080600;
	Thu, 14 Jan 2010 01:50:01 GMT (envelope-from gnats)
Resent-Date: Thu, 14 Jan 2010 01:50:01 GMT
Resent-Message-Id: <201001140150.o0E1o1Z6080600@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,
	"Steven G. Kargl" <kargl@troutmask.apl.washington.edu>
Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:4f8:fff6::34])
	by hub.freebsd.org (Postfix) with ESMTP id 1D1F3106568B
	for <FreeBSD-gnats-submit@freebsd.org>;
	Thu, 14 Jan 2010 01:43:18 +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 D507B8FC08
	for <FreeBSD-gnats-submit@freebsd.org>;
	Thu, 14 Jan 2010 01:43:17 +0000 (UTC)
Received: from troutmask.apl.washington.edu (localhost.apl.washington.edu
	[127.0.0.1])
	by troutmask.apl.washington.edu (8.14.3/8.14.3) with ESMTP id
	o0E1hOlt040084 for <FreeBSD-gnats-submit@freebsd.org>;
	Wed, 13 Jan 2010 17:43:24 -0800 (PST)
	(envelope-from kargl@troutmask.apl.washington.edu)
Received: (from kargl@localhost)
	by troutmask.apl.washington.edu (8.14.3/8.14.3/Submit) id
	o0E1hOm8040083; Wed, 13 Jan 2010 17:43:24 -0800 (PST)
	(envelope-from kargl)
Message-Id: <201001140143.o0E1hOm8040083@troutmask.apl.washington.edu>
Date: Wed, 13 Jan 2010 17:43:24 -0800 (PST)
From: "Steven G. Kargl" <kargl@troutmask.apl.washington.edu>
To: FreeBSD-gnats-submit@FreeBSD.org
X-Send-Pr-Version: 3.113
Cc: 
Subject: 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" <kargl@troutmask.apl.washington.edu>
List-Id: Standards compliance <freebsd-standards.freebsd.org>
List-Unsubscribe: <http://lists.freebsd.org/mailman/listinfo/freebsd-standards>, 
	<mailto:freebsd-standards-request@freebsd.org?subject=unsubscribe>
List-Archive: <http://lists.freebsd.org/pipermail/freebsd-standards>
List-Post: <mailto:freebsd-standards@freebsd.org>
List-Help: <mailto:freebsd-standards-request@freebsd.org?subject=help>
List-Subscribe: <http://lists.freebsd.org/mailman/listinfo/freebsd-standards>, 
	<mailto:freebsd-standards-request@freebsd.org?subject=subscribe>
X-List-Received-Date: Thu, 14 Jan 2010 01:50:02 -0000


>Number:         142803
>Category:       standards
>Synopsis:       j0 Bessel function inaccurate near zeros of the function
>Confidential:   no
>Severity:       serious
>Priority:       low
>Responsible:    freebsd-standards
>State:          open
>Quarter:        
>Keywords:       
>Date-Required:
>Class:          sw-bug
>Submitter-Id:   current-users
>Arrival-Date:   Thu Jan 14 01:50:01 UTC 2010
>Closed-Date:
>Last-Modified:
>Originator:     Steven G. Kargl
>Release:        FreeBSD 9.0-CURRENT amd64
>Organization:
APL-UW
>Environment:
System: FreeBSD troutmask.apl.washington.edu 9.0-CURRENT FreeBSD 9.0-CURRENT #3 r201269M: Mon Jan 4 13:52:03 PST 2010 kargl@troutmask.apl.washington.edu:/usr/obj/usr/src/sys/SPEW amd64


	
>Description:

The j0 bessel function supplied by libm is fairly inaccurate at
arguments at and near a zero of the function.  Here's the first
30 zeros computed by j0f, my implementation of j0f, a 4000-bit
significand computed via MPFR (the assumed exact value), and the
relative absolute error. 

    x        my j0f(x)     libm j0f(x)    MPFR j0        my err  libm err
 2.404825  5.6434434E-08  5.9634296E-08  5.6434400E-08      1.64 152824.59
 5.520078  2.4476664E-08  2.4153294E-08  2.4476659E-08      0.31  18878.52
 8.653728  1.0355323E-07  1.0359805E-07  1.0355306E-07      6.36   1694.47
11.791534 -2.4511966E-09 -3.5193941E-09 -3.5301714E-09  78243.14    781.53
14.930918 -6.6019510E-09 -6.3911618E-09 -6.4815052E-09   8962.94   6722.88
18.071064  5.1734359E-09  5.3149818E-09  5.1532318E-09   1362.83  10910.50
21.211637 -1.5023797E-07 -1.5002509E-07 -1.5023348E-07   1213.07  56347.01
24.352472  1.2524691E-07  1.2516310E-07  1.2524570E-07     41.65   2834.53
27.493479  5.4330716E-08  5.4263626E-08  5.4331104E-08     18.77   3261.75
30.634607  1.2205560E-07  1.2203689E-07  1.2205546E-07      4.85    645.39
33.775822 -2.0213101E-07 -2.0206903E-07 -2.0213095E-07      5.48   6263.95
36.917099  8.4751605E-08  8.4749573E-08  8.4751581E-08      0.99     82.59
40.058426 -1.7484851E-08 -1.7475532E-08 -1.7484840E-08      0.91    767.56
43.199791 -9.2091398E-08 -9.2135146E-08 -9.2091406E-08      2.47  13530.51
46.341187  2.1663259E-07  2.1664336E-07  2.1663259E-07      0.16    268.90
49.482609 -1.2502527E-07 -1.2504157E-07 -1.2502526E-07      2.69  23512.60
52.624050  1.8706569E-07  1.8707487E-07  1.8706569E-07      0.01    251.43
55.765511 -2.0935554E-08 -2.0932896E-08 -2.0935556E-08      0.19    227.03
58.906982  1.5637660E-07  1.5634730E-07  1.5637661E-07      0.26    892.22
62.048470  3.5779891E-08  3.5787338E-08  3.5779888E-08      0.14    403.17

Note, my j0f(x) currently uses double precision to accumulate intermediate
values.  Below x = 10, I use the ascending series to compute the value.
Above x = 10, I'm using an asymptotic approximation.  I haven't investigated
whether additional terms in an asymptotic approximation would pull 'my err'
for x = 11, 14, 18, and 21 closer to the exact value.  

>How-To-Repeat:

>Fix:

I'll get to it eventually.


>Release-Note:
>Audit-Trail:
>Unformatted:

From owner-freebsd-standards@FreeBSD.ORG  Thu Jan 14 02:59:20 2010
Return-Path: <owner-freebsd-standards@FreeBSD.ORG>
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 8A8C6106566B;
	Thu, 14 Jan 2010 02:59:20 +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 52BF98FC0C;
	Thu, 14 Jan 2010 02:59:20 +0000 (UTC)
Received: from troutmask.apl.washington.edu (localhost.apl.washington.edu
	[127.0.0.1])
	by troutmask.apl.washington.edu (8.14.3/8.14.3) with ESMTP id
	o0E2xQs8040432; Wed, 13 Jan 2010 18:59:26 -0800 (PST)
	(envelope-from kargl@troutmask.apl.washington.edu)
Received: (from kargl@localhost)
	by troutmask.apl.washington.edu (8.14.3/8.14.3/Submit) id
	o0E2xQaU040431; Wed, 13 Jan 2010 18:59:26 -0800 (PST)
	(envelope-from kargl)
From: "Steven G. Kargl" <kargl@troutmask.apl.washington.edu>
Message-Id: <201001140259.o0E2xQaU040431@troutmask.apl.washington.edu>
In-Reply-To: <201001140150.o0E1o1Zg080595@freefall.freebsd.org>
To: FreeBSD-gnats-submit@FreeBSD.org, freebsd-standards@FreeBSD.org
Date: Wed, 13 Jan 2010 18:59:26 -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: 
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 <freebsd-standards.freebsd.org>
List-Unsubscribe: <http://lists.freebsd.org/mailman/listinfo/freebsd-standards>, 
	<mailto:freebsd-standards-request@freebsd.org?subject=unsubscribe>
List-Archive: <http://lists.freebsd.org/pipermail/freebsd-standards>
List-Post: <mailto:freebsd-standards@freebsd.org>
List-Help: <mailto:freebsd-standards-request@freebsd.org?subject=help>
List-Subscribe: <http://lists.freebsd.org/mailman/listinfo/freebsd-standards>, 
	<mailto:freebsd-standards-request@freebsd.org?subject=subscribe>
X-List-Received-Date: Thu, 14 Jan 2010 02:59:20 -0000

I forgot to mention that I believe this
issue will be found in j1(3) as well.
I haven't tried any tests with y0, y1, 
jn or yn.  I'd be might surprised if 
these were robust.

-- 
Steve
http://troutmask.apl.washington.edu/~kargl/

From owner-freebsd-standards@FreeBSD.ORG  Thu Jan 14 03:00:11 2010
Return-Path: <owner-freebsd-standards@FreeBSD.ORG>
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 AC1821065695
	for <freebsd-standards@hub.freebsd.org>;
	Thu, 14 Jan 2010 03:00:11 +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 582BC8FC17
	for <freebsd-standards@hub.freebsd.org>;
	Thu, 14 Jan 2010 03:00:11 +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 o0E30BJv040984
	for <freebsd-standards@freefall.freebsd.org>;
	Thu, 14 Jan 2010 03:00:11 GMT
	(envelope-from gnats@freefall.freebsd.org)
Received: (from gnats@localhost)
	by freefall.freebsd.org (8.14.3/8.14.3/Submit) id o0E30B8V040981;
	Thu, 14 Jan 2010 03:00:11 GMT (envelope-from gnats)
Date: Thu, 14 Jan 2010 03:00:11 GMT
Message-Id: <201001140300.o0E30B8V040981@freefall.freebsd.org>
To: freebsd-standards@FreeBSD.org
From: "Steven G. Kargl" <kargl@troutmask.apl.washington.edu>
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" <kargl@troutmask.apl.washington.edu>
List-Id: Standards compliance <freebsd-standards.freebsd.org>
List-Unsubscribe: <http://lists.freebsd.org/mailman/listinfo/freebsd-standards>, 
	<mailto:freebsd-standards-request@freebsd.org?subject=unsubscribe>
List-Archive: <http://lists.freebsd.org/pipermail/freebsd-standards>
List-Post: <mailto:freebsd-standards@freebsd.org>
List-Help: <mailto:freebsd-standards-request@freebsd.org?subject=help>
List-Subscribe: <http://lists.freebsd.org/mailman/listinfo/freebsd-standards>, 
	<mailto:freebsd-standards-request@freebsd.org?subject=subscribe>
X-List-Received-Date: Thu, 14 Jan 2010 03:00:11 -0000

The following reply was made to PR standards/142803; it has been noted by GNATS.

From: "Steven G. Kargl" <kargl@troutmask.apl.washington.edu>
To: FreeBSD-gnats-submit@FreeBSD.org, freebsd-standards@FreeBSD.org
Cc:  
Subject: Re: standards/142803: j0 Bessel function inaccurate near
 zeros of the function
Date: Wed, 13 Jan 2010 18:59:26 -0800 (PST)

 I forgot to mention that I believe this
 issue will be found in j1(3) as well.
 I haven't tried any tests with y0, y1, 
 jn or yn.  I'd be might surprised if 
 these were robust.
 
 -- 
 Steve
 http://troutmask.apl.washington.edu/~kargl/

From owner-freebsd-standards@FreeBSD.ORG  Thu Jan 14 03:18:00 2010
Return-Path: <owner-freebsd-standards@FreeBSD.ORG>
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 65A6D1065676;
	Thu, 14 Jan 2010 03:18:00 +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 2E76C8FC14;
	Thu, 14 Jan 2010 03:18:00 +0000 (UTC)
Received: from troutmask.apl.washington.edu (localhost.apl.washington.edu
	[127.0.0.1])
	by troutmask.apl.washington.edu (8.14.3/8.14.3) with ESMTP id
	o0E3I6in040548; Wed, 13 Jan 2010 19:18:06 -0800 (PST)
	(envelope-from kargl@troutmask.apl.washington.edu)
Received: (from kargl@localhost)
	by troutmask.apl.washington.edu (8.14.3/8.14.3/Submit) id
	o0E3I689040547; Wed, 13 Jan 2010 19:18:06 -0800 (PST)
	(envelope-from kargl)
From: "Steven G. Kargl" <kargl@troutmask.apl.washington.edu>
Message-Id: <201001140318.o0E3I689040547@troutmask.apl.washington.edu>
In-Reply-To: <201001140150.o0E1o1Zg080595@freefall.freebsd.org>
To: FreeBSD-gnats-submit@FreeBSD.org, freebsd-standards@FreeBSD.org
Date: Wed, 13 Jan 2010 19:18:06 -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: 
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 <freebsd-standards.freebsd.org>
List-Unsubscribe: <http://lists.freebsd.org/mailman/listinfo/freebsd-standards>, 
	<mailto:freebsd-standards-request@freebsd.org?subject=unsubscribe>
List-Archive: <http://lists.freebsd.org/pipermail/freebsd-standards>
List-Post: <mailto:freebsd-standards@freebsd.org>
List-Help: <mailto:freebsd-standards-request@freebsd.org?subject=help>
List-Subscribe: <http://lists.freebsd.org/mailman/listinfo/freebsd-standards>, 
	<mailto:freebsd-standards-request@freebsd.org?subject=subscribe>
X-List-Received-Date: Thu, 14 Jan 2010 03:18:00 -0000

Well, I just checked the bug report on the
web, and the nicely formatted table looks 
like garbage due to the stupidity of crushing
consecutive spaces down to a single space.

-- 
Steve
http://troutmask.apl.washington.edu/~kargl/

From owner-freebsd-standards@FreeBSD.ORG  Thu Jan 14 03:20:08 2010
Return-Path: <owner-freebsd-standards@FreeBSD.ORG>
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 97A81106566C
	for <freebsd-standards@hub.freebsd.org>;
	Thu, 14 Jan 2010 03:20:08 +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 6BADC8FC15
	for <freebsd-standards@hub.freebsd.org>;
	Thu, 14 Jan 2010 03:20:08 +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 o0E3K8sx059169
	for <freebsd-standards@freefall.freebsd.org>;
	Thu, 14 Jan 2010 03:20:08 GMT
	(envelope-from gnats@freefall.freebsd.org)
Received: (from gnats@localhost)
	by freefall.freebsd.org (8.14.3/8.14.3/Submit) id o0E3K8IV059166;
	Thu, 14 Jan 2010 03:20:08 GMT (envelope-from gnats)
Date: Thu, 14 Jan 2010 03:20:08 GMT
Message-Id: <201001140320.o0E3K8IV059166@freefall.freebsd.org>
To: freebsd-standards@FreeBSD.org
From: "Steven G. Kargl" <kargl@troutmask.apl.washington.edu>
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" <kargl@troutmask.apl.washington.edu>
List-Id: Standards compliance <freebsd-standards.freebsd.org>
List-Unsubscribe: <http://lists.freebsd.org/mailman/listinfo/freebsd-standards>, 
	<mailto:freebsd-standards-request@freebsd.org?subject=unsubscribe>
List-Archive: <http://lists.freebsd.org/pipermail/freebsd-standards>
List-Post: <mailto:freebsd-standards@freebsd.org>
List-Help: <mailto:freebsd-standards-request@freebsd.org?subject=help>
List-Subscribe: <http://lists.freebsd.org/mailman/listinfo/freebsd-standards>, 
	<mailto:freebsd-standards-request@freebsd.org?subject=subscribe>
X-List-Received-Date: Thu, 14 Jan 2010 03:20:08 -0000

The following reply was made to PR standards/142803; it has been noted by GNATS.

From: "Steven G. Kargl" <kargl@troutmask.apl.washington.edu>
To: FreeBSD-gnats-submit@FreeBSD.org, freebsd-standards@FreeBSD.org
Cc:  
Subject: Re: standards/142803: j0 Bessel function inaccurate near
 zeros of the function
Date: Wed, 13 Jan 2010 19:18:06 -0800 (PST)

 Well, I just checked the bug report on the
 web, and the nicely formatted table looks 
 like garbage due to the stupidity of crushing
 consecutive spaces down to a single space.
 
 -- 
 Steve
 http://troutmask.apl.washington.edu/~kargl/

From owner-freebsd-standards@FreeBSD.ORG  Thu Jan 14 10:26:57 2010
Return-Path: <owner-freebsd-standards@FreeBSD.ORG>
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 19BBF106566C;
	Thu, 14 Jan 2010 10:26:57 +0000 (UTC)
	(envelope-from brde@optusnet.com.au)
Received: from mail05.syd.optusnet.com.au (mail05.syd.optusnet.com.au
	[211.29.132.186])
	by mx1.freebsd.org (Postfix) with ESMTP id A47098FC13;
	Thu, 14 Jan 2010 10:26:56 +0000 (UTC)
Received: from c122-106-161-214.carlnfd1.nsw.optusnet.com.au
	(c122-106-161-214.carlnfd1.nsw.optusnet.com.au [122.106.161.214])
	by mail05.syd.optusnet.com.au (8.13.1/8.13.1) with ESMTP id
	o0EAQnMt017970
	(version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO);
	Thu, 14 Jan 2010 21:26:53 +1100
Date: Thu, 14 Jan 2010 21:26:49 +1100 (EST)
From: Bruce Evans <brde@optusnet.com.au>
X-X-Sender: bde@delplex.bde.org
To: "Steven G. Kargl" <kargl@troutmask.apl.washington.edu>
In-Reply-To: <201001140143.o0E1hOm8040083@troutmask.apl.washington.edu>
Message-ID: <20100114200057.C62558@delplex.bde.org>
References: <201001140143.o0E1hOm8040083@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 <freebsd-standards.freebsd.org>
List-Unsubscribe: <http://lists.freebsd.org/mailman/listinfo/freebsd-standards>, 
	<mailto:freebsd-standards-request@freebsd.org?subject=unsubscribe>
List-Archive: <http://lists.freebsd.org/pipermail/freebsd-standards>
List-Post: <mailto:freebsd-standards@freebsd.org>
List-Help: <mailto:freebsd-standards-request@freebsd.org?subject=help>
List-Subscribe: <http://lists.freebsd.org/mailman/listinfo/freebsd-standards>, 
	<mailto:freebsd-standards-request@freebsd.org?subject=subscribe>
X-List-Received-Date: Thu, 14 Jan 2010 10:26:57 -0000

On Wed, 13 Jan 2010, Steven G. Kargl wrote:

>> Description:
>
> The j0 bessel function supplied by libm is fairly inaccurate at
> arguments at and near a zero of the function.  Here's the first
> 30 zeros computed by j0f, my implementation of j0f, a 4000-bit
> significand computed via MPFR (the assumed exact value), and the
> relative absolute error.

This is a very hard and relatively unimportant problem.

The corresponding problem for trigonometric functions is fairly hard
and relatively very important.  It is solved in fdlibm (and this in
FreeBSD libm) essentially by using a table of all the relevant zeros,
with the necessary several thousands of binary digits of accuracy for
many of the zeros.  Since trigonometric functions are periodic, it is
not necessary to have a table of approximations (e.g., polynomials)
for each of the zeros of interest.  There are O(2**100) zeros of
interest for 128-bit long doubles, so simple tables wouldn't work even
for the zeros.

The corresponding problem for lgamma() is relatively simple and
relatively important.  It is solved in places like the Intel ia64 libm
by using a table of all the relevant zeros and a table of approximations
for each of the zeros.  There are only about 70 relevant zeros (since
zeroes less than about -70 are so close to the (pseudo-)poles that
they cannot be represented in double precision (unlike the poles which,
since they are at the negative integers, can be represented down to
about -10**53 in double precision)).

For the corresponding problem for at least the simplest bessel function
j0(), the zeros are distributed like the zeros of sin(), but they
aren't exactly at multiples of Pi like those of sin(), so just finding
all the relevant ones to the necessary thousands of binary digits of
accuracy is a hard problem.  At least j0() is bounded like sin() (and
unlike lgamma()); thus we know that there are O(2**100) relevant zeros,
so it is impossible to find them all in advance.  Since bessel functions
aren't periodic, it is also unclear how they could be calculated
accurately near the zeros without using a different special approximation
near every zero.  In general, such calculations need to be done in extra
precision even for the linear term, giving relatively minor extra
complications.

I once tried to caclulate lgamma(z) near just the first zero of lgamma(),
using only extra precision applied to standard formalas (Stirling
approximation and not reflection since extra precision was only easy to
obtain than for the former), to see how far I could get using only
extra precision.  It took approximately doubled double precision (double
precision with most intermediate calculations done in extra precision)
to get near float precision for lgamma().

It might be worth making bessel functions accurate near just the first
few zeros.

>    x        my j0f(x)     libm j0f(x)    MPFR j0        my err  libm err
> 2.404825  5.6434434E-08  5.9634296E-08  5.6434400E-08      1.64 152824.59
> 5.520078  2.4476664E-08  2.4153294E-08  2.4476659E-08      0.31  18878.52
> 8.653728  1.0355323E-07  1.0359805E-07  1.0355306E-07      6.36   1694.47
> 11.791534 -2.4511966E-09 -3.5193941E-09 -3.5301714E-09  78243.14    781.53

Hmm.

> ...
> 62.048470  3.5779891E-08  3.5787338E-08  3.5779888E-08      0.14    403.17

The largest errors across all 2**32 float values for the 1-parameter
float precision bessel functions in libm is about 3.7 Mulps for j0()
and about 17 Gulps for y1():

%%%
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
%%%


> Note, my j0f(x) currently uses double precision to accumulate intermediate
> values.  Below x = 10, I use the ascending series to compute the value.
> Above x = 10, I'm using an asymptotic approximation.  I haven't investigated
> whether additional terms in an asymptotic approximation would pull 'my err'
> for x = 11, 14, 18, and 21 closer to the exact value.

Stirling asymptotic is what barely worked (with significant extra
precision) for the first zero of lgamma().  lgamma() is especially bad
for higher zeros of lgamma(), since there is a pole very near the zero.
I guess the behaviour for j0() is not so bad since there are no poles,
but the asymptotic method also seems to work not so well (compared with
a power series method) on lgamma(x) for x > 1 where there are no poles.

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?).

Bruce

From owner-freebsd-standards@FreeBSD.ORG  Thu Jan 14 10:30:05 2010
Return-Path: <owner-freebsd-standards@FreeBSD.ORG>
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 360C21065670
	for <freebsd-standards@hub.freebsd.org>;
	Thu, 14 Jan 2010 10:30:05 +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 244BF8FC0C
	for <freebsd-standards@hub.freebsd.org>;
	Thu, 14 Jan 2010 10:30:05 +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 o0EAU4hj066637
	for <freebsd-standards@freefall.freebsd.org>;
	Thu, 14 Jan 2010 10: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 o0EAU4ON066631;
	Thu, 14 Jan 2010 10:30:04 GMT (envelope-from gnats)
Date: Thu, 14 Jan 2010 10:30:04 GMT
Message-Id: <201001141030.o0EAU4ON066631@freefall.freebsd.org>
To: freebsd-standards@FreeBSD.org
From: Bruce Evans <brde@optusnet.com.au>
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 <brde@optusnet.com.au>
List-Id: Standards compliance <freebsd-standards.freebsd.org>
List-Unsubscribe: <http://lists.freebsd.org/mailman/listinfo/freebsd-standards>, 
	<mailto:freebsd-standards-request@freebsd.org?subject=unsubscribe>
List-Archive: <http://lists.freebsd.org/pipermail/freebsd-standards>
List-Post: <mailto:freebsd-standards@freebsd.org>
List-Help: <mailto:freebsd-standards-request@freebsd.org?subject=help>
List-Subscribe: <http://lists.freebsd.org/mailman/listinfo/freebsd-standards>, 
	<mailto:freebsd-standards-request@freebsd.org?subject=subscribe>
X-List-Received-Date: Thu, 14 Jan 2010 10:30:05 -0000

The following reply was made to PR standards/142803; it has been noted by GNATS.

From: Bruce Evans <brde@optusnet.com.au>
To: "Steven G. Kargl" <kargl@troutmask.apl.washington.edu>
Cc: FreeBSD-gnats-submit@FreeBSD.org, freebsd-standards@FreeBSD.org
Subject: Re: standards/142803: j0 Bessel function inaccurate near zeros of
 the function
Date: Thu, 14 Jan 2010 21:26:49 +1100 (EST)

 On Wed, 13 Jan 2010, Steven G. Kargl wrote:
 
 >> Description:
 >
 > The j0 bessel function supplied by libm is fairly inaccurate at
 > arguments at and near a zero of the function.  Here's the first
 > 30 zeros computed by j0f, my implementation of j0f, a 4000-bit
 > significand computed via MPFR (the assumed exact value), and the
 > relative absolute error.
 
 This is a very hard and relatively unimportant problem.
 
 The corresponding problem for trigonometric functions is fairly hard
 and relatively very important.  It is solved in fdlibm (and this in
 FreeBSD libm) essentially by using a table of all the relevant zeros,
 with the necessary several thousands of binary digits of accuracy for
 many of the zeros.  Since trigonometric functions are periodic, it is
 not necessary to have a table of approximations (e.g., polynomials)
 for each of the zeros of interest.  There are O(2**100) zeros of
 interest for 128-bit long doubles, so simple tables wouldn't work even
 for the zeros.
 
 The corresponding problem for lgamma() is relatively simple and
 relatively important.  It is solved in places like the Intel ia64 libm
 by using a table of all the relevant zeros and a table of approximations
 for each of the zeros.  There are only about 70 relevant zeros (since
 zeroes less than about -70 are so close to the (pseudo-)poles that
 they cannot be represented in double precision (unlike the poles which,
 since they are at the negative integers, can be represented down to
 about -10**53 in double precision)).
 
 For the corresponding problem for at least the simplest bessel function
 j0(), the zeros are distributed like the zeros of sin(), but they
 aren't exactly at multiples of Pi like those of sin(), so just finding
 all the relevant ones to the necessary thousands of binary digits of
 accuracy is a hard problem.  At least j0() is bounded like sin() (and
 unlike lgamma()); thus we know that there are O(2**100) relevant zeros,
 so it is impossible to find them all in advance.  Since bessel functions
 aren't periodic, it is also unclear how they could be calculated
 accurately near the zeros without using a different special approximation
 near every zero.  In general, such calculations need to be done in extra
 precision even for the linear term, giving relatively minor extra
 complications.
 
 I once tried to caclulate lgamma(z) near just the first zero of lgamma(),
 using only extra precision applied to standard formalas (Stirling
 approximation and not reflection since extra precision was only easy to
 obtain than for the former), to see how far I could get using only
 extra precision.  It took approximately doubled double precision (double
 precision with most intermediate calculations done in extra precision)
 to get near float precision for lgamma().
 
 It might be worth making bessel functions accurate near just the first
 few zeros.
 
 >    x        my j0f(x)     libm j0f(x)    MPFR j0        my err  libm err
 > 2.404825  5.6434434E-08  5.9634296E-08  5.6434400E-08      1.64 152824.59
 > 5.520078  2.4476664E-08  2.4153294E-08  2.4476659E-08      0.31  18878.52
 > 8.653728  1.0355323E-07  1.0359805E-07  1.0355306E-07      6.36   1694.47
 > 11.791534 -2.4511966E-09 -3.5193941E-09 -3.5301714E-09  78243.14    781.53
 
 Hmm.
 
 > ...
 > 62.048470  3.5779891E-08  3.5787338E-08  3.5779888E-08      0.14    403.17
 
 The largest errors across all 2**32 float values for the 1-parameter
 float precision bessel functions in libm is about 3.7 Mulps for j0()
 and about 17 Gulps for y1():
 
 %%%
 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
 %%%
 
 
 > Note, my j0f(x) currently uses double precision to accumulate intermediate
 > values.  Below x = 10, I use the ascending series to compute the value.
 > Above x = 10, I'm using an asymptotic approximation.  I haven't investigated
 > whether additional terms in an asymptotic approximation would pull 'my err'
 > for x = 11, 14, 18, and 21 closer to the exact value.
 
 Stirling asymptotic is what barely worked (with significant extra
 precision) for the first zero of lgamma().  lgamma() is especially bad
 for higher zeros of lgamma(), since there is a pole very near the zero.
 I guess the behaviour for j0() is not so bad since there are no poles,
 but the asymptotic method also seems to work not so well (compared with
 a power series method) on lgamma(x) for x > 1 where there are no poles.
 
 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?).
 
 Bruce

From owner-freebsd-standards@FreeBSD.ORG  Thu Jan 14 22:49:39 2010
Return-Path: <owner-freebsd-standards@FreeBSD.ORG>
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 CA601106566B;
	Thu, 14 Jan 2010 22:49:39 +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 ABC7C8FC16;
	Thu, 14 Jan 2010 22:49:39 +0000 (UTC)
Received: from troutmask.apl.washington.edu (localhost.apl.washington.edu
	[127.0.0.1])
	by troutmask.apl.washington.edu (8.14.3/8.14.3) with ESMTP id
	o0EMnk1t047429; Thu, 14 Jan 2010 14:49:46 -0800 (PST)
	(envelope-from kargl@troutmask.apl.washington.edu)
Received: (from kargl@localhost)
	by troutmask.apl.washington.edu (8.14.3/8.14.3/Submit) id
	o0EMnjLS047428; Thu, 14 Jan 2010 14:49:45 -0800 (PST)
	(envelope-from kargl)
From: "Steven G. Kargl" <kargl@troutmask.apl.washington.edu>
Message-Id: <201001142249.o0EMnjLS047428@troutmask.apl.washington.edu>
In-Reply-To: <20100114200057.C62558@delplex.bde.org>
To: Bruce Evans <brde@optusnet.com.au>
Date: Thu, 14 Jan 2010 14:49:45 -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 <freebsd-standards.freebsd.org>
List-Unsubscribe: <http://lists.freebsd.org/mailman/listinfo/freebsd-standards>, 
	<mailto:freebsd-standards-request@freebsd.org?subject=unsubscribe>
List-Archive: <http://lists.freebsd.org/pipermail/freebsd-standards>
List-Post: <mailto:freebsd-standards@freebsd.org>
List-Help: <mailto:freebsd-standards-request@freebsd.org?subject=help>
List-Subscribe: <http://lists.freebsd.org/mailman/listinfo/freebsd-standards>, 
	<mailto:freebsd-standards-request@freebsd.org?subject=subscribe>
X-List-Received-Date: Thu, 14 Jan 2010 22:49:39 -0000

Bruce Evans wrote:
> On Wed, 13 Jan 2010, Steven G. Kargl wrote:
> 
> >> Description:
> >
> > The j0 bessel function supplied by libm is fairly inaccurate at
> > arguments at and near a zero of the function.  Here's the first
> > 30 zeros computed by j0f, my implementation of j0f, a 4000-bit
> > significand computed via MPFR (the assumed exact value), and the
> > relative absolute error.
> 
> 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.

> >    x        my j0f(x)     libm j0f(x)    MPFR j0        my err  libm err
> > 2.404825  5.6434434E-08  5.9634296E-08  5.6434400E-08      1.64 152824.59
> > 5.520078  2.4476664E-08  2.4153294E-08  2.4476659E-08      0.31  18878.52
> > 8.653728  1.0355323E-07  1.0359805E-07  1.0355306E-07      6.36   1694.47
> > 11.791534 -2.4511966E-09 -3.5193941E-09 -3.5301714E-09  78243.14    781.53
> 
> Hmm.

I forgot to mention that 'my err' and 'libm err' are
in units of epsilon (ie, FLT_EPSILON for j0f).


> > Note, my j0f(x) currently uses double precision to accumulate intermediate
> > values.  Below x = 10, I use the ascending series to compute the value.
> > Above x = 10, I'm using an asymptotic approximation.  I haven't investigated
> > whether additional terms in an asymptotic approximation would pull 'my err'
> > for x = 11, 14, 18, and 21 closer to the exact value.
> 

(snip)

> 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:

    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
   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
   24.352472  1.2524569E-07  1.2516310E-07  1.2524570E-07      0.28   2834.53
   27.493479  5.4331110E-08  5.4263626E-08  5.4331104E-08      0.29   3261.75
   30.634607  1.2205545E-07  1.2203689E-07  1.2205546E-07      0.09    645.39
   33.775822 -2.0213095E-07 -2.0206903E-07 -2.0213095E-07      0.27   6263.95
   36.917099  8.4751576E-08  8.4749573E-08  8.4751581E-08      0.18     82.59
   40.058426 -1.7484838E-08 -1.7475532E-08 -1.7484840E-08      0.12    767.56
   43.199791 -9.2091398E-08 -9.2135146E-08 -9.2091406E-08      2.47  13530.51
   46.341187  2.1663259E-07  2.1664336E-07  2.1663259E-07      0.16    268.90
   49.482609 -1.2502527E-07 -1.2504157E-07 -1.2502526E-07      2.69  23512.60
   52.624050  1.8706569E-07  1.8707487E-07  1.8706569E-07      0.01    251.43
   55.765511 -2.0935557E-08 -2.0932896E-08 -2.0935556E-08      0.10    227.04
   58.906982  1.5637660E-07  1.5634730E-07  1.5637661E-07      0.28    892.23
   62.048470  3.5779891E-08  3.5787338E-08  3.5779899E-08      0.42    402.61

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.

-- 
Steve
http://troutmask.apl.washington.edu/~kargl/

From owner-freebsd-standards@FreeBSD.ORG  Thu Jan 14 22:50:03 2010
Return-Path: <owner-freebsd-standards@FreeBSD.ORG>
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 3C2411065670
	for <freebsd-standards@hub.freebsd.org>;
	Thu, 14 Jan 2010 22:50:03 +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 101CB8FC0A
	for <freebsd-standards@hub.freebsd.org>;
	Thu, 14 Jan 2010 22:50:03 +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 o0EMo21q011723
	for <freebsd-standards@freefall.freebsd.org>;
	Thu, 14 Jan 2010 22:50:02 GMT
	(envelope-from gnats@freefall.freebsd.org)
Received: (from gnats@localhost)
	by freefall.freebsd.org (8.14.3/8.14.3/Submit) id o0EMo2Pq011721;
	Thu, 14 Jan 2010 22:50:02 GMT (envelope-from gnats)
Date: Thu, 14 Jan 2010 22:50:02 GMT
Message-Id: <201001142250.o0EMo2Pq011721@freefall.freebsd.org>
To: freebsd-standards@FreeBSD.org
From: "Steven G. Kargl" <kargl@troutmask.apl.washington.edu>
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" <kargl@troutmask.apl.washington.edu>
List-Id: Standards compliance <freebsd-standards.freebsd.org>
List-Unsubscribe: <http://lists.freebsd.org/mailman/listinfo/freebsd-standards>, 
	<mailto:freebsd-standards-request@freebsd.org?subject=unsubscribe>
List-Archive: <http://lists.freebsd.org/pipermail/freebsd-standards>
List-Post: <mailto:freebsd-standards@freebsd.org>
List-Help: <mailto:freebsd-standards-request@freebsd.org?subject=help>
List-Subscribe: <http://lists.freebsd.org/mailman/listinfo/freebsd-standards>, 
	<mailto:freebsd-standards-request@freebsd.org?subject=subscribe>
X-List-Received-Date: Thu, 14 Jan 2010 22:50:03 -0000

The following reply was made to PR standards/142803; it has been noted by GNATS.

From: "Steven G. Kargl" <kargl@troutmask.apl.washington.edu>
To: Bruce Evans <brde@optusnet.com.au>
Cc: FreeBSD-gnats-submit@FreeBSD.org, freebsd-standards@FreeBSD.org
Subject: Re: standards/142803: j0 Bessel function inaccurate near
 zeros of the function
Date: Thu, 14 Jan 2010 14:49:45 -0800 (PST)

 Bruce Evans wrote:
 > On Wed, 13 Jan 2010, Steven G. Kargl wrote:
 > 
 > >> Description:
 > >
 > > The j0 bessel function supplied by libm is fairly inaccurate at
 > > arguments at and near a zero of the function.  Here's the first
 > > 30 zeros computed by j0f, my implementation of j0f, a 4000-bit
 > > significand computed via MPFR (the assumed exact value), and the
 > > relative absolute error.
 > 
 > 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.
 
 > >    x        my j0f(x)     libm j0f(x)    MPFR j0        my err  libm err
 > > 2.404825  5.6434434E-08  5.9634296E-08  5.6434400E-08      1.64 152824.59
 > > 5.520078  2.4476664E-08  2.4153294E-08  2.4476659E-08      0.31  18878.52
 > > 8.653728  1.0355323E-07  1.0359805E-07  1.0355306E-07      6.36   1694.47
 > > 11.791534 -2.4511966E-09 -3.5193941E-09 -3.5301714E-09  78243.14    781.53
 > 
 > Hmm.
 
 I forgot to mention that 'my err' and 'libm err' are
 in units of epsilon (ie, FLT_EPSILON for j0f).
 
 
 > > Note, my j0f(x) currently uses double precision to accumulate intermediate
 > > values.  Below x = 10, I use the ascending series to compute the value.
 > > Above x = 10, I'm using an asymptotic approximation.  I haven't investigated
 > > whether additional terms in an asymptotic approximation would pull 'my err'
 > > for x = 11, 14, 18, and 21 closer to the exact value.
 > 
 
 (snip)
 
 > 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:
 
     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
    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
    24.352472  1.2524569E-07  1.2516310E-07  1.2524570E-07      0.28   2834.53
    27.493479  5.4331110E-08  5.4263626E-08  5.4331104E-08      0.29   3261.75
    30.634607  1.2205545E-07  1.2203689E-07  1.2205546E-07      0.09    645.39
    33.775822 -2.0213095E-07 -2.0206903E-07 -2.0213095E-07      0.27   6263.95
    36.917099  8.4751576E-08  8.4749573E-08  8.4751581E-08      0.18     82.59
    40.058426 -1.7484838E-08 -1.7475532E-08 -1.7484840E-08      0.12    767.56
    43.199791 -9.2091398E-08 -9.2135146E-08 -9.2091406E-08      2.47  13530.51
    46.341187  2.1663259E-07  2.1664336E-07  2.1663259E-07      0.16    268.90
    49.482609 -1.2502527E-07 -1.2504157E-07 -1.2502526E-07      2.69  23512.60
    52.624050  1.8706569E-07  1.8707487E-07  1.8706569E-07      0.01    251.43
    55.765511 -2.0935557E-08 -2.0932896E-08 -2.0935556E-08      0.10    227.04
    58.906982  1.5637660E-07  1.5634730E-07  1.5637661E-07      0.28    892.23
    62.048470  3.5779891E-08  3.5787338E-08  3.5779899E-08      0.42    402.61
 
 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.
 
 -- 
 Steve
 http://troutmask.apl.washington.edu/~kargl/

From owner-freebsd-standards@FreeBSD.ORG  Fri Jan 15 00:38:46 2010
Return-Path: <owner-freebsd-standards@FreeBSD.ORG>
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 0BED1106566B;
	Fri, 15 Jan 2010 00:38:46 +0000 (UTC)
	(envelope-from brde@optusnet.com.au)
Received: from mail04.syd.optusnet.com.au (mail04.syd.optusnet.com.au
	[211.29.132.185])
	by mx1.freebsd.org (Postfix) with ESMTP id 90A028FC16;
	Fri, 15 Jan 2010 00:38:45 +0000 (UTC)
Received: from c220-239-227-214.carlnfd1.nsw.optusnet.com.au
	(c220-239-227-214.carlnfd1.nsw.optusnet.com.au [220.239.227.214])
	by mail04.syd.optusnet.com.au (8.13.1/8.13.1) with ESMTP id
	o0F0cdGw006782
	(version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO);
	Fri, 15 Jan 2010 11:38:42 +1100
Date: Fri, 15 Jan 2010 11:38:39 +1100 (EST)
From: Bruce Evans <brde@optusnet.com.au>
X-X-Sender: bde@delplex.bde.org
To: "Steven G. Kargl" <kargl@troutmask.apl.washington.edu>
In-Reply-To: <201001142249.o0EMnjLS047428@troutmask.apl.washington.edu>
Message-ID: <20100115110443.Y63344@delplex.bde.org>
References: <201001142249.o0EMnjLS047428@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 <freebsd-standards.freebsd.org>
List-Unsubscribe: <http://lists.freebsd.org/mailman/listinfo/freebsd-standards>, 
	<mailto:freebsd-standards-request@freebsd.org?subject=unsubscribe>
List-Archive: <http://lists.freebsd.org/pipermail/freebsd-standards>
List-Post: <mailto:freebsd-standards@freebsd.org>
List-Help: <mailto:freebsd-standards-request@freebsd.org?subject=help>
List-Subscribe: <http://lists.freebsd.org/mailman/listinfo/freebsd-standards>, 
	<mailto:freebsd-standards-request@freebsd.org?subject=subscribe>
X-List-Received-Date: Fri, 15 Jan 2010 00:38:46 -0000

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

From owner-freebsd-standards@FreeBSD.ORG  Fri Jan 15 00:40:03 2010
Return-Path: <owner-freebsd-standards@FreeBSD.ORG>
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 B63CE106566B
	for <freebsd-standards@hub.freebsd.org>;
	Fri, 15 Jan 2010 00:40:03 +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 A39E78FC16
	for <freebsd-standards@hub.freebsd.org>;
	Fri, 15 Jan 2010 00:40:03 +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 o0F0e3KN006896
	for <freebsd-standards@freefall.freebsd.org>;
	Fri, 15 Jan 2010 00:40:03 GMT
	(envelope-from gnats@freefall.freebsd.org)
Received: (from gnats@localhost)
	by freefall.freebsd.org (8.14.3/8.14.3/Submit) id o0F0e3JF006895;
	Fri, 15 Jan 2010 00:40:03 GMT (envelope-from gnats)
Date: Fri, 15 Jan 2010 00:40:03 GMT
Message-Id: <201001150040.o0F0e3JF006895@freefall.freebsd.org>
To: freebsd-standards@FreeBSD.org
From: Bruce Evans <brde@optusnet.com.au>
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 <brde@optusnet.com.au>
List-Id: Standards compliance <freebsd-standards.freebsd.org>
List-Unsubscribe: <http://lists.freebsd.org/mailman/listinfo/freebsd-standards>, 
	<mailto:freebsd-standards-request@freebsd.org?subject=unsubscribe>
List-Archive: <http://lists.freebsd.org/pipermail/freebsd-standards>
List-Post: <mailto:freebsd-standards@freebsd.org>
List-Help: <mailto:freebsd-standards-request@freebsd.org?subject=help>
List-Subscribe: <http://lists.freebsd.org/mailman/listinfo/freebsd-standards>, 
	<mailto:freebsd-standards-request@freebsd.org?subject=subscribe>
X-List-Received-Date: Fri, 15 Jan 2010 00:40:03 -0000

The following reply was made to PR standards/142803; it has been noted by GNATS.

From: Bruce Evans <brde@optusnet.com.au>
To: "Steven G. Kargl" <kargl@troutmask.apl.washington.edu>
Cc: Bruce Evans <brde@optusnet.com.au>, FreeBSD-gnats-submit@FreeBSD.org,
        freebsd-standards@FreeBSD.org
Subject: Re: standards/142803: j0 Bessel function inaccurate near zeros of
 the function
Date: Fri, 15 Jan 2010 11:38:39 +1100 (EST)

 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

From owner-freebsd-standards@FreeBSD.ORG  Fri Jan 15 01:53:11 2010
Return-Path: <owner-freebsd-standards@FreeBSD.ORG>
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 C2A9A106568D;
	Fri, 15 Jan 2010 01:53:11 +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 97AC98FC14;
	Fri, 15 Jan 2010 01:53:11 +0000 (UTC)
Received: from troutmask.apl.washington.edu (localhost.apl.washington.edu
	[127.0.0.1])
	by troutmask.apl.washington.edu (8.14.3/8.14.3) with ESMTP id
	o0F1rHXb048196; Thu, 14 Jan 2010 17:53:17 -0800 (PST)
	(envelope-from kargl@troutmask.apl.washington.edu)
Received: (from kargl@localhost)
	by troutmask.apl.washington.edu (8.14.3/8.14.3/Submit) id
	o0F1rHu7048195; Thu, 14 Jan 2010 17:53:17 -0800 (PST)
	(envelope-from kargl)
From: "Steven G. Kargl" <kargl@troutmask.apl.washington.edu>
Message-Id: <201001150153.o0F1rHu7048195@troutmask.apl.washington.edu>
In-Reply-To: <20100115110443.Y63344@delplex.bde.org>
To: Bruce Evans <brde@optusnet.com.au>
Date: Thu, 14 Jan 2010 17:53:17 -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 <freebsd-standards.freebsd.org>
List-Unsubscribe: <http://lists.freebsd.org/mailman/listinfo/freebsd-standards>, 
	<mailto:freebsd-standards-request@freebsd.org?subject=unsubscribe>
List-Archive: <http://lists.freebsd.org/pipermail/freebsd-standards>
List-Post: <mailto:freebsd-standards@freebsd.org>
List-Help: <mailto:freebsd-standards-request@freebsd.org?subject=help>
List-Subscribe: <http://lists.freebsd.org/mailman/listinfo/freebsd-standards>, 
	<mailto:freebsd-standards-request@freebsd.org?subject=subscribe>
X-List-Received-Date: Fri, 15 Jan 2010 01:53:11 -0000

Bruce Evans wrote:
> On Thu, 14 Jan 2010, Steven G. Kargl wrote:
> > Bruce Evans wrote:
> 
> >> 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.

Below x=10, I use the ascending series.  Above x=0, I use an asymptotic
approximation (AA).   x = 11.79... is the first zero I hit with the AA.

> >   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.

The AA is j0(x) = (P(x) * cos(x+pi/4) + Q(x) * sin(x+pi/4)) where 
P(x) and Q(x) are infinite, divergent sums.  I use the first 11 
terms in P(x) and Q(x) to achieve the 'my err' = 75.9.  Additional
terms cause an increase in 'my err' as does fewer terms.  This is
probably the limit of double precision.

I haven't investigated the intervals around the values I listed.  So,
there may be larger errors that are yet to be found.

BTW, the MPFR website has a document that describes all their algorithms.
They claim that the AA can be used for |x| > p*log(2)/2 where p is
the precision of x; however, in the mpfr code the criterion is |x| > p/2.

http://www.mpfr.org/algo.html

> 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:

My values are also computed on amd64.  I seldomly use i386 for numerical
work.  A quick change to my test code to use the double precision j0()
suggests that it has sufficient precision for the comparison you've
done.

-- 
Steve
http://troutmask.apl.washington.edu/~kargl/

From owner-freebsd-standards@FreeBSD.ORG  Fri Jan 15 02:00:04 2010
Return-Path: <owner-freebsd-standards@FreeBSD.ORG>
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 2D55C106566C
	for <freebsd-standards@hub.freebsd.org>;
	Fri, 15 Jan 2010 02:00: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 ED7A28FC16
	for <freebsd-standards@hub.freebsd.org>;
	Fri, 15 Jan 2010 02:00:03 +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 o0F203fn075062
	for <freebsd-standards@freefall.freebsd.org>;
	Fri, 15 Jan 2010 02:00:03 GMT
	(envelope-from gnats@freefall.freebsd.org)
Received: (from gnats@localhost)
	by freefall.freebsd.org (8.14.3/8.14.3/Submit) id o0F203QL075061;
	Fri, 15 Jan 2010 02:00:03 GMT (envelope-from gnats)
Date: Fri, 15 Jan 2010 02:00:03 GMT
Message-Id: <201001150200.o0F203QL075061@freefall.freebsd.org>
To: freebsd-standards@FreeBSD.org
From: "Steven G. Kargl" <kargl@troutmask.apl.washington.edu>
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" <kargl@troutmask.apl.washington.edu>
List-Id: Standards compliance <freebsd-standards.freebsd.org>
List-Unsubscribe: <http://lists.freebsd.org/mailman/listinfo/freebsd-standards>, 
	<mailto:freebsd-standards-request@freebsd.org?subject=unsubscribe>
List-Archive: <http://lists.freebsd.org/pipermail/freebsd-standards>
List-Post: <mailto:freebsd-standards@freebsd.org>
List-Help: <mailto:freebsd-standards-request@freebsd.org?subject=help>
List-Subscribe: <http://lists.freebsd.org/mailman/listinfo/freebsd-standards>, 
	<mailto:freebsd-standards-request@freebsd.org?subject=subscribe>
X-List-Received-Date: Fri, 15 Jan 2010 02:00:04 -0000

The following reply was made to PR standards/142803; it has been noted by GNATS.

From: "Steven G. Kargl" <kargl@troutmask.apl.washington.edu>
To: Bruce Evans <brde@optusnet.com.au>
Cc: FreeBSD-gnats-submit@FreeBSD.org, freebsd-standards@FreeBSD.org
Subject: Re: standards/142803: j0 Bessel function inaccurate near
 zeros of the function
Date: Thu, 14 Jan 2010 17:53:17 -0800 (PST)

 Bruce Evans wrote:
 > On Thu, 14 Jan 2010, Steven G. Kargl wrote:
 > > Bruce Evans wrote:
 > 
 > >> 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.
 
 Below x=10, I use the ascending series.  Above x=0, I use an asymptotic
 approximation (AA).   x = 11.79... is the first zero I hit with the AA.
 
 > >   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.
 
 The AA is j0(x) = (P(x) * cos(x+pi/4) + Q(x) * sin(x+pi/4)) where 
 P(x) and Q(x) are infinite, divergent sums.  I use the first 11 
 terms in P(x) and Q(x) to achieve the 'my err' = 75.9.  Additional
 terms cause an increase in 'my err' as does fewer terms.  This is
 probably the limit of double precision.
 
 I haven't investigated the intervals around the values I listed.  So,
 there may be larger errors that are yet to be found.
 
 BTW, the MPFR website has a document that describes all their algorithms.
 They claim that the AA can be used for |x| > p*log(2)/2 where p is
 the precision of x; however, in the mpfr code the criterion is |x| > p/2.
 
 http://www.mpfr.org/algo.html
 
 > 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:
 
 My values are also computed on amd64.  I seldomly use i386 for numerical
 work.  A quick change to my test code to use the double precision j0()
 suggests that it has sufficient precision for the comparison you've
 done.
 
 -- 
 Steve
 http://troutmask.apl.washington.edu/~kargl/