From mboxrd@z Thu Jan 1 00:00:00 1970 X-Spam-Checker-Version: SpamAssassin 3.4.4 (2020-01-24) on inbox.vuxu.org X-Spam-Level: X-Spam-Status: No, score=-3.3 required=5.0 tests=MAILING_LIST_MULTI, RCVD_IN_DNSWL_MED,RCVD_IN_MSPIKE_H3,RCVD_IN_MSPIKE_WL autolearn=ham autolearn_force=no version=3.4.4 Received: (qmail 8659 invoked from network); 25 Jan 2021 22:32:10 -0000 Received: from mother.openwall.net (195.42.179.200) by inbox.vuxu.org with ESMTPUTF8; 25 Jan 2021 22:32:10 -0000 Received: (qmail 3858 invoked by uid 550); 25 Jan 2021 22:32:08 -0000 Mailing-List: contact musl-help@lists.openwall.com; run by ezmlm Precedence: bulk List-Post: List-Help: List-Unsubscribe: List-Subscribe: List-ID: Reply-To: musl@lists.openwall.com Received: (qmail 3840 invoked from network); 25 Jan 2021 22:32:07 -0000 Date: Mon, 25 Jan 2021 17:31:55 -0500 From: Rich Felker To: Bruno Haible Cc: musl@lists.openwall.com Message-ID: <20210125223155.GY23432@brightrain.aerifal.cx> References: <1726645.2jRaziT4Rs@omega> MIME-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: <1726645.2jRaziT4Rs@omega> User-Agent: Mutt/1.5.21 (2010-09-15) Subject: Re: [musl] expl function imprecise on arm64 and s390x On Mon, Jan 25, 2021 at 09:55:40PM +0100, Bruno Haible wrote: > The expl() function in musl libc 1.2.2, on arm64 and s390x, returns results > that are far from as precise as the 'long double' type. > > Test program: > =============================================================================== > #ifndef __NO_MATH_INLINES > # define __NO_MATH_INLINES 1 /* for glibc */ > #endif > #include > #include > #include > #undef expl > extern > #ifdef __cplusplus > "C" > #endif > long double expl (long double); > static long double dummy (long double x) { return 0; } > int main (int argc, char *argv[]) > { > long double (* volatile my_expl) (long double) = argc ? expl : dummy; > int result = 0; > { > const long double TWO_LDBL_MANT_DIG = /* 2^LDBL_MANT_DIG */ > (long double) (1U << ((LDBL_MANT_DIG - 1) / 5)) > * (long double) (1U << ((LDBL_MANT_DIG - 1 + 1) / 5)) > * (long double) (1U << ((LDBL_MANT_DIG - 1 + 2) / 5)) > * (long double) (1U << ((LDBL_MANT_DIG - 1 + 3) / 5)) > * (long double) (1U << ((LDBL_MANT_DIG - 1 + 4) / 5)); > long double x = 11.358L; > long double y = my_expl (x); > long double z = my_expl (- x); > long double err = (y * z - 1.0L) * TWO_LDBL_MANT_DIG; > printf ("LDBL_MANT_DIG = %d\ny = %.*Lg\nz = %.*Lg\nerr = %g\n", > LDBL_MANT_DIG, > (int)(LDBL_MANT_DIG * 0.30103), y, > (int)(LDBL_MANT_DIG * 0.30103), z, > (double)err); > if (!(err >= -100.0L && err <= 100.0L)) > result |= 4; > } > return result; > } > =============================================================================== > $ gcc -Wall foo.c -lm > > Precise values: > y = 85647.901279331790624928290655026607849952741287... > z = 1.16757093292759527899833485690336857485297224766...e-5 > > Output on glibc/x86_64: > LDBL_MANT_DIG = 64 > y = 85647.90127933179059 > z = 1.167570932927595279e-05 > err = -0.5 > > Output on musl libc/arm64: > LDBL_MANT_DIG = 113 > y = 85647.90127933183975983411073684692 > z = 1.167570932927594569211357522497963e-05 > err = -1.77747e+17 > > Output on musl libc/s390x: > LDBL_MANT_DIG = 113 > y = 85647.90127933183975983411073684692 > z = 1.167570932927594569211357522497963e-05 > err = -1.77747e+17 > > You can see that only the first ca. 15 decimal digits are correct. > > Numerical approximation of transcendental functions is a long-studied > domain of numerical analysis. A recent introduction to the field, for > the exponential function, is [1]. > > Bruno > > [1] https://justinwillmert.com/articles/2020/numerically-computing-the-exponential-function-with-polynomial-approximations/ Thanks. I suspect you may find others; it's a known issue that not all of the long double functions are appropriately accurate for archs where long double is quad. sqrtl (a serious error because it's required to be correctly rounded) was fixed in the last release cycle. Hopefully we'll get the rest in better shape soon. Rich