arseniiv/continued_fraction.py

## continued_fraction.py
"""
Reference implementation of continued fraction routines
by 0.5°.
"""

from __future__ import annotations
from dataclasses import InitVar, dataclass
from itertools import pairwise
from math import floor, inf, isinf, isnan, isqrt, lcm, sqrt
from typing import Iterable, Iterator, Literal, Sequence, overload
from fractions import Fraction
from more_itertools import last

__all__ = ('Infinity', 'Num', 'FracLin', 'Quadratic',
           'from_cf', 'cf', 'convergents', 'from_periodic_cf',)


@dataclass(frozen=True, slots=True)
class Infinity:
    """Projective infinity ∞. Useful for dealing with continued fractions."""
    def __mul__(self, other: Num) -> Infinity:
        if isinstance(other, int | Fraction | float):
            if other == 0 or isnan(other):
                raise ZeroDivisionError(f'can’t multiply ∞ and {other}')
            return self
        if isinstance(other, Infinity):
            return self
        return NotImplemented

    def __rmul__(self, other: int | Fraction | float) -> Infinity:
        if isinstance(other, int | Fraction | float):
            if other == 0 or isnan(other):
                raise ZeroDivisionError(f'can’t multiply {other} and ∞')
            return self
        return NotImplemented

    def __truediv__(self, other: Num) -> Infinity:
        if isinstance(other, int | Fraction | float):
            if isinf(other) or isnan(other):
                raise ZeroDivisionError(f'can’t divide ∞ by {other}')
            return self
        if isinstance(other, Infinity):
            raise ZeroDivisionError('can’t divide ∞ by ∞')
        return NotImplemented

    def __rtruediv__(self, other: int | Fraction | float) -> int:
        if isinstance(other, int | Fraction | float):
            if isinf(other) or isnan(other):
                raise ZeroDivisionError(f'can’t divide {other} by ∞')
            return 0
        return NotImplemented

    def __float__(self) -> float:
        return inf


type Num = int | float | Fraction | Infinity


def from_cf(f: Iterable[int]) -> Fraction | Infinity:
    return last(convergents(f), Infinity())


def cf(x: Num | Quadratic) -> Iterator[int]:
    if isinstance(x, Infinity) or isinf(x):
        return
    if isinstance(x, int | float):
        x = Fraction(x)
    if isinstance(x, Quadratic):
        head, period = x.periodic_cf()
        yield from head
        while True:
            yield from period

    while True:
        a, frac = divmod(x, 1)
        yield a
        if not frac:
            break
        x = 1 / frac


def convergents(f: Iterable[int], semi: bool = False) -> Iterator[Fraction]:
    m_prev, m, n_prev, n = 0, 1, 1, 0
    for a in f:
        if semi:
            for ai in range(1, a):  # 1 to a-1
                yield Fraction(ai * m + m_prev, ai * n + n_prev)
        m_prev, m = m, (a * m + m_prev)
        n_prev, n = n, (a * n + n_prev)
        yield Fraction(m, n)


@dataclass(frozen=True, slots=True)
class FracLin:
    """Integer fractional linear transformation x ↦ (a x + b) / (c x + d)."""
    a: int
    b: int
    c: int
    d: int

    def __post_init__(self) -> None:
        if abs(self.determinant) != 1:
            raise ValueError('a * d - b * c should be ±1')

    @property
    def determinant(self) -> Literal[1, -1]:
        # asserted in __post_init__
        return self.a * self.d - self.b * self.c # type: ignore

    @staticmethod
    def from_cf(f: Iterable[int]) -> FracLin:
        a, b, c, d = 1, 0, 0, 1
        for n in f:
            # for a single cf term: n + 1/x = (n * x + 1) / (1 * x + 0)
            # so (a, b, c, d) = (n, 1, 1, 0); composing with this gives:
            a, b, c, d = (a * n + b), a, (c * n + d), c
        return FracLin(a, b, c, d)

    @staticmethod
    def id() -> FracLin:
        return FracLin(1, 0, 0, 1)

    def __mul__(self, other: FracLin) -> FracLin:
        if isinstance(other, FracLin):
            a, b, c, d = self.a, self.b, self.c, self.d
            a2, b2, c2, d2 = other.a, other.b, other.c, other.d
            return FracLin(a * a2 + b * c2,
                           a * b2 + b * d2,
                           c * a2 + d * c2,
                           c * b2 + d * d2)
        return NotImplemented

    @property
    def inverse(self) -> FracLin:
        if self.determinant == 1:
            return FracLin(self.d, -self.b, -self.c, self.a)
        return FracLin(-self.d, self.b, self.c, -self.a)

    @property
    def final(self) -> Fraction | Infinity:
        """This transformation as a value of cf; equivalently f.final = f(inf)."""
        if not self.c:
            return Infinity()
        return Fraction(self.a, self.c)

    @overload
    def __call__(self, arg: float) -> float | Infinity: ...
    @overload
    def __call__(self, arg: int | Fraction | Infinity) -> Fraction | Infinity: ...
    @overload
    def __call__(self, arg: Quadratic) -> Quadratic | Infinity: ...
    def __call__(self, arg):
        if isinstance(arg, Quadratic):
            return arg.apply_fraclin(self)
        if isinstance(arg, Fraction | int):
            arg_n, arg_d = arg.numerator, arg.denominator
        elif isinstance(arg, Infinity) or isinf(arg):
            arg_n, arg_d = 1, 0
        else: # isinstance(arg, float)
            arg_n, arg_d = arg, 1
        numer = self.a * arg_n + self.b * arg_d
        denom = self.c * arg_n + self.d * arg_d
        if isinstance(numer, int):
            numer = Fraction(numer)
        try:
            return numer / denom
        except ZeroDivisionError:
            return Infinity()

    def __str__(self) -> str:
        return f'(({self.a} {self.b}) ({self.c} {self.d}))'


def _rational_sqrt(x: Fraction) -> Fraction | None:
    n, d = x.as_integer_ratio()
    if n == isqrt(n) ** 2 and d == isqrt(d) ** 2:
        return Fraction(n, d)
    return None

def _compare_sqrt_with_rat(s: Fraction, x: Fraction) -> Literal[-1, 0, 1]:
    # sqrt(s) > x
    #   [when x ≥ 0]  s > x**2
    #   [when x < 0]  True
    # sqrt(s) < x
    #   [when x ≥ 0]  s < x**2
    #   [when x < 0]  False
    if s < 0:
        raise ValueError('s should be nonnegative')
    if x < 0 or s > x ** 2:
        return 1
    if s < x ** 2:
        return -1
    return 0


@dataclass(frozen=True, slots=True)
class Quadratic:
    """Quadratic irrationality q + (1 if plus else -1) * sqrt(s)."""
    q: Fraction
    plus: bool
    s: Fraction
    _nocheck: InitVar[bool] = False

    @staticmethod
    def make(q: Fraction, plus: bool, s: Fraction) -> Fraction | Quadratic:
        if s < 0:
            raise ValueError('sqrt-part should be nonnegative')
        sqrt_s = _rational_sqrt(s)
        if sqrt_s is not None:
            return q + sqrt_s
        return Quadratic(q, plus, s, _nocheck=True)

    def __str__(self) -> str:
        addition = '+' if self.plus else '-'
        if self.s.is_integer():
            radicand = str(self.s)
        else:
            radicand = f'({self.s})'
        if not self.q:
            return f'{addition}√{radicand}'
        return f'{self.q} {addition} √{radicand}'

    def __post_init__(self, _nocheck: bool) -> None:
        if _nocheck:
            return
        if self.s <= 0:
            raise ValueError('sqrt-part should be positive')
        if _rational_sqrt(self.s) is not None:
            raise ValueError('sqrt-part should not be a square')

    def __float__(self) -> float:
        return self.q + (1 if self.plus else -1) * sqrt(self.s)

    def __gt__(self, other: int | Fraction) -> bool:
        if not isinstance(other, int | Fraction):
            return NotImplemented
        # q + sqrt(s) > x <=> sqrt(s) > x - q
        # q - sqrt(s) > x <=> sqrt(s) < q - x
        sign = 1 if self.plus else -1
        return _compare_sqrt_with_rat(self.s, (other - self.q) * sign) == sign

    def __lt__(self, other: int | Fraction) -> bool:
        if not isinstance(other, int | Fraction):
            return NotImplemented
        # q + sqrt(s) < x <=> sqrt(s) < x - q
        # q - sqrt(s) < x <=> sqrt(s) > q - x
        sign = 1 if self.plus else -1
        return _compare_sqrt_with_rat(self.s, (other - self.q) * sign) == -sign

    def divmod1(self) -> tuple[int, Quadratic]:
        m = lcm(self.q.denominator, self.s.denominator)
        q_ = m * self.q
        s_ = m**2 * self.s
        assert q_.denominator == 1 == s_.denominator
        floor_sqrt = isqrt(s_.numerator)
        if not self.plus:
            # note that floor(-x) = -ceiling(x) and isqrt(n) = floor(sqrt(n))
            if floor_sqrt ** 2 < s_.numerator:
                floor_sqrt += 1  # make it a ceiling
            floor_sqrt = -floor_sqrt
        floor_ = floor((q_ + floor_sqrt) / m)
        frac_ = Quadratic(self.q - floor_, self.plus, self.s, _nocheck=True)
        assert 0 < frac_ < 1
        return floor_, frac_

    @property
    def inverse(self) -> Quadratic:
        # (q ± sqrt(s)) (q ∓ sqrt(s)) = q**2 - s
        # (q ± sqrt(s))⁻¹ = (q ∓ sqrt(s)) / (q**2 - s) =
        #                 = (q / d) ∓ (sgn d) sqrt(s / d**2)
        d = self.q ** 2 - self.s
        new_plus = self.plus ^ (d > 0)
        assert d != 0, 'you shouldn’t have used unchecked construction'
        return Quadratic(self.q / d, new_plus, self.s / (d ** 2),
                         _nocheck=True)

    def periodic_cf(self) -> tuple[tuple[int, ...], tuple[int, ...]]:
        x = self
        terms: list[int] = []
        visited: list[Quadratic] = []
        while x not in visited:
            visited.append(x)
            a, frac = x.divmod1()
            terms.append(a)
            x = frac.inverse

        #                   |    -- cycle --    |
        # terms:   [ t1, t2, t3, t4, t5, t6, t7] (t3)
        # visited: [a,  b,  c,  d,  e,  f,  g] (c)
        i = visited.index(x)
        return tuple(terms[:i]), tuple(terms[i:])

    def apply_fraclin(self, f: FracLin) -> Quadratic | Infinity:
        # (a q + a sqrt(s) + b) / (c q + c sqrt(s) + d)
        if f.c == 0 == f.d:
            # s isn’t a square => denominator == 0 only in this case
            return Infinity()
        # num = (aq+b + a sqrt(s)) (cq+d - c sqrt(s)) =
        #     = (aq+b cq+d - a c s) + (a cq+d - c aq+b) sqrt(s)
        # NB (a cq+d - c aq+b) = a c q + a d - a c q - b c = a d - b c ≠ 0
        # den = (cq+d + c sqrt(s)) (cq+d - c sqrt(s)) =
        #     = (cq+d)**2 - c**2 s
        aq_b = f.a * self.q + f.b
        cq_d = f.c * self.q + f.d
        denominator = cq_d ** 2 - self.s * f.c ** 2
        new_q = (aq_b * cq_d - f.a * f.c * self.s) / denominator
        s_multiplier = f.determinant / denominator
        new_plus = self.plus if s_multiplier > 0 else not self.plus
        new_s = self.s * abs(s_multiplier) ** 2
        return Quadratic(new_q, new_plus, new_s, _nocheck=True)


def from_periodic_cf(head: Sequence[int], period: Sequence[int]) -> Fraction | Quadratic | Infinity:
    if not period:
        return from_cf(head)
    if not any(period):  # only zeros
        raise ValueError('period should not have only zeros')
    if head:
        # CF[head, y...] = FracLin.from_cf(head)( CF[y...] )
        return FracLin.from_cf(head)(from_periodic_cf((), period))
    assert not head

    # y = [period, y] = FracLin.from_cf(period)(y) =: A(y)
    # [y 1] is a eigenvector of A =: [(a b) (c d)]
    # ------------------
    # [(a b) (c d)] [z w] = λ [z w]
    # { a z + b w = λ z
    # { c z + d w = λ w
    # ...let's use wolframalpha and add careful checks about infinity and det±1
    # λ_1 = 1/2 (-sqrt(a^2 - 2 a d + 4 b c + d^2) + a + d)
    # v_1 = (-(-a + d + sqrt(a^2 + 4 b c - 2 a d + d^2)) / (2 c), 1)
    # λ_2 = 1/2 (sqrt(a^2 - 2 a d + 4 b c + d^2) + a + d)
    # v_2 = (-(-a + d - sqrt(a^2 + 4 b c - 2 a d + d^2)) / (2 c), 1)
    # -------------
    # let's simplify it
    # D = sqrt(D'); D' = (a - d)^2 + 4 b c
    # λ1 = (a + d - D) / 2
    # v1 = ((a - d) - D) / 2 c
    # λ2 = (a + d + D) / 2
    # v2 = ((a - d) + D) / 2 c
    # ------------
    # :::none of λ may end up being zero:
    # λ1 = 0 <=> (a + d) ≥ 0 and det = 0 (shouldn't happen)
    # λ2 = 0 <=> (a + d) ≤ 0 and det = 0 (shouldn't happen)
    # -------------
    # > Galois proved that the regular continued fraction which represents a
    # > quadratic surd ζ is purely periodic if and only if ζ is a reduced surd.
    # > Surd v1 is reduced := v1 > 1 and -1 < v2 < 0
    # so we can choose them this way... again if we can compare them!!
    # -----------------
    # :::what if D' = 0, is the c=0 case included here?
    # then λ ~ (a + d) and v = F(a - d, 2 c) if c ≠ 0
    # :::but what if a + d = 0 ? can that happen?
    # then a d - b c = - a a - b c
    # but we already have D' = 0 so a a + b c = 0 so our det is zero!!
    # :::show by manual solution c = 0 ==> v = ∞  [surprise: not always!]
    # if c = 0: (d - λ) w = 0 so either w = 0 or λ = d
    # then (a - d) z = -b w and w is arbitrary, well but
    # z/w = -b / (a - d)
    # so we can return just a normal fraction and not Infinity if a ≠ d !!!!!

    f = FracLin.from_cf(period)
    if not f.c:
        if f.a == f.d:
            return Infinity()
        return Fraction(-f.b, f.a - f.d)
    disc = (f.a - f.d) ** 2 + 4 * f.b * f.c
    q = Fraction(f.a - f.d, 2 * f.c)
    if not disc:
        return q
    s = Fraction(disc, 4 * f.c ** 2)
    y1 = Quadratic(q, True, s)
    y2 = Quadratic(q, False, s)
    if y1 > 1:
        assert -1 < y2 < 0
        return y1
    if y2 > 1:
        assert -1 < y1 < 0
        return y2
    assert False, 'shouldn’t happen'


def _main() -> None:
    if False:
        x1 = Fraction('46_748/51_745')
        semis = list(convergents(cf(x1), semi=True))
        for c1, c2 in pairwise(semis):
            print(f'{c1} and {c2}:')
            n1, d1 = c1.numerator, c1.denominator
            n2, d2 = c2.numerator, c2.denominator
            print(f'{n1 * d2} - {n2 * d1} = {n1 * d2 - n2 * d1}')

    if False:
        A, B = [3, 1], [1, 2, 5, 2]
        trans_ab = FracLin.from_cf(A + B)
        assert trans_ab == FracLin.from_cf(A) * FracLin.from_cf(B)
        assert trans_ab.final == from_cf(A + B)
        assert FracLin.from_cf([]).final == from_cf([]) == Infinity()
        assert trans_ab.inverse * trans_ab == FracLin.id()
        assert ((trans_ab.final * trans_ab.inverse(0)) /
                (trans_ab(0) * trans_ab.inverse.final) == 1)
        # the last is not always true though because we calculate
        # (a/c) (−b/a) / (b/d) (−d/c), so if any of a b c d is zero then boom
        # BUT it’s true we can determine each of a b c d having all the others

    if True:
        from math import isclose
        from itertools import islice
        phi_float = (sqrt(5) + 1) / 2
        phi = Quadratic(Fraction(1, 2), True, Fraction(5, 4))
        assert isclose(phi, phi_float)
        assert 1 < phi < 2
        inv_phi = Quadratic(Fraction(-1, 2), True, Fraction(5, 4))
        assert 0 < inv_phi < 1
        assert str(inv_phi) == '-1/2 + √(5/4)'
        assert phi.inverse == inv_phi
        assert inv_phi.inverse == phi
        assert phi.periodic_cf() == ((), (1,))
        assert inv_phi.periodic_cf() == ((0,), (1,))
        x = Quadratic(Fraction(11, 5), False, Fraction(13, 7))
        x_cf = ((0, 1, 5),
                (6, 1, 25, 1, 1, 1, 3, 4, 1, 1, 1, 4, 3, 1, 1, 1, 25,
                 1, 6, 6, 2, 1, 1, 3, 2, 14, 68, 14, 2, 3, 1, 1, 2, 6))
        assert isclose(x, (11/5 - sqrt(13/7)))
        assert x.apply_fraclin(FracLin(0, 1, 1, 0)) == x.inverse
        assert x.apply_fraclin(FracLin.id()) == x
        assert x.inverse.inverse == x
        assert x.periodic_cf() == x_cf
        assert from_periodic_cf(*x_cf) == x
        assert isclose(x, from_cf(islice(cf(x), 15)))
        assert from_periodic_cf(*x.periodic_cf()) == x
        y_cf =        ((1, 5, 1), (3, 2, 1, 1))
        y_cf_normal = ((1, 5), (1, 3, 2, 1))
        y = from_periodic_cf(*y_cf)
        assert isinstance(y, Quadratic)
        assert y.periodic_cf() == y_cf_normal
        assert from_periodic_cf(*y_cf_normal) == y
        assert from_periodic_cf((), ()) == from_cf([])
        assert from_periodic_cf((1, 1, 1), ()) == from_cf([1, 1, 1])
        silver = from_periodic_cf((), (2,))
        assert isclose(silver, from_cf(islice(cf(silver), 15)))
        assert str(silver) == '1 + √2'
        assert from_periodic_cf((), (2, 2, 2)) == silver
        assert from_periodic_cf((2, 2), (2, 2)) == silver


if __name__ == '__main__':
    _main()
	"""
	Reference implementation of continued fraction routines
	by 0.5°.
	"""

	from __future__ import annotations
	from dataclasses import InitVar, dataclass
	from itertools import pairwise
	from math import floor, inf, isinf, isnan, isqrt, lcm, sqrt
	from typing import Iterable, Iterator, Literal, Sequence, overload
	from fractions import Fraction
	from more_itertools import last

	__all__ = ('Infinity', 'Num', 'FracLin', 'Quadratic',
	'from_cf', 'cf', 'convergents', 'from_periodic_cf',)


	@dataclass(frozen=True, slots=True)
	class Infinity:
	"""Projective infinity ∞. Useful for dealing with continued fractions."""
	def __mul__(self, other: Num) -> Infinity:
	if isinstance(other, int \| Fraction \| float):
	if other == 0 or isnan(other):
	raise ZeroDivisionError(f'can’t multiply ∞ and {other}')
	return self
	if isinstance(other, Infinity):
	return self
	return NotImplemented

	def __rmul__(self, other: int \| Fraction \| float) -> Infinity:
	if isinstance(other, int \| Fraction \| float):
	if other == 0 or isnan(other):
	raise ZeroDivisionError(f'can’t multiply {other} and ∞')
	return self
	return NotImplemented

	def __truediv__(self, other: Num) -> Infinity:
	if isinstance(other, int \| Fraction \| float):
	if isinf(other) or isnan(other):
	raise ZeroDivisionError(f'can’t divide ∞ by {other}')
	return self
	if isinstance(other, Infinity):
	raise ZeroDivisionError('can’t divide ∞ by ∞')
	return NotImplemented

	def __rtruediv__(self, other: int \| Fraction \| float) -> int:
	if isinstance(other, int \| Fraction \| float):
	if isinf(other) or isnan(other):
	raise ZeroDivisionError(f'can’t divide {other} by ∞')
	return 0
	return NotImplemented

	def __float__(self) -> float:
	return inf


	type Num = int \| float \| Fraction \| Infinity


	def from_cf(f: Iterable[int]) -> Fraction \| Infinity:
	return last(convergents(f), Infinity())


	def cf(x: Num \| Quadratic) -> Iterator[int]:
	if isinstance(x, Infinity) or isinf(x):
	return
	if isinstance(x, int \| float):
	x = Fraction(x)
	if isinstance(x, Quadratic):
	head, period = x.periodic_cf()
	yield from head
	while True:
	yield from period

	while True:
	a, frac = divmod(x, 1)
	yield a
	if not frac:
	break
	x = 1 / frac


	def convergents(f: Iterable[int], semi: bool = False) -> Iterator[Fraction]:
	m_prev, m, n_prev, n = 0, 1, 1, 0
	for a in f:
	if semi:
	for ai in range(1, a): # 1 to a-1
	yield Fraction(ai * m + m_prev, ai * n + n_prev)
	m_prev, m = m, (a * m + m_prev)
	n_prev, n = n, (a * n + n_prev)
	yield Fraction(m, n)


	@dataclass(frozen=True, slots=True)
	class FracLin:
	"""Integer fractional linear transformation x ↦ (a x + b) / (c x + d)."""
	a: int
	b: int
	c: int
	d: int

	def __post_init__(self) -> None:
	if abs(self.determinant) != 1:
	raise ValueError('a * d - b * c should be ±1')

	@property
	def determinant(self) -> Literal[1, -1]:
	# asserted in __post_init__
	return self.a * self.d - self.b * self.c # type: ignore

	@staticmethod
	def from_cf(f: Iterable[int]) -> FracLin:
	a, b, c, d = 1, 0, 0, 1
	for n in f:
	# for a single cf term: n + 1/x = (n * x + 1) / (1 * x + 0)
	# so (a, b, c, d) = (n, 1, 1, 0); composing with this gives:
	a, b, c, d = (a * n + b), a, (c * n + d), c
	return FracLin(a, b, c, d)

	@staticmethod
	def id() -> FracLin:
	return FracLin(1, 0, 0, 1)

	def __mul__(self, other: FracLin) -> FracLin:
	if isinstance(other, FracLin):
	a, b, c, d = self.a, self.b, self.c, self.d
	a2, b2, c2, d2 = other.a, other.b, other.c, other.d
	return FracLin(a * a2 + b * c2,
	a * b2 + b * d2,
	c * a2 + d * c2,
	c * b2 + d * d2)
	return NotImplemented

	@property
	def inverse(self) -> FracLin:
	if self.determinant == 1:
	return FracLin(self.d, -self.b, -self.c, self.a)
	return FracLin(-self.d, self.b, self.c, -self.a)

	@property
	def final(self) -> Fraction \| Infinity:
	"""This transformation as a value of cf; equivalently f.final = f(inf)."""
	if not self.c:
	return Infinity()
	return Fraction(self.a, self.c)

	@overload
	def __call__(self, arg: float) -> float \| Infinity: ...
	@overload
	def __call__(self, arg: int \| Fraction \| Infinity) -> Fraction \| Infinity: ...
	@overload
	def __call__(self, arg: Quadratic) -> Quadratic \| Infinity: ...
	def __call__(self, arg):
	if isinstance(arg, Quadratic):
	return arg.apply_fraclin(self)
	if isinstance(arg, Fraction \| int):
	arg_n, arg_d = arg.numerator, arg.denominator
	elif isinstance(arg, Infinity) or isinf(arg):
	arg_n, arg_d = 1, 0
	else: # isinstance(arg, float)
	arg_n, arg_d = arg, 1
	numer = self.a * arg_n + self.b * arg_d
	denom = self.c * arg_n + self.d * arg_d
	if isinstance(numer, int):
	numer = Fraction(numer)
	try:
	return numer / denom
	except ZeroDivisionError:
	return Infinity()

	def __str__(self) -> str:
	return f'(({self.a} {self.b}) ({self.c} {self.d}))'


	def _rational_sqrt(x: Fraction) -> Fraction \| None:
	n, d = x.as_integer_ratio()
	if n == isqrt(n) 2 and d == isqrt(d) 2:
	return Fraction(n, d)
	return None

	def _compare_sqrt_with_rat(s: Fraction, x: Fraction) -> Literal[-1, 0, 1]:
	# sqrt(s) > x
	# [when x ≥ 0] s > x**2
	# [when x < 0] True
	# sqrt(s) < x
	# [when x ≥ 0] s < x**2
	# [when x < 0] False
	if s < 0:
	raise ValueError('s should be nonnegative')
	if x < 0 or s > x ** 2:
	return 1
	if s < x ** 2:
	return -1
	return 0


	@dataclass(frozen=True, slots=True)
	class Quadratic:
	"""Quadratic irrationality q + (1 if plus else -1) * sqrt(s)."""
	q: Fraction
	plus: bool
	s: Fraction
	_nocheck: InitVar[bool] = False

	@staticmethod
	def make(q: Fraction, plus: bool, s: Fraction) -> Fraction \| Quadratic:
	if s < 0:
	raise ValueError('sqrt-part should be nonnegative')
	sqrt_s = _rational_sqrt(s)
	if sqrt_s is not None:
	return q + sqrt_s
	return Quadratic(q, plus, s, _nocheck=True)

	def __str__(self) -> str:
	addition = '+' if self.plus else '-'
	if self.s.is_integer():
	radicand = str(self.s)
	else:
	radicand = f'({self.s})'
	if not self.q:
	return f'{addition}√{radicand}'
	return f'{self.q} {addition} √{radicand}'

	def __post_init__(self, _nocheck: bool) -> None:
	if _nocheck:
	return
	if self.s <= 0:
	raise ValueError('sqrt-part should be positive')
	if _rational_sqrt(self.s) is not None:
	raise ValueError('sqrt-part should not be a square')

	def __float__(self) -> float:
	return self.q + (1 if self.plus else -1) * sqrt(self.s)

	def __gt__(self, other: int \| Fraction) -> bool:
	if not isinstance(other, int \| Fraction):
	return NotImplemented
	# q + sqrt(s) > x <=> sqrt(s) > x - q
	# q - sqrt(s) > x <=> sqrt(s) < q - x
	sign = 1 if self.plus else -1
	return _compare_sqrt_with_rat(self.s, (other - self.q) * sign) == sign

	def __lt__(self, other: int \| Fraction) -> bool:
	if not isinstance(other, int \| Fraction):
	return NotImplemented
	# q + sqrt(s) < x <=> sqrt(s) < x - q
	# q - sqrt(s) < x <=> sqrt(s) > q - x
	sign = 1 if self.plus else -1
	return _compare_sqrt_with_rat(self.s, (other - self.q) * sign) == -sign

	def divmod1(self) -> tuple[int, Quadratic]:
	m = lcm(self.q.denominator, self.s.denominator)
	q_ = m * self.q
	s_ = m*2 self.s
	assert q_.denominator == 1 == s_.denominator
	floor_sqrt = isqrt(s_.numerator)
	if not self.plus:
	# note that floor(-x) = -ceiling(x) and isqrt(n) = floor(sqrt(n))
	if floor_sqrt ** 2 < s_.numerator:
	floor_sqrt += 1 # make it a ceiling
	floor_sqrt = -floor_sqrt
	floor_ = floor((q_ + floor_sqrt) / m)
	frac_ = Quadratic(self.q - floor_, self.plus, self.s, _nocheck=True)
	assert 0 < frac_ < 1
	return floor_, frac_

	@property
	def inverse(self) -> Quadratic:
	# (q ± sqrt(s)) (q ∓ sqrt(s)) = q**2 - s
	# (q ± sqrt(s))⁻¹ = (q ∓ sqrt(s)) / (q**2 - s) =
	# = (q / d) ∓ (sgn d) sqrt(s / d**2)
	d = self.q ** 2 - self.s
	new_plus = self.plus ^ (d > 0)
	assert d != 0, 'you shouldn’t have used unchecked construction'
	return Quadratic(self.q / d, new_plus, self.s / (d ** 2),
	_nocheck=True)

	def periodic_cf(self) -> tuple[tuple[int, ...], tuple[int, ...]]:
	x = self
	terms: list[int] = []
	visited: list[Quadratic] = []
	while x not in visited:
	visited.append(x)
	a, frac = x.divmod1()
	terms.append(a)
	x = frac.inverse

	# \| -- cycle -- \|
	# terms: [ t1, t2, t3, t4, t5, t6, t7] (t3)
	# visited: [a, b, c, d, e, f, g] (c)
	i = visited.index(x)
	return tuple(terms[:i]), tuple(terms[i:])

	def apply_fraclin(self, f: FracLin) -> Quadratic \| Infinity:
	# (a q + a sqrt(s) + b) / (c q + c sqrt(s) + d)
	if f.c == 0 == f.d:
	# s isn’t a square => denominator == 0 only in this case
	return Infinity()
	# num = (aq+b + a sqrt(s)) (cq+d - c sqrt(s)) =
	# = (aq+b cq+d - a c s) + (a cq+d - c aq+b) sqrt(s)
	# NB (a cq+d - c aq+b) = a c q + a d - a c q - b c = a d - b c ≠ 0
	# den = (cq+d + c sqrt(s)) (cq+d - c sqrt(s)) =
	# = (cq+d)2 - c2 s
	aq_b = f.a * self.q + f.b
	cq_d = f.c * self.q + f.d
	denominator = cq_d ** 2 - self.s * f.c ** 2
	new_q = (aq_b * cq_d - f.a * f.c * self.s) / denominator
	s_multiplier = f.determinant / denominator
	new_plus = self.plus if s_multiplier > 0 else not self.plus
	new_s = self.s * abs(s_multiplier) ** 2
	return Quadratic(new_q, new_plus, new_s, _nocheck=True)


	def from_periodic_cf(head: Sequence[int], period: Sequence[int]) -> Fraction \| Quadratic \| Infinity:
	if not period:
	return from_cf(head)
	if not any(period): # only zeros
	raise ValueError('period should not have only zeros')
	if head:
	# CF[head, y...] = FracLin.from_cf(head)( CF[y...] )
	return FracLin.from_cf(head)(from_periodic_cf((), period))
	assert not head

	# y = [period, y] = FracLin.from_cf(period)(y) =: A(y)
	# [y 1] is a eigenvector of A =: [(a b) (c d)]
	# ------------------
	# [(a b) (c d)] [z w] = λ [z w]
	# { a z + b w = λ z
	# { c z + d w = λ w
	# ...let's use wolframalpha and add careful checks about infinity and det±1
	# λ_1 = 1/2 (-sqrt(a^2 - 2 a d + 4 b c + d^2) + a + d)
	# v_1 = (-(-a + d + sqrt(a^2 + 4 b c - 2 a d + d^2)) / (2 c), 1)
	# λ_2 = 1/2 (sqrt(a^2 - 2 a d + 4 b c + d^2) + a + d)
	# v_2 = (-(-a + d - sqrt(a^2 + 4 b c - 2 a d + d^2)) / (2 c), 1)
	# -------------
	# let's simplify it
	# D = sqrt(D'); D' = (a - d)^2 + 4 b c
	# λ1 = (a + d - D) / 2
	# v1 = ((a - d) - D) / 2 c
	# λ2 = (a + d + D) / 2
	# v2 = ((a - d) + D) / 2 c
	# ------------
	# :::none of λ may end up being zero:
	# λ1 = 0 <=> (a + d) ≥ 0 and det = 0 (shouldn't happen)
	# λ2 = 0 <=> (a + d) ≤ 0 and det = 0 (shouldn't happen)
	# -------------
	# > Galois proved that the regular continued fraction which represents a
	# > quadratic surd ζ is purely periodic if and only if ζ is a reduced surd.
	# > Surd v1 is reduced := v1 > 1 and -1 < v2 < 0
	# so we can choose them this way... again if we can compare them!!
	# -----------------
	# :::what if D' = 0, is the c=0 case included here?
	# then λ ~ (a + d) and v = F(a - d, 2 c) if c ≠ 0
	# :::but what if a + d = 0 ? can that happen?
	# then a d - b c = - a a - b c
	# but we already have D' = 0 so a a + b c = 0 so our det is zero!!
	# :::show by manual solution c = 0 ==> v = ∞ [surprise: not always!]
	# if c = 0: (d - λ) w = 0 so either w = 0 or λ = d
	# then (a - d) z = -b w and w is arbitrary, well but
	# z/w = -b / (a - d)
	# so we can return just a normal fraction and not Infinity if a ≠ d !!!!!

	f = FracLin.from_cf(period)
	if not f.c:
	if f.a == f.d:
	return Infinity()
	return Fraction(-f.b, f.a - f.d)
	disc = (f.a - f.d) ** 2 + 4 * f.b * f.c
	q = Fraction(f.a - f.d, 2 * f.c)
	if not disc:
	return q
	s = Fraction(disc, 4 * f.c ** 2)
	y1 = Quadratic(q, True, s)
	y2 = Quadratic(q, False, s)
	if y1 > 1:
	assert -1 < y2 < 0
	return y1
	if y2 > 1:
	assert -1 < y1 < 0
	return y2
	assert False, 'shouldn’t happen'


	def _main() -> None:
	if False:
	x1 = Fraction('46_748/51_745')
	semis = list(convergents(cf(x1), semi=True))
	for c1, c2 in pairwise(semis):
	print(f'{c1} and {c2}:')
	n1, d1 = c1.numerator, c1.denominator
	n2, d2 = c2.numerator, c2.denominator
	print(f'{n1 * d2} - {n2 * d1} = {n1 * d2 - n2 * d1}')

	if False:
	A, B = [3, 1], [1, 2, 5, 2]
	trans_ab = FracLin.from_cf(A + B)
	assert trans_ab == FracLin.from_cf(A) * FracLin.from_cf(B)
	assert trans_ab.final == from_cf(A + B)
	assert FracLin.from_cf([]).final == from_cf([]) == Infinity()
	assert trans_ab.inverse * trans_ab == FracLin.id()
	assert ((trans_ab.final * trans_ab.inverse(0)) /
	(trans_ab(0) * trans_ab.inverse.final) == 1)
	# the last is not always true though because we calculate
	# (a/c) (−b/a) / (b/d) (−d/c), so if any of a b c d is zero then boom
	# BUT it’s true we can determine each of a b c d having all the others

	if True:
	from math import isclose
	from itertools import islice
	phi_float = (sqrt(5) + 1) / 2
	phi = Quadratic(Fraction(1, 2), True, Fraction(5, 4))
	assert isclose(phi, phi_float)
	assert 1 < phi < 2
	inv_phi = Quadratic(Fraction(-1, 2), True, Fraction(5, 4))
	assert 0 < inv_phi < 1
	assert str(inv_phi) == '-1/2 + √(5/4)'
	assert phi.inverse == inv_phi
	assert inv_phi.inverse == phi
	assert phi.periodic_cf() == ((), (1,))
	assert inv_phi.periodic_cf() == ((0,), (1,))
	x = Quadratic(Fraction(11, 5), False, Fraction(13, 7))
	x_cf = ((0, 1, 5),
	(6, 1, 25, 1, 1, 1, 3, 4, 1, 1, 1, 4, 3, 1, 1, 1, 25,
	1, 6, 6, 2, 1, 1, 3, 2, 14, 68, 14, 2, 3, 1, 1, 2, 6))
	assert isclose(x, (11/5 - sqrt(13/7)))
	assert x.apply_fraclin(FracLin(0, 1, 1, 0)) == x.inverse
	assert x.apply_fraclin(FracLin.id()) == x
	assert x.inverse.inverse == x
	assert x.periodic_cf() == x_cf
	assert from_periodic_cf(*x_cf) == x
	assert isclose(x, from_cf(islice(cf(x), 15)))
	assert from_periodic_cf(*x.periodic_cf()) == x
	y_cf = ((1, 5, 1), (3, 2, 1, 1))
	y_cf_normal = ((1, 5), (1, 3, 2, 1))
	y = from_periodic_cf(*y_cf)
	assert isinstance(y, Quadratic)
	assert y.periodic_cf() == y_cf_normal
	assert from_periodic_cf(*y_cf_normal) == y
	assert from_periodic_cf((), ()) == from_cf([])
	assert from_periodic_cf((1, 1, 1), ()) == from_cf([1, 1, 1])
	silver = from_periodic_cf((), (2,))
	assert isclose(silver, from_cf(islice(cf(silver), 15)))
	assert str(silver) == '1 + √2'
	assert from_periodic_cf((), (2, 2, 2)) == silver
	assert from_periodic_cf((2, 2), (2, 2)) == silver


	if __name__ == '__main__':
	_main()