import rustymath as rm
import rustymath.mathish as mmMaths and statistics
Elementary and special functions, complex numbers, statistical reductions and tolerant comparisons.
Elementary functions
Functions accept scalars or arrays. Integer inputs to float functions compute at an appropriate float width. Float arrays retain their width:
rm.exp([0., 1, 2]), rm.log([1., 10, 100]), rm.sin([0., 0.5, 1])(Array([1.0, 2.7182818, 7.3890561], dtype=float64),
Array([0.0, 2.3025851, 4.6051702], dtype=float64),
Array([0.0, 0.47942554, 0.84147098], dtype=float64))root(n, y) gives the real nth root of y. logb(base, x) takes its base first:
rm.root(3, [8., 27., 1000.]), rm.logb(2, [8., 16.])(Array([2.0, 3.0, 10.0], dtype=float64), Array([3.0, 4.0], dtype=float64))expm1(x) and log1p(x) retain small changes that exp(x)-1 and log(1+x) lose to rounding. Recover a small continuously compounded growth rate from its proportional change:
rate = 1e-16
change = rm.expm1(rate)
rm.exp(rate) - 1, change, rm.log1p(change)(0.0, 1e-16, 1e-16)| Family | Functions |
|---|---|
| Exponentials and logarithms | exp, expm1, log, log1p, log2, log10, logb |
| Powers and roots | power, root, sqrt, cbrt |
| Circular and hyperbolic | sin, cos, tan, arcsin, arccos, arctan, arctan2; sinh, cosh, tanh, arcsinh, arccosh, arctanh |
| Integer bits | bitwise_and, bitwise_or, bitwise_xor, bitwise_not, left_shift, right_shift |
| Boolean operations and predicates | logical_and, logical_or, logical_xor, logical_xnor; isfinite, isinf, isnan |
hypot(x, y) computes a two-dimensional vector’s magnitude without intermediate overflow or underflow from squaring its components. The specialized roots sqrt1mx2, sqrt1px2 and sqrtx2m1 compute their combined expressions directly. For real |x| >= 1, sqrtx2m1(x) is sign(x) * sqrt(x*x-1).
round(x, decimals) rounds values to a decimal place, with ties to even. rint rounds to whole values. Neither changes how many digits the array’s repr displays. degrees and radians convert angular units.
rm.round([1.2345, -2.71828], 2), rm.degrees([0., 1.5707963267948966])(Array([1.23, -2.72], dtype=float64), Array([0.0, 90.0], dtype=float64))derivative(f, *args) gives a native elementwise function’s derivative at its arguments: an array for a function of one argument, and a tuple of its partial derivatives for a function of two. value_and_derivative(f, *args) also gives f’s result, which each slope reads where its formula needs it, so f runs once. A function with no derivative, such as floor, raises ValueError:
value, slope = rm.value_and_derivative(rm.tanh, [-1., 0, 1])
slope, 1 - value**2, rm.derivative(rm.power, 2., 3.)The Rust names follow Rust’s standard library instead. Python’s log, power, arcsin and arctan2 are Rust’s ln, pow, asin and atan2.
Complex numbers
Complex exp, log, power, sin, cos, sinh and cosh use vectorized real loops. Real inputs to sqrt and log give NaN outside their real domains. Supply complex values to request principal complex values:
rm.sqrt([-4., 4.]), rm.sqrt([-4+0j, 3+4j]), rm.log([-1+0j])(Array([nan, 2.0], dtype=float64),
Array([2j, (2+1j)], dtype=complex128),
Array([3.1415927j], dtype=complex128))cis(theta) constructs a unit complex number with angle theta in radians. Components, magnitude and angle are real-valued even for complex inputs. Predicates and equality give Booleans. conj reflects the imaginary component:
z = 2 * rm.cis([0., 0.5, 1.])
z, rm.real(z), rm.imag(z), rm.absolute(z), rm.angle(z), rm.conj(z)(Array([(2+0j), (1.7551651+0.95885108j), (1.0806046+1.682942j)], dtype=complex128),
Array([2.0, 1.7551651, 1.0806046], dtype=float64),
Array([0.0, 0.95885108, 1.682942], dtype=float64),
Array([2.0, 2.0, 2.0], dtype=float64),
Array([0.0, 0.5, 1.0], dtype=float64),
Array([(2-0j), (1.7551651-0.95885108j), (1.0806046-1.682942j)], dtype=complex128))Special functions
The special functions use SciPy’s names:
| Family | Functions |
|---|---|
| Gamma and log gamma | gamma(x), gammaln(x), loggamma(x), digamma(x), polygamma(n, x) |
| Error function and complement | erf(x), erfc(x), erfinv(p), erfcinv(q) |
| Regularized incomplete gamma | gammainc(a, x), gammaincc(a, x), gammaincinv(a, p) |
| Beta and regularized incomplete beta | beta(a, b), betainc(a, b, x), betaincinv(a, b, p) |
| Factorial and binomial | factorial(x), binom(n, k), binom_dn(n, k), binom_dk(n, k) |
| Lambert W | lambertw(x) |
gammaln(x) gives the real log of the magnitude of gamma. loggamma(x) gives the principal branch of log gamma as a complex value, including for real inputs. Root-level factorial extends factorial to floating-point and complex arguments. Use number_theory.factorial for exact integers:
rm.gamma([0.5, 5, 10]), rm.erfinv(0.5)(Array([1.7724539, 24.0, 362880.0], dtype=float64), 0.4769362762044699)digamma(x) is the derivative of gammaln(x), so ψ(1) is minus the Euler–Mascheroni constant. binom_dn(n, k) and binom_dk(n, k) give the binomial coefficient’s slopes in n and k. For a non-whole n with a whole k up to 64, binom and binom_dn are within 2k ulps, not 4. C(n, 2) is n(n - 1)/2, whose slope in n at 10 is 9.5:
rm.digamma(1.), rm.binom_dn(10., 2)(-0.5772156649015329, 9.5)Invert a regularized incomplete beta probability, then recover the probability from its quantile:
quantile = rm.betaincinv(2, 5, 0.9)
quantile, rm.betainc(2, 5, quantile)(0.5103163065514917, 0.9)lambertw supplies the principal branch W₀. It solves w * exp(w) = x:
w = rm.lambertw(10.)
w, w * rm.exp(w)(1.7455280027406994, 10.0)Reductions and statistics
With no axis, reductions summarize all items. Pass an axis position, name or tuple of axes to summarize those axes. Here each row holds one series. Clipping bounds the original values in one pass:
measurements = rm.array([[1., 2, 3], [4, 5, 6]])
rm.mean(measurements, axis=1), rm.std(measurements, axis=1), rm.clip(measurements, 2., 5.)(Array([2.0, 5.0], dtype=float64),
Array([0.81649658, 0.81649658], dtype=float64),
Array([[2.0, 2.0, 3.0], [4.0, 5.0, 5.0]], dtype=float64))diff takes adjacent differences along an axis; ptp gives maximum minus minimum, and trace sums the diagonal of each matrix in the last two axes. Integer differences wrap where they don’t fit, as NumPy’s do. Boolean differences report a change (XOR).
rm.diff(measurements), rm.ptp(measurements, axis=1), rm.trace(measurements)(Array([[1.0, 1.0], [1.0, 1.0]], dtype=float64),
Array([2.0, 2.0], dtype=float64),
6.0)var and std divide by n - ddof, with ddof=0 by default. average broadcasts weights by the leading-axis rule. A single row of observation weights applies to every series:
weights = [[1., 1, 2]]
rm.average(measurements, axis=1, weights=weights), rm.std(measurements, axis=1, ddof=1)(Array([2.25, 5.25], dtype=float64), Array([1.0, 1.0], dtype=float64))cov and corrcoef treat rows as variables unless rowvar=False. Float16 summaries accumulate and return float32. Complex variance and norms return real values:
rm.cov(measurements), rm.corrcoef(measurements)(Array([[1.0, 1.0], [1.0, 1.0]], dtype=float64),
Array([[1.0, 1.0], [1.0, 1.0]], dtype=float64))norm computes an elementwise p-norm over the selected axes. ord=1 sums magnitudes, ord=0 counts nonzeros, and ord=float('inf') selects the greatest magnitude. Two-dimensional input does not imply a matrix-induced norm:
rm.norm([3., 4.]), rm.isfinite([1., float('inf'), float('nan')])(5.0, Array([True, False, False], dtype=bool))histogram returns counts and bin edges. Give a bin count and range, or explicit increasing edges. The final bin includes its upper edge. Values outside the edges and NaNs are ignored:
rm.histogram(measurements, bins=3, range=(0, 6)), rm.histogram(measurements, bins=[0., 2, 5, 6])((Array([1, 2, 3], dtype=int64), Array([0.0, 2.0, 4.0, 6.0], dtype=float64)),
(Array([1, 3, 2], dtype=int64), Array([0.0, 2.0, 5.0, 6.0], dtype=float64)))Scans keep a running result. accumulate(f, x) applies a native binary function from the left, over the flattened input or one axis by position or name. A scalar start precedes each lane without adding an output item. Here subtraction spends a budget of 10 in order. logical_scan takes any Boolean dyad as a 2 by 2 truth table; this table computes running exclusive-or:
rm.accumulate(rm.subtract, [3., 2, 1], start=10), rm.logical_scan([False, True, True, False], [[False, True], [True, False]])(Array([7.0, 5.0, 4.0], dtype=float64),
Array([False, True, False, False], dtype=bool))sliding(f, x, width, axis=-1) folds each window of width neighbouring items by a native binary function with an associative reduction: add, multiply, maximum, minimum or a logical one. A lane of n items gives n - width + 1 folds. Windows of more than 8 items fold by van Herk-Gil-Werman’s method, at a cost for each item that doesn’t grow with the width. sliding_argmax and sliding_argmin give the position along the axis of each window’s first maximum or minimum: 5, where the 9 is, for each window that holds it. A fold along an axis of a view that windows gives runs the same sliding fold over the view’s storage, without building the windows: rm.sum(rm.windows(x, 3), axis=1) gives what sliding(add, x, 3) gives.
v = [3., 1, 4, 1, 5, 9, 2, 6]
rm.sliding(rm.add, [1., 2, 3, 4, 5], 3), rm.sum(rm.windows([1., 2, 3, 4, 5], 3), axis=1), rm.sliding(rm.maximum, v, 4), rm.sliding_argmax(v, 4)(Array([6.0, 9.0, 12.0], dtype=float64),
Array([4.0, 5.0, 9.0, 9.0, 9.0], dtype=float64),
Array([2, 4, 5, 5, 5], dtype=int64))Correlation and convolution
correlate(x, kernel, mode='valid') slides kernel along x and gives the sum of products of the window and the kernel at each position. convolve reverses the kernel first. The modes are SciPy’s: 'valid' keeps windows inside x, 'full' keeps every overlap, and 'same' gives one position for each of x’s, centred as NumPy centres them. The kernel isn’t conjugated.
rm.correlate([1., 2, 3, 4], [1., 10]), rm.convolve([1., 2, 3, 4], [1., 10], 'full')(Array([21.0, 32.0, 43.0], dtype=float64),
Array([1.0, 12.0, 23.0, 34.0, 40.0], dtype=float64))A kernel can have more axes than the axes it slides along. slide sets how many leading axes of x it slides along, and the kernel contracts x’s other axes, as a convolution layer contracts an image’s channels. The kernel’s remaining axes give the outputs at each position. Here a 3 by 3 kernel moves 2 positions at a time over a 9 by 8 image with 3 channels, and gives 4 outputs at each of 4 by 3 positions:
rm.correlate(rm.ones((9, 8, 3)), rm.ones((3, 3, 3, 4)), slide=2, stride=2).shape(4, 3, 4)f and g replace the sum and the product, as in linalg.inner. The maximum over each window of sums is a dilation:
rm.correlate([0., 5, 0, 0, 3], [0., 1, 0], f=rm.maximum, g=rm.add)Array([6.0, 5.0, 3.0], dtype=float64)overlap_add(frames, step=1, axis=-2) is the transpose of taking windows. It adds each frame back where its window was read. Items that several frames cover add up. An inverse short-time Fourier transform uses it:
rm.overlap_add([[1., 1, 1], [1., 1, 1], [1., 1, 1]], 2)Array([1.0, 1.0, 2.0, 1.0, 2.0, 1.0, 1.0], dtype=float64)Mathish mode
rustymath.mathish provides J-inspired conventions through the same function names as rustymath. For float64 and complex128, comparisons use relative tolerance 1e-14, zero divided by zero gives zero, and logb(1., 1.) gives zero. Narrower floats retain exact comparisons and raw division. Integer true division follows its floating computation type. Machine-integer floor division keeps its zero-divisor error.
The namespace does not change ordinary rm functions or array operators. Generic operations use their operand’s conventions, as in rm.accumulate(mm.divide, values). The namespace’s equal, not_equal, less, less_equal, greater and greater_equal select tolerant comparison explicitly:
mathish = mm.equal(0.1 + 0.2, 0.3), mm.divide(0., 0.), mm.logb(1., 1.)
ordinary = rm.equal(0.1 + 0.2, 0.3), rm.divide(0., 0.), rm.logb(1., 1.)
ordinary, mathish((False, nan, nan), (True, 0.0, 0.0))mm.floor, mm.ceil, mm.floor_divide, mm.remainder, mm.gcd and mm.lcm use the same tolerance. A value tolerantly equal to a whole number rounds to that number. mm.floor_divide(0., 0.) gives zero. Division of a nonzero real by zero still gives infinity. Algorithmic stopping tolerances are independent of these conventions. Ordinary rm.floor gives:
rm.floor([2.9999999999999996, 2.5])Array([2.0, 2.0], dtype=float64)floor_divide and remainder agree on near-integer quotients. Here the exact float quotient is below 3. remainder(x, 0) gives x. The arguments follow NumPy’s order, opposite to the Rust residue(y, x):
exact = rm.floor_divide(0.3, 0.1), rm.remainder(0.3, 0.1)
tolerant = mm.floor_divide(0.3, 0.1), mm.remainder(0.3, 0.1)
exact, tolerant((2.0, 0.09999999999999998), (3.0, 0.0))For complex inputs, floor rounds onto the Gaussian-integer lattice using APL’s rule. Remainder uses this complex floor of the quotient. gcd uses Euclid’s algorithm on Gaussian integers:
rm.floor(1.7 + 2.8j), rm.remainder(5+0j, 2+1j), rm.gcd(5+0j, 2+1j)((1+3j), 0j, (2+1j))Checking accuracy
import numpy as npTo check the accuracy yourself, compare with NumPy. np.spacing gives the size of one ulp at each value. Dividing the difference by it gives the distance between the two results in ulps:
xs = np.linspace(-10, 10, 100_001)
y = np.exp(xs)
(np.abs(rm.exp(xs) - y) / np.spacing(y)).max()np.float64(2.0)