standards/82654: C99 long double math functions are missing

David Schultz das at FreeBSD.ORG
Tue Jun 28 20:30:15 GMT 2005


The following reply was made to PR standards/82654; it has been noted by GNATS.

From: David Schultz <das at FreeBSD.ORG>
To: "Steven G. Kargl" <kargl at troutmask.apl.washington.edu>
Cc: FreeBSD-gnats-submit at FreeBSD.ORG
Subject: Re: standards/82654: C99 long double math functions are missing
Date: Tue, 28 Jun 2005 16:19:48 -0400

 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


More information about the freebsd-standards mailing list