From owner-freebsd-standards@FreeBSD.ORG Tue Sep 13 00:30:55 2005 Return-Path: X-Original-To: freebsd-standards@FreeBSD.org Delivered-To: freebsd-standards@FreeBSD.org Received: from mx1.FreeBSD.org (mx1.freebsd.org [216.136.204.125]) by hub.freebsd.org (Postfix) with ESMTP id 95C9016A420; Tue, 13 Sep 2005 00:30:55 +0000 (GMT) (envelope-from kargl@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 20D3A43D76; Tue, 13 Sep 2005 00:30:28 +0000 (GMT) (envelope-from kargl@troutmask.apl.washington.edu) Received: from troutmask.apl.washington.edu (localhost [127.0.0.1]) by troutmask.apl.washington.edu (8.13.4/8.13.4) with ESMTP id j8D0URma000278; Mon, 12 Sep 2005 17:30:27 -0700 (PDT) (envelope-from kargl@troutmask.apl.washington.edu) Received: (from kargl@localhost) by troutmask.apl.washington.edu (8.13.4/8.13.1/Submit) id j8D0URQm000277; Mon, 12 Sep 2005 17:30:27 -0700 (PDT) (envelope-from kargl) From: "Steven G. Kargl" Message-Id: <200509130030.j8D0URQm000277@troutmask.apl.washington.edu> In-Reply-To: <200506252340.j5PNeCQd005977@freefall.freebsd.org> To: FreeBSD-gnats-submit@FreeBSD.org, freebsd-standards@FreeBSD.org Date: Mon, 12 Sep 2005 17:30:27 -0700 (PDT) X-Mailer: ELM [version 2.4ME+ PL121h (25)] MIME-Version: 1.0 Content-Type: multipart/mixed; boundary=ELM1126571427-113-0_ Content-Transfer-Encoding: 7bit X-Content-Filtered-By: Mailman/MimeDel 2.1.5 Cc: Subject: Re: standards/82654: C99 long double math functions are missing 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: Tue, 13 Sep 2005 00:30:55 -0000 --ELM1126571427-113-0_ Content-Transfer-Encoding: 7bit Content-Type: text/plain; charset=US-ASCII FreeBSD-gnats-submit@FreeBSD.org wrote: > Thank you very much for your problem report. > It has the internal identification `standards/82654'. > The individual assigned to look at your > report is: freebsd-standards. > > You can access the state of your problem report at any time > via this link: > > http://www.freebsd.org/cgi/query-pr.cgi?pr=82654 > > >Category: standards > >Responsible: freebsd-standards > >Synopsis: C99 long double math functions are missing > >Arrival-Date: Sat Jun 25 23:40:12 GMT 2005 > Attached is a new implementation of logl(3) and diff to the exp.3 man page. The method to my implementation of the long double log should be acceptable for 64-bit precision and follow the following steps: 1) Deal with special cases 2) Split the argument into fraction, f, and exponent, n, with frexpl 3) Use a look-up table for 128 entries for log(f_n) where the interval [1/2,1) is split into 128 equally spaced intervals. 4) Use polynomial approximations for interpolation a) Start with Taylor's series where x is in [0,1/128] and truncate retaining the x**9 term, i.e., error of order 8e-22. b) Transform polynomial into range [-1,1] c) Rewrite transformed polynomial in an expansion of Chebyshev polynomials d) Reduce the Chebyshev polynomials to a polynomial of order x**8. 5) Run a boat load of tests. In the following test, "./log 1000000 10" means 1 million random arguments to logf, log, and logl are generated such that x = f*2**e and e is in [-10,10] and f is in [1/2,1). On i386 (53-bit precision), a comparison of logf, log, and logl to GMP/MPFR shows the following results. The legend is read as follows: R --> a difference of 1 in the last decimal digit 1 --> a difference of more than 1 in last decimal digit 2 --> a difference occurs in the 2nd to last decimal digit 3 --> a difference occurs in the 3rd to last decimal digit 4 --> a difference occurs in the 4th (or larger) to last decimal digit FLT --> 24-bit floating point DBL --> 53-bit double floating point LDBL -> 53-bit long double floating point on i386 -> 64-bit long double floating point on amd64 kargl[216] ./log 1000000 1 Checking logf, log, and logl ... R 1 2 3 4 Range of Failures FLT: 1128 85 18 0 10 9.90536e-01 1.09909e+00 DBL: 314483 185399 6972 7062 16245 3.67892680689692e-01 1.99999806284904e+00 LDBL: 316817 185429 6965 7071 16240 3.67892680689692e-01 1.99999806284904e+00 kargl[217] ./log 1000000 10 Checking logf, log, and logl ... R 1 2 3 4 Range of Failures FLT: 237 14 6 0 0 9.90536e-01 1.05196e+00 DBL: 141135 38632 1289 1269 2959 3.67997204419225e-01 2.71815394610167e+00 LDBL: 144160 38639 1289 1266 2957 3.67901860736310e-01 2.71815394610167e+00 kargl[218] ./log 1000000 120 Checking logf, log, and logl ... R 1 2 3 4 Range of Failures FLT: 32 0 0 0 0 DBL: 23673 3498 128 134 263 3.68287028279155e-01 2.71191326901317e+00 LDBL: 44026 3503 128 134 263 3.68287028279155e-01 2.71191326901317e+00 kargl[219] ./log 1000000 1020 Checking log and logl ... R 1 2 3 4 Range of Failures DBL: 3761 439 17 9 32 3.70691583026201e-01 2.67239280417562e+00 LDBL: 28123 437 17 10 31 3.70691583026201e-01 2.70551693812013e+00 kargl[220] ./log 1000000 16000 Checking logl ... R 1 2 3 4 Range of Failures LDBL: 30873 30 0 0 4 9.05885221436620e-01 2.32185203954577e+00 For i386 (53-bit precision), we see the implementation of logl is consistent with the implementation of log. Note, my comparison function converts the long double x to an ASCII string via scanf and the mpfr type is also converted to an ASCII string. Thus, a rounding error can occur in my numerical implementation, gdtoa, or mpfr_get_str and is the reason why I make a distinction between R and 1 in my tables. Note, the "Range of Failures" excludes the differences under R and is fairly well localized to a small range in x. On amd64 (64-bit precision), we find troutmask:sgk[211] ./log 1000000 1 Checking logf, log, and logl ... R 1 2 3 4 Range of Failures FLT: 4105 104 23 0 13 9.90094e-01 1.10419e+00 DBL: 5586 0 0 0 0 LDBL: 124664 631825 61534 8867 95289 2.50107496511191130e-01 1.99999793246388435e+00 troutmask:sgk[212] ./log 1000000 10 Checking logf, log, and logl ... R 1 2 3 4 Range of Failures FLT: 1583 17 2 0 1 9.90414e-01 1.09967e+00 DBL: 2511 0 0 0 0 LDBL: 480101 131205 12991 1654 17324 1.00236730759206694e-03 1.02341860389709473e+03 troutmask:sgk[213] ./log 1000000 120 Checking logf, log, and logl ... R 1 2 3 4 Range of Failures FLT: 254 1 0 0 0 1.05113e+00 1.05113e+00 DBL: 407 0 0 0 0 LDBL: 113720 12165 1176 151 1566 5.22837305538814689e-05 2.19174624633789062e+04 troutmask:sgk[214] ./log 1000000 1020 Checking log and logl ... R 1 2 3 4 Range of Failures DBL: 64 0 0 0 0 LDBL: 26339 1465 153 17 170 1.04225044424310909e-04 1.45063449859619141e+04 troutmask:sgk[215] ./log 1000000 16000 Checking logl ... R 1 2 3 4 Range of Failures LDBL: 15995 78 10 2 12 2.05640358899472631e-04 2.63248902320861816e+02 If I restrict the test program and logl to 53-bit precision for logl, we find all columns marked with 1, 2, 3, and 4 are zero for logl. If neither das or bde object, I would like to see logl.c committed to libm. With logl committed, I have implementations of log2f, log2, and log2l also ready for the tree. -- Steve http://troutmask.apl.washington.edu/~kargl/ --ELM1126571427-113-0_ Content-Transfer-Encoding: 7bit Content-Type: text/x-patch Content-Disposition: attachment; filename=exp.3.diff Content-Description: --- /mnt1/src/lib/msun/man/exp.3 Wed Aug 24 21:02:10 2005 +++ exp.3 Mon Sep 12 16:51:12 2005 @@ -45,8 +45,10 @@ .Nm expm1f , .Nm log , .Nm logf , +.Nm logl , .Nm log10 , .Nm log10f , +.Nm log10l , .Nm log1p , .Nm log1pf , .Nm pow , @@ -72,10 +74,14 @@ .Fn log "double x" .Ft float .Fn logf "float x" +.Ft "long double" +.Fn logl "long double x" .Ft double .Fn log10 "double x" .Ft float .Fn log10f "float x" +.Ft "long double" +.Fn log10l "long double x" .Ft double .Fn log1p "double x" .Ft float @@ -109,16 +115,18 @@ .Fa x . .Pp The -.Fn log -and the -.Fn logf +.Fn log , +.Fn logf , +and +.Fn logl functions compute the value of the natural logarithm of argument .Fa x . .Pp The -.Fn log10 -and the -.Fn log10f +.Fn log10 , +.Fn log10f , +and +.Fn log10l functions compute the value of the logarithm of argument .Fa x to base 10. --ELM1126571427-113-0_--