Table of contents

• Overview
  • Examples Modules
• Angle
• Arrays
• Combinatorics
• Complex
• Core
  • Precision
  • DD: double-double
• Field
• Fraction
  • Rational Numbers
  • Continued Fractions
• Gamma
• Primes
• Quaternion
• Root Finder
• Dependencies

Examples Modules

In addition to the modules above, the Commons Numbers source distribution contains example code demonstrating library functionality and/or providing useful development utilities. These modules are not part of the public API of the library and no guarantees are made concerning backwards compatibility. The example module parent page contains a listing of available modules.

Precision

The Precision class provides comparison of floating-point numbers using relative, absolute and unit of least precision (ULP) tolerances. Comparison using the relative tolerance \( \epsilon \) is symmetric as the relative delta is computed using the largest magnitude of the input pair of numbers:

\( \frac{|a - b|}{\operatorname{max}(|a|, |b|)} \le \epsilon \)

Comparisons using the ULP argument specifies the allowed distance between two numbers. Numbers are equal if they have (ulp - 1) representable numbers or fewer between a and b. Use ulp = 0 for binary equality.

Note: The default implementation for equality allows no numbers between a and b. It is equivalent to Precision.equals(a, b, 1). Thus using a floating-point delta of 0.0 may indicate numbers are equal when they are not binary equal (see example below).

// Default allows no numbers between
Precision.equals(1000.0, 1000.0)                              // true
Precision.equals(1000.0, 1000.0 + Math.ulp(1000.0))           // true
Precision.equals(1000.0, 1000.0 + 2 * Math.ulp(1000.0))       // false

// Absolute - tolerance is floating-point
Precision.equals(1000.0, 1001.0)                              // false
Precision.equals(1000.0, 1001.0, 1.0)                         // true
Precision.equals(1000.0, 1000.0 + Math.ulp(1000.0), 0.0)      // true  (but not binary equal)

// ULP - tolerance is integer
Precision.equals(1000.0, 1001.0)                              // false
Precision.equals(1000.0, 1001.0, 1)                           // false
Precision.equals(1000.0, 1000.0 + 2 * Math.ulp(1000.0), 1)    // false
Precision.equals(1000.0, 1000.0 + 2 * Math.ulp(1000.0), 2)    // true  (n - 1) numbers between
Precision.equals(1000.0, 1000.0 + 3 * Math.ulp(1000.0), 2)    // false

// Relative
Precision.equalsWithRelativeTolerance(1000.0, 1001.0, 1e-6)   // false
Precision.equalsWithRelativeTolerance(1000.0, 1001.0, 1e-3)   // true

Equality using NaN values will be false. To allow two NaN values to be equal requires the use of the method equalsIncludingNaN. Special handling of NaN values allows the default implementation to be optimized for the typical use case comparing finite numbers. Note that infinity values are equal.

Precision.equals(Double.NaN, 1000.0)                                   // false
Precision.equals(Double.NaN, Double.NaN)                               // false
Precision.equalsIncludingNaN(Double.NaN, Double.NaN)                   // true

Precision.equals(Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY)   // true
Precision.equals(Double.NEGATIVE_INFINITY, Double.NEGATIVE_INFINITY)   // true
Precision.equals(Double.POSITIVE_INFINITY, Double.NEGATIVE_INFINITY)   // false

Precision provides a compareTo method to provide ordering of numbers that are not considered equal using the tolerance:

Precision.compareTo(100, 100, 0.0)   // 0
Precision.compareTo(100, 101, 1.0)   // 0
Precision.compareTo(100, 102, 1.0)   // -1
Precision.compareTo(102, 100, 1.0)   // 1

The comparison uses the natural ordering and is symmetric with regard to NaN values. The comparison using <, ==, > is extended in the DoubleEquivalence interface that supports <=, >= and signum.

Precision.DoubleEquivalence eq = Precision.doubleEquivalenceOfEpsilon(1.0);
eq.lt(100, 100)    // false
eq.lte(100, 100)   // true
eq.eq(100, 100)    // true
eq.gte(100, 100)   // true
eq.gt(100, 100)    // false

Precision also provides methods for other floating-point operations such as rounding and computing a representable delta from a number and an ideal delta.

// Round to a number of digits to the right of the decimal point
Precision.round(678.125, 4)              // 678.125
Precision.round(678.125, 3)              // 678.125
Precision.round(678.125, 2)              // 678.13
Precision.round(678.125, 1)              // 678.1
Precision.round(678.125, 0)              // 678.0
Precision.round(678.125, -1)             // 680.0
Precision.round(678.125, -2)             // 700.0

Precision.representableDelta(1.0, 0.1)   // 0.10000000000000009

DD: double-double

The DD class provides an extended precision floating-point number. A double-double is an unevaluated sum of two IEEE double precision numbers capable of representing at least 106 bits of significand. A normalized double-double number \( (x, xx) \) satisfies the condition that the parts are non-overlapping in magnitude such that:

|x| > |xx|
x == x + xx

The implementation assumes a normalized representation during operations on a DD number and computes results as a normalized representation. The number is limited to the same exponent range as a double. Thus the number can be considered as a standard IEEE double with additional information that would be lost to rounding. This information can be used in arithmetic operations so that compound calculation can reduce the error present in the final result returned as a double.

The number \( (x, xx) \) may also be referred to using the labels high and low to indicate the magnitude of the parts as \( ( x_{hi}, x_{lo} ) \), or using a numerical suffix for the parts as \( ( x_0, x_1 ) \). The numerical suffix is typically used when the extended floating-point number has an arbitrary number of parts, for example a quad-double is \( ( x_0, x_1, x_2, x_3 ) \). The parts of the DD can be accessed using the hi() and lo() methods and the evaluated sum \( x + xx \) as doubleValue().

The DD class is immutable. Instances are constructed with factory methods. A DD can be constructed from a double, int or long primitive. In the case of a long this maintains the full 64-bits of information; this is not possible when using implicit type conversion to a double. These conversions can be reversed as the DD provides type conversions using the java.lang.Number methods.

double x = Math.PI;
int    y = 42;
long   z = -8564728970587006436L;
DD.of(x).doubleValue()  // = x
DD.of(y).intValue()     // = y
DD.of(z).longValue()    // = z
                        // (long) (double) z != z

The DD class can be created as the result of a double operation and the round-off error. In the case of addition and multiplication the result is exact. For division the result will have the closest 106-bit representation to the exact result. Note that use of the factory constructors for double arithmetic is preferred over constructing the arguments as a DD from a double and performing the equivalent operation in DD arithmetic:

double a = 1.2345678901234567;
double b = 123.45678901234567;
DD w = DD.ofProduct(a, b);     // (152.41578753238835,-1.0325951435749745E-14)
DD x = DD.ofSum(a, b);         // (124.69135690246912,-1.1102230246251565E-15)
DD y = DD.ofDifference(a, b);  // (-122.22222112222221,-1.1102230246251565E-15)
DD z = DD.fromQuotient(1, 3);  // (0.3333333333333333,1.850371707708594E-17)
// w.hi() = a * b
// x.hi() = a + b
// y.hi() = a - b
// z.hi() = 1.0 / 3

// Inefficient
DD zz = DD.of(1).divide(DD.of(3));
z.equals(zz)   // true

The DD can also be converted from and to a BigDecimal. The conversion from may lose precision. The conversion to is exact but only supported for finite numbers as BigDecimal does not support infinity and NaN. The method isFinite() can be used to verify that the evaluated sum \( x + xx \) is finite; this is only true when both parts are finite and conversion to BigDecimal is possible. Note that BigDecimal can be used for numerical equality and ranking operations as the DD class does not support java.lang.Comparable, and the equals(Object) method is true for numerical equality of the parts, not numerical equality of their evaluated sum. The following demonstrates a round trip of \( \pi \) to 50 decimal places. The DD is limited to approximately 34 decimal places and this is split into the double value for \( \pi \) and a double representing the rounding error.

BigDecimal pi = new BigDecimal("3.14159265358979323846264338327950288419716939937510");
DD x = DD.from(pi);         // (3.141592653589793,1.2246467991473532E-16)
pi.compareTo(x.bigDecimalValue())  // != 0
// x.hi() = Math.PI
// x.lo() = pi.subtract(new BigDecimal(Math.PI)).doubleValue()

DD nan = DD.of(Double.NaN);
nan.isFinite()          // false
nan.bigDecimalValue()   // throws NumberFormatException

Note: It is not possible to directly specify the two parts of the DD number. The two parts must be added using a sum to ensure the DD maintains a normalized representation. If the two parts already represent a number \( (x, xx) \) such that \( x == x + xx \) then the magnitudes of the parts will be unchanged; any signed zeros may be subject to a sign change.

The DD class provides the arithmetic operations add, subtract, multiply, and divide for a double or DD argument. int and long arguments can be used as a double due to implicit type conversion. However if the long value has more than 53-bits of precision then it should first be converted to a DD.

long x = -8564728970587006436L;
DD.ONE.add(x).longValue()         // != x + 1
DD.ONE.add(DD.of(x)).longValue()  // == x + 1

The arithmetic methods are not intended to provide an exactly rounded result to 106-bit precision. A compromise has been made between accuracy and performance so that the DD class is a viable alternative to using BigDecimal for high accuracy arithmetic when the ultimate result is required as a double. The approximate error bounds of operations are supplied in the class javadoc. Typical errors bounds are within a relative error of 1-4 eps of the exact result where eps is \( 2^{-106} \).

The following demonstrates a simple sum where the double representation of a fraction is inexact and results in a rounding error. This is removed by extending the fraction to a double-double representation:

1.0 / 2 + 1.0 / 3 + 1.0 / 6            // = 0.9999999999999999
DD z = DD.fromQuotient(1, 2)
         .add(DD.fromQuotient(1, 3))
         .add(DD.fromQuotient(1, 6));  // (1.0,-4.622231866529366E-33)
z.doubleValue()                        // = 1.0

Note that the error on the final result is influenced by the calculation. For summation of terms of the same sign a very large number of operations would be required to observe a 1 ULP error in the final double result. If terms of the opposite sign are summed then smaller magnitude intermediate terms can be lost due to the limited 106-bit precision. In this case almost total cancellation will produce the incorrect result.

double a = 1;
double b = Math.pow(2, 53);
double c = Math.pow(2, 106);
DD.of(a).add(b).add(c).subtract(c).subtract(b)  // (0, 0)
                                                // NOT (1, 0)

When using products the compound error will accumulate faster but is not influenced by the sign of the operands. A power function pow(int) is provided for compound multiplication of the same argument that is more efficient but does not reduce the error bound.

The DDMath class provides special support for operations using the DD class. The purpose is to supplement the arithmetic operations in the DD class providing greater accuracy at the cost of performance. This includes a power function that returns the result in a fractional representation to avoid over or underflow, with a lower error bound than the equivalent function in DD.

double a = 1.5 * Math.pow(2, 1023);
double b = 4 * Math.pow(2, -1022);
DD x = DD.of(a);
DD y = DD.of(b);

// Values too extreme for DD arithmetic!
x.multiply(y).isFinite()  // false

// Create fractional representation as [0.5, 1) * 2^b
int[] xb = {0};
int[] yb = {0};
x = x.frexp(xb);       // (0.75, 0) * 2^1024
y = y.frexp(yb);       // (0.5, 0)  * 2^-1019

DD z = x.multiply(y);  // (0.375, 0)
// Rescale by 2^5
z.scalb(xb[0] + yb[0]).doubleValue()  // = a * b

Math.fma(Double.MAX_VALUE, 2.0, -Double.MAX_VALUE)  // = Double.MAX_VALUE
Math.fma(Double.MIN_VALUE, 0.5,  Double.MIN_VALUE)  // = Double.MIN_VALUE * 1.5

double x = ...
double y = ...
DD z = DD.ofProduct(x, y);
z.hi()  // = x * y;
z.lo()  // = Math.fma(x, y, -x * y)

Rational Numbers

The Fraction class is a representation of a rational number p / q using int values for the numerator p and denominator q.

Instances are constructed with factory methods:

int p = 3;
int q = 4;
Fraction f = Fraction.of(p, q);
// f == 3 / 4

Fraction can be created using a double. The floating-point number may not be exactly representable; a continued fraction is used to converge on a rational representation within a specified absolute tolerance, or a maximum denominator. An ArithmeticException is thrown if the algorithm fails to converge or the number is not representable.

int maxIterations = 100;
//                               tolerance
Fraction a = Fraction.from(0.6152, 0.02,     maxIterations);
Fraction b = Fraction.from(0.6152, 1.0e-7,   maxIterations);
// a = 3 / 5      == 0.6
// b = 769 / 1250 == 0.6152 (exact)

//                               max denominator
Fraction c = Fraction.from(0.6152, 9);
Fraction d = Fraction.from(0.6152, 9999);
// c = 3 / 5      == 0.6
// d = 769 / 1250 == 0.6152 (exact)

Fraction e = Fraction.from(Math.pow(2, 31));
// e = Integer.MIN_VALUE / -1

// *** Throws an ArithmeticException ***
Fraction.from(Math.pow(2, 32));

The fraction is represented in reduced form using the greatest common divisor. For example 240 / 256 is reduced to 15 / 16. The sign can be either in the numerator or the denominator; the sign of the fraction is accessed using the signum() method. Fraction equality can be evaluated using equals(). The class implements Comparable<Fraction> to allow sorting using the natural order imposed by the fraction's sign and magnitude.

Fraction f1 = Fraction.of(-240, 256);
Fraction f2 = Fraction.of(15, -16);

f1.signum()                 == -1
f1.equals(f2)               == true
f1.compareTo(f2)            == 0
f1.compareTo(Fraction.ZERO) == -1

Fraction.of(Integer.MIN_VALUE, Integer.MIN_VALUE).compareTo(Fraction.ONE) == 0

Arithmetic operations are provided that return new instances of Fraction. Note that fraction operations are exact.

Fraction result = Fraction.of(1, 2).add(Fraction.of(1, 3))
                                   .add(Fraction.of(1, 6));
// result == 1 / 1
// Note: 1.0 / 2 + 1.0 / 3 + 1.0 / 6 == 0.9999999999999999 (inexact)

Fraction extends java.lang.Number and can be converted to other primitive numbers with possible loss of precision:

double a = Fraction.of(1, 8).doubleValue();
double b = Fraction.of(1, 3).doubleValue();
int    c = Fraction.of(8, 3).intValue();
// a == 0.125     (exact)
// b == 1.0 / 3   (inexact)
// c == 2         (inexact, whole number part of 2 2/3)

Fraction is limited to rational numbers representable using int values. Operations that cannot be represented will raise an ArithmeticException. Intermediate computations may use the long datatype and arithmetic is relatively fast.

The BigFraction class uses BigInteger for the numerator and denominator and can provide exact arithmetic up to the limitations of BigInteger. This is documented as at least +/- 2^Integer.MAX_VALUE (exclusive) and may support values outside of that range.

The BigFraction class provides all the arithmetic operations of Fraction and can be used when Fraction cannot represent the rational numbers; there will be a performance cost to using BigFraction due to the switch from primitive values to BigInteger for computations.

BigFraction can be constructed from int, long or BigInteger numerator and denominators.

BigFraction a = BigFraction.of(1, 3);
// a == 1 / 3

The BigFraction conversion of double values can be exact, or adjusted to be within a provided tolerance. Note that exact conversion may create an unexpected rational number as it is limited by the precision of the input double value which may be inexact:

BigFraction a = BigFraction.from(1.0 / 3);
// a == 6004799503160661 / 18014398509481984

//                                        max denominator
BigFraction b = BigFraction.from(1.0 / 3, 3);
// b == 1 / 3

Continued Fractions

A generalized continued fraction is an expression of the form:

\( x = b_0 + \cfrac{a_1}{b_1 + \cfrac{a_2}{b_2 + \cfrac{a_3}{b_3 + \cfrac{a_4}{b_4 + \ddots\,}}}} \)

Classes are provided to evaluate the expression through iteration to a provided error tolerance. The user must supply the coefficients \( a \) and \( b \).

The ContinuedFraction class requires a function to supply the nth coefficients; the API allows the fraction function to be evaluated for different arguments. This is useful when the coefficients can be computed using a formula based on the index and the argument. The class is abstract and must be extended to provide the coefficients \( a \) and \( b \). For example the following fraction:

\( \tan(z) = \cfrac{z}{1 - \cfrac{z^2}{3 - \cfrac{z^2}{5 - \ddots\,}}} \)

Can be computed using:

ContinuedFraction cf = new ContinuedFraction() {
    @Override
    public double getA(int n, double x) {
        return n < 2 ? x : -x * x;
    }

    @Override
    public double getB(int n, double x) {
        return n == 0 ? 0 : 2 * n - 1;
    }
};

double eps = 1e-8;
double z = 0.125;
double result = cf.evaluate(z, eps);
// result ~ Math.tan(z)

The GeneralizedContinuedFraction class requires a generator for the next pair of coefficients. This is useful when the coefficients can be computed using the previous values. The class evaluates the fraction using static methods and a supplied generator of coefficients. This can be done with the generator supplying the initial term \( b_0 \) where \( a_0 \) is ignored, or with the generator starting at \( (a_1, b_1) \) and the term \( b_0 \) provided separately.

For example the golden ratio:

\( x = 1 + \cfrac{1}{1 + \cfrac{1}{1 + \cfrac{1}{1 + \ddots\,}}} \)

Can be computed as:

double eps = 1e-10;
double gr1 = GeneralizedContinuedFraction.value(() -> Coefficient.of(1, 1), eps);
double gr2 = GeneralizedContinuedFraction.value(1, () -> Coefficient.of(1, 1), eps);
// gr1 ~ gr2 ~ 1.6180339887498948...

Whether to separate the term \( b_0 \) depends on the fraction and the algorithm convergence. See the GeneralizedContinuedFraction documentation for details. Note that it is also possible to evaluate part of the fraction using only later terms and then computing the final result by including the early terms.

An example using a recursive generator for the evaluation of \( \tan(z) \) is computed using:

static class Tan implements Supplier<Coefficient> {
    private double a;
    private double b = -1;

    /**
     * @param z Argument z
     */
    Tan(double z) {
        a = z;
    }

    @Override
    public Coefficient get() {
        b += 2;
        // Special first case
        if (b == 1) {
            double z = a;
            // Update the term for the rest of the fraction.
            // The continuant is subtracted from the b terms, thus all the
            // remaining a terms are negative.
            a *= -z;
            return Coefficient.of(z, b);
        }
        return Coefficient.of(a, b);
    }
}

double z = 0.125;
double tan1 = GeneralizedContinuedFraction.value(0, new Tan(z));

// Advance 1 term
Tan t = new Tan(z);
Coefficient c = t.get();
double tan2 = c.getA() / GeneralizedContinuedFraction.value(c.getB(), t);

// tan1 ~  Math.tan(z)
// tan2 == Math.tan(z)  (more accurate)

A simple continued fraction will have an \( a \) coefficient of 1. For example the following fraction:

\( \frac{415}{93} = 4 + \cfrac{1}{2 + \cfrac{1}{6 + \cfrac{1}{7}}} \)

can be evaluated as a simple continued fraction using:

private static Supplier<Coefficient> simpleContinuedFraction(double... coefficients) {
    return new Supplier<Coefficient>() {
        /* iteration. */
        private int n;
        /* denominator terms. */
        private double[] b = coefficients.clone();

        @Override
        public Coefficient get() {
            if (n != b.length) {
                // Return the next term
                return Coefficient.of(1, b[n++]);
            }
            return Coefficient.of(0, 1);
        }
    };
}

double result1 = GeneralizedContinuedFraction.value(simpleContinuedFraction(4, 2, 6, 7));
double result2 = GeneralizedContinuedFraction.value(4, simpleContinuedFraction(2, 6, 7));
// result1 ~ result2 ~ 415.0 / 93.0