spire.math

Algebraic provides an exact number type for algebraic numbers. Algebraic numbers are roots of polynomials with rational coefficients. With it, we can represent expressions involving addition, multiplication, division, n-roots (eg. sqrt or cbrt), and roots of rational polynomials. So, it is similar Rational, but adds roots as a valid, exact operation. The cost is that this will not be as fast as Rational for many operations.

In general, you can assume all operations on this number type are exact, except for those that explicitly construct approximations to an Algebraic number, such as toBigDecimal.

For an overview of the ideas, algorithms, and proofs of this number type, you can read the following papers:

• "On Guaranteed Accuracy Computation." C. K. Yap.
• "Recent Progress in Exact Geometric Computation." C. Li, S. Pion, and C. K. Yap.
• "A New Constructive Root Bound for Algebraic Expressions" by C. Li & C. K. Yap.
• "A Separation Bound for Real Algebraic Expressions." C. Burnikel, et al.

Abstract class that can be used to implement custom binary merges with e.g. special collision behavior or an ordering that is not defined via an Order[T] typeclass

Merge that uses binary search to reduce the number of comparisons

This can be orders of magnitude quicker than a linear merge for types that have a relatively expensive comparison operation (e.g. Rational, BigInt, tuples) and will not be much slower than linear merge even in the worst case for types that have a very fast comparison (e.g. Int)

Complex numbers. Depending on the underlying scalar T, can represent the Gaussian integers (T = BigInt/SafeLong), the Gaussian rationals (T = Rational) or the complex number field (T: Field).

Note that we require T to be at least CRing, a commutative ring, so the implementation below is slightly less general than the Cayley-Dickson construction.

FastComplex is an ugly, beautiful hack.

The basic idea is to encode two 32-bit Floats into a single 64-bit Long. The lower-32 bits are the "real" Float and the upper-32 are the "imaginary" Float.

Since we're overloading the meaning of Long, all the operations have to be defined on the FastComplex object, meaning the syntax for using this is a bit ugly. To add to the ugly beauty of the whole thing I could imagine defining implicit operators on Long like +@, -@, *@, /@, etc.

You might wonder why it's even worth doing this. The answer is that when you need to allocate an array of e.g. 10-20 million complex numbers, the GC overhead of using any object is HUGE. Since we can't build our own "pass-by-value" types on the JVM we are stuck doing an encoding like this.

Here are some profiling numbers for summing an array of complex numbers, timed against a concrete case class implementation using Float (in ms):

size | encoded |  class
 1M |     5.1 |    5.8
 5M |    28.5 |   91.7
 10M |    67.7 |  828.1
 20M |   228.0 | 2687.0

Not bad, eh?

Value class which encodes two floating point values in a Long.

We get (basically) unboxed complex numbers using this hack. The underlying implementation lives in the FastComplex object.

An implementation of insertion sort.

Works well for small arrays but due to quadratic complexity is not generally optimal.

Integral number types, where / is truncated division.

Interval represents a set of values, usually numbers.

Intervals have upper and lower bounds. Each bound can be one of four kinds:

• Closed: The boundary value is included in the interval. * Open: The boundary value is excluded from the interval. * Unbound: There is no boundary value. * Empty: The interval itself is empty.

When the underlying type of the interval supports it, intervals may be used in arithmetic. There are several possible interpretations of interval arithmetic: the interval can represent uncertainty about a single value (for instance, a quantity +/- tolerance in engineering) or it can represent all values in the interval simultaneously. In this implementation we have chosen to use the probabillistic interpretation.

One common pitfall with interval arithmetic is that many familiar algebraic relations do not hold. For instance, given two intervals a and b:

a == b does not imply a * a == a * b

Consider a = b = [-1, 1]. Since any number times itself is non-negative, a * a = [0, 1]. However, a * b = [-1, 1], since we may actually have a=1 and b=-1.

These situations will result in loss of precision (in the form of wider intervals). The result is not wrong per se, but less accurate than it could be.

These intervals should not be used with floating point bounds, as proper rounding is not implemented. Generally, the JVM is not an easy platform to perform robust arithmetic, as the IEEE 754 rounding modes cannot be set.

==Overview== A simple implementation of N-dimensional dual numbers, for automatically computing exact derivatives of functions. This code (and documentation) closely follow the one in Google's "Ceres" library of non-linear least-squares solvers (see Sameer Agarwal, Keir Mierle, and others: Ceres Solver.)

While a complete treatment of the mechanics of automatic differentiation is beyond the scope of this header (see http://en.wikipedia.org/wiki/Automatic_differentiation for details), the basic idea is to extend normal arithmetic with an extra element "h" such that h != 0, but h^2^ = 0. Dual numbers are extensions of the real numbers analogous to complex numbers: whereas complex numbers augment the reals by introducing an imaginary unit i such that i^2^ = -1, dual numbers introduce an "infinitesimal" unit h such that h^2^ = 0. Analogously to a complex number c = x + yi, a dual number d = x * yh has two components: the "real" component x, and an "infinitesimal" component y. Surprisingly, this leads to a convenient method for computing exact derivatives without needing to manipulate complicated symbolic expressions.

For example, consider the function

 f(x) = x * x ,

evaluated at 10. Using normal arithmetic,

f(10) = 100, and df/dx(10) = 20.

Next, augment 10 with an infinitesimal h to get:

 f(10 + h) = (10 + h) * (10 + h)
           = 100 + 2 * 10 * h + h * h
           = 100 + 20 * h       +---
                   +-----       |
                   |            +--- This is zero
                   |
                   +----------------- This is df/dx

Note that the derivative of f with respect to x is simply the infinitesimal component of the value of f(x + h). So, in order to take the derivative of any function, it is only necessary to replace the numeric "object" used in the function with one extended with infinitesimals. The class Jet, defined in this header, is one such example of this, where substitution is done with generics.

To handle derivatives of functions taking multiple arguments, different infinitesimals are used, one for each variable to take the derivative of. For example, consider a scalar function of two scalar parameters x and y:

 f(x, y) = x * x + x * y

Following the technique above, to compute the derivatives df/dx and df/dy for f(1, 3) involves doing two evaluations of f, the first time replacing x with x + h, the second time replacing y with y + h.

For df/dx:

 f(1 + h, y) = (1 + h) * (1 + h) + (1 + h) * 3
             = 1 + 2 * h + 3 + 3 * h
             = 4 + 5 * h

 Therefore df/dx = 5

For df/dy:

 f(1, 3 + h) = 1 * 1 + 1 * (3 + h)
             = 1 + 3 + h
             = 4 + h

 Therefore df/dy = 1

To take the gradient of f with the implementation of dual numbers ("jets") in this file, it is necessary to create a single jet type which has components for the derivative in x and y, and pass them to a routine computing function f. It is convenient to use a generic version of f, that can be called also with non-jet numbers for standard evaluation:

 def f[@specialized(Double) T : Field](x: T, y: T): T = x * x + x * y

 val xValue = 9.47892774
 val yValue = 0.287740

 // The "2" means there should be 2 dual number components.
 implicit val dimension = JetDim(2)
 val x: Jet[Double] = xValue + Jet.h[Double](0);  // Pick the 0th dual number for x.
 val y: Jet[Double] = yValue + Jet.h[Double](1);  // Pick the 1th dual number for y.

 val z: Jet[Double] = f(x, y);
 println("df/dx = " + z.infinitesimal(0) + ", df/dy = " + z.infinitesimal(1));

For the more mathematically inclined, this file implements first-order "jets". A 1st order jet is an element of the ring

 T[N] = T[t_1, ..., t_N] / (t_1, ..., t_N)^2

which essentially means that each jet consists of a "scalar" value 'a' from T and a 1st order perturbation vector 'v' of length N:

 x = a + \sum_i v[i] t_i

A shorthand is to write an element as x = a + u, where u is the perturbation. Then, the main point about the arithmetic of jets is that the product of perturbations is zero:

 (a + u) * (b + v) = ab + av + bu + uv
                   = ab + (av + bu) + 0

which is what operator* implements below. Addition is simpler:

 (a + u) + (b + v) = (a + b) + (u + v).

The only remaining question is how to evaluate the function of a jet, for which we use the chain rule:

 f(a + u) = f(a) + f'(a) u

where f'(a) is the (scalar) derivative of f at a.

By pushing these things through generics, we can write routines that at same time evaluate mathematical functions and compute their derivatives through automatic differentiation.

Used to implicitly define the dimensionality of the Jet space.

Simple linear merge

Interface for a merging strategy object.

In-place merge sort implementation. This sort is stable but does mutate the given array. It is an in-place sort but it does allocate a temporary array of the same size as the input. It uses InsertionSort for sorting very small arrays.

Convenient apply and implicits for Numbers

A NumberTag provides information about important implementations details of numbers. For instance, it includes information about whether we can expect arithmetic to overflow or produce invalid values, the bounds of the number if they exist, whether it is an approximate or exact number type, etc.

TODO

3. LiteralOps? Literal conversions?
4. Review operator symbols?
5. Support for more operators?
6. Start to worry about things like e.g. pow(BigInt, BigInt)

Polynomial A univariate polynomial class and EuclideanRing extension trait for arithmetic operations. Polynomials can be instantiated using any type C for which a Ring[C] and Eq[C] are in scope, with exponents given by Int values. Some operations require more precise algebraic structures, such as GCDRing, EuclideanRing or Field to be in scope.

Quaternions defined over a subset A of the real numbers.

In-place quicksort implementation. It is not stable, but does not allocate extra space (other than stack). Like MergeSort, it uses InsertionSort for sorting very small arrays.

Provides a type to do safe long arithmetic. This type will never overflow, but rather convert the underlying long to a BigInteger as need and back down to a Long when possible.

Given a function for finding approximate medians, this will create an exact median finder.

Interface for a sorting strategy object.

Object providing in-place sorting capability for arrays.

Sorting.sort() uses quickSort() by default (in-place, not stable, generally fastest but might hit bad cases where it is quadratic. Also provides mergeSort() (in-place, stable, uses extra memory, still pretty fast) and insertionSort(), which is slow except for small arrays.

Implementation of three-valued logic.

This type resembles Boolean, but has three values instead of two:

• Trilean.True (equivalent to true)
• Trilean.False (equivalent to false)
• Trilean.Unknown

Trilean supports the same operations that Boolean does, and as long as all values are True or False, the results will be the same. However, the truth tables have to be extended to work with unknown:

not:

-+-
T|F
U|U
F|T

and:

 |T U F
-+-----
T|T U F
U|U U F
F|F F F

or:

 |T U F
-+-----
T|T T T
U|T U U
F|T U F

Trilean is implemented as a value type, so in most cases it will only have the overhead of a single Int. However, in some situations it will be boxed.

abs

ceil

choose (binomial coefficient)

e

Integer Euclidean division, equotmod, equot, emod

exp

factorial

fibonacci

floor

gcd

lcm

log

max

min

An implementation of the shifting n-th root algorithm for BigDecimal. For the BigDecimal a, this is guaranteed to be accurate up to the precision specified in ctxt.

See http://en.wikipedia.org/wiki/Shifting_nth_root_algorithm

pi

pow

Exponentiation function, e.g. x^y

If base^ex doesn't fit in a Long, the result will overflow (unlike Math.pow which will return +/- Infinity).

round

signum

sqrt