standards/82654: C99 long double math functions are missing
Steven G. Kargl
kargl at troutmask.apl.washington.edu
Mon Sep 12 17:30:55 PDT 2005
FreeBSD-gnats-submit at 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/
-------------- next part --------------
A non-text attachment was scrubbed...
Name: exp.3.diff
Type: text/x-patch
Size: 861 bytes
Desc:
Url : http://lists.freebsd.org/pipermail/freebsd-standards/attachments/20050912/5309bd4b/exp.3.bin
More information about the freebsd-standards
mailing list