From owner-freebsd-standards@FreeBSD.ORG Fri Sep 28 15:24:50 2007 Return-Path: Delivered-To: freebsd-standards@freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:4f8:fff6::34]) by hub.freebsd.org (Postfix) with ESMTP id C6A7D16A4D4 for ; Fri, 28 Sep 2007 15:24:50 +0000 (UTC) (envelope-from sgk@troutmask.apl.washington.edu) Received: from troutmask.apl.washington.edu (troutmask.apl.washington.edu [128.208.78.105]) by mx1.freebsd.org (Postfix) with ESMTP id D1D9813C494 for ; Fri, 28 Sep 2007 15:24:48 +0000 (UTC) (envelope-from sgk@troutmask.apl.washington.edu) Received: from troutmask.apl.washington.edu (localhost.apl.washington.edu [127.0.0.1]) by troutmask.apl.washington.edu (8.14.1/8.14.1) with ESMTP id l8SFMRcp039284 for ; Fri, 28 Sep 2007 08:22:27 -0700 (PDT) (envelope-from sgk@troutmask.apl.washington.edu) Received: (from sgk@localhost) by troutmask.apl.washington.edu (8.14.1/8.14.1/Submit) id l8SFMRu4039283 for freebsd-standards@freebsd.org; Fri, 28 Sep 2007 08:22:27 -0700 (PDT) (envelope-from sgk) Date: Fri, 28 Sep 2007 08:22:27 -0700 From: Steve Kargl To: freebsd-standards@freebsd.org Message-ID: <20070928152227.GA39233@troutmask.apl.washington.edu> Mime-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline User-Agent: Mutt/1.4.2.3i Subject: long double broken on i386? X-BeenThere: freebsd-standards@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list List-Id: Standards compliance List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Fri, 28 Sep 2007 15:24:51 -0000 So, in my never ending crusade to implement missing C99 math functions, I decided to tackle sinl, cosl, and tanl. I have implementations that use argument reduction from k_rem_pio2.c and Chebyshev interpolation for the trig functions over the range [0,pi/4]. On amd64, these routines give <= 1/2 ULP for the various values I've tested in the reduced range. Now, if I test these on i386, I see between 400 and 1200 ULP. I spent a week or so trying to understand the descrepancy and "fix" the problem using a double-double precision arithemetic. I finally found the problem! /usr/include/float.h on i386 is lying about its numerical representation of long double. In particular, at least these 3 values appear to be wrong #define LDBL_MANT_DIG 64 #define LDBL_EPSILON 1.0842021724855044340E-19L #define LDBL_DIG 18 #include int main (void) { int i; long double a; a = 1.L; for (i = 1; i < 64; i++) { a /= 2; if (1.L - a == 1.L) break; } printf("%d %Le\n", i, a); return 0; } amd64:sgk[206] ./z 64 1.084202e-19 i386:kargl[205] ./z 54 5.551115e-17 The above results for i386 are consistent with DBL_MANT_DIG and DBL_EPSILON. Is this intentional, and should float.h be fixed? -- Steve