# 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