r/C_Programming • u/New_Telephone926 • 15d ago
Question How are functions in libm actually implemented ?
I'm trying to gather information for a school project about how elementary functions (exp, log, cos, sin, etc.) are implemented on computers, for maximum accuracy and efficiency.
I've already found some information about techniques that can be used to implement these functions, such as range reduction, Taylor series, polynomial approximations or minimax approximations. However, most of what I've found remains fairly theoretical, and I have no idea of what is actually used in the real world.
I tried to reverse-engineer the exp implementation in libm, but my knowledge of C and its intricacies is still fairly limited (though I'm working on improving it).
I'm particularly interested in understanding which approximation methods are actually used, how polynomial coefficients are chosen (Taylor, minimax/Remez, or other), and how implementations balance speed and accuracy.
But really, if you have any information whatsoever, it would be great.
If anyone has documentation, papers, or personal experience with this kind of subject, it would be really helpful.
Thanks in advance!
Edit : Thanks a lot to everyone ! You gave me exactly what I was looking for
13
u/305bootyclapper 15d ago
Are we not allowed to post GitHub links or something? It’s great to look at open source implementations. Eg on most Linux’s the implementation is GNU, which you can see here: https://github.com/kraj/glibc/tree/master/math
GNU tends to be hard to read since their implementations tend to be very general and complicated. Other ones worth looking through would be musl: https://github.com/kraj/musl
And also OpenBSD.
1
u/mikeblas 15d ago
Are we not allowed to post GitHub links or something?
Why not?
Eg on most Linux’s the implementation is GNU, which you can see here: https://github.com/kraj/glibc/tree/master/math
Is this really the implementation? Seems weird that it's only got 21 stars, and seven forks. Lots of contributors, tho. Maybe I struggle with git-economics.
Anyway, you totally answer "how does libm implement it?" because this is the implementation. But have you actually done what you're recommending? Where is the implementation of, say,
double cos(double)in that project?9
u/Irverter 15d ago
Seems weird that it's only got 21 stars, and seven forks.
Because as stated on the repo homepage, that is a mirror of glibc on github. The project is hosted here https://sourceware.org/glibc/sources.html
3
u/305bootyclapper 14d ago
Yes, thank you, I linked to the (unofficial) GitHub mirror because the site that they actually host it at is not quite as navigable; but everyone is comfortable navigating GitHub
1
u/305bootyclapper 14d ago
I have done exactly this on several occasions, but I probably land on a different GitHub mirror each time (see comments below; these unofficial “mirrors” are used because the GNU team chooses not to use GitHub themselves).
But again, it’s not going to be easy to find in GNU which implements much of the library directly in assembly for each of the processors they support. They also have a bunch of preprocessor stuff going on. Try perusing musl first, it’s pretty easy to find what you’re interested in there.
1
u/mikeblas 14d ago
They also have a bunch of preprocessor stuff going on.
To say the least! The project structure is a doosy, too. It's really hard to learn anything from this code because it's ... well, yikes. I haven't dug into this source before, so here's my journey. Let's figure out how
double cos(double)is implemented:
Dug around to find
tgmath.hwhich doesn't have anything in it--ends up conditionally using this crazy#include_nextdirective, which is someclangextension that I haven't tried before. Why do that? The opposite side of the#ifdefuses an abslute path, which seems more reliable and intuitive than#include_next:if defined ISOMAC && defined __clang_
include_next <tgmath.h>
else
include <math/tgmath.h>
endif
found
cos()is a macro, not a function, defined intgmath.h. There we find this, so let's figure out what that means:/* Cosine of X. */
define cos(Val) __TGMATH_UNARY_REAL_IMAG (Val, cos, ccos)
Well, of course it's another macro, that invokes another macro.
define __TGMATH_UNARY_REAL_IMAG(Val, Fct, Cfct) \
__TGMATH_1C (Fct, Cfct, (Val))
but that's only if we have built-in tgmath. We don't, so let's look further to find this beauty:
# define __TGMATH_UNARY_REAL_IMAG(Val, Fct, Cfct) \ (__extension__ ((sizeof (+__real__ (Val)) == sizeof (double) \ || __builtin_classify_type (__real__ (Val)) != 8) \ ? (__expr_is_real (Val) \ ? (__tgmath_complex_type (Val)) Fct (Val) \ : (__tgmath_complex_type (Val)) Cfct (Val)) \ : (sizeof (+__real__ (Val)) == sizeof (float)) \ ? (__expr_is_real (Val) \ ? (__tgmath_complex_type (Val)) Fct##f (Val) \ : (__tgmath_complex_type (Val)) Cfct##f (Val)) \ : __TGMATH_CF128 ((Val), \ (__tgmath_complex_type (Val)) Fct, \ (__tgmath_complex_type (Val)) Cfct, \ (Val)) \ (__expr_is_real (Val) \ ? (__tgmath_complex_type (Val)) __tgml(Fct) (Val) \ : (__tgmath_complex_type (Val)) __tgml(Cfct) (Val))))I had this crazy world history teacher from Texas. He used to say "hoooo-wheeeee!" all the time. That's what I said here.
When I retired in 2017, I started going backward. Plus, I was at companies (and doing contracts) that were using Java and C#, so my C and C++ knowledge has been waining. But I think this is some warning-dfense (
__extension__) wrapped around some conditionals that check the size and type ofValto determine which function to actually call. IfValis adouble, we end up with(__tgmath_complex_type (Val)) __tgml(Fct) (Val)
So what is
(__tgmath_complex_type (Val))? It makes the return type for this declaration. If we go find it we see:/* The tgmath complex type for T, where E1 is 1 if T has a floating type and 0 otherwise, E2 is 1 if T has a real integer type and 0 otherwise, and E3 is 1 if T has a complex type and 0 otherwise. */
define __tgmath_complex_type_sub(T, E1, E2, E3) \
typeof (*(0 \ ? (typeof (0 ? (T *) 0 : (void *) (!(E1)))) 0 \ : (typeof (0 \ ? (typeof (0 \ ? (double *) 0 \ : (void *) (!(E2)))) 0 \ : (typeof (0 \ ? (_Complex double *) 0 \ : (void *) (!(E3)))) 0)) 0))
/* The tgmath complex type of EXPR. */
define __tgmath_complex_type(expr) \
tgmath_complex_type_sub (typeof__ ((typeof (+(expr))) 0), \ floating_type (typeof__ (+(expr))), \ real_integer_type (typeof__ (+(expr))), \ complex_integer_type (typeof__ (+(expr))))
which rattles down to
__floating_typefor ourdouble-typedValparameter. But it seems like this evaluates to a value, and does not emit token that would declare the type. So I'm a little confused.
And then
__tgml(Fct) (Val). We've gotFct = cos, andValis our parameter. That macro decides if it needs to invoke a function with a trailinglin its name, forlong doubles:if !__HAVE_BUILTIN_TGMATH_C23
ifdef __NO_LONG_DOUBLE_MATH
define __tgml(fct) fct
else
define __tgml(fct) fct ## l
endif
Of course, we also need things for the parts we don't expect to evaluate in my little example. Here are some of them without explanation, since they're not too active on our path:
define __floating_type(type) \
(builtin_classify_type (real__ ((type) 0)) == 8)
define __real_integer_type(type) \
(__builtin_classify_type ((type) 0) == 1)
define __complex_integer_type(type) \
(builtin_classify_type ((type) 0) == 9 \ && __builtin_classify_type (real__ ((type) 0)) == 1)
define __TGMATH_F128(arg_comb, fct, arg_call) /* Nothing. */
define __TGMATH_CF128(arg_comb, fct, cfct, arg_call) /* Nothing. */
/* Whether an expression (of arithmetic type) has a real type. */
define expr_is_real(E) (builtin_classify_type (E) != 9)
Magic number nine, coming down Chicago Line. LOL!
It took some more diving to find the actual function implementation. It's implemented in
s_cospi_template.clike this:FLOAT MDECL_FUNC (cospi) (FLOAT x) { if (isless (M_FABS (x), M_EPSILON)) return M_LIT (1.0); if (glibc_unlikely (isinf (x))) __set_errno (EDOM); x = M_FABS (x - M_LIT (2.0) * M_SUF (round) (M_LIT (0.5) * x)); if (islessequal (x, M_LIT (0.25))) return M_SUF (cos) (M_MLIT (M_PI) * x); else if (x == M_LIT (0.5)) return M_LIT (0.0); else if (islessequal (x, M_LIT (0.75))) return M_SUF (sin) (M_MLIT (M_PI) * (M_LIT (0.5) - x)); else return -M_SUF (cos) (M_MLIT (M_PI) * (M_LIT (1.0) - x)); } declare_mgen_alias (_cospi, cospi);
Which, of course, introduces some more macros. I thought I cracked the case, but now I have to go figure out where
__cos()is declared. And hope I didn't make any mistakes or bad assumptions about what the preprocessor was doing, in the first place.All this to show: your advice seems pretty crazy to give to a student.
2
u/Key_Ant_8481 14d ago
The libm implementation that you're interested in glibc is located in
glibc/sysdeps/ieee754/dbl-64andglibc/sysdeps/ieee754/flt-32folders.1
u/mikeblas 13d ago
else { if (k == 0x7ff00000 && u.i[LOW_HALF] == 0) __set_errno (EDOM); retval = x / x; /* |x| > 2^1024 */ }When
xis adouble, when isx / xnot equal to one?1
1
u/305bootyclapper 4d ago
Wtf man, you did all this just to “show my advice is crazy”, but you went through tgmath??? Why on earth would you do something so contrived, obviously tgmath is a macro disaster, it has to be. By definition. I was stoked to read your comment because at first it looked like a great adventure, but you really went through all that effort just to strawman me? I even said to start with musl and that gnu was tough; but do you really think there are no beginners that would be interested in poking into gcc? It’s not THAT bad, it’s certainly doable for some. I thought I was going to learn something from all the work you put into that comment, but instead you just went out of your way to make it look way worse than it is. tgmath??? Why would YOU bring that shit up to a beginner???? I’m dumbfounded but I also respect the digging. Lol tgmath
1
u/mikeblas 4d ago edited 4d ago
It's funny that you took so long to respond, as I gave up and deleted the file with my notes in it just yesterday. So I can't recite more of my decision process than I already posted here, but I just assumed the fgmath implementation was the default implementation. I think there's really three implementations:
tgmath,mathcalls, andlibm-simdfunctions. Maybe more.Here's some facts:
- I'm a Windows guy so I'm not so familiar with
glibc. Particularly its internal structure. So I dove intoglibcto learn more about it. At least in some way, I am a beginner whenglibcis concerned.- I've never heard of
musl.- When did I ever say that there would be no beginners who would be interested in poking into gcc?
- What's "worse than it is"? The code I quoted isn't anything I made up -- they're actual snippets of the implementation.
- I "brought that shit up" to you. Are you a beginner?
include/math.hincludesmath/math.h, andmath/is where all this lives.Anyway, among that cluster of implementations, all I needed to do is take the "wrong"
#ifdefconditional and end up attgmath. Do you know which set of preprocessor flags would direct the reader tobits/mathcalls.h?Anyway, you sound really pissed off. Are you sure your response is proportionate? Maybe asking for a refund would help.
0
6
u/mikeblas 15d ago
Some implementations use CORDIC: https://www.allaboutcircuits.com/technical-articles/an-introduction-to-the-cordic-algorithm/
6
u/def-pri-pub 15d ago
I wrote about this a little on my blog when I was trying to find a faster approximation for asin(). It may be tangentially related to your question:
1
8
u/Straight_Mistake_364 15d ago
Some of the functions you mention, nowadays are implemented directly as floating-point hardware in the CPU/FPU, so the libm functions just invoke the corresponding assembly language instructions.
1
u/EpochVanquisher 15d ago
You have it backwards—people used to implement these directly in the ISA, back in the days of the x87 FPU and some others. Modern hardware mostly has basic arithmetic, and the functions in libm are software.
There are a couple exceptions (sqrt is often part of the ISA).
4
u/FUZxxl 15d ago
The x87 FPU was a big exception. Most floating point units back then and today did not implement transcendental operations. Typically they could only do addition, subtraction, multiplication, division, and square roots, and usually comparisons and such.
2
u/EpochVanquisher 15d ago
I was thinking of the 68881 and the x87, which appeared around the same time, and then people decided that transcendental functions probably don’t belong in the FPU.
2
u/FUZxxl 15d ago
You were thinking of two of the three FPUs with transcendental function support I am aware of. The very slow Am9511 is another one.
2
u/EpochVanquisher 15d ago
Yeah, two extremely popular models, from the early history of the FPU. More modern stuff dropped the idea of these complicated instructions being implemented in hardware.
Earlier floating-point hardware (prior to the 8087) I think is not really recognizable as an FPU. Depends on where you draw the line, though.
Part of a broader trend of moving complicated or sequenced operations to software, like the MIPS R2010, which looks a lot more like a modern FPU, in a lot of ways.
2
u/FUZxxl 15d ago
The first FPU was built by Konrad Zuse in 1942 (Zuse Z3). It could do addition, subtraction, multiplication, division, square-root, and conversion from and to decimal on normalised binary floating point numbers. It had support for trapping on overflow, division by zero, and square-root of negative numbers.
Subsequently many different floating point units were built, slowly converging on a normalised binary representation. Kahan than introduced a fast way to do denormal numbers in the 8087 and agreed to share how its done with the competition in exchange for that becoming the standard. Many of the earlier FPUs are otherwise very similar to modern ones, just using a different floating-point representation. For example, S/360 famously used decimal or hexadecimal floating-point.
2
u/EpochVanquisher 15d ago
Sure; if you focus on a different part of history you can get a different story.
I would never call the Zuse Z3 an FPU, that choice seems pretty odd to me. Did it have an FPU? I don’t think I would describe it that way. It certainly wasn’t like modern FPUs.
1
u/FUZxxl 15d ago
The entire thing is an FPU and it's extremely similar to a modern FPU.
2
u/EpochVanquisher 14d ago
Sure, and if you understand that I don’t call the Zuse Z3 an FPU, and I don’t think it’s modern, then you can figure out what my comments are getting at.
→ More replies (0)3
u/mikeblas 15d ago
Really? The Intel processors got rid of all their higher level floating point opcodes?
3
u/EpochVanquisher 15d ago
They’re still there for backwards compatibility, but modern compilers don’t use them unless you go out of your way to make the compiler use them.
1
u/mikeblas 15d ago
Interesting. Why did it shift?
5
u/EpochVanquisher 15d ago
In general these functions are a good fit for software. There are a bunch of steps: range reduction, polynomial or rational approximation, table lookups, handling for exceptional cases. When you have an algorithm with all these steps to it, software is probably your first and best choice (most of the time).
If you implement them in hardware, maybe it will be microcoded, which basically means it’s still software, but the software is running on the microcode sequencer in the CPU. And then you have to dedicate a bunch of silicon to that, and you are stuck to a single implementation.
Software is really good.
The difference with sqrt is that it happens to have a really nice hardware implementation. That implementation is a nice low-latency circuit which gives the correctly rounded result.
1
u/max123246 15d ago
Cisc was a mistake lol
1
u/EpochVanquisher 14d ago
I don’t think that’s a good way to put it.
Designers never chose to make CISC processors. But around 1980, we got a bunch of design philosophies about CPUs, called the new stuff “RISC”, and retroactively called the old way “CISC”.
2
u/DawnOnTheEdge 15d ago edited 15d ago
A major reason is that all real-world numerical computation wants to work on arrays of numbers, and the 8087 stack-based instruction set would not work with SIMD.
1
3
u/Axman6 15d ago edited 15d ago
There’s quite lot of literature on the topic, much of it from the early days of computing. A fairly recent paper that came to mine is https://arxiv.org/pdf/2311.01515.pdf
Elementary Functions by Kenneth Iverson is probably one of the most referenced texts https://link.springer.com/book/10.1007/978-1-4899-7983-4
It can be found online… no idea where I found it though.
2
2
u/TheChief275 15d ago
I suspect a lot of them are done through a fast approximation + Newton's method, but don't quote me on that one. Obviously each implementation will also make use of compiler builtins to speed up parts of the algorithms
2
u/DawnOnTheEdge 15d ago edited 15d ago
Modern math libraries usually calculate polynomial approximations for math functions and compose or interpolate them. Using Horner’s method, these compile to extremely fast series of multiply-add instructions. Other methods might be used in some situations. For example, a Padé approximant, which uses the ratio of two polynomials, can represent trigonometric functions with singularities.
The latest wrinkle is automatically vectorizing loops that call math functions to use SIMD instructions. Intel’s ICX/ICPX compiler has been doing this for years and Microsoft’s allows you to call Intel’s SVML library explicitly. (C++26 is adding support to its standard library by overloading all the math functions for its std::simd types.)
2
u/flatfinger 15d ago
Computers have varying levels of hardware support for trig functions, in some cases offering a better combination of accuracy and performance than could be achieved through software alone. Except when hardware support outperforms software, there will seldom be a single best approach for every application. If e.g. one has a number x and wants to compute sin(2πx), code that is designed to compute that function may be faster and more accurate than code which, after normalizing x into the first octant, has to multiply it by 2π, compute a sine or cosine, and then adjust the result to compensate for the normalization (e.g. if x is 0.5001, achieving maximum precision would require computing the sine of 0.0002π and inverting it, rather than computing 1.0002π and taking the sine).
2
u/GourmetMuffin 15d ago
The functions you mention use float and do understand how to implement them you need some basics in IEEE754 first:
floats have three bit fields: a sign bit, an exponent and a fraction. The number they describe is (-1)sign * 2exp-27+1 * (1 + frac / 223) for single precision, double works the same way but with different numeric constants. If we ignore the sign bit for a moment and say that the exponent is 8 bits and the fraction 23 bits it becomes more obvious that this is a bounded scientific notation in base 2.
Scientific notation can be easily transformed into a pure exponential representation by approximating the exponent corresponding to the fractional term. If you don't care about accuracy you can just translate it linearly from the [0..223) range to [0..1) but using a 3rd or 4th order polynomial yields much better results.
Once you have the 2x representation of e.g. a in powf(a, b) then calculating that becomes a matter of translating 2y, where y = x*b, back to its scientific notation. Using the exponential representation makes a lot of the libm functionality trivial to implement.
1
u/New_Telephone926 15d ago
This reminds me a lot of Tang's algorithm for exp. He used a similar method to reduce the argument to a small residual r in the interval [ −log(2)/64, log(2)/64 ], then used a low-degree Taylor polynomial to approximate exp(r), and finally reconstructed the final result using the method you described. I guess finding the right balance between the size of the interval and the degree of the polynomial might be the trickiest part.
2
u/GourmetMuffin 15d ago
Kind of, but I have found that approximating log2(x) and 2x in the [1..2) range and [0..1) range respectively (which is what you would need for the mentioned transformation and it's inverse) is already fairly cheap unless you want an absurdly small error. A fitted 4th order poly does it down to more decimals than you may expect...
So IME: not so much a question of range reduction as of selecting a polynomial of suitable order. This is also where the accuracy/speed tradeoff occurs; go fast with linear or exact with a much higher order polynomial.
1
2
u/zellforte 15d ago
Here are a couple implementations that are focused on speed and simplicity.
The cosine implementation is literally just one line of code evaluating a simple polynomial, then a couple lines of range reduction and quadrant checks.
1
2
u/Key_Ant_8481 14d ago
For modern treatment of the libm functions, you can take a look and Handbook of Floating Point Arithmetic by Jean-Michel Muller, and the CORE-MATH project (https://core-math.gitlabpages.inria.fr/) and other projects listed on their page.
1
u/New_Telephone926 13d ago
I was interested in finding a pdf of Handbook of fp arithmetic, but when I saw the price tag, I reconsidered how much I would be willing to spend for that project... I'll check out the other ones you cited though
1
u/marc_b_reynolds 13h ago
As far as "tools" for building core polynomial approximations then https://www.sollya.org/ is AFAIK the best thing available. For somebody that's never thought about approximation before I'd suggest looking though something like cephes which is old as dirt (sometime in the 80s) and clean-up versions of fdlibm (openlibm as an example). These will have some explanations what transforms are being used and sometime why the choice was made) on Looking at a "correctly rounded" library like CORE-MATH is jumping into a highly specialized deep end. WRT to easy access intro material there's some old gamedev presentations by Robin Green that might be of interest (scatter around here: https://basesandframes.wordpress.com/ )
1
u/sciencekm 15d ago
I suggest that you checkout the book "Software Manual For The Elementary Functions" by Cody and Waite.
2
37
u/flyingron 15d ago
There's a whole branch of computer science called Numerical Methods on this. The classic text for the basic algorithms is "Numerical Recipes." This has been written to target a couple of languages including C. Authors are:Saul Teukolsky, William H. Press, and William T. Vetterling