From owner-svn-src-stable@freebsd.org Tue Sep 12 00:26:58 2017 Return-Path: Delivered-To: svn-src-stable@mailman.ysv.freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:1900:2254:206a::19:1]) by mailman.ysv.freebsd.org (Postfix) with ESMTP id 2478AE0A509; Tue, 12 Sep 2017 00:26:58 +0000 (UTC) (envelope-from rlibby@FreeBSD.org) Received: from repo.freebsd.org (repo.freebsd.org [IPv6:2610:1c1:1:6068::e6a:0]) (using TLSv1.2 with cipher ECDHE-RSA-AES256-GCM-SHA384 (256/256 bits)) (Client did not present a certificate) by mx1.freebsd.org (Postfix) with ESMTPS id 0023A7DBA6; Tue, 12 Sep 2017 00:26:57 +0000 (UTC) (envelope-from rlibby@FreeBSD.org) Received: from repo.freebsd.org ([127.0.1.37]) by repo.freebsd.org (8.15.2/8.15.2) with ESMTP id v8C0QvoD014068; Tue, 12 Sep 2017 00:26:57 GMT (envelope-from rlibby@FreeBSD.org) Received: (from rlibby@localhost) by repo.freebsd.org (8.15.2/8.15.2/Submit) id v8C0Qvii014065; Tue, 12 Sep 2017 00:26:57 GMT (envelope-from rlibby@FreeBSD.org) Message-Id: <201709120026.v8C0Qvii014065@repo.freebsd.org> X-Authentication-Warning: repo.freebsd.org: rlibby set sender to rlibby@FreeBSD.org using -f From: Ryan Libby Date: Tue, 12 Sep 2017 00:26:56 +0000 (UTC) To: src-committers@freebsd.org, svn-src-all@freebsd.org, svn-src-stable@freebsd.org, svn-src-stable-11@freebsd.org Subject: svn commit: r323475 - in stable/11/lib/msun: src tests X-SVN-Group: stable-11 X-SVN-Commit-Author: rlibby X-SVN-Commit-Paths: in stable/11/lib/msun: src tests X-SVN-Commit-Revision: 323475 X-SVN-Commit-Repository: base MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit X-BeenThere: svn-src-stable@freebsd.org X-Mailman-Version: 2.1.23 Precedence: list List-Id: SVN commit messages for all the -stable branches of the src tree List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Tue, 12 Sep 2017 00:26:58 -0000 Author: rlibby Date: Tue Sep 12 00:26:56 2017 New Revision: 323475 URL: https://svnweb.freebsd.org/changeset/base/323475 Log: MFC r323003,r323004: r323003: lib/msun: avoid referring to broken LDBL_MAX r323004: lib/msun: add more csqrt unit tests for precision and overflow Modified: stable/11/lib/msun/src/catrig.c stable/11/lib/msun/src/s_csqrtl.c stable/11/lib/msun/tests/csqrt_test.c Directory Properties: stable/11/ (props changed) Modified: stable/11/lib/msun/src/catrig.c ============================================================================== --- stable/11/lib/msun/src/catrig.c Mon Sep 11 23:47:49 2017 (r323474) +++ stable/11/lib/msun/src/catrig.c Tue Sep 12 00:26:56 2017 (r323475) @@ -469,8 +469,13 @@ clog_for_large_values(double complex z) /* * Avoid overflow in hypot() when x and y are both very large. - * Divide x and y by E, and then add 1 to the logarithm. This depends - * on E being larger than sqrt(2). + * Divide x and y by E, and then add 1 to the logarithm. This + * depends on E being larger than sqrt(2), since the return value of + * hypot cannot overflow if neither argument is greater in magnitude + * than 1/sqrt(2) of the maximum value of the return type. Likewise + * this determines the necessary threshold for using this method + * (however, actually use 1/2 instead as it is simpler). + * * Dividing by E causes an insignificant loss of accuracy; however * this method is still poor since it is uneccessarily slow. */ Modified: stable/11/lib/msun/src/s_csqrtl.c ============================================================================== --- stable/11/lib/msun/src/s_csqrtl.c Mon Sep 11 23:47:49 2017 (r323474) +++ stable/11/lib/msun/src/s_csqrtl.c Tue Sep 12 00:26:56 2017 (r323475) @@ -42,8 +42,16 @@ __FBSDID("$FreeBSD$"); */ #pragma STDC CX_LIMITED_RANGE ON -/* We risk spurious overflow for components >= LDBL_MAX / (1 + sqrt(2)). */ -#define THRESH (LDBL_MAX / 2.414213562373095048801688724209698L) +/* + * We risk spurious overflow for components >= LDBL_MAX / (1 + sqrt(2)). + * Rather than determining the fully precise value at which we might + * overflow, just use a threshold of approximately LDBL_MAX / 4. + */ +#if LDBL_MAX_EXP != 0x4000 +#error "Unsupported long double format" +#else +#define THRESH 0x1p16382L +#endif long double complex csqrtl(long double complex z) Modified: stable/11/lib/msun/tests/csqrt_test.c ============================================================================== --- stable/11/lib/msun/tests/csqrt_test.c Mon Sep 11 23:47:49 2017 (r323474) +++ stable/11/lib/msun/tests/csqrt_test.c Tue Sep 12 00:26:56 2017 (r323475) @@ -214,28 +214,94 @@ test_nans(void) /* * Test whether csqrt(a + bi) works for inputs that are large enough to - * cause overflow in hypot(a, b) + a. In this case we are using - * csqrt(115 + 252*I) == 14 + 9*I - * scaled up to near MAX_EXP. + * cause overflow in hypot(a, b) + a. Each of the tests is scaled up to + * near MAX_EXP. */ static void test_overflow(int maxexp) { long double a, b; long double complex result; + int exp, i; - a = ldexpl(115 * 0x1p-8, maxexp); - b = ldexpl(252 * 0x1p-8, maxexp); - result = t_csqrt(CMPLXL(a, b)); - assert(creall(result) == ldexpl(14 * 0x1p-4, maxexp / 2)); - assert(cimagl(result) == ldexpl(9 * 0x1p-4, maxexp / 2)); + assert(maxexp > 0 && maxexp % 2 == 0); + + for (i = 0; i < 4; i++) { + exp = maxexp - 2 * i; + + /* csqrt(115 + 252*I) == 14 + 9*I */ + a = ldexpl(115 * 0x1p-8, exp); + b = ldexpl(252 * 0x1p-8, exp); + result = t_csqrt(CMPLXL(a, b)); + assert(creall(result) == ldexpl(14 * 0x1p-4, exp / 2)); + assert(cimagl(result) == ldexpl(9 * 0x1p-4, exp / 2)); + + /* csqrt(-11 + 60*I) = 5 + 6*I */ + a = ldexpl(-11 * 0x1p-6, exp); + b = ldexpl(60 * 0x1p-6, exp); + result = t_csqrt(CMPLXL(a, b)); + assert(creall(result) == ldexpl(5 * 0x1p-3, exp / 2)); + assert(cimagl(result) == ldexpl(6 * 0x1p-3, exp / 2)); + + /* csqrt(225 + 0*I) == 15 + 0*I */ + a = ldexpl(225 * 0x1p-8, exp); + b = 0; + result = t_csqrt(CMPLXL(a, b)); + assert(creall(result) == ldexpl(15 * 0x1p-4, exp / 2)); + assert(cimagl(result) == 0); + } } +/* + * Test that precision is maintained for some large squares. Set all or + * some bits in the lower mantdig/2 bits, square the number, and try to + * recover the sqrt. Note: + * (x + xI)**2 = 2xxI + */ +static void +test_precision(int maxexp, int mantdig) +{ + long double b, x; + long double complex result; + uint64_t mantbits, sq_mantbits; + int exp, i; + + assert(maxexp > 0 && maxexp % 2 == 0); + assert(mantdig <= 64); + mantdig = rounddown(mantdig, 2); + + for (exp = 0; exp <= maxexp; exp += 2) { + mantbits = ((uint64_t)1 << (mantdig / 2 )) - 1; + for (i = 0; + i < 100 && mantbits > ((uint64_t)1 << (mantdig / 2 - 1)); + i++, mantbits--) { + sq_mantbits = mantbits * mantbits; + /* + * sq_mantibts is a mantdig-bit number. Divide by + * 2**mantdig to normalize it to [0.5, 1), where, + * note, the binary power will be -1. Raise it by + * 2**exp for the test. exp is even. Lower it by + * one to reach a final binary power which is also + * even. The result should be exactly + * representable, given that mantdig is less than or + * equal to the available precision. + */ + b = ldexpl((long double)sq_mantbits, + exp - 1 - mantdig); + x = ldexpl(mantbits, (exp - 2 - mantdig) / 2); + assert(b == x * x * 2); + result = t_csqrt(CMPLXL(0, b)); + assert(creall(result) == x); + assert(cimagl(result) == x); + } + } +} + int main(void) { - printf("1..15\n"); + printf("1..18\n"); /* Test csqrt() */ t_csqrt = _csqrt; @@ -255,41 +321,56 @@ main(void) test_overflow(DBL_MAX_EXP); printf("ok 5 - csqrt\n"); + test_precision(DBL_MAX_EXP, DBL_MANT_DIG); + printf("ok 6 - csqrt\n"); + /* Now test csqrtf() */ t_csqrt = _csqrtf; test_finite(); - printf("ok 6 - csqrt\n"); + printf("ok 7 - csqrt\n"); test_zeros(); - printf("ok 7 - csqrt\n"); + printf("ok 8 - csqrt\n"); test_infinities(); - printf("ok 8 - csqrt\n"); + printf("ok 9 - csqrt\n"); test_nans(); - printf("ok 9 - csqrt\n"); + printf("ok 10 - csqrt\n"); test_overflow(FLT_MAX_EXP); - printf("ok 10 - csqrt\n"); + printf("ok 11 - csqrt\n"); + test_precision(FLT_MAX_EXP, FLT_MANT_DIG); + printf("ok 12 - csqrt\n"); + /* Now test csqrtl() */ t_csqrt = csqrtl; test_finite(); - printf("ok 11 - csqrt\n"); + printf("ok 13 - csqrt\n"); test_zeros(); - printf("ok 12 - csqrt\n"); + printf("ok 14 - csqrt\n"); test_infinities(); - printf("ok 13 - csqrt\n"); + printf("ok 15 - csqrt\n"); test_nans(); - printf("ok 14 - csqrt\n"); + printf("ok 16 - csqrt\n"); test_overflow(LDBL_MAX_EXP); - printf("ok 15 - csqrt\n"); + printf("ok 17 - csqrt\n"); + + test_precision(LDBL_MAX_EXP, +#ifndef __i386__ + LDBL_MANT_DIG +#else + DBL_MANT_DIG +#endif + ); + printf("ok 18 - csqrt\n"); return (0); }