Bug 82654 - C99 long double math functions are missing
Summary: C99 long double math functions are missing
Status: Closed FIXED
Alias: None
Product: Base System
Classification: Unclassified
Component: standards (show other bugs)
Version: 6.0-CURRENT
Hardware: Any Any
: Normal Affects Only Me
Assignee: freebsd-numerics (Nobody)
URL:
Keywords:
Depends on:
Blocks:
 
Reported: 2005-06-26 00:40 UTC by Steven G. Kargl
Modified: 2014-06-04 18:33 UTC (History)
1 user (show)

See Also:


Attachments
file.diff (14.72 KB, patch)
2005-06-26 00:40 UTC, Steven G. Kargl
no flags Details | Diff
hypotl.c (2.03 KB, text/x-c++src)
2005-06-26 18:00 UTC, Steven G. Kargl
no flags Details
cabsl.c (1.53 KB, text/x-c++src)
2005-06-26 18:00 UTC, Steven G. Kargl
no flags Details
hypot.3.diff (952 bytes, patch)
2005-06-26 18:00 UTC, Steven G. Kargl
no flags Details | Diff
logl.c (13.06 KB, text/x-c++src)
2005-09-13 01:30 UTC, Steven G. Kargl
no flags Details
exp.3.diff (860 bytes, patch)
2005-09-13 01:30 UTC, Steven G. Kargl
no flags Details | Diff

Note You need to log in before you can comment on or make changes to this bug.
Description Steven G. Kargl 2005-06-26 00:40:12 UTC
Most of the long double math functions as prescribed by C99 are
missing.

Fix: The enclosed patch implements logl(), log10l(), sqrtl(), and cbrtl().
I'm sure someone will want bit twiddling or assembly code, but the 
c code works on both i386 and amd64.
Comment 1 Steven G. Kargl 2005-06-26 18:00:57 UTC
FreeBSD-gnats-submit@FreeBSD.org wrote:

> >Category:       standards
> >Responsible:    freebsd-standards
> >Synopsis:       C99 long double math functions are missing
> >Arrival-Date:   Sat Jun 25 23:40:12 GMT 2005
> 

The attached files implement hypotl() and cabsl().  The
documentation has been updated.  Note hypot(3) manpage
refers to a non-existent cabs.c file.

-- 
Steve
http://troutmask.apl.washington.edu/~kargl/
Comment 2 David Schultz freebsd_committer freebsd_triage 2005-06-28 21:19:48 UTC
On Sat, Jun 25, 2005, Steven G. Kargl wrote:
> The enclosed patch implements logl(), log10l(), sqrtl(), and cbrtl().
> I'm sure someone will want bit twiddling or assembly code, but the 
> c code works on both i386 and amd64.

Cool.  I don't have much time to look at this right now, but I
definitely will want to go through this later.  I have some
general comments about accuracy, though.

IEEE-754 says that algebraic functions (sqrt and cbrt in this
case) should always produce correctly rounded results.  The sqrt()
and sqrtf() implementations handle this by computing an extra bit
of the result, and using this bit and whether the remainder is 0
or not to determine which way to round.  This can be a bit tricky
to get right, but it can probably be done straightforwardly by
following fdlibm's example.  Simply using native floating-point
arithmetic as you have done will probably not suffice, unfortunately,
so this is something that definitely needs to be fixed.

The transcendental functions (e.g. logl() and log10l()) are not
required to be correctly rounded because it is not known how to
ensure correct rounding in a bounded amount of time.  However, the
guarantee made by fdlibm and most other math libraries is that it
will always be correctly rounded, except for a small percentage of
cases that are very close to halfway between two representable
numbers.  For illustration, this might mean that 0.125000000000001
gets rounded to 0.12 instead of 0.13 if we had two decimal digits
of accuracy.

Now, technically speaking, there's no *requirement* that these
transcendental functions be reasonably accurate.  The old BSD math
library often gave errors of several ulps or worse on particular
``bad'' inputs.  But it is certainly desirable that they work at
least as well as their double and float counterparts.  One could
argue that most people don't care about the last few bits of
accuracy, but some people do (think Intel Pentium bug), and I worry
that adding routines with mediocre accuracy now will mean that
nobody will bother writing better ones later.  Consider, for
instance, that glibc's implementations of fma() and most of the
complex math functions have been broken for years because they
were implemented by people who wanted to claim standards conformance
without fully understanding what they were doing.  Then again, I
can't argue too much against your implementations given that nobody
seems to have implemented more accurate BSD-licensed routines yet.
If you'd like, I can point you to what are considered the cutting-
edge papers on how to implement these functions in software, but I
don't have time to work on it myself in the forseeable future.

Two other minor points:

- Looking briefly at your logl() and log10l() implementations,
  I'm concerned about accuracy at inputs very close to 0.

- From my notes, the ``lowest-hanging fruit'', i.e. the unimplemented
  long double functions that would be the easiest to implement
  accurately, are fmodl(), remainderl(), and remquol().  These are
  easy mainly because they can be implemented as modified versions
  of their double counterparts, with a minimal amount of special-
  casing for various long double implementations.

By the way, it's really great that someone has taken an interest
in this.  One of these days, I should have more time to work on it...

--David
Comment 3 Steven G. Kargl 2005-09-13 01:30:27 UTC
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/
Comment 4 Steven G. Kargl 2005-09-14 20:54:40 UTC
My latest effort to implement logl() has some problems.
In particular, bde has shown me that catastrophic 
cancellation can lead to enormous errors for value
near 1.  So, back to the drawing broad.

-- 
Steve
http://troutmask.apl.washington.edu/~kargl/
Comment 5 Mark Linimon freebsd_committer freebsd_triage 2007-07-10 16:48:38 UTC
State Changed
From-To: open->closed

Closed at submitter's request.
Comment 6 Mark Linimon freebsd_committer freebsd_triage 2007-07-11 23:19:21 UTC
State Changed
From-To: closed->open

Although submitter has lost interest in this one, on second examination 
it sounds as though it might still be a problem, so reopen.
Comment 7 Gavin Atkinson freebsd_committer freebsd_triage 2013-06-13 12:09:50 UTC
Responsible Changed
From-To: freebsd-standards->freebsd-numerics

Over to -numerics.  It looks like at least half of this has 
already made it into releases, and the other half is in head, so 
it may be that this PR can be put into "patched" state now.
Comment 8 Steve Kargl freebsd_committer freebsd_triage 2014-06-04 18:33:46 UTC
All of the long double functions mentioned in this PR have been implemented,
and have been available in libm for over a year.