From owner-freebsd-numerics@freebsd.org Tue Sep 4 00:15:00 2018 Return-Path: Delivered-To: freebsd-numerics@mailman.ysv.freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2610:1c1:1:606c::19:1]) by mailman.ysv.freebsd.org (Postfix) with ESMTP id 482B4FFA82E for ; Tue, 4 Sep 2018 00:15:00 +0000 (UTC) (envelope-from sgk@troutmask.apl.washington.edu) Received: from troutmask.apl.washington.edu (troutmask.apl.washington.edu [128.95.76.21]) (using TLSv1.2 with cipher ECDHE-RSA-AES256-GCM-SHA384 (256/256 bits)) (Client CN "troutmask", Issuer "troutmask" (not verified)) by mx1.freebsd.org (Postfix) with ESMTPS id C52EA7BFBF for ; Tue, 4 Sep 2018 00:14:56 +0000 (UTC) (envelope-from sgk@troutmask.apl.washington.edu) Received: from troutmask.apl.washington.edu (localhost [127.0.0.1]) by troutmask.apl.washington.edu (8.15.2/8.15.2) with ESMTPS id w83NvOFm095479 (version=TLSv1.2 cipher=DHE-RSA-AES256-GCM-SHA384 bits=256 verify=NO) for ; Mon, 3 Sep 2018 16:57:24 -0700 (PDT) (envelope-from sgk@troutmask.apl.washington.edu) Received: (from sgk@localhost) by troutmask.apl.washington.edu (8.15.2/8.15.2/Submit) id w83NvO86095478 for freebsd-numerics@freebsd.org; Mon, 3 Sep 2018 16:57:24 -0700 (PDT) (envelope-from sgk) Date: Mon, 3 Sep 2018 16:57:24 -0700 From: Steve Kargl To: freebsd-numerics@freebsd.org Subject: j0 (and y0) in the range 2 <= x < (p/2)*log(2) Message-ID: <20180903235724.GA95333@troutmask.apl.washington.edu> Reply-To: sgk@troutmask.apl.washington.edu MIME-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline User-Agent: Mutt/1.10.1 (2018-07-13) X-BeenThere: freebsd-numerics@freebsd.org X-Mailman-Version: 2.1.27 Precedence: list List-Id: "Discussions of high quality implementation of libm functions." List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Tue, 04 Sep 2018 00:15:00 -0000 Anyone know where the approximations for j0 (and y0) come from? msun/src/e_j0.c states * for x in (2,inf) * j0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)-q0(x)*sin(x0)) * where x0 = x-pi/4. It is better to compute sin(x0),cos(x0) * as follow: * cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4) * = 1/sqrt(2) * (cos(x) + sin(x)) * sin(x0) = sin(x)cos(pi/4)-cos(x)sin(pi/4) * = 1/sqrt(2) * (sin(x) - cos(x)) * (To avoid cancellation, use * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) * to compute the worse one.) p0(x) and q0(x) are divergent asymptotic series. If I extract pzero() and qzero() from e_j0.c and compare the results against summing truncated versions of p0(x) and q0(x), there are no obvious connections. Reading the documentation for the algorithms used in MPFR suggests that x >= p/2*log(2), where p is precision of x, is required for use of the large argument approximation for j0(x). In double precision, p = 53, so we have x >= 18.368... Consider x=18.4 and sum up to N = 31 in the asymptotic series: % ./pq 30 18.4 p = 9.997932830701132e-01, q = -6.781826311540553e-03 <-- series pp = 9.997932830701132e-01, qq = -6.781826311540509e-03 <-- pzero,qzero ulp(p, pp) = 0.000000e+00 ulp(q, qq) = 2.550000e+01 This is almost reasonable if 25.5 ULP is acceptable in q0(x). Note the series are computed in long double with 64 bits of precision. Now, comparing x = 2 and summing N = 4 (best results). % ./pq 4 2 p = 9.894313812255859e-01, q = -5.334472656250000e-02 pp = 9.862158212188928e-01, qq = -5.647769967932505e-02 ulp(p, pp) = 1.448159e+13 ulp(q, qq) = 2.257545e+14 For values of N > 4, the series start to diverge! So, how does msun use the large argument approximation for j0(x)? -- Steve