From nobody Sat Dec 18 17:51:51 2021 X-Original-To: freebsd-current@mlmmj.nyi.freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2610:1c1:1:606c::19:1]) by mlmmj.nyi.freebsd.org (Postfix) with ESMTP id 2796418E7825; Sat, 18 Dec 2021 17:51:53 +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.3 with cipher TLS_AES_256_GCM_SHA384 (256/256 bits) key-exchange X25519 server-signature RSA-PSS (4096 bits) server-digest SHA256 client-signature RSA-PSS (2048 bits) client-digest SHA256) (Client CN "troutmask", Issuer "troutmask" (not verified)) by mx1.freebsd.org (Postfix) with ESMTPS id 4JGYMF03RXz4cfv; Sat, 18 Dec 2021 17:51:52 +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.16.1/8.16.1) with ESMTPS id 1BIHppYb071383 (version=TLSv1.3 cipher=TLS_AES_256_GCM_SHA384 bits=256 verify=NO); Sat, 18 Dec 2021 09:51:51 -0800 (PST) (envelope-from sgk@troutmask.apl.washington.edu) Received: (from sgk@localhost) by troutmask.apl.washington.edu (8.16.1/8.16.1/Submit) id 1BIHpp8v071382; Sat, 18 Dec 2021 09:51:51 -0800 (PST) (envelope-from sgk) Date: Sat, 18 Dec 2021 09:51:51 -0800 From: Steve Kargl To: Mark Murray Cc: freebsd-hackers@freebsd.org, freebsd-current@freebsd.org Subject: Re: What to do about tgammal? Message-ID: <20211218175151.GA71197@troutmask.apl.washington.edu> References: <20211204185352.GA20452@troutmask.apl.washington.edu> <20211213022223.GA41440@troutmask.apl.washington.edu> <813F29E3-8478-4282-9518-5943DE7B5492@FreeBSD.org> <20211214215106.GA50381@troutmask.apl.washington.edu> <20211218035222.GA68916@troutmask.apl.washington.edu> <6C888EBF-1734-4EDC-8DBF-D2BA2454C37D@FreeBSD.org> List-Id: Discussions about the use of FreeBSD-current List-Archive: https://lists.freebsd.org/archives/freebsd-current List-Help: List-Post: List-Subscribe: List-Unsubscribe: Sender: owner-freebsd-current@freebsd.org MIME-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: <6C888EBF-1734-4EDC-8DBF-D2BA2454C37D@FreeBSD.org> X-Rspamd-Queue-Id: 4JGYMF03RXz4cfv X-Spamd-Bar: ---- Authentication-Results: mx1.freebsd.org; none X-Spamd-Result: default: False [-4.00 / 15.00]; REPLY(-4.00)[] X-ThisMailContainsUnwantedMimeParts: N On Sat, Dec 18, 2021 at 10:41:14AM +0000, Mark Murray wrote: > > Hmm. I think my understanding of ULP is missing something? > > I thought that ULP could not be greater than the mantissa size > in bits? > > I.e., I thought it represents average rounding error (compared with > "perfect rounding"), not truncation error, as the above very large > ULPs suggest. > The definition of ULP differs according which expert you choose to follow. :-) For me (a non-expert), ULP is measured in the system of the "accurate answer", which is assumed to have many more bits of precision than the "approximate answer". From a very old das@ email and for long double I have /* From a das email: * * ulps = err * (2**(LDBL_MANT_DIG - e)), where e = ilogb(accurate). * * I use: * * mpfr_sub(err, accurate, approximate, RNDN); * mpfr_abs(err, err, RNDN); * mpfr_mul_2si(ulpx, err, -mpfr_get_exp(accurate) + LDBL_MANT_DIG, RNDN); * ulp = mpfr_get_d(ulpx, RNDN); * if (ulp 0.506 && mpfr_cmpabs(err, ldbl_minx) 0) { * print warning...; * } */ Here, "err = |accurate - approximate|" is done in the precision of the "accurate answer", and RNDN is round-to-nearest. The line 'mpfr_mul_2si(...)' is doing the scaling to manipulate the exponent of the radix 2. This is agreement with Goldberg. IIRC, Jean-Michel Muller et al have much more detailed discussion about error and ULPs in their book "Handbook of Floating-Point Arithmetic". https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html So, on ld80 and tgamma, I have % tlibm -l -x 8 -X 1755 -n 100000 -P tgamma Interval tested for tgammal: [8,1755] ulp <= 0.5: 90.207% 90207 | 90.207% 90207 0.5 < ulp < 0.6: 7.085% 7085 | 97.292% 97292 0.6 < ulp < 0.7: 2.383% 2383 | 99.675% 99675 0.7 < ulp < 0.8: 0.314% 314 | 99.989% 99989 0.8 < ulp < 0.9: 0.011% 11 | 100.000% 100000 Max ulp: 0.853190 at 9.236118561185611856579e+02 The "ulp <= 0.5" line indicates the number of correctly rounded case. If the value of ULP exceeds 1, then you are starting to lose more than 1 bit in the significant. Consider individual points. Here's a correctly rounded case: % tlibm -l -a 9.236118561185000000 tgamma xu = LD80C(0x93c72441a44a57a9, 3, 9.23611856118499999994e+00L), libmu = LD80C(0x82f85b3fe4b36f9b, 16, 6.70567128873706962722e+04L), mpfru = LD80C(0x82f85b3fe4b36f9b, 16, 6.70567128873706962722e+04L), ULP = 0.42318 Notice, ULP < 0.5 and the bits in libm and mpfr answer rounded to long double are all the same as is the decimal representation. Now consider the max ULP case: % tlibm -l -a 9.236118561185611856579e+02 tgamma xu = LD80C(0xe6e728a690c4fa87, 9, 9.23611856118561185658e+02L), libmu = LD80C(0xba6bea661a7eda3c, 7762, 5.72944398777756327202e+2336L), mpfru = LD80C(0xba6bea661a7eda3d, 7762, 5.72944398777756327245e+2336L), ULP = 0.85319 ULP is < 1, which is desireable. The bits agree except for the last place where we're off by 1, ie, 0xc vs 0xd. The decimal representation is of course off. So, now on grimoirei without my patch, I have % ./tlibm_libm -l -a 1.59255925592559255925592559255925594e+01 tgamma x = 1.59255925592559255925592559255925594e+01L libm = 1.06660136733166064453125000000000000e+12L mpfr = 1.06660136733166211967839038834033459e+12L, ULP = 13932370826141802496.00000 You have 16 correct decimal digits out of 36, which is the best you can expect from mapping tgammal to tgamma. ULP is significant larger than 1. The ulp is almost guaranteed to be larger than 2**(113-53). Now with my patch and worse case for 10000 x in [10:16] % ./tlibm_lmath -l -a 1.59255925592559255925592559255925594e+01 tgamma x = 1.59255925592559255925592559255925594e+01L libm = 1.06660136733166211967839038834033897e+12 mpfr = 1.06660136733166211967839038834033459e+12L, ULP = 41.40877 I don't print out the hex representation in ld128, but you see the number of correct decimal digits is 33 digits compared to 36. -- Steve