Math libraries

If there’s one thing we associate with computing more than the concept of Turing Completeness (a.k.a. The Church-Turing Thesis) it is mathematical computation! CPUs have dedicated hardware to make this rapid & memory-efficient for some finite number of base2-digits (a.k.a. bits). This page documents facilities refining & exposing digital computers’ arithmatic expertise.

bc Infix Maths

Sometimes you want to perform math on the commandline for which there’s the bc & dc commands, bc being discussed here.

After initiliazing its I/O bc parses its commandline flags & trailing args both given directly & via BC_ENV_ARGS envvar. POSIXLY_CORRECT & BC_LINE_LENGTH are also checked. Some “collection” (including some common numbers & “load” func), nametree, & labels, are initialized.

Uses backwards bug-compatible emulated floating point.


Then bc registers a SIGINT signal handler, opens the first sourcefile or standardin, initializes possibly LibEdit or LibReadline if available & in REPL mode, runs a Bison/Flex immediately-evaluating (conditionally via a generate or run_code function) parser consulting those collections, outputs a newline if that was silenced, & successfully exits.

The execution is a straightforward opcode stackmachine interpretor, for the sake of function calls. Has a standard prelude.

Flex is configured to concatenate multiple files & possibly read from libReadline.

P.S. I gather I’d find the older implementation more interesting…

dc Postfix Maths

Above I described the bc command, so here I’ll study the corresponding dc command. Both evaluating math.

dc starts by initializers some common emulated numbers, virtual registers, & stubs for strings & arrays before parsing commandline flags & iterating over args carefully opened as files or stdin. For each it lexes (via a handwritten lexer) the file treating the tokens as immediately-evaluated opcodes. With seperate interpretors for strings vs numbers.

Outputs flushed on exit.

There’s support submodules in dc for interpreting array, “system” & I/O, numbers (wrapping bc’s implementation), stack (which gets incorporated into all the opcodes) & registers, & strings. Amongst a couple core utilities.

LibGNU MultiPrecision

Sometimes you want to avoid bugs caused by the finite number of bits a CPU-native number consists of. GNU provides a handful of libraries for dealing with this, & today I’ll start studying GNU MultiPrecision (LibGMP).

Amongst the main codebase (I’ll go over the subsystems starting tomorrow), I see:

Some of those macros wraps Assembly code.

MiniGMP

Within GMP there’s a submodule named “Mini GMP” which initializes, normalizes, clears, conanicalizes, set fields, compare, & finalizes mpz_t & mpq_t structures. And ofcourse there’s mathematical operators. These are all straightforward with possibly a couple branches.

You can swap fields between 2 mpq_ts.

Convert to or from doubles involves a loop taking an exponent of a power-of-2. Or it can convert to/from strings.

The actual logic with most of the loops is implemented as macros.


There’s another implementation file containing macros implemented similarly. As well as utilities for crashing, managing memory & managing those callbacks, memory-copying & -comparison loops, loop to compute normalized sizes. And ofcourse loops over the words in multiprecision number(s) to perform mathematical operators on them whilst propagating carries, Mini GMP implements these in C.

Divide is more complex involving multiple fastpath routines. Which are involved in stringifying.


There’s routines for computing greatest-common-divisor on single words or multiprecision numbers through repeated refinement. With divide-by-2 fastpaths. Similarly there’s Lowest-Common-Multiple routines.

There’s exponentiation routines which may involve repeated multiplication, or divide&conquer algorithms via division.

There’s a binary-search for squareroots with a fastpath.

There’s checks for perfect squares.

There’s repeated multiplication factorials


There’s repeated bitwise-operators computing Jacobi numbers.

Repeated addition/subtraction/multiplication/division Lucas Modulo operation, supporting its own stringifier.

Repeated exponentiation labelled “Miller-Rabin”.

Relatively simple prime-number testers.

Bittesting, absolutes, bit-subtraction, bitsetting, bitclearing, & more logical operators. Hamming distance via ops on each bit. Array-scanning. Sizing per base. Stringifying abstractions.

LibMPF

GNU MultiPrecision includes a “MPF” sublibrary”, which I’ll summarize today.

This operates upon an mpf_t structure with the following utilities:

LibMPN

GNU MultiPrecision includes a LibMPN subsystem implemented in Assembly for a wide variety of CPUs. Normally I’d discuss the x86_64 implementation since that’s what I’m running, but instead I’ll study the ARM64 implementation since that’s simpler!

This includes routines for:

There’s dataheaders determining which opcodes to use on Cora53, Cora57, Cora72, Cora73, XGene1, & other ARM CPUs.

On Cora53 it LibMPN implements special variations of it’s comparison routines.

LibMPQ

GNU Multiprecision has a LibMPQ sublibrary implementing fractions, consisting of several mostly-trivial utilities around a struct holding a numerator & denominator as MPN & MPZ numbers. These utilities include:

0/d is normalized to 0/1.


Also I see utilities for:

LibMPZ

The sublibrary around the MPZ numeric datastructure includes:

LibMPQ

Today I’ll continue studying MPQ, which largely serves to wrap various MPN routines with memory management & normalization. These include:

Also there’s converters:

As well as:

Regarding prime numbers:

LibMPZ

Reading through the rest of GNU Multiprecisions MPZ sublibrary, I see:

Text serialization

Here I’m looking over GNU MultiPrecision’s textual serialization library.

This has several public API wrappers around the core logic for different function signatures. And methodtables of underlying utilities, specified by those public APIs.

There’s a replacement printf implementation repeatedly looking for “%” & examining subsequent character(s) calling the core logic where appropriate.

The core logic works similarly parsing the core logic into a “param” struct & calling the appropriate serialization function, some of which are defined here. With macros handling the serialization & functions handling the formatting.

Random-number generation

GNU MultiPrecision includes a random-number-generating sublibrary.

This includes a handful of wrapper functions converting the seed value from different types. Or calling methods. There’s some global variables.

There’s a routine which repeatedly generates numbers until its less than n, avoiding infinite loops. With a bitmask to improve odds.

There’s an abstraction around selecting a randomness-formula, or another to set it.


There’s an implementation of the Mersenne Twister PRNG formula upon MPZ. With a seperate methodtable for seedable MT-random. Includes a default-state lookuptable & plenty of bittwiddling I don’t understand.

Or there’s a more traditional Linear-Congruential random-number generator implemented upon MPN. Again, plenty of bittwiddling!

Includes a variant dealing in base2-exponentials.

Text Deserialization

GNU MultiPrecision includes a text deserialization sublibrary. This works basically the same way as the serialization sublibrary, if not simpler.

Handling public APIs & methodtables similarly.

There’s a core number scanner deferring to existing parsers (from MPF, MPQ, or MPZ subparsers), to be wrapped in a scanf replacement.

LibGMP Tuning

GNU Multiprecision cares a lot about performance (and correctness). As such there’s a subsystem which uses C macros to tweak various APIs to capture performance metrics. With the core “cycle-counting” code implemented in Assembly for a variety of CPUs to avoid incurring any overhead & generating misleading numbers.

There’s some noops implemented, presumably to test this infrastructure.

There’s utilities implemented around this to eg. try out some magic numbers to find the most optimal one.


There’s also utility commands to perform statistical analysis on the results & render charts via GNU PlotLib.

There’s some shared utility libraries for these commands, including one to tweak CPU settings. And a benchmark suite, with its own testrunner.

There’s perl-script to generate profiled versions of the MPN library.

And I’m not even going to dig through the testsuite…

##LibMPFR floatingpoint arithmatic LibMPFR implements floatingpoint numbers to a configurable amount of precision similar to how the Floating-Point-Unit in your CPU works. In short mostly this implements various mathematical formulas I barely understand often with a refinement loop & rounding. Whilst adding special-cases for NaN or infinity, frequently fastpaths, & often wrapper functions. Ocassionally lookuptables, some of which are CPU-specific, are consulted for some of these functions. The most complex logic in LibMPFR is involved in parsing & serializing strings.


In long, here’s my commentary as I skimmed through the entire toolbox:






LibMPC

LibMPC is a surprisingly simple wrapper around LibMPFR (which I described the other day) implementing multiprecision complex numbers.

Functions wrapper MPFR’s on a single property includes:

There’s a constant version number. There’s headerfiles ofcourse.

There’s memory-management utilities.


MPFR wrappers upon both properties include:


Going from least complex to most of the remaining APIs: