from math import hypot, acos, cos, sin, pi, atan2

class Segment:
    def __init__(self, p1, p2):
        self.p1 = p1
        self.p2 = p2

    def length(self):
        return hypot(self.p1.x - self.p2.x, self.p1.y - self.p2.y)

    def __repr__(self):
        return f"Segment {self.p1} --- {self.p2}"

    def intersect(self, circle):
        """
        Alternative implementation using discriminant method for line-circle intersection
        More numerically stable for edge cases
        """
        # Translate so circle is at origin
        p1_translated = self.p1 - circle.center
        p2_translated = self.p2 - circle.center

        # Line segment parameterization: P(t) = p1 + t * (p2 - p1), t in [0,1]
        dx = p2_translated.x - p1_translated.x
        dy = p2_translated.y - p1_translated.y

        # Quadratic equation coefficients for |P(t)|² = r²
        a = dx * dx + dy * dy
        b = 2 * (p1_translated.x * dx + p1_translated.y * dy)
        c = p1_translated.x * p1_translated.x + p1_translated.y * p1_translated.y - circle.radius * circle.radius

        # Handle degenerate case: segment is a point
        if a < EPS * EPS:
            return abs(c) <= EPS

        # Calculate discriminant
        discriminant = b * b - 4 * a * c

        # No intersection if discriminant is negative
        if discriminant < -EPS:
            return False

        # Check if intersection points are within segment bounds [0, 1]
        if discriminant < EPS:  # One solution (tangent)
            t = -b / (2 * a)
            return -EPS <= t <= 1 + EPS
        else:  # Two solutions
            sqrt_discriminant = sqrt(discriminant)
            t1 = (-b - sqrt_discriminant) / (2 * a)
            t2 = (-b + sqrt_discriminant) / (2 * a)

            # Check if either intersection is within segment bounds
            return ((-EPS <= t1 <= 1 + EPS) or (-EPS <= t2 <= 1 + EPS) or 
                    (t1 < -EPS and t2 > 1 + EPS))


class Point:
    def __init__(self, x, y):
        self.x = x
        self.y = y

    def __str__(self):
        return f"Point({{{self.x}, {self.y}}})"
        return f"({self.x}, {self.y})"

    def __add__(self, other):
        return Point(self.x + other.x, self.y + other.y)

    def __sub__(self, other):
        return Point(self.x - other.x, self.y - other.y)

    def __neg__(self):
        return Point(-self.x, -self.y)

    def __rmul__(self, k):
        return Point(k * self.x, k * self.y)

    def __eq__(self, other):
        return self.x == other.x and self.y == other.y

    def __lt__(self, other):
        return (self.x, self.y) <= (other.x, other.y)


class Circle:
    def __init__(self, center, radius):
        self.center = center
        self.radius = radius

    def __repr__(self):
        return f"Circle(Point({{{self.center.x}, {self.center.y}}}), {self.radius})"
        return f"Circle of center {self.center} and radius {self.radius}"

    # finds the two points (if they exist) of the circle that have tangents passing through point p
    def tangents_toPoint(self, p):
        c, r = self.center, self.radius
        d = Segment(c, p).length()
        if r > d:
            return (False, [])
        u = max(-1, min(r / d, 1))
        t = acos(u)
        q = c + r / d * (p - c)
        p1 = Point(
            c.x + cos(t) * (q.x - c.x) - sin(t) * (q.y - c.y),
            c.y + sin(t) * (q.x - c.x) + cos(t) * (q.y - c.y),
        )
        p2 = Point(
            c.x + cos(t) * (q.x - c.x) + sin(t) * (q.y - c.y),
            c.y - sin(t) * (q.x - c.x) + cos(t) * (q.y - c.y),
        )
        return (True, [p1, p2])

    # Given two circles c1 and c2, finds (if they exist) the points p1, p2 in c1, q1 and q2 in c2 such that (p1, q2) and (p2, q1) are internal bi tangents
    def internal_bitangents(self, other):
        r1, r2 = self.radius, other.radius
        x1, y1 = self.center.x, self.center.y
        x2, y2 = other.center.x, other.center.y
        dis = hypot(x2 - x1, y2 - y1)
        if r1 + r2 >= dis:
            return (False, [])
        u = max(-1, min((r1 + r2) / dis, 1))
        t = acos(u)

        p = self.center + r1 / dis * (other.center - self.center)

        p1 = Point(
            x1 + cos(t) * (p.x - x1) - sin(t) * (p.y - y1),
            y1 + sin(t) * (p.x - x1) + cos(t) * (p.y - y1),
        )
        p2 = Point(
            x1 + cos(t) * (p.x - x1) + sin(t) * (p.y - y1),
            y1 - sin(t) * (p.x - x1) + cos(t) * (p.y - y1),
        )

        #########
        q = other.center + r2 / dis * (self.center - other.center)

        q1 = Point(
            x2 + cos(t) * (q.x - x2) - sin(t) * (q.y - y2),
            y2 + sin(t) * (q.x - x2) + cos(t) * (q.y - y2),
        )
        q2 = Point(
            x2 + cos(t) * (q.x - x2) + sin(t) * (q.y - y2),
            y2 - sin(t) * (q.x - x2) + cos(t) * (q.y - y2),
        )

        return (True, [p1, q1, p2, q2])

    # Given two circles c1 and c2, finds (if they exist) the points p1, p2 in c1, q1 and q2 in c2 such that (p1, q2) and (p2, q1) are external bi tangents

    def external_bitangents(self, other):
        r1, r2 = self.radius, other.radius
        x1, y1 = self.center.x, self.center.y
        x2, y2 = other.center.x, other.center.y
        dis = hypot(x2 - x1, y2 - y1)
        if abs(r1 - r2) > dis:
            return (False, [])
        u = max(-1, min((r1 - r2) / dis, 1))
        t = acos(u)

        p = self.center + r1 / dis * (other.center - self.center)

        p1 = Point(
            x1 + cos(t) * (p.x - x1) - sin(t) * (p.y - y1),
            y1 + sin(t) * (p.x - x1) + cos(t) * (p.y - y1),
        )
        p2 = Point(
            x1 + cos(t) * (p.x - x1) + sin(t) * (p.y - y1),
            y1 - sin(t) * (p.x - x1) + cos(t) * (p.y - y1),
        )

        #########
        q = self.center + (dis + r2) / dis * (other.center - self.center)
        q1 = Point(
            x2 + cos(t) * (q.x - x2) - sin(t) * (q.y - y2),
            y2 + sin(t) * (q.x - x2) + cos(t) * (q.y - y2),
        )
        q2 = Point(
            x2 + cos(t) * (q.x - x2) + sin(t) * (q.y - y2),
            y2 - sin(t) * (q.x - x2) + cos(t) * (q.y - y2),
        )

        return (True, [p1, q1, p2, q2])

    # gets the intersection with another circle
    def intersect(self, other):
        assert type(other) == Circle, "Element is not a circle"
        ca, cb = self.center, other.center
        ra, rb = self.radius, other.radius
        d = Segment(ca, cb).length()
        # No intersection
        if ra + rb < d or abs(ra - rb) > d:
            return (False, [])
        a = (ra**2 + d**2 - rb**2) / (2 * d)
        # 1 intersection point
        if abs(a - ra) <= 1e-12:
            c = ca + a / d * (cb - ca)
            return (False, [c])
        # 2 intersection points
        d = ca + ra / d * (cb - ca)
        u = max(-1, min(a / ra, 1))
        t = acos(u)
        p = Point(
            ca.x + cos(t) * (d.x - ca.x) - sin(t) * (d.y - ca.y),
            ca.y + sin(t) * (d.x - ca.x) + cos(t) * (d.y - ca.y),
        )
        q = Point(
            ca.x + cos(t) * (d.x - ca.x) + sin(t) * (d.y - ca.y),
            ca.y - sin(t) * (d.x - ca.x) + cos(t) * (d.y - ca.y),
        )
        return (True, [p, q])
        
# a, b = Point(-3, 1), Point(4.25, 0)
# cs = [Circle(Point(0, 0), 2.5), Circle(Point(1.5, 2), 0.5),
#       Circle(Point(3.5, 1), 1), Circle(Point(3.5, -1.7), 1.2)]
        
# for i in range(4):
#     for j in range(i + 1, 4):
#         b, r = cs[i].external_bitangents(cs[j])
#         print(i, j, b)
#         print('\n'.join(str(p) for p in r))

a = Point(3.1777024751372713,-1.0431819935673712)
c = Circle(Point(1.8478065823942857,-0.43572320663297526), 1.2977733144010024)
b, r = c.tangents_toPoint(a)
for p in r:
    print(p)

Embed on website

To embed this program on your website, copy the following code and paste it into your website's HTML: