Fundamental numerical types in Piranha

Any computer algebra system needs to be able to represent fundamental numerical types such as integers, rationals, (approximations of) reals, etc. C++ and Python provide basic numerical types that can be used with Piranha’s generic algorithms and data structures. However, especially in C++, the numerical types provided by the language - while having the advantage of providing very high performance - can be unsuitable for certain applications (e.g., all the standard integral types provided by C++ have a finite range).

On the C++ side Piranha provides a set of additional fundamental numerical types that can interoperate with the C++ fundamental types and that can be used with Piranha’s generic algorithms and data structures (e.g., as coefficient types in a polynomial). On the Python side, Piranha’s fundamental types are automatically translated to/from corresponding Python types whenever Pyranha’s routine and data structures are used.

The integer type

Piranha’s integer type is used to represent signed integers of arbitrary size (that is, the range of the type is limited only by the available memory). The type is based on the mpz_t type from the GMP library. As an optimisation, for small numbers integer avoids using dynamic memory allocation and can thus have better performance than a straightforward mpz_t wrapper.

integer behaves much like a standard C++ integral type:

  • a default-constructed integer object is initialised to zero;
  • an integer object can be constructed from all the fundamental C++ types (construction from floating-point types results in the truncation of the original value);
  • an integer object can be converted to all the the fundamental C++ types (if the conversion to a C++ integral type overflows the target type, an error is thrown);
  • the division operator performs truncated division;
  • in mixed-mode operations, the rank of integer is higher than any other C++ integral type but still lower than floating-point types (e.g., integer + long results in an integer, integer + float results in a float).

The following C++ code illustrates some features of the integer type:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
#include <iostream>
#include <stdexcept>

#include <piranha/piranha.hpp>

using namespace piranha;

int main()
{
    init();

    // Various ways of constructing an integer.
    std::cout << integer{42} << '\n';
    std::cout << integer{"42"} << '\n';
    std::cout << integer{42.123} << '\n';

    // You can construct integers of arbitrary size with the string constructor.
    std::cout << integer{"12345678987654321234567898765432"} << '\n';

    // Arithmetics and interoperability with fundamental C++ types.
    std::cout << (integer{21} * 2) << '\n';
    std::cout << (1u + integer{41}) << '\n';
    std::cout << (43. - integer{1}) << '\n';
    std::cout << (integer{85} / 2) << '\n';
    std::cout << (integer{42} == 42) << '\n';

    // Exponentiation is provided by the math::pow() function.
    std::cout << math::pow(integer{42}, 2) << '\n';

    // Conversion to a C++ integral type can fail.
    try {
        (void)static_cast<unsigned char>(-integer{"12345678987654321"});
    } catch (const std::overflow_error &) {
        std::cout << "Overflow!\n";
    }

    // The "_z" suffix can be used to create integer literals.
    auto n = 42_z;
    std::cout << (n == 42) << '\n';

    // Calculate (42 choose 21) via the math::binomial() function.
    std::cout << math::binomial(42_z, 21_z) << '\n';
}

In the first lines, a few possible ways of constructing integer objects are shown, including a constructor from string that allows to initialise an integer with arbitrarily large values. In the following lines, some examples of arithmetics with basic C++ types are illustrated.

On line 28, we can see how the math::pow() function can be invoked in order to compute the exponentiation of an integer. C++ lacks a dedicated exponentiation operator and thus the functionality is provided by a function instead. The math::pow() function is part of Piranha’s mathematical library.

Line 38 highlights a couple of handy C++11 features. The snippet 42_z is a user-defined literal: it results in the construction of an integer object via the string "42". The constructed temporary object is then assigned to a variable n in the statement auto n = 42_z. The type of n is automatically deduced via the keyword auto from the right-hand side of the assignment, thus n is an integer object. This is arguably the C++ syntax that most closely matches Python’s syntax.

On line 42, we use the math::binomial() function from the math library to compute the binomial coefficient \({42 \choose 21}\), passing the two arguments as integer objects created via the user-defined literal _z.

As expected, on the Python side things look simpler:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
# Import the pyranha.math submodule.
from pyranha import math

# Various ways of constructing an int.
print(42)
print(int("42"))
print(int(42.123))

# Arbitrarily large ints are supported directly,
# no need to go through a string conversion.
print(12345678987654321234567898765432)

# Interoperability with float.
print(43. - 1)
# Truncated division in Python 2.x, floating-point division in Python 3.x.
print(85 / 2)

# Exponentiation via the '**' operator.
print(42**2)

# Calculate (42 choose 21) using the binomial() function from the math
# submodule.
print(math.binomial(42, 21))

Python indeed provides multiprecision integers as basic types. Some notable differences exist between Python 2.x and Python 3.x in this regard:

  • Python 2.x has two different integral types, int and long, representing small and large integers respectively. Python 3.x instead has a single integral type named int. Piranha is able to deal with both cases;
  • in Python 2.x, integer division is truncated and returns an integral; in Python 3.x, integral division is a true division returning a float.

In the first lines of the Python code, a few ways of constructing an integral object are shown. Contrary to C++, big integers can be constructed directly without passing through a string constructor. Additionally, Python has a dedicated exponentiation operator called **, thus it is not necessary to resort to a separate function for this operation.

On the last line, we see the first real example of using a Piranha function from Python. The pyranha.math.binomial() function will take Python integers as arguments, convert them to C++ integer objects and pass them to the C++ function math::binomial(). The output value of the C++ function will be converted back to a Python integer which will then be returned by the pyranha.math.binomial() function.

The rational type

A second fundamental type provided by Piranha is the rational class, which represents arbitrary-precision rational numbers as integer numerator-denominator pairs. The rational type in Piranha extends smoothly the standard C++ numerical hierarchy, and it obeys the following basic rules:

  • a default-constructed rational object is initialised to zero;
  • a rational object is always kept in the usual canonical form consisting of coprime numerator and denominator, with the denominator always strictly positive. Zero is always represented as 0/1;
  • a rational object can be converted to/from all the basic C++ numerical types and integer (the conversion to an integral type computes the truncated division of numerator and denominator);
  • in mixed-mode operations, the rank of rational is higher than that of integer but lower than that of floating-point types.

The following C++ code showcases a few features of the rational class:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
#include <iostream>
#include <stdexcept>

#include <piranha/piranha.hpp>

using namespace piranha;

int main()
{
    init();

    // Various ways of constructing a rational.
    std::cout << rational{42} << '\n';
    std::cout << rational{"42"} << '\n';
    std::cout << rational{1.5} << '\n';
    std::cout << rational{42, 12} << '\n';
    std::cout << rational{"42/12"} << '\n';
    // The integer class can also be used to construct a rational.
    std::cout << rational{42_z} << '\n';
    std::cout << rational{42_z, 12_z} << '\n';
    // Mixed integral types in the numerator-denominator constructor can be used.
    std::cout << rational{42ull, 12_z} << '\n';

    // Arithmetics and interoperability with fundamental C++ types and integer.
    std::cout << (rational{42, 13} * 2) << '\n';
    std::cout << (1u + rational{42, 13}) << '\n';
    std::cout << (43. - rational{1, 2}) << '\n';
    std::cout << (rational{85} / 13) << '\n';
    std::cout << (rational{84, 2} == 42) << '\n';
    std::cout << (42_z >= rational{84, 3}) << '\n';

    // Exponentiation is provided by the math::pow() function.
    std::cout << math::pow(rational{42, 13}, 2) << '\n';
    std::cout << math::pow(rational{42, 13}, -3_z) << '\n';

    // Conversion to a C++ integral type can fail.
    try {
        (void)static_cast<unsigned char>(-rational{42, 5});
    } catch (const std::overflow_error &) {
        std::cout << "Overflow!\n";
    }
    // Conversion to integer cannot fail, and yields
    // the truncated value.
    std::cout << static_cast<integer>(rational{10, 3}) << '\n';

    // The "_q" suffix can be used to create rational literals.
    auto r = 42_q;
    std::cout << (r == 42) << '\n';
    // The literal can also be used to initialise a rational from numerator and
    // denominator without using any explicit constructor.
    r = 42 / 13_q;
    std::cout << r << '\n';

    // Calculate (42/13 choose 21) via the math::binomial() function.
    std::cout << math::binomial(42 / 13_q, 21_z) << '\n';
}

In the first code block, a few ways of constructing a rational are shown. In addition to the constructors from interoperable types, a rational can also be constructed from an integral numerator-denominator pair. The string representation for a rational consists, intuitively, of the representation of numerator and denominator in base 10, separated by a / sign. The denominator can be omitted if it is unitary.

In the second code block, some examples of arithmetic and logical operations involving rational and interoperable types are displayed.

On lines 33-34, we use the math::pow() function to compute integral powers of a rational. math::pow() is able to accept many different combinations of argument types. In this specific case, with a rational base and an integral exponent, the result will be exact and its type will be rational.

In the fourth code block, a couple of examples of conversion from rational are shown. Conversion to C++ integral types can fail in case of overflow, whereas conversion to integer will never fail.

In the fifth code block, a few usages of the user-defined literal _q are displayed. _q can be appended to an integer literal to signal that a rational object is to be constructed using that literal. This can be combined with the division operator to yield a user-friendly syntax to initialise rational objects, as shown on line 51:

r = 42/13_q;

This line is effectively interpreted by the compiler as:

r = 42/rational{"13"};

In the last code block, we can see another invocation of the math::binomial() function. This time the top argument is a rational, wheras the bottom argument is an integer. The specialisation of the binomial function for these two types will yield the exact result as a rational.

On the Python side things are again simpler:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
# Import the pyranha.math submodule.
from pyranha import math
# Import the standard Fraction class.
from fractions import Fraction as F

# Various ways of constructing a rational.
print(F(42))
print(F("42"))
print(F(1.5))
print(F(42, 12))
print(F("42/12"))

# Arithmetics and interoperability with float and int.
print(F(42, 13) * 2)
print(1 + F(42, 13))
print(43. - F(1, 2))
print(F(85) / 13)
print(F(84, 2) == 42)
print(42 >= F(84, 3))

# Exponentiation is provided by the standard '**' operator.
print(F(42, 13)**2)
print(F(42, 13)**-3)

# Conversion to int yields the truncated value.
print(int(F(10, 3)))

# Calculate (42/13 choose 21) via the math::binomial() function.
print(math.binomial(F(42, 13), 21))

Although Python does not provide a rational type as a builtin, a rational class named Fraction is available in the standard fractions module since Python 2.6. A few simple examples of usage of the Fraction class are shown in the Python code.

Fraction instances are automatically converted to/from rational by Pyranha as needed. For instance, on the last line of the Python code we see another usage of the pyranha.math.binomial() function. This time the arguments, of type Fraction and int, are automatically converted to rational and integer before being passed to the math::binomial() C++ function. The rational result of the C++ function is then converted back to Fraction and returned.

The real type

The third basic numerical type provided by Piranha is called real, and it represents arbitrary-precision floating-point numbers. It consists of a thin wrapper around the mpfr_t type from the GNU MPFR library. real is essentially a floating-point type whose number of significant digits can be selected at runtime.

The real type obeys the following basic rules:

  • a default-constructed real object is initialised to zero;
  • the precision (i.e., the number of significant digits) is measured in bits, and it can be set at different values for different instances of real. The default value is 113 bits (similar to IEEE 754 quadruple-precision);
  • a real object can be converted to/from all the basic C++ numerical types, integer and rational (the conversion to an integral type truncates the original value);
  • in mixed-mode operations, the rank of real is higher than that of integer, rational and any numeric C++ type;
  • operations involving multiple real instances with different precisions will produce a result with the highest precision among the operands.

The following C++ code showcases a few features of the real class:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
#include <iostream>
#include <stdexcept>

#include <piranha/piranha.hpp>

using namespace piranha;

int main()
{
    init();

    // Various ways of constructing a real.
    std::cout << real{42} << '\n';
    std::cout << real{42.3} << '\n';
    std::cout << real{"1.2345"} << '\n';
    std::cout << real{43_z} << '\n';
    std::cout << real{43 / 45_q} << '\n';
    // Construct with a custom precision of 500 bits.
    std::cout << real{42, 500} << '\n';
    std::cout << real{42.3, 500} << '\n';
    std::cout << real{"1.2345", 500} << '\n';
    // Construct non-finite values.
    std::cout << real{"inf"} << '\n';
    std::cout << real{"-inf"} << '\n';
    std::cout << real{"nan"} << '\n';

    // Basic arithmetic operations.
    std::cout << real{42} + 1 << '\n';
    std::cout << real{42} * 2 << '\n';
    std::cout << 1.5 / real{42} << '\n';
    std::cout << 1_z + real{42} << '\n';
    std::cout << 1 / 2_q - real{42} << '\n';
    // The precision of the result is the highest precision among the operands.
    std::cout << (real{42} / real{7.1, 300}) << '\n';
    std::cout << (real{42} / real{7.1, 300}).get_prec() << '\n';

    // Conversion to a C++ integral type can fail.
    try {
        (void)static_cast<unsigned char>(-real{42.5});
    } catch (const std::overflow_error &) {
        std::cout << "Overflow!\n";
    }
    // Conversion to integer cannot fail for finite values, and yields
    // the truncated value.
    std::cout << static_cast<integer>(real{10.3}) << '\n';

    // The "_r" suffix can be used to create real literals.
    auto r = 42.123_r;
    std::cout << (r == real{"42.123"}) << '\n';

    // Calculate (-42.123 choose 21) via the math::binomial() function.
    std::cout << math::binomial(-42.123_r, 21) << '\n';
}

In the first block, a few ways of constructing real objects are shown: from C++ arithmetic types, from string, from integer and rational. The two-arguments constructors allow construction with custom precision. It is important to note that construction from a C++ floating-point type is different from string construction. That is,

real{1.23};

is in general different from

real{"1.23"};

In the first case, the 1.23 literal is first converted to a double by the compiler, and the double is then used to initialise the real object. In the second case, the real object is initialised directly with the string "1.23". On a typical desktop machine, in the first case printing to screen the constructed real object will produce

1.22999999999999998223643160599749535

In the second case the output will be:

1.22999999999999999999999999999999998

In the second code block, a few examples of arithmetic operations involving real are shown. Operations involving multiple real objects with different precision will produce a result with the highest precision among the operands.

In the third code block, a few usages of the explicit conversion operator are shown.

In the fourth code block, the _r string literal is introduced. Similarly to _z and _q, it can be used to construct a real object with default precision from a floating-point literal. Construction via _r is equivalent to construction from string.

Finally, in the last block, we can see another invocation of the math::binomial() function. This time the top argument is a real, wheras the bottom argument is an int. The specialisation of the binomial function for these two types will yield the approximate result as a real, computed with the default precision of 113 bits.

On the Python side, C++ real objects are automatically converted to/from mpf objects from the mpmath Python library.

Note

The mpmath module needs not to be installed when compiling Pyranha. Its presence is detected by Pyranha at runtime.

Whereas in the case of integers and rationals the conversion to/from Python is straightforward and unambiguous, in the case of real things are slightly more complicated. In mpmath, there is a global precision setting which is user-configurable (either via mpmath.mp.prec or mpmath.mp.dps) and one works without specifying different precisions for different mpf objects.

The rules of conversion between real and mpf are the following:

  • when a real object is returned to Python by some C++ function exposed in Pyranha, it will be converted to an mpf object with the current mpmath global precision;
  • when an mpf object is passed from Python into a C++ function exposed in Pyranha accepting a real argument, the mpf will be converted to a real which will have the same precision as the global mpmath precision.

Potential pitfalls

TLDR