From owner-freebsd-numerics@freebsd.org Tue Sep 4 03:57:48 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 541C8FDE9C3 for ; Tue, 4 Sep 2018 03:57:48 +0000 (UTC) (envelope-from stephen@missouri.edu) Received: from um-tip2-missouri-out.um.umsystem.edu (um-tip2-missouri-out.um.umsystem.edu [198.209.49.149]) (using TLSv1.2 with cipher ECDHE-RSA-AES256-GCM-SHA384 (256/256 bits)) (Client CN "um-tip1.um.umsystem.edu", Issuer "InCommon RSA Server CA" (verified OK)) by mx1.freebsd.org (Postfix) with ESMTPS id CCD6784BC8 for ; Tue, 4 Sep 2018 03:57:47 +0000 (UTC) (envelope-from stephen@missouri.edu) X-IronPort-Anti-Spam-Filtered: true X-IronPort-Anti-Spam-Result: =?us-ascii?q?A2E+BgCVAY5b/yI40cYXAUIOCAUBAQEBA?= =?us-ascii?q?wEBAQkBAQGCVwFMgQ9tEigKjFmNXZgnAQoYAQoJAoEZgl9GAoNpOBQBAgEBAgE?= =?us-ascii?q?BAgICaRwMgmhLagEBAQEBASMCDQdcAQEBAgIBASshIBsCAQgRBAEBLycBCQEdC?= =?us-ascii?q?AIEAQcHBAEHFQSDAYEdZA+GLJ4Gh3MBB4F+BYpvggCBEoIUfoMbAQICgT4BAYV?= =?us-ascii?q?2Aod6hSWONwkChjKFTINxH45ZiyeIFwICAgIJAhSBWCI0gSFyEzuCbIM2AQiHV?= =?us-ascii?q?oJmgh46bwEBDoEGigoECxeBCAGBGwEB?= X-IPAS-Result: =?us-ascii?q?A2E+BgCVAY5b/yI40cYXAUIOCAUBAQEBAwEBAQkBAQGCVwF?= =?us-ascii?q?MgQ9tEigKjFmNXZgnAQoYAQoJAoEZgl9GAoNpOBQBAgEBAgEBAgICaRwMgmhLa?= =?us-ascii?q?gEBAQEBASMCDQdcAQEBAgIBASshIBsCAQgRBAEBLycBCQEdCAIEAQcHBAEHFQS?= =?us-ascii?q?DAYEdZA+GLJ4Gh3MBB4F+BYpvggCBEoIUfoMbAQICgT4BAYV2Aod6hSWONwkCh?= =?us-ascii?q?jKFTINxH45ZiyeIFwICAgIJAhSBWCI0gSFyEzuCbIM2AQiHVoJmgh46bwEBDoE?= =?us-ascii?q?GigoECxeBCAGBGwEB?= Received: from ex-n22.um.umsystem.edu (HELO EX2-N22.um.umsystem.edu) ([198.209.56.34]) by um-tip2-exch-relay.um.umsystem.edu with ESMTP; 03 Sep 2018 22:56:29 -0500 Received: from EX2-N14.um.umsystem.edu (198.209.56.22) by EX2-N22.um.umsystem.edu (198.209.56.34) with Microsoft SMTP Server (version=TLS1_2, cipher=TLS_ECDHE_RSA_WITH_AES_128_CBC_SHA256_P256) id 15.1.1531.3; Mon, 3 Sep 2018 22:56:28 -0500 Received: from EX2-N14.um.umsystem.edu ([198.209.56.22]) by EX2-N14.um.umsystem.edu ([198.209.56.22]) with mapi id 15.01.1531.003; Mon, 3 Sep 2018 22:56:28 -0500 From: "Montgomery-Smith, Stephen" To: "freebsd-numerics@freebsd.org" , "sgk@troutmask.apl.washington.edu" Subject: Re: j0 (and y0) in the range 2 <= x < (p/2)*log(2) Thread-Topic: j0 (and y0) in the range 2 <= x < (p/2)*log(2) Thread-Index: AQHUQ+RdF989zhJjQEG1rlb5E3D8e6TffLMg Date: Tue, 4 Sep 2018 03:56:28 +0000 Message-ID: References: <20180903235724.GA95333@troutmask.apl.washington.edu> In-Reply-To: <20180903235724.GA95333@troutmask.apl.washington.edu> Accept-Language: en-US Content-Language: en-US X-MS-Has-Attach: X-MS-TNEF-Correlator: x-originating-ip: [72.161.255.145] MIME-Version: 1.0 Content-Type: text/plain; charset="us-ascii" Content-Transfer-Encoding: quoted-printable X-Content-Filtered-By: Mailman/MimeDel 2.1.27 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 03:57:48 -0000 A quick google search turned up this https://www.cl.cam.ac.uk/~jrh13/papers/bessel.pdf which has the functions p0 and q0. Maybe this was the basis of this code. ________________________________ From: owner-freebsd-numerics@freebsd.org on behalf of Steve Kargl Sent: Monday, September 3, 2018 6:57:24 PM To: freebsd-numerics@freebsd.org Subject: j0 (and y0) in the range 2 <=3D x < (p/2)*log(2) Anyone know where the approximations for j0 (and y0) come from? msun/src/e_j0.c states * for x in (2,inf) * j0(x) =3D sqrt(2/(pi*x))*(p0(x)*cos(x0)-q0(x)*sin(x0)) * where x0 =3D x-pi/4. It is better to compute sin(x0),cos(x0) * as follow: * cos(x0) =3D cos(x)cos(pi/4)+sin(x)sin(pi/4) * =3D 1/sqrt(2) * (cos(x) + sin(x)) * sin(x0) =3D sin(x)cos(pi/4)-cos(x)sin(pi/4) * =3D 1/sqrt(2) * (sin(x) - cos(x)) * (To avoid cancellation, use * sin(x) +- cos(x) =3D -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 >=3D 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 =3D 53, so we have x >=3D 18.368... Consider x=3D18.4 and sum up to N =3D 31 in the asymptotic series: % ./pq 30 18.4 p =3D 9.997932830701132e-01, q =3D -6.781826311540553e-03 <-- series pp =3D 9.997932830701132e-01, qq =3D -6.781826311540509e-03 <-- pzero,qze= ro ulp(p, pp) =3D 0.000000e+00 ulp(q, qq) =3D 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 =3D 2 and summing N =3D 4 (best results). % ./pq 4 2 p =3D 9.894313812255859e-01, q =3D -5.334472656250000e-02 pp =3D 9.862158212188928e-01, qq =3D -5.647769967932505e-02 ulp(p, pp) =3D 1.448159e+13 ulp(q, qq) =3D 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 _______________________________________________ freebsd-numerics@freebsd.org mailing list https://lists.freebsd.org/mailman/listinfo/freebsd-numerics To unsubscribe, send any mail to "freebsd-numerics-unsubscribe@freebsd.org"