From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: from us-smtp-delivery-124.mimecast.com (us-smtp-delivery-124.mimecast.com [170.10.129.124]) by sourceware.org (Postfix) with ESMTPS id CFE0C396D816 for ; Wed, 16 Nov 2022 20:33:07 +0000 (GMT) DMARC-Filter: OpenDMARC Filter v1.4.1 sourceware.org CFE0C396D816 Authentication-Results: sourceware.org; dmarc=pass (p=none dis=none) header.from=redhat.com Authentication-Results: sourceware.org; spf=pass smtp.mailfrom=redhat.com DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=redhat.com; s=mimecast20190719; t=1668630787; h=from:from:reply-to:reply-to:subject:subject:date:date: message-id:message-id:to:to:cc:cc:mime-version:mime-version: content-type:content-type:in-reply-to:in-reply-to: references:references; bh=Rm9UbgwFPR7rmh+VSl56/Mfqg8MNt9Umqca8CyZL75o=; b=dJoJzWNFRAVF7eUES03uus9DJ7GBt8GkiLQN+NNYnof/DP9814yAlVBh3mzpTLd+KWyplJ lsocRMKi991YnNbM39KMynjQGBlJyBWSOJaBwSJepxI6e39lHHXhJZA79unEhy9kOpdODU jMV1H8ICCIEZz3+8VrXEK0vvz36g3AY= Received: from mimecast-mx02.redhat.com (mimecast-mx02.redhat.com [66.187.233.88]) by relay.mimecast.com with ESMTP with STARTTLS (version=TLSv1.2, cipher=TLS_ECDHE_RSA_WITH_AES_256_GCM_SHA384) id us-mta-659-S_XpQNoMPKapDvJ1LHgjow-1; Wed, 16 Nov 2022 15:33:04 -0500 X-MC-Unique: S_XpQNoMPKapDvJ1LHgjow-1 Received: from smtp.corp.redhat.com (int-mx10.intmail.prod.int.rdu2.redhat.com [10.11.54.10]) (using TLSv1.2 with cipher AECDH-AES256-SHA (256/256 bits)) (No client certificate requested) by mimecast-mx02.redhat.com (Postfix) with ESMTPS id 936A2833AEC; Wed, 16 Nov 2022 20:33:03 +0000 (UTC) Received: from tucnak.zalov.cz (unknown [10.39.192.21]) by smtp.corp.redhat.com (Postfix) with ESMTPS id 4F409492B04; Wed, 16 Nov 2022 20:33:03 +0000 (UTC) Received: from tucnak.zalov.cz (localhost [127.0.0.1]) by tucnak.zalov.cz (8.17.1/8.17.1) with ESMTPS id 2AGKWwxx002590 (version=TLSv1.3 cipher=TLS_AES_256_GCM_SHA384 bits=256 verify=NOT); Wed, 16 Nov 2022 21:32:58 +0100 Received: (from jakub@localhost) by tucnak.zalov.cz (8.17.1/8.17.1/Submit) id 2AGKWuEE002589; Wed, 16 Nov 2022 21:32:56 +0100 Date: Wed, 16 Nov 2022 21:32:56 +0100 From: Jakub Jelinek To: Joseph Myers , Aldy Hernandez Cc: GCC patches , Andrew MacLeod Subject: Re: [PATCH] [range-ops] Implement sqrt. Message-ID: Reply-To: Jakub Jelinek References: <20221113200553.440728-1-aldyh@redhat.com> <6150f7fd-5a57-c138-f65e-8dc3bf13d11a@codesourcery.com> MIME-Version: 1.0 In-Reply-To: <6150f7fd-5a57-c138-f65e-8dc3bf13d11a@codesourcery.com> X-Scanned-By: MIMEDefang 3.1 on 10.11.54.10 X-Mimecast-Spam-Score: 0 X-Mimecast-Originator: redhat.com Content-Type: multipart/mixed; boundary="ABBrJ4vX7X83UTxn" Content-Disposition: inline X-Spam-Status: No, score=-3.5 required=5.0 tests=BAYES_00,DKIMWL_WL_HIGH,DKIM_SIGNED,DKIM_VALID,DKIM_VALID_AU,DKIM_VALID_EF,RCVD_IN_DNSWL_NONE,RCVD_IN_MSPIKE_H2,SPF_HELO_NONE,SPF_NONE,TXREP autolearn=ham autolearn_force=no version=3.4.6 X-Spam-Checker-Version: SpamAssassin 3.4.6 (2021-04-09) on server2.sourceware.org List-Id: --ABBrJ4vX7X83UTxn Content-Type: text/plain; charset=us-ascii Content-Disposition: inline On Mon, Nov 14, 2022 at 09:55:29PM +0000, Joseph Myers wrote: > On Sun, 13 Nov 2022, Jakub Jelinek via Gcc-patches wrote: > > > So, I wonder if we don't need to add a target hook where targets will be > > able to provide upper bound on error for floating point functions for > > different floating point modes and some way to signal unknown accuracy/can't > > be trusted, in which case we would give up or return just the range for > > VARYING. > > Note that the figures given in the glibc manual are purely empirical > (largest errors observed for inputs in the glibc testsuite on a system > that was then used to update the libm-test-ulps files); they don't > constitute any kind of guarantee about either the current implementation > or the API, nor are they formally verified, nor do they come from > exhaustive testing (though worst cases from exhaustive testing for float > may have been added to the glibc testsuite in some cases). (I think the > only functions known to give huge errors for some inputs, outside of any > IBM long double issues, are the Bessel functions and cpow functions. But > even if other functions don't have huge errors, and some > architecture-specific implementations might have issues, there are > certainly some cases where errors can exceed the 9ulp threshold on what > the libm tests will accept in libm-test-ulps files, which are thus > considered glibc bugs. (That's 9ulp from the correctly rounded value, > computed in ulp of that value. For IBM long double it's 16ulp instead, > treating the format as having a fixed 106 bits of precision. Both figures > are empirical ones chosen based on what bounds sufficed for most libm > functions some years ago; ideally, with better implementations of some > functions we could probably bring those numbers down.)) I know I can't get guarantees without formal proofs and even ulps from reported errors are better than randomized testing. But I think at least for non-glibc we want to be able to get a rough idea of the usual error range in ulps. This is what I came up with so far (link with gcc -o ulp-tester{,.c} -O2 -lmpfr -lm ), it still doesn't verify that functions are always within the mathematical range of results ([-0.0, Inf] for sqrt, [-1.0, 1.0] for sin/cos etc.), guess that would be useful and verify the program actually does what is intended. One can supply just one argument (number of tests, first 46 aren't really random) or two, in the latter case the second should be upward, downward or towardzero to use non-default rounding mode. The idea is that we'd collect ballpark estimates for roundtonearest and then estimates for the other 3 rounding modes, the former would be used without -frounding-math, max over all 4 rounding modes for -frounding-math as gcc will compute using mpfr always in round to nearest. Jakub --ABBrJ4vX7X83UTxn Content-Type: text/plain; charset=us-ascii Content-Disposition: attachment; filename="ulp-tester.c" #ifdef THIS_TYPE static THIS_TYPE THIS_FUNC (ulp) (THIS_TYPE val) { if (__builtin_isnormal (val)) return THIS_FUNC (ldexp) (THIS_LIT (1.0), THIS_FUNC (ilogb) (val) - THIS_MANT_DIG + 1); else return THIS_FUNC (ldexp) (THIS_LIT (1.0), THIS_MIN_EXP - THIS_MANT_DIG); } static void THIS_FUNC (test) (THIS_TYPE (*fn) (THIS_TYPE), int (*mpfr_fn) (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t), const char *name, unsigned long count) { char buf1[256], buf2[256], buf3[256]; mpfr_set_default_prec (THIS_MANT_DIG); mpfr_t v; mpfr_init2 (v, THIS_MANT_DIG); THIS_TYPE max_ulp = THIS_LIT (0.0); volatile THIS_TYPE val = THIS_LIT (0.0); int m = 0; for (unsigned long i = 0; i < count; ++i) { if (m == 0) m = 1; else if (m <= 10) { val = THIS_FUNC (nextafter) (val, THIS_LIT (1.0)); if ((m == 1 && val == THIS_MIN) || m > 1) ++m; } else if (m == 11) { val = THIS_LIT (1.0); for (int j = 0; j < 10; j++) val = THIS_FUNC (nextafter) (val, THIS_LIT (0.0)); ++m; } else if (m <= 32) { val = THIS_FUNC (nextafter) (val, THIS_LIT (1.0)); ++m; } else if (m == 33) { val = THIS_MAX; for (int j = 0; j < 10; j++) val = THIS_FUNC (nextafter) (val, THIS_LIT (0.0)); ++m; } else if (m <= 45) { val = THIS_FUNC (nextafter) (val, THIS_LIT (1.0)); ++m; } else val = THIS_FUNC (get_rand) (); if (__builtin_isnan (val)) continue; THIS_TYPE given = fn (val); THIS_MPFR_SET (v, val, MPFR_RNDN); mpfr_fn (v, v, MPFR_RNDN); THIS_TYPE expected = THIS_MPFR_GET (v, MPFR_RNDN); if ((!__builtin_isnan (given)) != (!__builtin_isnan (expected)) || __builtin_isinf (given) != __builtin_isinf (expected)) { THIS_SNPRINTF (buf1, val); THIS_SNPRINTF (buf2, given); THIS_SNPRINTF (buf3, expected); printf ("%s (%s) = %s rather than %s\n", name, buf1, buf2, buf3); } else if (!__builtin_isnan (given) && !__builtin_isinf (given)) { THIS_TYPE this_ulp = THIS_FUNC (fabs) (given - expected) / THIS_FUNC (ulp) (expected); if (this_ulp > max_ulp) max_ulp = this_ulp; } } printf ("%s max error %.1fulp\n", name, (float) max_ulp); } #undef THIS_TYPE #undef THIS_LIT #undef THIS_FUNC #undef THIS_MIN_EXP #undef THIS_MANT_DIG #undef THIS_MIN #undef THIS_MAX #undef THIS_MPFR_SET #undef THIS_MPFR_GET #undef THIS_SNPRINTF #else #define _GNU_SOURCE 1 #include #include #include #include #include #include #if defined(__FLT128_DIG__) && defined(__GLIBC_PREREQ) #if __GLIBC_PREREQ (2, 26) #define TEST_FLT128 #define MPFR_WANT_FLOAT128 #endif #endif #include #include static long rand_n; static int rand_c; static uint32_t get_rand32 (void) { uint32_t ret = 0; if (rand_c == 0) { ret = random () & 0x7fffffff; rand_c = 31; } else ret = rand_n & (((uint32_t) 1 << rand_c) - 1); ret <<= (32 - rand_c); rand_n = random (); ret |= rand_n & (((uint32_t) 1 << (32 - rand_c)) - 1); rand_n >>= (32 - rand_c); rand_c = 31 - (32 - rand_c); return ret; } static uint64_t get_rand64 (void) { return (((uint64_t) get_rand32 ()) << 32) | get_rand32 (); } static float get_randf (void) { uint32_t i = get_rand32 (); float f; memcpy (&f, &i, sizeof (f)); return f; } static double get_rand (void) { uint64_t i = get_rand64 (); double d; memcpy (&d, &i, sizeof (d)); return d; } static long double get_randl (void) { long double ld; uint64_t i = get_rand64 (); memcpy (&ld, &i, sizeof (i)); if (sizeof (long double) == 12) { uint32_t j = get_rand32 (); memcpy ((char *) &ld + 8, &j, sizeof (j)); } else if (sizeof (long double) == 16) { i = get_rand64 (); memcpy ((char *) &ld + 8, &i, sizeof (i)); } return ld; } #ifdef TEST_FLT128 static long double get_randf128 (void) { _Float128 f128; uint64_t i = get_rand64 (); memcpy (&f128, &i, sizeof (i)); i = get_rand64 (); memcpy ((char *) &f128 + 8, &i, sizeof (i)); return f128; } #endif #define THIS_TYPE float #define THIS_LIT(v) v##f #define THIS_FUNC(v) v##f #define THIS_MIN_EXP __FLT_MIN_EXP__ #define THIS_MANT_DIG __FLT_MANT_DIG__ #define THIS_MIN __FLT_MIN__ #define THIS_MAX __FLT_MAX__ #define THIS_MPFR_SET mpfr_set_flt #define THIS_MPFR_GET mpfr_get_flt #define THIS_SNPRINTF(buf, x) snprintf ((buf), sizeof (buf), "%a", (x)); #include "ulp-tester.c" #define THIS_TYPE double #define THIS_LIT(v) v #define THIS_FUNC(v) v #define THIS_MIN_EXP __DBL_MIN_EXP__ #define THIS_MANT_DIG __DBL_MANT_DIG__ #define THIS_MIN __DBL_MIN__ #define THIS_MAX __DBL_MAX__ #define THIS_MPFR_SET mpfr_set_d #define THIS_MPFR_GET mpfr_get_d #define THIS_SNPRINTF(buf, x) snprintf ((buf), sizeof (buf), "%a", (x)); #include "ulp-tester.c" #define THIS_TYPE long double #define THIS_LIT(v) v##L #define THIS_FUNC(v) v##l #define THIS_MIN_EXP __LDBL_MIN_EXP__ #define THIS_MANT_DIG __LDBL_MANT_DIG__ #define THIS_MIN __LDBL_MIN__ #define THIS_MAX __LDBL_MAX__ #define THIS_MPFR_SET mpfr_set_ld #define THIS_MPFR_GET mpfr_get_ld #define THIS_SNPRINTF(buf, x) snprintf ((buf), sizeof (buf), "%La", (x)); #include "ulp-tester.c" #ifdef TEST_FLT128 #define THIS_TYPE _Float128 #define THIS_LIT(v) v##F128 #define THIS_FUNC(v) v##f128 #define THIS_MIN_EXP __FLT128_MIN_EXP__ #define THIS_MANT_DIG __FLT128_MANT_DIG__ #define THIS_MIN __FLT128_MIN__ #define THIS_MAX __FLT128_MAX__ #define THIS_MPFR_SET mpfr_set_float128 #define THIS_MPFR_GET mpfr_get_float128 #define THIS_SNPRINTF(buf, x) strfromf128 ((buf), sizeof (buf), "%a", (x)); #include "ulp-tester.c" #else #define testf128(fn, mpfr_fn, name, count) do { } while (0) #endif int main (int argc, const char **argv) { const char *arg; char *endptr; (void) argc; if (argc <= 1) arg = ""; else arg = argv[1]; unsigned long count = strtoul (arg, &endptr, 10); if (endptr == arg) { fprintf (stderr, "ulp-tester number_of_iterations rnd\n"); return 1; } const char *rnd = "tonearest"; if (argc >= 3) rnd = argv[2]; if (strcmp (rnd, "upward") == 0) fesetround (FE_UPWARD); else if (strcmp (rnd, "downward") == 0) fesetround (FE_DOWNWARD); else if (strcmp (rnd, "towardzero") == 0) fesetround (FE_TOWARDZERO); #define TESTS(fn) \ testf (fn##f, mpfr_##fn, #fn "f", count); \ test (fn, mpfr_##fn, #fn, count); \ testl (fn##l, mpfr_##fn, #fn "l", count); \ testf128 (fn##f128, mpfr_##fn, #fn "f128", count) TESTS (sqrt); TESTS (sin); TESTS (cos); TESTS (exp10); } #endif --ABBrJ4vX7X83UTxn--