From owner-freebsd-numerics@FreeBSD.ORG Wed Mar 26 00:22:17 2014 Return-Path: Delivered-To: freebsd-numerics@freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:1900:2254:206a::19:1]) (using TLSv1 with cipher ADH-AES256-SHA (256/256 bits)) (No client certificate requested) by hub.freebsd.org (Postfix) with ESMTPS id 136BBFEE; Wed, 26 Mar 2014 00:22:17 +0000 (UTC) Received: from troutmask.apl.washington.edu (troutmask.apl.washington.edu [128.95.76.21]) (using TLSv1.2 with cipher ECDHE-RSA-AES256-GCM-SHA384 (256/256 bits)) (No client certificate requested) by mx1.freebsd.org (Postfix) with ESMTPS id CD7F3998; Wed, 26 Mar 2014 00:22:16 +0000 (UTC) Received: from troutmask.apl.washington.edu (localhost.apl.washington.edu [127.0.0.1]) by troutmask.apl.washington.edu (8.14.8/8.14.8) with ESMTP id s2Q0M5jP010046 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-GCM-SHA384 bits=256 verify=NO); Tue, 25 Mar 2014 17:22:05 -0700 (PDT) (envelope-from sgk@troutmask.apl.washington.edu) Received: (from sgk@localhost) by troutmask.apl.washington.edu (8.14.8/8.14.8/Submit) id s2Q0M5j5010045; Tue, 25 Mar 2014 17:22:05 -0700 (PDT) (envelope-from sgk) Date: Tue, 25 Mar 2014 17:22:05 -0700 From: Steve Kargl To: freebsd-numerics@freebsd.org, freebsd-toolchain@freebsd.org Subject: clang is almost useless for complex arithmetic Message-ID: <20140326002205.GA9940@troutmask.apl.washington.edu> MIME-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline User-Agent: Mutt/1.5.23 (2014-03-12) X-BeenThere: freebsd-numerics@freebsd.org X-Mailman-Version: 2.1.17 Precedence: list List-Id: "Discussions of high quality implementation of libm functions." List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Wed, 26 Mar 2014 00:22:17 -0000 It appears that clang developers have chosen the naive complex division algorithm, and it does not matter whether one turns CX_LIMITED_RANGE on or off. This means that if one uses clang with complex types, one must be careful with the range of values allowed in complex division. In other words, implementation of complex libm routines cannot use complex data types and must fallback to a decomposition into real and imaginary components. For the record, I am using % cc --version FreeBSD clang version 3.4 (tags/RELEASE_34/final 197956) 20140216 Target: x86_64-unknown-freebsd11.0 With the program that follows my sig, if I do % cc -o z -O -DOFF cmplx.c -lm && ./z > pragma_off % cc -o z -O cmplx.c -lm && ./z > pragma_on % diff -u pragma_off pragma_on | wc -l 0 % cat pragma_on Expected: 0x1p-1023 -0x1p-1023 Compiler: 0x0p+0 -0x0p+0 Improved: 0x1p-1023 -0x1p-1023 Expected: 0x1p+1023 0x0p+0 Compiler: inf nan Improved: 0x1p+1023 0x0p+0 Expected: 0x1p+346 -0x1p-1008 Compiler: nan -0x0p+0 Improved: 0x1p+346 -0x1p-1008 Expected: 0x1p+1023 0x0p+0 Compiler: inf 0x0p+0 Improved: 0x1p+1023 0x0p+0 Expected: 0x1p+364 -0x1p-1072 Compiler: nan -0x0p+0 Improved: 0x1p+364 -0x1p-1072 Expected: 0x1p-1072 0x1p+20 Compiler: 0x0p+0 nan Improved: 0x1p-1072 0x1p+20 Expected: 0x1.ffffffffff8p+961 0x1.ffffffffff8p+982 Compiler: nan nan Improved: 0x1.ffffffffff8p+961 0x1.ffffffffff8p+982 Expected: 0x1.3333333333333p-1 0x1.999999999999ap-3 Compiler: nan nan Improved: 0x1.3333333333333p-1 0x1.999999999999ap-3 Expected: 0x1p-9 -0x1p-9 Compiler: nan nan Improved: 0x1p-9 -0x1p-9 Expected: 0x1p-279 0x1.f8p-729 Compiler: 0x1p-279 0x0p+0 Improved: 0x1p-279 0x1.f8p-729 Here, "Expected" is the expected result. "Compiler" is what clang produces. "Improved" implements the algorithm suggested in Michael Baudin, Robert L. Smith, "A Robust Complex Division in Scilab," http://arxiv.org/abs/1210.4539 -- Steve #include #include #include #include #include #ifdef OFF #pragma CX_LIMITED_RANGE off #else #pragma CX_LIMITED_RANGE on #endif typedef union { double complex f; double a[2]; } double_complex; #define REALPART(z) ((z).a[0]) #define IMAGPART(z) ((z).a[1]) static __inline double complex cpack(double x, double y) { double_complex z; REALPART(z) = x; IMAGPART(z) = y; return (z.f); } /* * For details of internal() and divi() see: * Michael Baudin, Robert L. Smith, "A Robust Complex Division in Scilab," * http://arxiv.org/abs/1210.4539 */ double complex internal(double a, double b, double c, double d, int l) { static const double ti = 0x1.p-1022, rti = 0x1.p1022, hu = 0x1.p1023, rhu = 0x1.p-1023; double e, f, r, t, ma, mi; r = d / c; if (r != 0) { if (fabs(c) < rhu) { t = c * hu + d * hu * r; e = (a * hu + b * hu * r) / t; f = (b * hu - a * hu * r) / t; } else if (fabs(c) > rti) { t = 1 / (c * ti + d * ti * r); e = (a * ti + b * ti * r) * t; f = (b * ti - a * ti * r) * t; } else { t = 1 / (c + d * r); ma = fmax(fabs(a), fabs(b)); mi = fmin(fabs(a), fabs(b)); if (ma < rti && mi * r > rhu) { e = (a + b * r) * t; f = (b - a * r) * t; } else { r = r * t; e = a * t + b * r; f = b * t - a * r; } } } else { t = 1 / c; e = (a + d * (b / c)) * t; f = (b - d * (a / c)) * t; } if (l == 1) f = -f; return (cpack(e, f)); } double complex divi(double complex z1, double complex z2) { double complex r; double a, b, c, d; a = creal(z1); b = cimag(z1); c = creal(z2); d = cimag(z2); if (fabs(d) <= fabs(c)) r = internal(a, b, c, d, 0); else r = internal(b, a, d, c, 1); return (r); } int main(void) { /* Numerators. */ double complex x[10] = { cpack(1, 1), cpack(1, 1), cpack(0x1.p1023, 0x1.p-1023), cpack(0x1.p1023, 0x1.p1023), cpack(0x1.p1020, 0x1.p-844), cpack(0x1.p-71, 0x1.p1021), cpack(0x1.p-347, 0x1.p-54), cpack(0x1.p-1074, 0x1.p-1074), cpack(0x1.p1015, 0x1.p-989), cpack(0x1.p-622, 0x1.p-1071) }; /* Denominators. */ double complex y[10] = { cpack(1, 0x1.p1023), cpack(0x1.p-1023, 0x1.p-1023), cpack(0x1.p677, 0x1.p-677), cpack(1, 1), cpack(0x1.p656, 0x1.p-780), cpack(0x1.p1001, 0x1.p-323), cpack(0x1.p-1037, 0x1.p-1058), cpack(0x1.p-1073, 0x1.p-1074), cpack(0x1.p1023, 0x1.p1023), cpack(0x1.p-343, 0x1.p-798) }; double complex r[10] = { cpack(0x1.p-1023, -0x1.p-1023), cpack(0x1.p1023, 0), cpack(0x1.p346, -0x1.p-1008), cpack(0x1.p1023, 0), cpack(0x1.p364, -0x1.p-1072), cpack(0x1.p-1072, 0x1.p20), cpack(3.898125604559113300E289, 8.174961907852353577E295), cpack(0.6, 0.2), cpack(0.001953125, -0.001953125), cpack(1.02951151789360578E-84, 6.97145987515076231E-220) }; double complex a, z[10]; int j; for (j = 0; j < 10; j++) { z[j] = x[j] / y[j]; a = divi(x[j], y[j]); printf("Expected: %la %la\n", creal(r[j]), cimag(r[j])); printf("Compiler: %la %la\n", creal(z[j]), cimag(z[j])); printf("Improved: %la %la\n", creal(a), cimag(a)); printf("\n"); } return 0; }