matrico’s Floating-Point Addon-Module fpmath

To extend the flonum functionality of CHICKEN Scheme and particularly the (chicken flonum) module, an addon-module named fpmath is implemented, that supplies various flonum specific functions for matrico specific, but also general, use. And the (chicken flonum) naming of prefixing functions with fp is utilized.

As the fpmath module is not very useful without the imported (chicken flonum) module, it is reexported. And as the fpmath module is also practically useful for high-level use in the overarching matrico module, matrico in turn reexports fpmath. You can take a look at the src/fpmath.scm source file here.

Converters

To bridge between the exact fixnum numbers and the inexact flonum numbers, converters are needed. R5RS provides a fixnum to flonum conversion, but its name inexact->exact is too unwieldy, thus a shorter alias fp is introduced. Another very useful converter is, given two fixnums representing numerator and denominator, a conversion to a flonum from such rational components, which is used below for constant provisioning.

  • (fp n) - Alias for exact->inexact, returns inexact flonum from exact fixnum argument.
  • (fp% n d) - Returns flonum from numerator fixnum and denominator fixnum arguments.

Predicates

A tricky issue is comparing flonums for equality, due to their inexactness. However, any equality comparison boils down to a comparison against zero, in example via the difference of two flonums. To handle this, two functions are implemented, one to compare against exact zero, which is emphasized by two question-mark predicate indicators, and one to compare against a tolerance:

  • (fpzero?? x) - Returns true for zero, else false.
  • (fpzero? x tol) - Returns true if first argument is less than second argument, else false.

Operators

The (chicken flonum) module only provides one unary elementary arithmetic operations, namely fpneg, the unary shortcut for (fp- 0.0 x). For division, addition, and multiplication no unary operations are defined. These are added here: In terms of division, the multiplicative inverse (fp/ 1.0 x) is defined, for addition the double, and for multiplication the square of the argument.

Also, an operation central to numerical computation, for example in the dot-product, and thus matrix multiplication, matrix decompositions, like the QR decomposition, or integrators such as Runge-Kutta methods, is included. This three operand operation of a combined flonum multiplication and addition is classically known as a Flop (short for FLoating-point OPeration) or “axpy” in BLAS, and technically as FMA (Fused Multiply-Add) in CPUs, or MAD (Multiply-Add) in GPUs. Also, this operation has improved accuracy, as only the final result is rounded to the closest flonum. This operation is important enough to be provided in the C library, and also in the flonum modules of various Schemes such as Gauche and MIT Scheme.

  • (fp*2 x) - Returns the double of the flonum argument.
  • (fp^2 x) - Returns the square of the flonum argument.
  • (fprec x) - Returns the reciprocal of the flonum argument.
  • (fp*+ x y z) - Returns the fused multiply-add of the flonum arguments.

In future versions of CHICKEN Scheme (maybe version 5.4), (chicken flonum) will include the fp*+ function, hence it is wrapped in a CHICKEN version dependent cond-expand.

Constants

The (chicken flonum) module does not provide any of the typical mathematical constants. So, the fpmath module fills this gap by defining the three most important real-valued constants (in my opinion, outside of: 0, 1, √2) as thunks (zero argument functions):

  • A circle constant: tau (exclusively), and (particularly) not pi.
  • Euler’s number: e
  • Golden ratio: phi

I am using thunks, because I am not storing the constants’ value to precision in-source, but computing them as the module loads. This could be done, in example, for tau via (8.0 * (atan 1.0)), for e via (exp 1.0), or for phi (* (+ 1.0 (sqrt 5.0)) 0.5). However, there are two issues with this: First, how would I know that these are correctly computed, so the functions are correct and the formulas are not suffering from numerical effects like annihilation. And second, I would need to call functions whose actual complexity I don’t know. Also, I want to use only a single elementary arithmetic function to approximate these constants. Since floating-point numbers are of finite precision, I will use rational approximation (like 355/113 for pi) and thus a mere single division operation. Luckily, this is facilitated by the “On-Line Encyclopedia of Integer Sequences” which provides numerators, denominators, and number of correct decimal digits, such that a rational approximation with accuracy matching flonums can be selected, and tested (verified) against the closed forms mentioned above:

I have not seen use of rational approximations in the context of approximate floating-point constants in numerics before, so maybe somebody has some pointers to other projects using this? Anyway, I think this is a nice exploitation of flonum inexactness.

Generalized Functions

Beyond the absolute value function fpabs in (chicken flonum) there are some further “generalized functions” (in simple terms, these mathematical mappings are generalized because they require conditionals, due to their piecewise defintion) I consider useful:

  • (fpdelta x) - Returns zero except for argument zero it returns one, see Delta.
  • (fpheaviside x) - Returns one for positive arguments zero otherwise, see Heaviside.
  • (fpsign x) - Returns minus one for negative, one for positive, and zero for zero arguments, see Sign.

Logarithms

The (chicken flonum) module only provides a natural logarithm fplog, which is the minimal necessary, but not too comfortable. So I am adding typical calculator functionality (in their naming), based on the logarithm rules:

  • (fpln x) - Alias for fplog, returns natural logarothm of flonum argument.
  • (fplb x) - Returns binary logarithm (base 2) of flonum argument.
  • (fplg x) - Returns decimal logarithm (base 10) of flonum argument.
  • (fplogb b x) - Returns logarithm to positive flonum base argument of flonum argument.

Note, that the utilized values ln(2) ((fplog 2.0)) and ln(10) ((fplog 10.0)) are neither hard-coded as constants nor computed during a function call to fplb or fplg, but pre-computed once on load of the module.

Hyperbolic Functions and Inverse Hyperbolic Functions

As I stated in an earlier post, a shortcoming of R5RS and also (chicken flonum) is the absence of hyperbolic functions and inverse hyperbolic functions. So, these are provided in this addon module via the exponential definitions of the hyperbolic functions and the logarithmic definitions of the inverse hyperbolic functions:

  • (fpsinh x) - Returns hyperbolic sine for flonum argument.
  • (fpcosh x) - Returns hyperbolic cosine for flonum argument.
  • (fptanh x) - Returns hyperbolic tangent for flonum argument.
  • (fpasinh x) - Returns area hyperbolic sine (inverse hyperbolic sine) for flonum argument.
  • (fpacosh x) - Returns area hyperbolic cosine (inverse hyperbolic cosine) for flonum argument.
  • (fpatanh x) - Returns area hyperbolic tangent (inverse hyperbolic tangent) flor flonum argument.

Haversed Trigonometric Functions

A lesser used variant of the trigonometric functions are the haversed sine (haversine) and haversed cosine (havercosine). These scaled and shifted sine and cosine are sometimes useful in systems and control theory, my field of mathematical research.

  • (fphsin x) - Returns haversine for flonum argument.
  • (fphcos x) - Returns havercosine for flonum argument.

Logarithmic Hyperbolic Functions

Another pair of lesser known functions are logarithmic hyperbolic functions, which are implemented in the GNU Scientific Library. While log-sinh is a differentiable function over the positive real numbers that behaves like a logarithm for small (positive) values and like a linear function for large (positive) values, the log-cosh is a differentiable function over all real numbers, that behaves like a quadratic function for small values and like the absolute value for large values.

  • (fplnsinh x) - Returns log-sinh for flonum argument.
  • (fplncosh x) - Returns log-cosh for flonum argument.

Special Functions

Some useful extra functions are implemented:

  • (fpsignsqrt x) - Returns the square root of the absolute value of the flonum argument times the sign of the argument.
  • (fpsinc x) - Returns the (unnormalized) cardinale sine of the flonum argument.
  • (fpsigm x) - Returns the sigmoid (standard logistic function) of the flonum argument.
  • (fpgauss x) - Returns the standard Gauss function (zero expectation, unit deviation) of the flonum argument.
  • (fpstirling x) - Returns the Stirling’s approximation for factorials of the flonum argument.

Note, that fpsinc can be normalized by scaling its argument, fpsigm and fpgauss can be generalized by shifting and scaling its argument and also scaling its result, while fpstirling can be tested with the utils module’s factorial function.

Utilities

To print matrices on the terminal a fixed element print width is necessary. Given the flonum-print-precision from (chicken flonum), this is only workable with a precision of 17, which in turn would mean on a standard 80x25 terminal window that only a matrix with 4 columns could be printed without additional breaks. Now, with a maximum width of 8 already a 9-column matrix becomes printable. This tapered flonum to string conversion is implemented by:

  • (fptaper x) - Returns 8-character string approximately describing flonum argument.

Practically this means the integer and fractional parts can occupy maximally 6 characters, as one character is reserved for the sign and another for the decimal point. So the smallest fully printable number is ±0.00001, where as the largest is ±99999.9. Any smaller number is partially printed as ±0.000… and larger numbers are printed as ±100000. or ±100000… respectively. To highlight an exact zero, it is printed as ____0___. Special values are printed as: ___±∞___ and ___NaN__.

However, as this function is supposed only for internal use, it is excluded from reexport in matrico.