Re: Support of log10?

[Nelson H. F. Beebe]

Peng Yu <pengyu.ut@gmail.com> reports this output:

>> ...
>> $ seq 3 | awk -v 'OFS=\t' -v OFMT=%.20g -e 'function log10(x) { return
>> log(x)/log(10); } { print "log10(" 10^$1 ")", log10(10^$1) }'
>> log10(10)    1
>> log10(100)   2
>> log10(1000)  2.9999999999999995559
>> 
>> I don't find a native support of log10(). But if I implement log10
>> with log, I got the above problem. The result of log10(1000) should be
>> exact 3. But it is not.
>> ...

Peng, you need to understand more about computer arithmetic: because
it has finite precision and range, it is NOT the same as mathematical
arithmetic.  Consider this simple C program:

% cat foo.c
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

int
main()
{
    double x, y;
    const double XMAX = 1.e15;

    (void)printf("Test of log10(x)\n");
    for (x = 1; x < XMAX; x*= 10)
    {
        y = log10(x);
        (void)printf("%10.1e\t%.17g\t%a\n", x, y, y);
    }

    (void)printf("\n");
    (void)printf("Test of log(x) / log(10)\n");

    for (x = 1; x < XMAX; x*= 10)
    {
        y = log(x) / log(10);
        (void)printf("%10.1e\t%.17g\t%a\n", x, y, y);
    }
    
    return (EXIT_SUCCESS);
}

Here is what it's output is on a CentOS 7.9 system on x86_64:

% cc foo.2 -lm && ./a.out
Test of log10(x)
   1.0e+00      0       0x0p+0
   1.0e+01      1       0x1p+0
   1.0e+02      2       0x1p+1
   1.0e+03      3       0x1.8p+1
   1.0e+04      4       0x1p+2
   1.0e+05      5       0x1.4p+2
   1.0e+06      6       0x1.8p+2
   1.0e+07      7       0x1.cp+2
   1.0e+08      8       0x1p+3
   1.0e+09      9       0x1.2p+3
   1.0e+10      10      0x1.4p+3
   1.0e+11      11      0x1.6p+3
   1.0e+12      12      0x1.8p+3
   1.0e+13      13      0x1.ap+3
   1.0e+14      14      0x1.cp+3

Test of log(x) / log(10)
   1.0e+00      0       0x0p+0
   1.0e+01      1       0x1p+0
   1.0e+02      2       0x1p+1
   1.0e+03      2.9999999999999996      0x1.7ffffffffffffp+1
   1.0e+04      4       0x1p+2
   1.0e+05      5       0x1.4p+2
   1.0e+06      5.9999999999999991      0x1.7ffffffffffffp+2
   1.0e+07      7       0x1.cp+2
   1.0e+08      8       0x1p+3
   1.0e+09      8.9999999999999982      0x1.1ffffffffffffp+3
   1.0e+10      10      0x1.4p+3
   1.0e+11      11      0x1.6p+3
   1.0e+12      11.999999999999998      0x1.7ffffffffffffp+3
   1.0e+13      12.999999999999998      0x1.9ffffffffffffp+3
   1.0e+14      14      0x1.cp+3

A careful implementation of log10(x) inside this system's -lm library
gets the expected answers for the EXACTLY-REPRESENTABLE argument; the
simple-minded division with log(x)/log(10) does not.

Please consult any of several books and papers on computer arithmetic,
recorded in

        http://www.math.utah.edu/pub/tex/bib/fparith.html
        http://www.math.utah.edu/pub/tex/bib/fparith.bib

That file has more than 6800 references, which should give some idea
that the subject is far from trivial.

A second bibliography on elementary and special functions

        http://www.math.utah.edu/pub/tex/bib/elefunt.html
        http://www.math.utah.edu/pub/tex/bib/elefunt.bib

collects almost 2400 references.

If you want just a couple of freely-accessible documents on the
subject, look at these two files:

        http://www.math.utah.edu/~beebe/fp.pdf
        http://www.math.utah.edu/~beebe/arith.pdf

The second is an extension of the first.

-------------------------------------------------------------------------------
- Nelson H. F. Beebe                    Tel: +1 801 581 5254                  -
- University of Utah                    FAX: +1 801 581 4148                  -
- Department of Mathematics, 110 LCB    Internet e-mail: beebe@math.utah.edu  -
- 155 S 1400 E RM 233                       beebe@acm.org  beebe@computer.org -
- Salt Lake City, UT 84112-0090, USA    URL: http://www.math.utah.edu/~beebe/ -
-------------------------------------------------------------------------------