#include "algebra/algebra.h"
#include <print>
using namespace algebra;
using namespace algebra::literals;
int main(int argc, char* argv[]) {
int e = 2;
integer i = 7_i + e;
rational r = 5/6_q;
rational a = r * i;
std::print("{} | {:.2f}\n", a, a); // prints 15/2 | 7.50
std::print("{:.20}\n", sqrt(2_q, 8)); // prints 1.41421356237309504880
decimal d = 1.1_d; // stored exactly (unlike float and double which can't represent this value)
std::print("{}\n", d); // prints 1.1
return 0;
}
- header-only and no dependencies
- full
constexprandstd::format()support - arbitrary precision and compact algebraic data types
integer/rational/real<>/decimalclasses behave similarly to built-inintandfloattypes (except for overflow)- no heap allocation for integer values in
[-UINT64, UINT64]range - all types cast to any built-in integer and floating point type, and construct from any built-in
integer;
rationalandreal<>also construct fromfloatanddoubleexactly, while anintegeris built from a floating point value withround_to_zero() - no silent overflow / failures (std::runtime_error is thrown)
- output using
std::format()/std::print()/std::ostream/.str() realallows more compact and efficient representation thanrational, but requires roundingreal<2>is similar to built-infloatanddouble, but with arbitrary long mantissa, and 32-bit exponentdecimalalias forreal<10>sizeof(integer)is 16 bytes andsizeof(rational)is 32 bytes, whilestd::vector<>is 24 bytes
- multiplication and division currently use
O(N^2)algorithms where N is number of 64-bit words used (mul_karatsuba()anddivide_bz()are available, but are not used by the operators yet) - the boolean and buffer operations on 2d regions are quadratic in the number of edges
real<Base>division is not exact: it rounds to a fixed number of digits
| header | contents |
|---|---|
algebra/algebra.h |
includes everything below |
algebra/integer.h |
integer and functions on it (also pulls in integer_class.h) |
algebra/rational.h |
rational and functions on it (also pulls in rational_class.h) |
algebra/real.h |
real<Base>, decimal (also pulls in real_class.h) |
algebra/xrational.h |
xrational: rational * sqrt(integer) |
algebra/expr.h |
symbolic expressions (expr, expr_ptr) |
algebra/vector.h |
Vec<D, T> with Vec2 / Vec3 / Vec4 aliases |
algebra/rational_vector.h |
qvec2/3/4 and xvec2/3/4 aliases plus mixed-type vector operators |
algebra/solve_linear.h |
small linear systems and determinants |
algebra/geometry.h |
Line3, Plane3 and their intersections; pulls in the distance and intersection headers below |
algebra/point_segment_squared_distance.h |
point to segment squared distance in 3d |
algebra/segment_segment_squared_distance.h |
segment to segment squared distance in 3d |
algebra/segment_segment_intersection.h |
segment intersection in 2d, as points or as parameters |
algebra/polygon2.h |
MultiPolygon2<T>: a 2d region as rings plus a complement flag |
algebra/polygon2_boolean.h |
union, intersection, difference and symmetric difference of regions |
algebra/polygon2_buffer.h |
dilate, erode and buffer by a convex structuring element |
algebra/polygon2_arc.h |
ArcPolygon2<T>: the same, with circular arc edges |
algebra/polygon2_arc_boolean.h |
ArcRegion<T>: boolean combinations of arc regions, as a tree |
algebra/dual.h |
dual<T> dual numbers for forward mode automatic differentiation |
algebra/kernels.h, algebra/util.h, algebra/types.h |
low level word array kernels and 128-bit helpers (names starting with __ are internal) |
Arbitrary precision signed integer. The magnitude lives in words and the sign in the sign of its
word count, so a value of one word or less needs no heap allocation at all.
Overloaded operators:
- arithmetic
+-*/%+=-=*=/=%=++-- - relational
<><=>===!= - shift
<<>><<=>>= - bitwise
~
Division by zero throws std::runtime_error, and so does a conversion to a built-in type that would
not fit.
- Allows low level access to the vector of individual words of this number.
- The words, least significant first.
- Accepts a leading
-, and'as a digit separator. Bases 2, 8, 10 and 16.
- Negative, zero or positive; the magnitude of the returned value is the word count.
- The least significant word, which is only meaningful when the value is not zero.
- Same as
a = -a, but in place and without memory allocation.
- Whether the value fits into that built-in type, which is what the corresponding cast requires.
- Writes into a caller provided buffer and returns the number of characters written.
- Number of set bits in the two's complement representation.
- Number of bytes used by the words of this number.
- Bit
iof the magnitude. Note that this is not the two's complement bit thatpopcount()counts.
- Remainder modulo a small constant, without a division. Non-negative, also for a negative value.
Always kept in lowest terms with a positive denominator.
- Initializes rational as a/b, and simplifies by removing common divisor.
- Same as
rational(num, den), but assuming they are already simplified.
- Exact conversion.
- Accepts
123,-1/2,1.25and1e-3forms.
- You can use
.simplify()after directly modifying.numand.denfields, to remove common factors from them. - It throws exception if
denis zero. - Note that
rationalis automatically simplified after all arithmetic operations.
- Swap
numanddenin-place. Throws exception ifnumis zero.
- Same as
a = -a, but performed in-place without memory allocation.
std::format() supports {:.N} and {:.Nf}, which round to N digits after the decimal point.
For integer, std::format() supports a fill and alignment ({:*>10}, {:>10}, {:^10}), a width,
and a base: b/B for binary, o for octal, d for decimal, x/X for hexadecimal.
num * Base**exp, with decimal as an alias for real<10>.
- Exact conversion; throws if the denominator is not a power of
Base.
- The nearest value with
digitsdigits after the point, with halves going away from zero.
- Moves trailing factors of
Basefromnumintoexp.
rational * sqrt(integer). Closed under multiplication and division; addition requires
compatible roots.
- Must be positive. Not fully simplified: it can still contain square factors.
expr_ptris an alias forstd::shared_ptr<expr>.- Node types:
expr_integer,expr_rational,expr_power,expr_sum,expr_product,expr_negation,expr_sin,expr_cos,expr_pi,expr_e,expr_var. - Constants:
ZERO_EXPR,ONE_EXPR,PI_EXPR,E_EXPR.
Overloaded operators:
- arithmetic
+-*/ - relational
<><=>===!=(these compare values, by determining the sign of the difference)
A region of the plane, as a set of rings plus a complement flag. Ring2<T> is a
std::vector<Vec2<T>>; the closing edge from back() to front() is implicit, so a ring never
repeats its first vertex.
Membership uses the nonzero winding rule, flipped by complement. An outer boundary winds counter
clockwise and a hole winds clockwise, so a hole cancels the shell containing it and nesting to any
depth works without tracking which ring is whose hole. The flag is what makes inversion exact and
free: the complement of a bounded region is unbounded and has no finite ring representation.
T needs exact arithmetic and division for the predicates to hold, which is why the default is
rational rather than integer.
- The empty region.
~MultiPolygon2<T>()is the whole plane.
Overloaded operators: ~ (complement), | & - ^ (boolean operations, in
polygon2_boolean.h), and ==, which compares the rings structurally rather than as point sets.
The same, with edges that are line segments or circular arcs. A ring is a std::vector<ArcVertex<T>>
(ArcRing2<T>), where each vertex carries the bulge of the edge leaving it:
tan(theta/4)for the arc's included angletheta, and 0 for a straight edge. A positive bulge puts the arc on the left of the edge, a negative one on the right.
That choice is what keeps everything rational: a rational bulge with rational endpoints gives a
rational centre and squared radius, so no coordinate is ever irrational. A bulge cannot describe a
full circle, since theta == 2*pi needs tan(pi/2); circle_ring() uses two half circle edges.
Members and operators mirror MultiPolygon2: rings, complement, is_empty(),
is_whole_plane(), ~ and ==.
A boolean combination of arc regions, kept as a tree and evaluated on demand. Arc regions have no
explicit boolean result: two arcs meet at cx +- sqrt(r*r - dy*dy), which is not rational, and
cutting a further arc at such a point nests the radicals. Membership is exactly computable, so
contains() on a combination is the combination of contains() on its operands, and every leaf test
is the exact rational predicate from polygon2_arc.h.
That gives exact union, intersection, difference, symmetric difference and complement over arc regions, closed under further combination, with no tolerance anywhere. What it does not give is a ring list or an exact area for a combination -- writing those down is the step that needs the irrational points.
- A leaf.
- How many arc regions the combination rests on.
Overloaded operators: | & - ^ and ~. Note that ~region is a strict negation of
contains(), unlike ~polygon, which flips the complement flag and so leaves the boundary belonging
to both sides.
- Throws
std::runtime_errorwith the source location whenvalueis false.
- All three assume the operands are already in
[0, m).
2**e.exp2()throws for a negative exponent.
a = a * a, using half the multiplications ofmul(a, a, out).
- Sub-quadratic multiplication.
qmust not aliasaorb.
acc += a * b/acc -= a * bwithout memory allocation.
quotandremhave to be different objects; either may aliasaorb.
- Returns the remainder, for a signed or unsigned built-in
bof any width.
- Recursive (Burnikel-Ziegler) division; same result as
div().
- Truncates towards zero, so the remainder carries the sign of the dividend. The result type is the
widest one that holds it: signed for a signed divisor, and
integerfor an unsigned one, where a negative remainder fits neither the divisor's type nor its signed counterpart. Divisors of every width up to 128 bits have an overload.
- All
mod()overloads return a value in[0, abs(b)), unlikeoperator%which truncates towards zero. Theinteger&overload replaces its argument in place, so it is chosen for a non const lvalue: spell the operandconst(or use the return value) to get the value form.
abs(a) > abs(b), minimizing memory allocation.
- Bitwise complement of the magnitude, and its two's complement. Both reject a negative value.
- uniformly sample from
[0, (2**n)-1]
- uniformly sample from
[0, count-1];counthas to be positive
- uniformly sample from
[min, max]
- returns
result * (base ** exp)
- Of the magnitudes, so the sign of either argument does not matter.
abs(a * b) / gcd(a, b), with the sign ofa * b.
- Largest
qwithq * q <= x.
- Alternative
isqrt()implementations, kept for benchmarking.
- Very fast, but only approximate for large values.
- Largest
qwithq**n <= a. Throws forn == 0, which is not a root.
- Sets
btosqrt(a)and returns true whenais a perfect square.
- Factors
sqrt(a)intowhole * sqrt(root), accumulating into already initialized arguments.
- Cheap filter that rejects ~98% of non-squares.
a % 63anda % 65, in one pass and without a division.
- The value truncated towards zero. Throws for nan and infinity.
- Deterministic Miller-Rabin.
- Miller-Rabin with the first
roundsprimes as bases (at most 40). - It returns false if n is composite and returns true if n is probably prime.
- Higher value of
roundsindicates more accuracy.
- Prime factorization as (factor, exponent) pairs.
- A divisor of
nstrictly between 1 andn, or 0 when Fermat's method does not find one quickly. A primenreports 0, since the difference of squares it factors into isn = a*a - b*bwitha - b == 1.
- assume that the operands are in
[0, m-1]range
- returns
(a**b) % m
- returns x such that
(a * x) mod m == 1, or false if such number doesn't exist
- Binomial coefficient (n over k).
- The same coefficient reduced modulo
m.
- The base has to be at least two.
- Divides all arguments by their common divisor.
- returns
a * b < c(cheaper than naive multiplication)
- returns
a < b * c(cheaper than naive multiplication)
- returns
a * b < c * d(cheaper than naive multiplication)
- returns
static_cast<unsigned __int128>(a >> e)without memory allocation
- returns
static_cast<uint64_t>(a >> e)without memory allocation integerconverts tocwords, a read only view of its words.
- Newton iteration; the number of correct digits roughly doubles per iteration.
- The part beyond the integer part, so
trunc(a) + fract(a) == a. The sign followsa, the same waystd::modfsplits a floating point value.
- The nearest integer, with halves going away from zero.
- The nearest multiple of
base**-digits, with halves going away from zero.
- round towards 0 to integer
- https://en.wikipedia.org/wiki/Chudnovsky_algorithm for computing PI
nis both the number of series terms and the number of square root iterations.
- Taylor series with
nterms.
- Scales all arguments by the same factor, preserving their ratios.
a * B**exp, forexp >= 0.
- returns
result * (base ** exp)
- Structural equality, unlike
operator==which compares values.
a->sign(), ornulloptwhen the sign cannot be determined.
- Lower and upper bound of the value, when they can be determined.
integer_value / rational_value / power_base / power_exp / sum_values / product_values / negation_value
- Accessors for the corresponding node type.
- Arithmetic operators,
==,dot(),dot2(),cross(),lerp(),abs(),is_zero(),argmax_abs(),div_colinear(),same_sign(),order()/strict_order()/loose_order(), and swizzles such asxy(),yz(),xzy().
- solves
A + B*x = 0, returningNonewhen there is no solution andAnywhen everyxis one
- solves
A + sB + tC = 0, false when there is no unique solution
- solves
A + sB + tC + rD = 0
- Plane equation is
f(x) = (n * x + d) / sqrt(den).
std::variant<None, Vec3<T>, Line3<T>, Plane3<T>> plane_intersection(const Plane3<T>&, const Plane3<T>&, const Plane3<T>&)
T segment_segment_squared_distance(const Vec3<T>& pa, const Vec3<T>& pb, const Vec3<T>& qa, const Vec3<T>& qb)
- Returns
None, a point, or a segment.segment_segment_intersection_param()returns the parameters instead of the points, andsegment_segment_intersects()returns 0, 1 or 2.
- Twice the signed area, which stays integral when
Tis. Positive is counter clockwise.
- Throws for an unbounded region, which has no finite area.
- Reverses the orientation of a ring, turning a shell into a hole and back.
- How many times the region wraps counter clockwise around
p. Undefined on the boundary, so callers teston_boundary()first.
- Closed region membership: the boundary belongs to the region on either side of a complement, so it
is contained by both
aand~a.
- The smallest axis aligned box containing every vertex. Throws for an unbounded or empty region.
- Drops repeated and collinear vertices, and rings that enclose no area (fewer than three vertices, or all vertices on one line). The interior is unchanged, but the boundary of a dropped sliver goes with it, so a point that was only on such a sliver stops being contained.
BoolOpisUnion,Intersection,DifferenceorSymmetricDifference. The operators|,&,-and^call it.- Exact for an exact
T, with no epsilon anywhere: every edge is cut at every crossing, each fragment is classified by sampling a point a provably safe step off its midpoint, and the surviving fragments are stitched into rings oriented with the interior on their left. - Cost is quadratic in the number of edges, for the cutting and for the classification.
- The result is not
simplify()ed: cutting leaves collinear vertices where the inputs met.
MultiPolygon2<T> buffer(const MultiPolygon2<T>& a, const T& size, Element element = square_element<T>)
- Positive size grows, negative shrinks, zero is the identity.
elementmaps a size to a convex structuring element and is any callable, so one that takes more than a size goes in as a lambda:buffer(a, r, [](const T& s) { return polygon_element(s, 16); }).- Note what is not offered: buffering by a Euclidean distance. Moving an edge out by
rneeds the unit normal, i.e.sqrt(dx*dx + dy*dy), and a round join needs a circular arc whose intersections with its neighbours are irrational as well. Neither is representable inrational, so the shape to buffer by is given explicitly instead and the result stays exact.
- Minkowski sum and its dual with a convex
bthat contains the origin.
max(|dx|, |dy|) <= r, i.e. buffering in the Chebyshev metric.
|dx| + |dy| <= r, i.e. the Manhattan metric.
- A convex polygon with
2*sidesvertices inscribed in the circle of radiusr, with rational coordinates from the Pythagorean parametrisation. It is a subset of the disk, so it under-approximates a round buffer as closely as wanted.
- Negates every vertex, giving the
-Bthat erosion needs. Note that the ring comes back with the opposite orientation;dilate()builds convex hulls from the element's vertices and tests the origin against a winding number, so neither cares.
- The counter clockwise convex hull, by monotone chain. Exact for an exact
T.
- The midpoint, centre and squared radius of an arc. All rational; the radius itself is not.
arc_center()andarc_radius2()require a non-zero bulge.
- Whether
plies on the arc (or straight edge) fromatob, endpoints included.
- The same vertices with every arc replaced by its chord.
- The chord polygon's winding number, corrected by one per arc whose circular segment contains
p. Undefined on the boundary and on the chord of an arc.
- Closed membership, as for
MultiPolygon2. A point on a chord is stepped off first, by a step that provably crosses nothing.
- The area of the chord polygon. The true area differs by the circular segments, whose area is
r*r*(theta - sin theta)/2and so is not representable, which is why there is no exactsigned_area()for an arc region.
- A closed circle, as the two half circle edges that a single bulge cannot express.
- Evaluates the tree at
p.
void area_bounds(const ArcRegion<T>& r, const Vec2<T>& min, const Vec2<T>& max, int depth, T& lower, T& undecided)
- Area by subdividing the given box, since the exact area involves
r*r*(theta - sin theta)/2. Boxes that sample as fully inside or fully outside settle, the rest are split untildepthruns out, and what still straddles the boundary is reported asundecided. - The test per box is five sample points, not an exact containment test, so
lowerandundecidedare a good estimate rather than a proven bound. The box has to contain the region.
- Dual numbers for forward mode automatic differentiation, with
+-*/andsqrtpowexplogsincostanatanabs.
using namespace algebra::literals;
| literal | type |
|---|---|
123_i |
integer |
1/2_q |
rational |
1.5_f |
real<2> |
1.1_d |
decimal |
2_x |
xrational |
3_e |
expr_ptr |