Skip to main content

Polynomials, exact arithmetic and number theory

Polynomial calculations, rational matrices and arbitrary-size integers. 
import rustymath as rm

Polynomials

rm.polynomial evaluates polynomials, finds roots, constructs coefficients from roots, differentiates and fits data, under the names of numpy.polynomial.polynomial. Coefficients start from the constant term. Evaluation, construction and derivatives take real or complex arrays. polyvalfromroots(x, roots, m=1) evaluates directly from roots without forming coefficients. The roots of x^3 - 6x^2 + 11x - 6 are 1, 2 and 3:

p = rm.polynomial
roots = rm.array([1., 2, 3])
c = p.polyfromroots(roots)
c, p.polyroots(c), p.polyval([0, 4], c), p.polyder(c), p.polyvalfromroots([0, 4], roots)
(Array([-6.0, 11.0, -6.0, 1.0], dtype=float64),
 Array([(1+0j), (2-0j), (3+0j)], dtype=complex128),
 Array([-6.0, 6.0], dtype=float64),
 Array([11.0, -12.0, 3.0], dtype=float64),
 Array([-6.0, 6.0], dtype=float64))

polyfit(x, y, degree) fits coefficients by least squares. Evaluate the fitted line at new points with the same polyval operation:

fit = p.polyfit([0., 1, 2, 3], [1.1, 2.9, 5.2, 6.8], 1)
fit, p.polyval([4., 5], fit)
(Array([1.09, 1.94], dtype=float64), Array([8.85, 10.79], dtype=float64))

Exact arithmetic

rm.exact computes with integers and fractions.Fraction. It provides rational gcd, lcm, remainder and binomial coefficients, matrix solves, pseudoinverses, determinants, permanents and row reduction, and polynomial construction, evaluation and derivatives. Arguments can be integers, Fractions, nested lists of them, or integer arrays. Results are Fractions, in nested lists for arrays. Matrix calls take one matrix.

from fractions import Fraction
rm.exact.solve([[2, 1], [1, 2]], [Fraction(1, 2), 0]), rm.exact.binom(Fraction(3, 2), 2)
([Fraction(1, 3), Fraction(-1, 6)], Fraction(3, 8))

An underdetermined system has many solutions. exact.solve chooses the one with least Euclidean norm. pinv gives a reusable matrix for the same solve:

a = [[1, 1]]
rm.exact.solve(a, [1]), rm.exact.pinv(a)
([Fraction(1, 2), Fraction(1, 2)], [[Fraction(1, 2)], [Fraction(1, 2)]])

Reduce an augmented system with columns set to its number of unknowns. .rows contains the nonzero reduced rows and .pivots their pivot columns. .determinant is the signed product of the unscaled pivots. It is the coefficient determinant when those columns form a square, full-rank matrix:

e = rm.exact.echelon([[2, 1, 1], [1, 2, 0]], columns=2)
e.rows, e.pivots, e.determinant
([[Fraction(1, 1), Fraction(0, 1), Fraction(2, 3)],
  [Fraction(0, 1), Fraction(1, 1), Fraction(-1, 3)]],
 [0, 1],
 Fraction(3, 1))

The permanent sums the same products as the determinant without permutation signs. It also exists in floating-point rm.linalg:

square = [[1, 2], [3, 4]]
rm.exact.det(square), rm.exact.permanent(square)
(Fraction(-2, 1), Fraction(10, 1))

Rational gcd and remainder compute without rounding. Exact polynomial operations use the same constant-first coefficients:

r = [Fraction(1, 2), Fraction(3, 4)]
rm.exact.gcd(r, Fraction(1, 3)), rm.exact.remainder(r, Fraction(1, 3)), rm.exact.polyval(r, [1, 0, Fraction(1, 2)])
([Fraction(1, 6), Fraction(1, 12)],
 [Fraction(1, 6), Fraction(1, 12)],
 [Fraction(9, 8), Fraction(41, 32)])

Number theory

rm.number_theory tests primality, factorises, finds neighbouring primes and gives Euler’s totient for arbitrary-size nonnegative Python integers. Primality uses Baillie–PSW: no composite below 2^64 passes, and larger results are probable primes with no known counterexample. Long calls stop on Ctrl-C:

nt = rm.number_theory
nt.factor(2**64 + 1), nt.prime_count(10**8), nt.next_prime(10**20)
([(274177, 1), (67280421310721, 1)], 5761455, 100000000000000000039)

primes(start, stop) lists a half-open interval; primes(stop) starts from zero. Enumeration, prime_count and factorial take unsigned 64-bit arguments. Listed primes must fit int64. number_theory.factorial returns an exact Python integer; root-level rm.factorial is floating-point:

nt.primes(100, 120), nt.is_prime(101), nt.totient(36), nt.factorial(30)
(Array([101, 103, 107, 109, 113], dtype=int64),
 True,
 12,
 265252859812191058636308480000000)

Mixed radix

rm.encode and rm.decode alias rm.integer.encode and decode. Digits occupy the last axis, most significant first. Bases [0, 60, 60] split seconds into hours, minutes and seconds. Zero leaves the hours unbounded:

seconds = rm.array([3661, 86399])
hms = rm.encode([0, 60, 60], seconds)
hms, rm.decode([0, 60, 60], hms)
(Array([[1, 1, 1], [23, 59, 59]], dtype=int64),
 Array([3661, 86399], dtype=int64))

Integer decoding wraps at int64’s range, as BPL’s ⊥ does, and integer encoding raises OverflowError when a quotient doesn’t fit int64. polynomial.decode and exact.decode accept field-valued bases and digits. Fractional seconds stay exact:

rm.exact.decode([0, 60, 60], [1, 1, Fraction(3, 2)])
Fraction(7323, 2)