Where can I find the mathematical routines in GCC's source? How do the math functions work?

Question

Good night. I'm a Math bachelor and I'm studying log() and series. And I do want to see how GCC compute all those stuffs, it will help me a lot, there's nothing inside math.h I alread read it. I crazy trying to find how GCC compute logarithms and square roots using the fastest method ever. I donwloaded the source but I can't find where are the math routines.

https://github.com/gcc-mirror/gcc

I just want to see it, I'm not a good programmer at all, my thing is math.

Answer 1

Mathematical functions are part of the C standard library, and GCC just uses those. If you want to look at the source code, you can either download the source code from the official glibc website (for the GNU C Library version, which is one of the most used), or use an online code browser. Here's the code for log() for example.

Since you are saying you're not that much of a programmer though, I doubt you'll find the GNU C Standard Library comprehensible. It is the result of decades of optimizations and compatibility adjustments, and the code is very complex. I would suggest to take a look at the musl C Library instead. The source code is much cleaner and more commented. Here's the log() function, and here's all the files regarding mathematical functions.

Finally, nor GCC or the C library have "the fastest method ever" to compute such functions. The goal of the C library is not to provide the fastest possible implementation of each mathematical function, but to provide a good enough implementation while still being portable enough to be used on multiple architectures, so those are still really fast, but most likely not "the fastest ever". In the best case, some mathematical function could even be reduced to a single CPU instruction if the CPU supports fast built-in hardware mathematical operations (like for example Intel x86 with fsqrt for the square root).

Answer 2

Take a look at this log implementation.

This is from fdlibm that has the implementations (following the IEEE-754) of a lot of math functions in C for humans.

From the implementation:

Method

1. Argument Reduction: find k and f such that

    x = 2^k * (1+f),
    where sqrt(2)/2 < 1+f < sqrt(2) .

2. Approximation of log(1+f).

Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
     = 2s + 2/3 s**3 + 2/5 s**5 + .....,
         = 2s + s*R

• We use a special Reme algorithm on [0,0.1716] to generate polynomial of degree 14 to approximate R The maximum error of this polynomial approximation is bounded by 2**-58.45. In other words,

                 2      4      6      8      10      12      14
    R(z) ~ Lg1*s +Lg2*s +Lg3*s +Lg4*s +Lg5*s  +Lg6*s  +Lg7*s
    (the values of Lg1 to Lg7 are listed in the program)

and

    |      2          14          |     -58.45
    | Lg1*s +...+Lg7*s    -  R(z) | <= 2 
    |                             |

Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2. In order to guarantee error in log below 1ulp, we compute log by

    log(1+f) = f - s*(f - R)    (if f is not too large)
    log(1+f) = f - (hfsq - s*(hfsq+R)). (better accuracy)

3. Finally,

     log(x) = k*ln2 + log(1+f).  
            = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo)))

• Here ln2 is split into two floating point number:

        ln2_hi + ln2_lo,

where n*ln2_hi is always exact for |n| < 2000.

Real implementation and special cases of the explanation you can check in this link.

Answer 3

Functions like log are part of the math library that's commonly called "libm". Implementations of the standard C library, typically come with an implementation of libm, so what you are looking for is most likely in glibc. You can find the implementation of log in glibc here: https://code.woboq.org/userspace/glibc/sysdeps/ieee754/dbl-64/e_log.c.html

There are some comments in the source code that give you hints about the algorithm used, but not a detailed explanation.

Of course there are different implementations of libm - for example there is openlibm and netlib fdlibm. The documentation of both explain the algorithm used. Here is how log is implemented in openlibm: https://github.com/JuliaMath/openlibm/blob/master/src/e_log.c

_{(Interesting - looks like log in openlibm and fdlibm have come from the same source)}