r"""
A loose collection of tools for Euclidean geometry in the plane.
.. SEEALSO::
:mod:`flatsurf.geometry.circle` for everything specific to circles in the plane
"""
######################################################################
# This file is part of sage-flatsurf.
#
# Copyright (C) 2016-2020 Vincent Delecroix
# 2020-2023 Julian Rüth
#
# sage-flatsurf is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 2 of the License, or
# (at your option) any later version.
#
# sage-flatsurf is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with sage-flatsurf. If not, see <https://www.gnu.org/licenses/>.
######################################################################
[docs]
def is_cosine_sine_of_rational(cos, sin, scaled=False):
r"""
Check whether the given pair is a cosine and sine of a same rational angle.
INPUT:
- ``cos`` -- a number
- ``sin`` -- a number
- ``scaled`` -- a boolean (default: ``False``); whether to allow ``cos``
and ``sin`` to be scaled by the same positive algebraic number
EXAMPLES::
sage: from flatsurf.geometry.euclidean import is_cosine_sine_of_rational
sage: c = s = AA(sqrt(2))/2
sage: is_cosine_sine_of_rational(c, s)
True
sage: c = AA(sqrt(3))/2
sage: s = AA(1/2)
sage: is_cosine_sine_of_rational(c, s)
True
sage: c = AA(sqrt(5)/2)
sage: s = (1 - c**2).sqrt()
sage: c**2 + s**2
1.000000000000000?
sage: is_cosine_sine_of_rational(c, s)
False
sage: c = (AA(sqrt(5)) + 1)/4
sage: s = (1 - c**2).sqrt()
sage: is_cosine_sine_of_rational(c, s)
True
sage: K.<sqrt2> = NumberField(x**2 - 2, embedding=1.414)
sage: is_cosine_sine_of_rational(K.zero(), -K.one())
True
TESTS::
sage: from pyexactreal import ExactReals # optional: pyexactreal # random output due to matplotlib warnings with some combinations of setuptools and matplotlib
sage: R = ExactReals() # optional: pyexactreal
sage: is_cosine_sine_of_rational(R.one(), R.zero()) # optional: pyexactreal
True
"""
from sage.all import AA
# We cannot check in AA due to https://github.com/flatsurf/exact-real/issues/172
# We just trust that non-algebraic elements won't allow conversion to AA.
# if cos not in AA:
# return False
# if sin not in AA:
# return False
if not scaled:
if cos**2 + sin**2 != 1:
return False
try:
cos = AA(cos)
except ValueError:
# This is a replacement for the "in AA" checked disabled above.
return False
cos = cos.as_number_field_element(embedded=True)
# We need an explicit conversion to the number field due to https://github.com/sagemath/sage/issues/35613
cos = cos[0](cos[1])
try:
sin = AA(sin)
except ValueError:
# This is a replacement for the "in AA" checked disabled above.
return False
sin = sin.as_number_field_element(embedded=True)
# We need an explicit conversion to the number field due to https://github.com/sagemath/sage/issues/35613
sin = sin[0](sin[1])
from sage.all import ComplexBallField
CBF = ComplexBallField(53)
x = CBF(cos) + CBF.gen(0) * CBF(sin)
xN = x
# Suppose that (cos, sin) are indeed sine and cosine of a rational angle.
# Then x = cos + I*sin generates a cyclotomic field C and for some N we
# have x^N = ±1. Since C is contained in the compositum of K=Q(cos) and
# L=Q(i*sin) and Q(cos) and Q(sin) are both contained in C, the degree of C
# is bounded from above by twice (accounting for the imaginary unit) the
# degrees of K and L. The degree of C is the totient of N which is bounded
# from below by n / (e^γ loglog n + 3 / loglog n) [cf. wikipedia].
degree_bound = 2 * cos.minpoly().degree() * sin.minpoly().degree()
from itertools import count
for n in count(2):
xN *= x
c = xN.real()
s = xN.imag()
if xN.real().contains_zero() or xN.imag().contains_zero():
c, s = cos, sin
for i in range(n - 1):
c, s = c * cos - s * sin, s * cos + c * sin
if c == 0 or s == 0:
return True
CBF = ComplexBallField(CBF.precision() * 2)
x = CBF(cos) + CBF.gen(0) * CBF(sin)
xN = x**n
from math import log
if n / (2.0 * log(log(n)) + 3 / log(log(n))) > 2 * degree_bound:
return False
[docs]
def acos(cos_angle, numerical=False):
r"""
Return the arccosine of ``cos_angle`` as a multiple of 2π, i.e., as a value
between 0 and 1/2.
INPUT:
- ``cos_angle`` -- a floating point number, the cosine of an angle
- ``numerical`` -- a boolean (default: ``False``); whether to return a
numerical approximation of the arccosine or try to reconstruct an exact
rational value for the arccosine (in radians.)
EXAMPLES::
sage: from flatsurf.geometry.euclidean import acos
sage: acos(1)
0
sage: acos(.5)
1/6
sage: acos(0)
1/4
sage: acos(-.5)
1/3
sage: acos(-1)
1/2
sage: acos(.25)
Traceback (most recent call last):
...
NotImplementedError: cannot recover a rational angle from these numerical results
sage: acos(.25, numerical=True)
0.2097846883724169
"""
import math
angle = math.acos(cos_angle) / (2 * math.pi)
assert 0 <= angle <= 0.5
if numerical:
return angle
# fast and dirty way using floating point approximation
from sage.all import RR
angle_rat = RR(angle).nearby_rational(0.00000001)
if angle_rat.denominator() > 256:
raise NotImplementedError(
"cannot recover a rational angle from these numerical results"
)
return angle_rat
[docs]
def angle(u, v, numerical=False):
r"""
Return the angle between the vectors ``u`` and ``v`` divided by `2 \pi`.
INPUT:
- ``u``, ``v`` - vectors
- ``numerical`` - boolean (default: ``False``), whether to return floating
point numbers
EXAMPLES::
sage: from flatsurf.geometry.euclidean import angle
As the implementation is dirty, we at least check that it works for all
denominator up to 20::
sage: u = vector((AA(1),AA(0)))
sage: for n in xsrange(1,20): # long time (1.5s)
....: for k in xsrange(1,n):
....: v = vector((AA(cos(2*k*pi/n)), AA(sin(2*k*pi/n))))
....: assert angle(u,v) == k/n
The numerical version (working over floating point numbers)::
sage: import math
sage: u = (1, 0)
sage: for n in xsrange(1,20):
....: for k in xsrange(1,n):
....: a = 2 * k * math.pi / n
....: v = (math.cos(a), math.sin(a))
....: assert abs(angle(u,v,numerical=True) * 2 * math.pi - a) < 1.e-10
If the angle is not rational, then the method returns an element in the real
lazy field::
sage: v = vector((AA(sqrt(2)), AA(sqrt(3))))
sage: a = angle(u, v)
Traceback (most recent call last):
...
NotImplementedError: cannot recover a rational angle from these numerical results
sage: a = angle(u, v, numerical=True)
sage: a # abs tol 1e-14
0.14102355421224375
sage: exp(2*pi.n()*CC(0,1)*a)
0.632455532033676 + 0.774596669241483*I
sage: v / v.norm()
(0.6324555320336758?, 0.774596669241484?)
"""
import math
u0 = float(u[0])
u1 = float(u[1])
v0 = float(v[0])
v1 = float(v[1])
cos_uv = (u0 * v0 + u1 * v1) / math.sqrt((u0 * u0 + u1 * u1) * (v0 * v0 + v1 * v1))
if cos_uv < -1.0:
assert cos_uv > -1.0000001
cos_uv = -1.0
elif cos_uv > 1.0:
assert cos_uv < 1.0000001
cos_uv = 1.0
angle = acos(cos_uv, numerical=numerical)
return 1 - angle if u0 * v1 - u1 * v0 < 0 else angle
# a neater way is provided below by working only with number fields
# but this method is slower...
# sqnorm_u = u[0]*u[0] + u[1]*u[1]
# sqnorm_v = v[0]*v[0] + v[1]*v[1]
#
# if sqnorm_u != sqnorm_v:
# # we need to take a square root in order that u and v have the
# # same norm
# u = (1 / AA(sqnorm_u)).sqrt() * u.change_ring(AA)
# v = (1 / AA(sqnorm_v)).sqrt() * v.change_ring(AA)
# sqnorm_u = AA.one()
# sqnorm_v = AA.one()
#
# cos_uv = (u[0]*v[0] + u[1]*v[1]) / sqnorm_u
# sin_uv = (u[0]*v[1] - u[1]*v[0]) / sqnorm_u
[docs]
def ccw(v, w):
r"""
Return a positive number if the turn from ``v`` to ``w`` is
counterclockwise, a negative number if it is clockwise, and zero if the two
vectors are collinear.
.. NOTE::
This function is sometimes also referred to as the wedge product or
simply the determinant. We chose the more customary name ``ccw`` from
computational geometry here.
EXAMPLES::
sage: from flatsurf.geometry.euclidean import ccw
sage: ccw((1, 0), (0, 1))
1
sage: ccw((1, 0), (-1, 0))
0
sage: ccw((1, 0), (0, -1))
-1
sage: ccw((1, 0), (1, 0))
0
"""
return v[0] * w[1] - v[1] * w[0]
[docs]
def is_parallel(v, w):
r"""
Return whether the vectors ``v`` and ``w`` are parallel (but not
anti-parallel.)
EXAMPLES::
sage: from flatsurf.geometry.euclidean import is_parallel
sage: is_parallel((0, 1), (0, 1))
True
sage: is_parallel((0, 1), (0, 2))
True
sage: is_parallel((0, 1), (0, -2))
False
sage: is_parallel((0, 1), (0, 0))
False
sage: is_parallel((0, 1), (1, 0))
False
TESTS::
sage: V = QQ**2
sage: is_parallel(V((0,1)), V((0,2)))
True
sage: is_parallel(V((1,-1)), V((2,-2)))
True
sage: is_parallel(V((4,-2)), V((2,-1)))
True
sage: is_parallel(V((1,2)), V((2,4)))
True
sage: is_parallel(V((0,2)), V((0,1)))
True
sage: is_parallel(V((1,1)), V((1,2)))
False
sage: is_parallel(V((1,2)), V((2,1)))
False
sage: is_parallel(V((1,2)), V((1,-2)))
False
sage: is_parallel(V((1,2)), V((-1,-2)))
False
sage: is_parallel(V((2,-1)), V((-2,1)))
False
"""
if ccw(v, w) != 0:
# vectors are not collinear
return False
return v[0] * w[0] + v[1] * w[1] > 0
[docs]
def is_anti_parallel(v, w):
r"""
Return whether the vectors ``v`` and ``w`` are anti-parallel, i.e., whether
``v`` and ``-w`` are parallel.
EXAMPLES::
sage: from flatsurf.geometry.euclidean import is_anti_parallel
sage: V = QQ**2
sage: is_anti_parallel(V((0,1)), V((0,-2)))
True
sage: is_anti_parallel(V((1,-1)), V((-2,2)))
True
sage: is_anti_parallel(V((4,-2)), V((-2,1)))
True
sage: is_anti_parallel(V((-1,-2)), V((2,4)))
True
sage: is_anti_parallel(V((1,1)), V((1,2)))
False
sage: is_anti_parallel(V((1,2)), V((2,1)))
False
sage: is_anti_parallel(V((0,2)), V((0,1)))
False
sage: is_anti_parallel(V((1,2)), V((1,-2)))
False
sage: is_anti_parallel(V((1,2)), V((-1,2)))
False
sage: is_anti_parallel(V((2,-1)), V((-2,-1)))
False
"""
return is_parallel(v, -w)
[docs]
def line_intersection(p1, p2, q1, q2):
r"""
Return the point of intersection between the line joining p1 to p2
and the line joining q1 to q2. If the lines do not have a single point of
intersection, we return None. Here p1, p2, q1 and q2 should be vectors in
the plane.
"""
if ccw(p2 - p1, q2 - q1) == 0:
return None
# Since the wedge product is non-zero, the following is invertible:
from sage.all import matrix
m = matrix([[p2[0] - p1[0], q1[0] - q2[0]], [p2[1] - p1[1], q1[1] - q2[1]]])
return p1 + (m.inverse() * (q1 - p1))[0] * (p2 - p1)
[docs]
def is_segment_intersecting(e1, e2):
r"""
Return whether the segments ``e1`` and ``e2`` intersect.
OUTPUT:
- ``0`` - do not intersect
- ``1`` - one endpoint in common
- ``2`` - non-trivial intersection
EXAMPLES::
sage: from flatsurf.geometry.euclidean import is_segment_intersecting
sage: is_segment_intersecting(((0,0),(1,0)),((0,1),(0,3)))
0
sage: is_segment_intersecting(((0,0),(1,0)),((0,0),(0,3)))
1
sage: is_segment_intersecting(((0,0),(1,0)),((0,-1),(0,3)))
2
sage: is_segment_intersecting(((-1,-1),(1,1)),((0,0),(2,2)))
2
sage: is_segment_intersecting(((-1,-1),(1,1)),((1,1),(2,2)))
1
"""
if e1[0] == e1[1] or e2[0] == e2[1]:
raise ValueError("degenerate segments")
elts = [e[i][j] for e in (e1, e2) for i in (0, 1) for j in (0, 1)]
from sage.structure.element import get_coercion_model
cm = get_coercion_model()
base_ring = cm.common_parent(*elts)
if isinstance(base_ring, type):
from sage.structure.coerce import py_scalar_parent
base_ring = py_scalar_parent(base_ring)
from sage.all import matrix
m = matrix(base_ring, 3)
xs1, ys1 = map(base_ring, e1[0])
xt1, yt1 = map(base_ring, e1[1])
xs2, ys2 = map(base_ring, e2[0])
xt2, yt2 = map(base_ring, e2[1])
m[0] = [xs1, ys1, 1]
m[1] = [xt1, yt1, 1]
m[2] = [xs2, ys2, 1]
s0 = m.det()
m[2] = [xt2, yt2, 1]
s1 = m.det()
if (s0 > 0 and s1 > 0) or (s0 < 0 and s1 < 0):
# e2 stands on one side of the line generated by e1
return 0
m[0] = [xs2, ys2, 1]
m[1] = [xt2, yt2, 1]
m[2] = [xs1, ys1, 1]
s2 = m.det()
m[2] = [xt1, yt1, 1]
s3 = m.det()
if (s2 > 0 and s3 > 0) or (s2 < 0 and s3 < 0):
# e1 stands on one side of the line generated by e2
return 0
if s0 == 0 and s1 == 0:
assert s2 == 0 and s3 == 0
if xt1 < xs1 or (xt1 == xs1 and yt1 < ys1):
xs1, xt1 = xt1, xs1
ys1, yt1 = yt1, ys1
if xt2 < xs2 or (xt2 == xs2 and yt2 < ys2):
xs2, xt2 = xt2, xs2
ys2, yt2 = yt2, ys2
if xs1 == xt1 == xs2 == xt2:
xs1, xt1, xs2, xt2 = ys1, yt1, ys2, yt2
assert xs1 < xt1 and xs2 < xt2, (xs1, xt1, xs2, xt2)
if (xs2 > xt1) or (xt2 < xs1):
return 0 # no intersection
elif (xs2 == xt1) or (xt2 == xs1):
return 1 # one endpoint in common
else:
assert (
xs1 <= xs2 < xt1
or xs1 < xt2 <= xt1
or (xs2 < xs1 and xt2 > xt1)
or (xs2 > xs1 and xt2 < xt1)
), (xs1, xt1, xs2, xt2)
return 2 # one dimensional
elif s0 == 0 or s1 == 0:
# treat alignment here
if s2 == 0 or s3 == 0:
return 1 # one endpoint in common
else:
return 2 # intersection in the middle
return 2 # middle intersection
[docs]
def is_between(e0, e1, f):
r"""
Check whether the vector ``f`` is strictly in the sector formed by the vectors
``e0`` and ``e1`` (in counter-clockwise order).
EXAMPLES::
sage: from flatsurf.geometry.euclidean import is_between
sage: is_between((1, 0), (1, 1), (2, 1))
True
sage: from itertools import product
sage: vecs = [(1, 0), (1, 1), (0, 1), (-1, 1), (-1, 0), (-1, -1), (0, -1), (1, -1)]
sage: for (i, vi), (j, vj), (k, vk) in product(enumerate(vecs), repeat=3):
....: assert is_between(vi, vj, vk) == ((i == j and i != k) or i < k < j or k < j < i or j < i < k), ((i, vi), (j, vj), (k, vk))
"""
if e0[0] * e1[1] > e1[0] * e0[1]:
# positive determinant
# [ e0[0] e1[0] ]^-1 = [ e1[1] -e1[0] ]
# [ e0[1] e1[1] ] [-e0[1] e0[0] ]
# f[0] * e1[1] - e1[0] * f[1] > 0
# - f[0] * e0[1] + e0[0] * f[1] > 0
return e1[1] * f[0] > e1[0] * f[1] and e0[0] * f[1] > e0[1] * f[0]
elif e0[0] * e1[1] == e1[0] * e0[1]:
# aligned vector
if e0[0] * e1[0] >= 0 and e0[1] * e1[1] >= 0:
return f[0] * e0[1] != f[1] * e0[0] or f[0] * e0[0] < 0 or f[1] * e0[1] < 0
else:
return e0[0] * f[1] > e0[1] * f[0]
else:
# negative determinant
# [ e1[0] e0[0] ]^-1 = [ e0[1] -e0[0] ]
# [ e1[1] e0[1] ] [-e1[1] e1[0] ]
# f[0] * e0[1] - e0[0] * f[1] > 0
# - f[0] * e1[1] + e1[0] * f[1] > 0
return e0[1] * f[0] < e0[0] * f[1] or e1[0] * f[1] < e1[1] * f[0]
[docs]
def solve(x, u, y, v):
r"""
Return (a,b) so that: x + au = y + bv
INPUT:
- ``x``, ``u``, ``y``, ``v`` -- two dimensional vectors
EXAMPLES::
sage: from flatsurf.geometry.euclidean import solve
sage: K.<sqrt2> = NumberField(x^2 - 2, embedding=AA(2).sqrt())
sage: V = VectorSpace(K,2)
sage: x = V((1,-sqrt2))
sage: y = V((1,1))
sage: a = V((0,1))
sage: b = V((-sqrt2, sqrt2+1))
sage: u = V((0,1))
sage: v = V((-sqrt2, sqrt2+1))
sage: a, b = solve(x,u,y,v)
sage: x + a*u == y + b*v
True
sage: u = V((1,1))
sage: v = V((1,sqrt2))
sage: a, b = solve(x,u,y,v)
sage: x + a*u == y + b*v
True
"""
d = -u[0] * v[1] + u[1] * v[0]
if d.is_zero():
raise ValueError("parallel vectors")
a = v[1] * (x[0] - y[0]) + v[0] * (y[1] - x[1])
b = u[1] * (x[0] - y[0]) + u[0] * (y[1] - x[1])
return (a / d, b / d)
[docs]
def projectivization(x, y, signed=True, denominator=None):
r"""
Return a simplified version of the projective coordinate [x: y].
If ``signed`` (the default), the second coordinate is made non-negative;
otherwise the coordinates keep their signs.
If ``denominator`` is ``False``, returns [x/y: 1] up to sign. Otherwise,
the returned coordinates have no denominator and no non-unit gcd.
TESTS::
sage: from flatsurf.geometry.euclidean import projectivization
sage: projectivization(2/3, -3/5, signed=True, denominator=True)
(10, -9)
sage: projectivization(2/3, -3/5, signed=False, denominator=True)
(-10, 9)
sage: projectivization(2/3, -3/5, signed=True, denominator=False)
(10/9, -1)
sage: projectivization(2/3, -3/5, signed=False, denominator=False)
(-10/9, 1)
sage: projectivization(-1/2, 0, signed=True, denominator=True)
(-1, 0)
sage: projectivization(-1/2, 0, signed=False, denominator=True)
(1, 0)
sage: projectivization(-1/2, 0, signed=True, denominator=False)
(-1, 0)
sage: projectivization(-1/2, 0, signed=False, denominator=False)
(1, 0)
"""
from sage.all import Sequence
parent = Sequence([x, y]).universe()
if y:
z = x / y
if denominator is True or (denominator is None and hasattr(z, "denominator")):
d = parent(z.denominator())
else:
d = parent(1)
if signed and y < 0:
d *= -1
return (z * d, d)
elif signed and x < 0:
return (parent(-1), parent(0))
else:
return (parent(1), parent(0))
[docs]
def slope(a, rotate=1):
r"""
Return either ``1`` (positive slope) or ``-1`` (negative slope).
If ``rotate`` is set to 1 then consider the edge as if it was rotated counterclockwise
infinitesimally.
EXAMPLES::
sage: from flatsurf.geometry.euclidean import slope
sage: slope((1, 1))
1
sage: slope((-1, 1))
-1
sage: slope((-1, -1))
1
sage: slope((1, -1))
-1
sage: slope((1, 0))
1
sage: slope((0, 1))
-1
sage: slope((-1, 0))
1
sage: slope((0, -1))
-1
sage: slope((1, 0), rotate=-1)
-1
sage: slope((0, 1), rotate=-1)
1
sage: slope((-1, 0), rotate=-1)
-1
sage: slope((0, -1), rotate=-1)
1
sage: slope((1, 0), rotate=0)
0
sage: slope((0, 1), rotate=0)
0
sage: slope((-1, 0), rotate=0)
0
sage: slope((0, -1), rotate=0)
0
sage: slope((0, 0))
Traceback (most recent call last):
...
ValueError: zero vector
"""
x, y = a
if not x and not y:
raise ValueError("zero vector")
if (x > 0 and y > 0) or (x < 0 and y < 0):
return 1
elif (x > 0 and y < 0) or (x < 0 and y > 0):
return -1
if rotate == 0:
return 0
if rotate == 1:
return 1 if x else -1
if rotate == -1:
return 1 if y else -1
raise ValueError("invalid argument rotate={}".format(rotate))