Files
rpmroute/src/freeroute/geometry/vector.py
T
rsp2k 33cb196bc8 Add exact planar points and vectors
Ports the point/vector foundation of geometry/planar: Side and Signum
(three-valued signs), Limits, FloatPoint (approximate), and the exact
IntPoint/IntVector plus projective RationalPoint/RationalVector.

Arithmetic model per type is documented in each module. The key
simplification over the Java source: Python's unbounded int makes the
exact orientation determinants and rational (x,y,z) coordinates trivial,
so side_of is kept exact (upstream uses a double for speed) and no
BigInteger or CRIT_INT overflow promotion is needed. Rational points use
the projective triple with z=0 denoting the point at infinity.

Tests cover collinearity, determinants beyond Java long range, rational
equality/reduction, and integer/rational promotion.
2026-07-11 18:01:21 -06:00

275 lines
8.8 KiB
Python

"""Vectors: ``Vector`` base, ``IntVector``, ``RationalVector``.
Ports ``geometry/planar/{Vector,IntVector,RationalVector}.java``.
**Arithmetic model:**
* :class:`IntVector` — exact, Python ``int`` ``(x, y)``. Determinants and
orientation are exact; Python's unbounded ``int`` removes the ``long`` /
``double`` overflow concern that forces FreeRouting to cap coordinates at
``CRIT_INT`` and use a ``double`` determinant in ``side_of`` (this port keeps
``side_of`` exact).
* :class:`RationalVector` — exact, projective ``(x, y, z)`` Python ``int``
representing the affine vector ``(x/z, y/z)``; ``z`` is kept non-negative.
Mirrors the upstream ``BigInteger`` triple. ``BigIntDirection`` is not needed
because ``int`` is unbounded, so normalized directions are always
:class:`~freeroute.geometry.direction.IntDirection`.
"""
from __future__ import annotations
import math
from .float_point import FloatPoint
from .side import Side, Signum
__all__ = ["Vector", "IntVector", "RationalVector"]
def _add_rational(a: tuple[int, int, int], b: tuple[int, int, int]) -> tuple[int, int, int]:
"""Add two projective coordinate triples (``BigIntAux.add_rational_coordinates``)."""
if a[2] == b[2]:
return (a[0] + b[0], a[1] + b[1], a[2])
return (
a[0] * b[2] + b[0] * a[2],
a[1] * b[2] + b[1] * a[2],
a[2] * b[2],
)
class Vector:
"""Abstract base for plane vectors (translations of points)."""
def _xyz(self) -> tuple[int, int, int]: # pragma: no cover - overridden
raise NotImplementedError
def is_zero(self) -> bool: # pragma: no cover - overridden
raise NotImplementedError
def is_multiple_of_45_degree(self) -> bool:
return self.is_orthogonal() or self.is_diagonal()
def is_orthogonal(self) -> bool: # pragma: no cover - overridden
raise NotImplementedError
def is_diagonal(self) -> bool: # pragma: no cover - overridden
raise NotImplementedError
def side_of(self, other: Vector) -> Side:
"""Side of ``self`` relative to the directed line from zero to ``other``.
Exact: uses the numerators of the projective coordinates (the positive
denominators ``z`` do not change the sign of the determinant).
"""
sx, sy, _ = self._xyz()
ox, oy, _ = other._xyz()
return Side.of(ox * sy - oy * sx)
def projection(self, other: Vector) -> Signum:
"""Sign of the scalar product of ``self`` and ``other`` (exact)."""
sx, sy, _ = self._xyz()
ox, oy, _ = other._xyz()
return Signum.of(sx * ox + sy * oy)
def scalar_product(self, other: Vector) -> float:
v1 = self.to_float()
v2 = other.to_float()
return v1.x * v2.x + v1.y * v2.y
def add(self, other: Vector) -> Vector:
if isinstance(self, IntVector) and isinstance(other, IntVector):
return IntVector(self.x + other.x, self.y + other.y)
x, y, z = _add_rational(self._xyz(), other._xyz())
return RationalVector(x, y, z)
def negate(self) -> Vector: # pragma: no cover - overridden
raise NotImplementedError
def to_float(self) -> FloatPoint: # pragma: no cover - overridden
raise NotImplementedError
def to_normalized_direction(self):
from .direction import IntDirection
x, y, _ = self._xyz()
g = math.gcd(abs(x), abs(y))
if g > 1:
x //= g
y //= g
return IntDirection(x, y)
def length_approx(self) -> float:
return self.to_float().size()
def cos_angle(self, other: Vector) -> float:
result = self.scalar_product(other)
return result / (self.to_float().size() * other.to_float().size())
def angle_approx(self, other: Vector | None = None) -> float:
if other is None:
return IntVector(1, 0).angle_approx(self)
result = math.acos(max(-1.0, min(1.0, self.cos_angle(other))))
if self.side_of(other) == Side.ON_THE_LEFT:
result = -result
return result
def add_to_point(self, point): # pragma: no cover - overridden
raise NotImplementedError
# Standard zero vector.
def _zero() -> IntVector:
return IntVector(0, 0)
class IntVector(Vector):
"""A vector with exact integer coordinates."""
__slots__ = ("x", "y")
def __init__(self, x: int, y: int) -> None:
self.x = int(x)
self.y = int(y)
def _xyz(self) -> tuple[int, int, int]:
return (self.x, self.y, 1)
def __eq__(self, other: object) -> bool:
if not isinstance(other, IntVector):
return NotImplemented
return self.x == other.x and self.y == other.y
def __hash__(self) -> int:
return hash((self.x, self.y))
def __repr__(self) -> str:
return f"IntVector({self.x}, {self.y})"
def is_zero(self) -> bool:
return self.x == 0 and self.y == 0
def negate(self) -> IntVector:
return IntVector(-self.x, -self.y)
def is_orthogonal(self) -> bool:
return self.x == 0 or self.y == 0
def is_diagonal(self) -> bool:
return abs(self.x) == abs(self.y)
def determinant(self, other: IntVector) -> int:
"""Exact determinant of ``self`` and ``other`` (``x*oy - y*ox``)."""
return self.x * other.y - self.y * other.x
def turn_90_degree(self, factor: int) -> IntVector:
n = factor % 4
if n == 0:
return IntVector(self.x, self.y)
if n == 1:
return IntVector(-self.y, self.x)
if n == 2:
return IntVector(-self.x, -self.y)
return IntVector(self.y, -self.x)
def mirror_at_y_axis(self) -> IntVector:
return IntVector(-self.x, self.y)
def mirror_at_x_axis(self) -> IntVector:
return IntVector(self.x, -self.y)
def to_float(self) -> FloatPoint:
return FloatPoint(self.x, self.y)
def add_to_point(self, point):
from .point import IntPoint
if isinstance(point, IntPoint):
return IntPoint(point.x + self.x, point.y + self.y)
return point.translate_by(self)
def change_length_approx(self, length: float) -> Vector:
from .point import ZERO_POINT
new_point = self.to_float().change_size(length)
return new_point.round().difference_by(ZERO_POINT)
class RationalVector(Vector):
"""A vector with exact rational coordinates ``(x/z, y/z)``, ``z`` >= 0."""
__slots__ = ("x", "y", "z")
def __init__(self, x: int, y: int, z: int) -> None:
if z >= 0:
self.x, self.y, self.z = int(x), int(y), int(z)
else:
self.x, self.y, self.z = -int(x), -int(y), -int(z)
@staticmethod
def from_int_vector(v: IntVector) -> RationalVector:
return RationalVector(v.x, v.y, 1)
def _xyz(self) -> tuple[int, int, int]:
return (self.x, self.y, self.z)
def __eq__(self, other: object) -> bool:
if not isinstance(other, RationalVector):
return NotImplemented
# (x1,y1,z1) == (x2,y2,z2) iff cross products vanish.
return self.x * other.z - other.x * self.z == 0 and self.y * other.z - other.y * self.z == 0
def __hash__(self) -> int:
if self.z != 0:
g = math.gcd(math.gcd(abs(self.x), abs(self.y)), self.z)
g = g or 1
return hash((self.x // g, self.y // g, self.z // g))
return hash((self.x, self.y, 0))
def __repr__(self) -> str:
return f"RationalVector({self.x}, {self.y}, {self.z})"
def is_zero(self) -> bool:
return self.x == 0 and self.y == 0
def negate(self) -> RationalVector:
return RationalVector(-self.x, -self.y, self.z)
def is_orthogonal(self) -> bool:
return self.x == 0 or self.y == 0
def is_diagonal(self) -> bool:
return abs(self.x) == abs(self.y)
def turn_90_degree(self, factor: int) -> RationalVector:
n = factor % 4
if n == 0:
return RationalVector(self.x, self.y, self.z)
if n == 1:
return RationalVector(-self.y, self.x, self.z)
if n == 2:
return RationalVector(-self.x, -self.y, self.z)
return RationalVector(self.y, -self.x, self.z)
def mirror_at_y_axis(self) -> RationalVector:
return RationalVector(-self.x, self.y, self.z)
def mirror_at_x_axis(self) -> RationalVector:
return RationalVector(self.x, -self.y, self.z)
def to_float(self) -> FloatPoint:
if self.z == 0:
big = 1.0e38
return FloatPoint(big, big)
return FloatPoint(self.x / self.z, self.y / self.z)
def add_to_point(self, point):
from .point import IntPoint
if isinstance(point, IntPoint):
new_x = self.z * point.x + self.x
new_y = self.z * point.y + self.y
from .point import RationalPoint
return RationalPoint(new_x, new_y, self.z)
return point.translate_by(self)