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.1 required=5.0 tests=DKIM_INVALID,DKIM_SIGNED, 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 5980 invoked from network); 31 Jan 2021 12:37:18 -0000 Received: from mother.openwall.net (195.42.179.200) by inbox.vuxu.org with ESMTPUTF8; 31 Jan 2021 12:37:18 -0000 Received: (qmail 29956 invoked by uid 550); 31 Jan 2021 12:37:14 -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 29926 invoked from network); 31 Jan 2021 12:37:13 -0000 DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; t=1612096622; s=strato-dkim-0002; d=clisp.org; h=References:In-Reply-To:Message-ID:Date:Subject:Cc:To:From:Cc:Date: From:Subject:Sender; bh=+C3k7+rofdlS0E9bOPMeFWJgZd+B/JNwmKfR8qWDKoU=; b=Wue2jDwEZZND65BTbRSnX6w5BxSRWJZksO2om404Vwfl8sCVaV7W0XbCRBKlQ2BNgz nEXjDNYFN+S8DI2h5sFXI9nXy0GQREEySpwUaIky5wliz2oUJkcxlbvA3bh43q3g+qmO 9UZU7LXkDY7x3q3KWhfDBbx8B8ey9lt2dVIcrGO5oQw7GwRkhiFv/AbkIohYTPIpXxkG NuthK56d+i5bDfZxqBMrghbVNcwWc421vTlCf7XLlhZRBIpTMkPyMYCmlvWGSuZVaTvT LBH6St+r+6W6cG3QTiaB0NZhomgQYfi2VTPUlWDRDtEEQ1QFn1kZPkoi9s165Q2ILLKp yN4g== X-RZG-AUTH: ":Ln4Re0+Ic/6oZXR1YgKryK8brlshOcZlIWs+iCP5vnk6shH+AHjwLuWOHqfyyvs=" X-RZG-CLASS-ID: mo00 From: Bruno Haible To: Rich Felker Cc: musl@lists.openwall.com Date: Sun, 31 Jan 2021 13:36:59 +0100 Message-ID: <3805114.vVB5bnNuiK@omega> User-Agent: KMail/5.1.3 (Linux/4.4.0-197-generic; KDE/5.18.0; x86_64; ; ) In-Reply-To: <20210125223155.GY23432@brightrain.aerifal.cx> References: <1726645.2jRaziT4Rs@omega> <20210125223155.GY23432@brightrain.aerifal.cx> MIME-Version: 1.0 Content-Transfer-Encoding: 7Bit Content-Type: text/plain; charset="us-ascii" Subject: [musl] expm1l, logl, log1pl, log2l, log10l, remainderl functions imprecise on arm64 and s390x Rich Felker wrote: > 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. ... > Hopefully we'll get the rest in better shape soon. Indeed. Here are more test cases. Hope you can use them. I computed the "correct values" using the arbitrary-precision arithmetic of GNU clisp. Test case for expm1l: =============================================================================== #ifndef __NO_MATH_INLINES # define __NO_MATH_INLINES 1 /* for glibc */ #endif #include #include #include #undef expm1l extern #ifdef __cplusplus "C" #endif long double expm1l (long double); static long double dummy (long double x) { return 0; } int main (int argc, char *argv[]) { long double (* volatile my_expm1l) (long double) = argc ? expm1l : 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_expm1l (x); long double z = my_expm1l (- x); volatile long double t = (1.0L + y) * z; long double err = (y + t) * 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 |= 1; } return result; } =============================================================================== Output on musl/arm64 and musl/s390x: LDBL_MANT_DIG = 113 y = 85646.90127933183975983411073684692 z = -0.9999883242906707492281270788225811 err = -1.11952e+22 Correct values: y = 85646.9012793317906249282906550266078499527412871 z = -0.999988324290670724047210016651430966314251470277 Test case for logl: =============================================================================== #ifndef __NO_MATH_INLINES # define __NO_MATH_INLINES 1 /* for glibc */ #endif #include #include #include #undef logl extern #ifdef __cplusplus "C" #endif long double logl (long double); static long double dummy (long double x) { return 0; } int main (int argc, char *argv[]) { long double (* volatile my_logl) (long double) = argc ? logl : 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 = 16.981137113807045L; long double y = my_logl (x); long double z = my_logl (1.0L / x); long double err = (y + z) * 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 |= 1; } return result; } =============================================================================== Output on musl/arm64 and musl/s390x: LDBL_MANT_DIG = 113 y = 2.832103146474115096253854062524624 z = -2.83210314647411554034306391258724 err = -2.30584e+18 Correct values: y = 2.83210314647411532843886058554627196598416729434 z = -2.83210314647411532843886058554627196598416729434 Test case for log1pl: =============================================================================== #ifndef __NO_MATH_INLINES # define __NO_MATH_INLINES 1 /* for glibc */ #endif #include #include #include #undef log1pl extern #ifdef __cplusplus "C" #endif long double log1pl (long double); static long double dummy (long double x) { return 0; } int main (int argc, char *argv[]) { long double (* volatile my_log1pl) (long double) = argc ? log1pl : 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_log1pl (x); long double z = my_log1pl (- x / (1.0L + x)); long double err = (y + z) * 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 >= -900.0L && err <= 900.0L)) result |= 1; } return result; } =============================================================================== Output on musl/arm64 and musl/s390x: LDBL_MANT_DIG = 113 y = 2.514303626638787925173801340861246 z = -2.51430362663878748108459149079863 err = 2.30584e+18 Correct values: y = 2.514303626638787808991104205181310474969595858254 z = -2.514303626638787808991104205181310474969595858254 Test case for log2l: =============================================================================== #ifndef __NO_MATH_INLINES # define __NO_MATH_INLINES 1 /* for glibc */ #endif #include #include #include #undef log2l extern #ifdef __cplusplus "C" #endif long double log2l (long double); static long double dummy (long double x) { return 0; } int main (int argc, char *argv[]) { long double (* volatile my_log2l) (long double) = argc ? log2l : 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_log2l (x); long double z = my_log2l (1.0L / x); long double err = (y + z) * 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 >= -10000.0L && err <= 10000.0L)) result |= 1; } if (DBL_MAX_EXP < LDBL_MAX_EXP) { long double x = ldexpl (1.0L, DBL_MAX_EXP); long double y = my_log2l (x); if (y > 0 && y + y == y) /* infinite? */ result |= 2; } return result; } =============================================================================== Output and exit code on musl/arm64 and musl/s390x: LDBL_MANT_DIG = 113 y = 3.50563691176277458794174890499562 z = -3.505636911762774143852539054933004 err = 2.30584e+18 3 Correct values: y = 3.505636911762774324741549890160677088025417867858 z = -3.505636911762774324741549890160677088025417867858 Test case for log10l: =============================================================================== #ifndef __NO_MATH_INLINES # define __NO_MATH_INLINES 1 /* for glibc */ #endif #include #include #include #undef log10l extern #ifdef __cplusplus "C" #endif long double log10l (long double); static long double dummy (long double x) { return 0; } int main (int argc, char *argv[]) { long double (* volatile my_log10l) (long double) = argc ? log10l : 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 = 7.90097792256024576L; long double y = my_log10l (x); long double z = my_log10l (1.0L / x); long double err = (y + z) * 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 |= 1; } return result; } =============================================================================== Output on musl/arm64 and musl/s390x: LDBL_MANT_DIG = 113 y = 0.8976808482634930363985859003150836 z = -0.8976808482634929253762834377994295 err = 5.76461e+17 Correct values: y = 0.897680848263492989548770123475961808070106674916 z = -0.897680848263492989548770123475961808070106674916 Test case for remainderl: =============================================================================== #ifndef __NO_MATH_INLINES # define __NO_MATH_INLINES 1 /* for glibc */ #endif #include #include #include #undef remainderl extern #ifdef __cplusplus "C" #endif long double remainderl (long double, long double); static long double dummy (long double x, long double y) { return 0; } int main (int argc, char *argv[]) { long double (* volatile my_remainderl) (long double, long double) = argc ? remainderl : 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.357996098166760874793827872740874983168L; long double y = 0.486497838502717923110029188864352615388L; long double z = my_remainderl (x, y) - (x - 23 * y); long double err = z * TWO_LDBL_MANT_DIG; printf ("LDBL_MANT_DIG = %d\nz = %.*Lg\nerr = %g\n", LDBL_MANT_DIG, (int)(LDBL_MANT_DIG * 0.30103), z, (double)err); if (!(err >= -32.0L && err <= 32.0L)) result |= 1; } return result; } =============================================================================== Output on musl/arm64 and musl/s390x: LDBL_MANT_DIG = 113 z = -6.067178037652900601980525052256813e-15 err = -3.15026e+19 Correct value: z = 0