Problem:
Cross Points of a Circe and a Line For given a circle $c$ and a line $l$, print the coordinates of the cross points of them. Input The input is given in the following format. $cx\; cy\; r$ $q$ $Line_1$ $Line_2$ : $Line_q$  In the first line, the center coordinate of the circle and its radius are given by $cx$, $cy$ and $r$. In the second line, the number of queries $q$ is given. In the following $q$ lines, as queries, $Line_i$ are given ($1 \leq i \leq q$) in the following format. $x_1\; y_1\; x_2\; y_2$ Each line is represented by two points $p1$ and $p2$ which the line crosses. The coordinate of $p1$ and $p2$ are given by ($x1$, $y1$) and ($x2$, $y2$) respectively. All input values are given in integers. Output For each query, print the coordinates of the cross points in the following rules.  If there is one cross point, print two coordinates with the same values.  Print the coordinate with smaller $x$ first. In case of a tie, print the coordinate with smaller $y$ first.  The output values should be in a decimal fraction with an error less than 0.000001. Constraints  $p1$ and $p2$ are different  The circle and line have at least one cross point  $1 \leq q \leq 1,000$  $-10,000 \leq cx, cy, x1, y1, x2, y2 \leq 10,000$  $1 \leq r \leq 10,000$  Sample Input and Output Sample Input 2 1 1
2
0 1 4 1
3 0 3 3 Sample Output 1.00000000 1.00000000 3.00000000 1.00000000
3.00000000 1.00000000 3.00000000 1.00000000

Solution (Python):
import math
def dot(a, b):
  return a.real * b.real + a.imag * b.imag

def projection(p0,p1,p2):
  a=p1-p0
  b=p2-p0
  pro = a * dot(a, b) / (abs(a) ** 2)
  t=p0+pro
  return t

def get_cross_point(p0,p1,p2):
  pro=projection(p0,p1,p2)
  e  =(p1-p0)/abs(p1-p0)
  base = (r**2 - abs(pro-p2)**2)**0.5
  if (pro-e*base).real == (pro+e*base).real:
    if (pro-e*base).imag < (pro+e*base).imag:
      ans1=pro-e*base
      ans2=pro+e*base
    else:
      ans1=pro+e*base
      ans2=pro-e*base
  elif (pro-e*base).real < (pro+e*base).real:
    ans1=pro-e*base
    ans2=pro+e*base
  else:
    ans1=pro+e*base
    ans2=pro-e*base
  return ans1.real, ans1.imag,ans2.real,ans2.imag

cx, cy, r = map(int, input().split())
p2=complex(cx,cy)
q = int(input())
for _ in range(q):
  x0,y0,x1,y1=map(int,input().split())
  p0=complex(x0,y0)
  p1=complex(x1,y1)
  print("{:.8f} {:.8f} {:.8f} {:.8f}".format(*get_cross_point(p0,p1,p2)))

Solution (Python):
from sys import stdin
import math

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

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

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

    def __mul__(self, k):
        return Vector(self.x * k, self.y * k)

    def __truediv__(self, k):
        return Vector(self.x / k, self.y / k)

    def __gt__(self, other):
        return self.x > other.x and self.y > other.yb

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

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

    def dot(self, other):
        return self.x * other.x + self.y * other.y

    
    def cross(self, other):
        return self.x * other.y - self.y * other.x

    def norm(self):
        return self.x * self.x + self.y * self.y

    def abs(self):
        return math.sqrt(self.norm())

class Point(Vector):
    def __init__(self, *args, **kargs):
        return super().__init__(*args, **kargs)

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

class Line(Segment):
    def __init__(self, *args, **kargs):
        return super().__init__(*args, **kargs)

class Circle:
    def __init__(self, c=Point(0, 0), r=1):
        self.c = c
        self.r = r

def project(s, p):
    base = s.p2 - s.p1
    hypo = p - s.p1
    r = hypo.dot(base) / base.norm()
    return s.p1 + base * r

def get_cross_point(c, l):
    pr = project(l, c.c)
    e = (l.p1- l.p2) / (l.p2 - l.p1).abs()
    base = math.sqrt( c.r * c.r - (pr - c.c).norm() )
    return pr + e * base, pr - e * base

def read_and_print_results(n, c):
    for _ in range(n):
        line = stdin.readline().strip().split()
        p0 = Vector(int(line[0]), int(line[1]))
        p1 = Vector(int(line[2]), int(line[3]))
        l = Line(p0, p1)
        p2, p3 = get_cross_point(c, l)
        x1, y1 = min((p2.x, p2.y), (p3.x, p3.y))
        x2, y2 = max((p2.x, p2.y), (p3.x, p3.y))
        print("{0:0.8f} {1:0.8f} {2:0.8f} {3:0.8f}".format(x1, y1, x2, y2))

x, y, r = [ int(i) for i in input().split() ]
c = Circle(Point(x, y), r)
n = int(input())
read_and_print_results(n, c)

Solution (Python):
from sys import stdin
readline = stdin.readline


def dot(a, b):
    return a.real * b.real + a.imag * b.imag


def unit(a):
    return a / abs(a)


def near_pos_online(p, a, b):
    unit_a_b = unit(b - a)
    dist_from_a = dot(unit_a_b, p - a)
    return a + unit_a_b * dist_from_a


def intersection_of_circle_and_line(center, r, a, b):
    near = near_pos_online(center, a, b)
    d = sqrt(r ** 2 - abs(center - near) ** 2) * unit(b - a)
    return [near + d, near - d]


from math import sqrt
from operator import attrgetter
def main():
    cx, cy, r = map(int, readline().split())
    center = complex(cx, cy)
    q = int(readline())
    for i in range(q):
        x1, y1, x2, y2 = map(int, readline().split())
        a, b = complex(x1, y1), complex(x2, y2)
        intersection = intersection_of_circle_and_line(center, r, a, b)
        intersection.sort(key=attrgetter('real', 'imag'))
        pre, post = intersection

        print(*list(map('{:.6f}'.format, [pre.real, pre.imag, post.real, post.imag])))

main()

Solution (Python):



"""
input:
2 1 1
2
0 1 4 1
3 0 3 3

output:
1.00000000 1.00000000 3.00000000 1.00000000
3.00000000 1.00000000 3.00000000 1.00000000
"""

import sys
import math


class Segment(object):
    __slots__ = ('source', 'target')

    def __init__(self, source, target):
        self.source = complex(source)
        self.target = complex(target)


class Circle(object):
    __slots__ = ('centre', 'radius')

    def __init__(self, x, y, r):
        self.centre = x + y * 1j
        self.radius = r


def dot(a, b):
    return a.real * b.real + a.imag * b.imag


def project(s, p):
    base_vector = s.target - s.source
    prj_ratio = dot(p - s.source, base_vector) / pow(abs(base_vector), 2)
    return s.source + base_vector * prj_ratio


def get_cross_point(c, l):
    prj_vector = project(l, c.centre)
    line_unit_vector = (l.target - l.source) / (abs(l.target - l.source))
    base = math.sqrt(pow(c.radius, 2) - pow(abs(prj_vector - c.centre), 2))
    p1, p2 = prj_vector + line_unit_vector * base, prj_vector - line_unit_vector * base

    if p1.real < p2.real:
        ans = (p1, p2)
    elif p1.real == p2.real:
        if p1.imag < p2.imag:
            ans = (p1, p2)
        else:
            ans = (p2, p1)
    else:
        ans = (p2, p1)
    return ans


def solve(_lines):
    for line in _lines:
        line_axis = tuple(map(int, line))
        p0, p1 = (x + y * 1j for x, y in zip(line_axis[::2], line_axis[1::2]))
        l = Segment(p0, p1)
        cp1, cp2 = get_cross_point(circle, l)
        print('{0:.8f} {1:.8f} {2:.8f} {3:.8f}'.format(cp1.real, cp1.imag, cp2.real, cp2.imag))

    return None


if __name__ == '__main__':
    _input = sys.stdin.readlines()

    cx, cy, radius = map(int, _input[0].split())
    q_num = int(_input[1])
    lines = map(lambda x: x.split(), _input[2:])

    circle = Circle(cx, cy, radius)
    solve(lines)

Solution (Python):



"""
input:
2 1 1
2
0 1 4 1
3 0 3 3

output:
1.00000000 1.00000000 3.00000000 1.00000000
3.00000000 1.00000000 3.00000000 1.00000000
"""

import sys
import math


class Segment(object):
    __slots__ = ('source', 'target')

    def __init__(self, source, target):
        self.source = complex(source)
        self.target = complex(target)


class Circle(object):
    __slots__ = ('centre', 'radius')

    def __init__(self, centre, radius):
        self.centre = complex(centre)
        self.radius = float(radius)


def dot(a, b):
    return a.real * b.real + a.imag * b.imag


def project(s, p):
    base_vector = s.target - s.source
    prj_ratio = dot(p - s.source, base_vector) / pow(abs(base_vector), 2)
    return s.source + base_vector * prj_ratio


def get_cross_point(c, l):
    prj_vector = project(l, c.centre)
    line_unit_vector = (l.target - l.source) / (abs(l.target - l.source))
    base = math.sqrt(pow(c.radius, 2) - pow(abs(prj_vector - c.centre), 2))
    p1, p2 = prj_vector + line_unit_vector * base, prj_vector - line_unit_vector * base

    if p1.real < p2.real:
        ans = (p1, p2)
    elif p1.real == p2.real:
        if p1.imag < p2.imag:
            ans = (p1, p2)
        else:
            ans = (p2, p1)
    else:
        ans = (p2, p1)
    return ans


def solve(_lines):
    for line in _lines:
        line_axis = tuple(map(int, line))
        p0, p1 = (x + y * 1j for x, y in zip(line_axis[::2], line_axis[1::2]))
        l = Segment(p0, p1)
        cp1, cp2 = get_cross_point(circle, l)
        print('{0:.8f} {1:.8f} {2:.8f} {3:.8f}'.format(cp1.real, cp1.imag, cp2.real, cp2.imag))

    return None


if __name__ == '__main__':
    _input = sys.stdin.readlines()

    cx, cy, r = map(int, _input[0].split())
    q_num = int(_input[1])
    lines = map(lambda x: x.split(), _input[2:])

    circle = Circle(cx + cy * 1j, r)
    solve(lines)

Solution (Python):



"""
input:
2 1 1
2
0 1 4 1
3 0 3 3

output:
1.00000000 1.00000000 3.00000000 1.00000000
3.00000000 1.00000000 3.00000000 1.00000000
"""

import math
import sys
from operator import attrgetter


class Segment(object):
    __slots__ = ('source', 'target')

    def __init__(self, source, target):
        self.source = complex(source)
        self.target = complex(target)


class Circle(object):
    __slots__ = ('centre', 'radius')

    def __init__(self, centre, radius):
        self.centre = complex(centre)
        self.radius = float(radius)


def dot(a, b):
    return a.real * b.real + a.imag * b.imag


def project(s, p):
    base_vector = s.target - s.source
    prj_ratio = dot(p - s.source, base_vector) / pow(abs(base_vector), 2)
    return s.source + base_vector * prj_ratio


def get_cross_point(c, l):
    prj_vector = project(l, c.centre)
    line_unit_vector = (l.target - l.source) / (abs(l.target - l.source))
    base = math.sqrt(pow(c.radius, 2) - pow(abs(prj_vector - c.centre), 2))
    ans = [prj_vector + line_unit_vector * base, prj_vector - line_unit_vector * base]

    ans.sort(key=attrgetter('real', 'imag'))
    return ans


def solve(_lines):
    for line in _lines:
        line_axis = tuple(map(int, line))
        p0, p1 = (x + y * 1j for x, y in zip(line_axis[::2], line_axis[1::2]))
        l = Segment(p0, p1)
        cp1, cp2 = get_cross_point(circle, l)
        print('{0:.8f} {1:.8f} {2:.8f} {3:.8f}'.format(cp1.real, cp1.imag, cp2.real, cp2.imag))

    return None


if __name__ == '__main__':
    _input = sys.stdin.readlines()

    cx, cy, r = map(int, _input[0].split())
    q_num = int(_input[1])
    lines = map(lambda x: x.split(), _input[2:])

    circle = Circle(cx + cy * 1j, r)
    solve(lines)

Solution (Python):



"""
input:
2 1 1
2
0 1 4 1
3 0 3 3

output:
1.00000000 1.00000000 3.00000000 1.00000000
3.00000000 1.00000000 3.00000000 1.00000000
"""

import math
import sys
from operator import attrgetter
from collections import namedtuple


def dot(a, b):
    return a.real * b.real + a.imag * b.imag


def project(s, p):
    base_vector = s.target - s.source
    prj_ratio = dot(p - s.source, base_vector) / pow(abs(base_vector), 2)
    return s.source + base_vector * prj_ratio


def get_cross_point(c, l):
    prj_vector = project(l, c.centre)
    line_unit_vector = (l.target - l.source) / (abs(l.target - l.source))
    base = math.sqrt(pow(c.radius, 2) - pow(abs(prj_vector - c.centre), 2))
    ans = [prj_vector + line_unit_vector * base, prj_vector - line_unit_vector * base]

    ans.sort(key=attrgetter('real', 'imag'))
    return ans


def solve(_lines):
    for line in _lines:
        line_axis = tuple(map(int, line))
        p0, p1 = (x + y * 1j for x, y in zip(line_axis[::2], line_axis[1::2]))
        l = Segment(source=p0, target=p1)
        cp1, cp2 = get_cross_point(circle, l)
        print('{0:.8f} {1:.8f} {2:.8f} {3:.8f}'.format(cp1.real, cp1.imag, cp2.real, cp2.imag))

    return None


if __name__ == '__main__':
    _input = sys.stdin.readlines()

    cx, cy, r = map(int, _input[0].split())
    q_num = int(_input[1])
    lines = map(lambda x: x.split(), _input[2:])

    Circle = namedtuple('Circle', ('centre', 'radius'))
    Segment = namedtuple('Segment', ('source', 'target'))

    circle = Circle(centre=cx + cy * 1j, radius=r)
    solve(lines)

Solution (Python):
def solve_eq(a,b,c):
    return (-b-(b**2-4*a*c)**0.5)/(2*a),(-b+(b**2-4*a*c)**0.5)/(2*a)

def ret_line(x1,y1,x2,y2):
    k = (y1-y2)/(x1-x2)
    d = y1-k*x1
    return k,d
    
x,y,r = map(float,input().split())
q = int(input())
x1,y1,x2,y2 = [],[],[],[]
for i in range(q):
    arr = list(map(float,input().split()))
    x1.append(arr[0])
    y1.append(arr[1])
    x2.append(arr[2])
    y2.append(arr[3])

for i in range(q):
    if x1[i] != x2[i]:
        k,d = ret_line(x1[i],y1[i],x2[i],y2[i])
        p1,p2 = solve_eq(1+k**2,2*k*d-2*k*y-2*x,x**2+y**2+d**2-2*y*d-r**2)
        print(p1,k*p1+d,p2,k*p2+d)
    if x1[i] == x2[i]:
        p1,p2 = solve_eq(1,-2*y,y**2+(x1[i]-x)**2-r**2)
        print(x1[i],p1,x1[i],p2)

Solution (Python):
def dot(c1, c2):
    return c1.real * c2.real + c1.imag * c2.imag

def cross(c1, c2):
    return c1.real * c2.imag - c1.imag * c2.real

def project(p1, p2, p3):
    
    
    base = p2 - p1
    return p1 + base * dot(p3 - p1, base) / (abs(base) ** 2)

def print_cross_point(c, r, p1, p2):
    
    
    pr = project(p1, p2, c)
    l = p2 - p1
    e = l / abs(l)
    base = (r ** 2 - abs(pr - c) ** 2) ** 0.5
    c_p1 = pr - e * base
    c_p2 = pr + e * base
    print("{0:.10f} {1:.10f} {2:.10f} {3:.10f}".format(c_p1.real, c_p1.imag,
                                                       c_p2.real, c_p2.imag))


import sys

file_input = sys.stdin

cx, cy, r = map(int, file_input.readline().split())
c = cx + cy * 1j

q = file_input.readline()

for line in file_input:
    x1, y1, x2, y2 = map(int, line.split())
    p1 = x1 + y1 * 1j
    p2 = x2 + y2 * 1j
    if x1 < x2:
        print_cross_point(c, r, p1, p2)
    elif x1 > x2 or y1 > y2:
        print_cross_point(c, r, p2, p1)
    else:
        print_cross_point(c, r, p1, p2)

Solution (Python):
cx, cy, r = map(int, input().split())
q = int(input())
for _ in range(q):
    x1, y1, x2, y2 = map(int, input().split())
    dx, dy = x2 - x1, y2 - y1
    px, py = cx - x1, cy - y1
    s2 = dx ** 2 + dy ** 2
    dot = dx * px + dy * py
    crs = dx * py - dy * px
    t1 = (dot + (s2 * r ** 2 - crs ** 2) ** 0.5) / s2
    t2 = (dot - (s2 * r ** 2 - crs ** 2) ** 0.5) / s2
    if (dx < 0 or (dx == 0 and dy < 0)) ^ (t1 > t2):
        t1, t2 = t2, t1
    print(x1 + dx * t1, y1 + dy * t1, x1 + dx * t2, y1 + dy * t2)

Solution (Python):
from math import sqrt
def solve(cx, cy, r, x1, y1, x2, y2):
    xd = x2 - x1; yd = y2 - y1
    X = x1 - cx; Y = y1 - cy
    a = xd**2 + yd**2
    b = xd * X + yd * Y
    c = X**2 + Y**2 - r**2
    D = b**2 - a*c
    s1 = (-b + sqrt(D)) / a
    s2 = (-b - sqrt(D)) / a
    return (x1 + xd*s1, y1 + yd*s1), (x1 + xd*s2, y1 + yd*s2)
cx, cy, r = map(int, raw_input().split())
q = input()
for i in xrange(q):
    x1, y1, x2, y2 = map(int, raw_input().split())

    p0, p1 = sorted(solve(cx, cy, r, x1, y1, x2, y2))
    print "%.08f %.08f %.08f %.08f" % (p0 + p1)

Solution (Python):
from math import sqrt
def intersection_cl(S0, S1, C0, r):
    x0, y0 = S0; x1, y1 = S1
    cx, cy = C0

    xd = x1 - x0; yd = y1 - y0
    X = x0 - cx; Y = y0 - cy
    A = xd**2 + yd**2
    B = xd*X + yd*Y
    C = X**2 + Y**2 - r**2
    if B**2 < A*C:
        
        return tuple()
    """
    if B**2 == A*C:
        
        s = -B / float(A)
        return tuple((x0 + xd*s, y0 + yd*s))
    """
    
    D = B**2 - A*C
    s1 = (-B + sqrt(D)) / A
    s2 = (-B - sqrt(D)) / A
    return ((x0 + xd*s1, y0 + yd*s1), (x0 + xd*s2, y0 + yd*s2))

cx, cy, r = map(int, input().split())
q = int(input())
for i in range(q):
    x1, y1, x2, y2 = map(int, input().split())

    p0, p1 = sorted(intersection_cl((x1, y1), (x2, y2), (cx, cy), r))
    print("%.08f %.08f %.08f %.08f" % (p0 + p1))

Solution (Python):
import sys
from itertools import starmap
import math
readline = sys.stdin.readline
EPS = 1e-9
ONLINE_FRONT = -2
CLOCKWISE = -1
ON_SEGMENT = 0
COUNTER_CLOCKWISE = 1
ONLINE_BACK = 2
class Circle(object):
    __slots__ = ('c', 'r')
    def __init__(self, c, r):
        self.c = c
        self.r = r
class Segment(object):
    __slots__ = ('fi', 'se')
    def __init__(self, fi, se):
        self.fi = fi
        self.se = se
Line = Segment
def cross(a, b):
    return a.real * b.imag - a.imag * b.real
def dot(a, b):
    return a.real * b.real + a.imag * b.imag
def norm(base):
    return abs(base) ** 2
def project(s, p2):
    base = s.fi - s.se
    r = dot(p2 - s.fi, base) / norm(base)
    return s.fi + base * r
def reflect(s, p):
    return p + (project(s, p) - p) * 2.0
def ccw(p1, p2, p3):
    a = p2 - p1
    b = p3 - p1
    if cross(a, b) > EPS: return 1
    if cross(a, b) < -EPS: return -1
    if dot(a, b) < -EPS: return 2
    if norm(a) < norm(b): return -2
    return 0
def intersect4(p1, p2, p3, p4):
    return (ccw(p1, p2, p3) * ccw(p1, p2, p4) <= 0 and
			ccw(p3, p4, p1) * ccw(p3, p4, p2) <= 0)
def intersect2(s1, s2):
    return intersect4(s1.fi, s1.se, s2.fi, s2.se)
def getDistance(a, b):  return abs(a - b)
def getDistanceLP(l, p):
    return abs(cross(l.se - l.fi, p - l.fi) / abs(l.se - l.fi))
def getDistanceSP(s, p):
    if dot(s.se - s.fi, p - s.fi) < 0.0: return abs(p - s.fi)
    if dot(s.fi - s.se, p - s.se) < 0.0: return abs(p - s.se)
    return getDistanceLP(s, p)
def getDistances(s1, s2):
    if intersect2(s1, s2): return 0.0
    return min(getDistanceSP(s1, s2.fi), getDistanceSP(s1, s2.se),
               getDistanceSP(s2, s1.fi), getDistanceSP(s2, s1.se))
def getCrossPoint(s1, s2):
	base = s2.se - s2.fi
	d1 = abs(cross(base, s1.fi - s2.fi))
	d2 = abs(cross(base, s1.se - s2.fi))
	t = d1 / (d1 + d2)
	return s1.fi + (s1.se - s1.fi) * t
def getCrossPoints(c, l):
	pr = project(l, c.c)
	e = (l.se - l.fi) / abs(l.se - l.fi)
	base = math.sqrt(c.r * c.r - norm(pr - c.c))
	return Segment(pr + e * base, pr - e * base)
x, y, r = map(int, readline().split())
ci = Circle(complex(x, y), r)
n = int(readline())
for _ in [0] * n:
    p0, p1 = starmap(complex, zip(*[map(int, readline().split())] * 2))
    ppp = getCrossPoints(ci, Line(p0, p1))
    if ppp.se.real < ppp.fi.real or (ppp.se.real == ppp.fi.real and ppp.se.imag < ppp.fi.imag):
        ppp.fi, ppp.se = ppp.se, ppp.fi
    print("{:.10f} {:.10f} {:.10f} {:.10f}".format(ppp.fi.real, ppp.fi.imag, ppp.se.real, ppp.se.imag))

Solution (Python):
import sys
from itertools import starmap
import math
readline = sys.stdin.readline
EPS = 1e-9
ONLINE_FRONT = -2
CLOCKWISE = -1
ON_SEGMENT = 0
COUNTER_CLOCKWISE = 1
ONLINE_BACK = 2
class Circle(object):
    __slots__ = ('c', 'r')
    def __init__(self, c, r):
        self.c = c
        self.r = r
class Segment(object):
    __slots__ = ('fi', 'se')
    def __init__(self, fi, se):
        self.fi = fi
        self.se = se
Line = Segment
def cross(a, b):
    return a.real * b.imag - a.imag * b.real
def dot(a, b):
    return a.real * b.real + a.imag * b.imag
def norm(base):
    return abs(base) ** 2
def project(s, p2):
    base = s.fi - s.se
    r = dot(p2 - s.fi, base) / norm(base)
    return s.fi + base * r
def reflect(s, p):
    return p + (project(s, p) - p) * 2.0
def ccw(p1, p2, p3):
    a = p2 - p1
    b = p3 - p1
    if cross(a, b) > EPS: return 1
    if cross(a, b) < -EPS: return -1
    if dot(a, b) < -EPS: return 2
    if norm(a) < norm(b): return -2
    return 0
def intersect4(p1, p2, p3, p4):
    return (ccw(p1, p2, p3) * ccw(p1, p2, p4) <= 0 and
			ccw(p3, p4, p1) * ccw(p3, p4, p2) <= 0)
def intersect2(s1, s2):
    return intersect4(s1.fi, s1.se, s2.fi, s2.se)
def getDistance(a, b):  return abs(a - b)
def getDistanceLP(l, p):
    return abs(cross(l.se - l.fi, p - l.fi) / abs(l.se - l.fi))
def getDistanceSP(s, p):
    if dot(s.se - s.fi, p - s.fi) < 0.0: return abs(p - s.fi)
    if dot(s.fi - s.se, p - s.se) < 0.0: return abs(p - s.se)
    return getDistanceLP(s, p)
def getDistances(s1, s2):
    if intersect2(s1, s2): return 0.0
    return min(getDistanceSP(s1, s2.fi), getDistanceSP(s1, s2.se),
               getDistanceSP(s2, s1.fi), getDistanceSP(s2, s1.se))
def getCrossPoint(s1, s2):
	base = s2.se - s2.fi
	d1 = abs(cross(base, s1.fi - s2.fi))
	d2 = abs(cross(base, s1.se - s2.fi))
	t = d1 / (d1 + d2)
	return s1.fi + (s1.se - s1.fi) * t
def getCrossPoints(c, l):
	pr = project(l, c.c)
	e = (l.se - l.fi) / abs(l.se - l.fi)
	base = math.sqrt(c.r * c.r - norm(pr - c.c))
	return Segment(pr + e * base, pr - e * base)
x, y, r = map(int, readline().split())
ci = Circle(complex(x, y), r)
n = int(readline())
for _ in [0] * n:
    p0, p1 = starmap(complex, zip(*[map(int, readline().split())] * 2))
    ppp = getCrossPoints(ci, Line(p0, p1))
    if ppp.se.real < ppp.fi.real or (ppp.se.real == ppp.fi.real and ppp.se.imag < ppp.fi.imag):
        ppp.fi, ppp.se = ppp.se, ppp.fi
    print(*map('{:.10f}'.format, [ppp.fi.real, ppp.fi.imag, ppp.se.real, ppp.se.imag]))

Solution (Python):
import sys
from operator import itemgetter, attrgetter
from itertools import starmap
import math
readline = sys.stdin.readline
EPS = 1e-9
ONLINE_FRONT = -2
CLOCKWISE = -1
ON_SEGMENT = 0
COUNTER_CLOCKWISE = 1
ONLINE_BACK = 2
class Circle(object):
    __slots__ = ('c', 'r')
    def __init__(self, c, r):
        self.c = c
        self.r = r
class Segment(object):
    __slots__ = ('fi', 'se')
    def __init__(self, fi, se):
        self.fi = fi
        self.se = se
Line = Segment
def cross(a, b):
    return a.real * b.imag - a.imag * b.real
def dot(a, b):
    return a.real * b.real + a.imag * b.imag
def norm(base):
    return abs(base) ** 2
def project(s, p2):
    base = s.fi - s.se
    r = dot(p2 - s.fi, base) / norm(base)
    return s.fi + base * r
def reflect(s, p):
    return p + (project(s, p) - p) * 2.0
def ccw(p1, p2, p3):
    a = p2 - p1
    b = p3 - p1
    if cross(a, b) > EPS: return 1
    if cross(a, b) < -EPS: return -1
    if dot(a, b) < -EPS: return 2
    if norm(a) < norm(b): return -2
    return 0
def intersect4(p1, p2, p3, p4):
    return (ccw(p1, p2, p3) * ccw(p1, p2, p4) <= 0 and
			ccw(p3, p4, p1) * ccw(p3, p4, p2) <= 0)
def intersect2(s1, s2):
    return intersect4(s1.fi, s1.se, s2.fi, s2.se)
def getDistance(a, b):  return abs(a - b)
def getDistanceLP(l, p):
    return abs(cross(l.se - l.fi, p - l.fi) / abs(l.se - l.fi))
def getDistanceSP(s, p):
    if dot(s.se - s.fi, p - s.fi) < 0.0: return abs(p - s.fi)
    if dot(s.fi - s.se, p - s.se) < 0.0: return abs(p - s.se)
    return getDistanceLP(s, p)
def getDistances(s1, s2):
    if intersect2(s1, s2): return 0.0
    return min(getDistanceSP(s1, s2.fi), getDistanceSP(s1, s2.se),
               getDistanceSP(s2, s1.fi), getDistanceSP(s2, s1.se))
def getCrossPoint(s1, s2):
	base = s2.se - s2.fi
	d1 = abs(cross(base, s1.fi - s2.fi))
	d2 = abs(cross(base, s1.se - s2.fi))
	t = d1 / (d1 + d2)
	return s1.fi + (s1.se - s1.fi) * t
def getCrossPoints(c, l):
    pr = project(l, c.c)
    e = (l.se - l.fi) / abs(l.se - l.fi)
    base = math.sqrt(c.r * c.r - norm(pr - c.c))
    return Segment(*sorted((pr + e * base, pr - e * base), key=attrgetter('real', 'imag')))
x, y, r = map(int, readline().split())
ci = Circle(complex(x, y), r)
n = int(readline())
for _ in [0] * n:
    p0, p1 = starmap(complex, zip(*[map(int, readline().split())] * 2))
    ppp = getCrossPoints(ci, Line(p0, p1))
    print(*map('{:.10f}'.format, (ppp.fi.real, ppp.fi.imag, ppp.se.real, ppp.se.imag)))

Solution (Python):

"""
http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_7_D

"""
import sys
from math import sqrt
from sys import stdin
input = stdin.readline


class Point(object):
    epsilon = 1e-10

    def __init__(self, x=0.0, y=0.0):
        if isinstance(x, tuple):
            self.x = x[0]
            self.y = x[1]
        else:
            self.x = x
            self.y = 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 __mul__(self, other):
        return Point(other * self.x, other * self.y)

    def __truediv__(self, other):
        return Point(self.x / other, self.y / other)

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

    def __eq__(self, other):
        from math import fabs
        if fabs(self.x - other.x) < Point.epsilon and fabs(self.y - other.y) < Point.epsilon:
            return True
        else:
            return False

    def norm(self):
        return self.x * self.x + self.y * self.y

    def __abs__(self):
        return sqrt(self.norm())

    def ccw(self, p0, p1):
        
        a = Vector(p1 - p0)
        b = Vector(self - p0)
        if Vector.cross(a, b) > Point.epsilon:
            return 1 
        elif Vector.cross(a, b) < -Point.epsilon:
            return -1 
        elif Vector.dot(a, b) < -Point.epsilon:
            return 2 
        elif a.norm() < b.norm():
            return -2 
        else:
            return 0 

    def project(self, s):
        
        base = Vector(s.p2 - s.p1)
        a = Vector(self - s.p1)
        r = Vector.dot(a, base)
        r /= base.norm()
        return s.p1 + base * r

    def reflect(self, s):
        
        proj = self.project(s)
        return self + (proj - self)*2

    def distance(self, s):
        
        if Vector.dot(s.p2-s.p1, self-s.p1) < 0.0:
            return abs(self - s.p1)
        if Vector.dot(s.p1-s.p2, self-s.p2) < 0.0:
            return abs(self - s.p2)
        return abs(Vector.cross(s.p2-s.p1, self-s.p1) / abs(s.p2-s.p1))


class Vector(Point):
    def __init__(self, x=0.0, y=0.0):
        if isinstance(x, tuple):
            self.x = x[0]
            self.y = x[1]
        elif isinstance(x, Point):
            self.x = x.x
            self.y = x.y
        else:
            self.x = x
            self.y = y

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

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

    def __mul__(self, other):
        return Vector(other * self.x, other * self.y)

    def __truediv__(self, other):
        return Vector(other / self.x, other / self.y)

    @classmethod
    def dot(cls, a, b):
        return a.x * b.x + a.y * b.y

    @classmethod
    def cross(cls, a, b):
        return a.x * b.y - a.y * b.x

    @classmethod
    def is_orthogonal(cls, a, b):
        return Vector.dot(a, b) == 0.0

    @classmethod
    def is_parallel(cls, a, b):
        return Vector.cross(a, b) == 0.0


class Segment(object):
    def __init__(self, p1=Point(), p2=Point()):
        if isinstance(p1, Point):
            self.p1 = p1
            self.p2 = p2
        elif isinstance(p1, tuple):
            self.p1 = Point(p1[0], p1[1])
            self.p2 = Point(p2[0], p2[1])

    def intersect(self, s):
        
        ans1 = s.p1.ccw(self.p1, self.p2) * s.p2.ccw(self.p1, self.p2)
        ans2 = self.p1.ccw(s.p1, s.p2) * self.p2.ccw(s.p1, s.p2)
        return ans1 <= 0 and ans2 <= 0

    def cross_point(self, s):
        
        base = s.p2 - s.p1
        d1 = abs(Vector.cross(base, self.p1-s.p1))
        d2 = abs(Vector.cross(base, self.p2-s.p1))
        t = d1 / (d1 + d2)
        return self.p1 + (self.p2 - self.p1) * t

    def distance(self, s):
        
        if self.intersect(s):
            return 0.0
        d1 = s.p1.distance(self)
        d2 = s.p2.distance(self)
        d3 = self.p1.distance(s)
        d4 = self.p2.distance(s)
        return min(d1, d2, d3, d4)


    @classmethod
    def is_orthogonal(cls, s1, s2):
        a = Vector(s1.p2 - s1.p1)
        b = Vector(s2.p2 - s2.p1)
        return Vector.is_orthogonal(a, b)

    @classmethod
    def is_parallel(cls, s1, s2):
        a = Vector(s1.p2 - s1.p1)
        b = Vector(s2.p2 - s2.p1)
        return Vector.is_parallel(a, b)


class Line(Segment):
    pass


class Cirle(object):
    def __init__(self, x, y=Point(), r=1.0):
        if isinstance(x, Point):
            self.c = x
            self.r = y
        elif isinstance(x, tuple):
            self.c = Point(x[0], x[1])
            self.r = r

    def cross_point(self, s):
        pr = self.c.project(s)
        e = (s.p2 - s.p1) / abs(s.p2 - s.p1)
        base = sqrt(self.r * self.r - (pr - self.c).norm())
        return pr + e * base, pr - e * base


def main(args):
    x, y, r = map(int, input().split())
    q = int(input())
    c = Cirle(Point(x, y), r)
    for _ in range(q):
        x1, y1, x2, y2 = map(int, input().split())
        p1 = Point(x1, y1)
        p2 = Point(x2, y2)
        s = Segment(p1, p2)
        ans1, ans2 = c.cross_point(s)
        if ans1.x < ans2.x:
            print('{:.8f} {:.8f} {:.8f} {:.8f}'.format(ans1.x, ans1.y, ans2.x, ans2.y))
        elif ans1.x == ans2.x and ans1.y <= ans2.y:
            print('{:.8f} {:.8f} {:.8f} {:.8f}'.format(ans1.x, ans1.y, ans2.x, ans2.y))
        else:
            print('{:.8f} {:.8f} {:.8f} {:.8f}'.format(ans2.x, ans2.y, ans1.x, ans1.y))


if __name__ == '__main__':
    main(sys.argv[1:])

Solution (Python):
from math import sqrt
cx,cy,r = map(int,input().split())
c = [cx, cy]
q = int(input())

class Line:pass

def dot(a, b):
    return sum([i * j for i,j in zip(a, b)])
def add(a,b):
    return [a[0] + b[0],a[1] + b[1]]
def sub(a, b):
    return [a[0] - b[0],a[1] - b[1]]
def multi(k, a):
    return [k*a[0],k*a[1]]
def cross(a, b):
    return  a[0] * b[1] - a[1] * b[0]
def norm(a):
    return sqrt(a[0] ** 2 + a[1] ** 2)    
def proj(p, l):
    base = sub(l.p2, l.p1)
    hypo = sub(p, l.p1)
    t = dot(hypo,base) / norm(base)
    return add(l.p1, multi(t/norm(base), base))
    
for i in range(q):
    x1,y1,x2,y2 = map(int, input().split())
    l = Line()
    l.p1 = [x1,y1]
    l.p2 = [x2,y2]
    pr = proj(c, l)
    e = multi(1/norm(sub(l.p2, l.p1)),sub(l.p2, l.p1))
    e2 = multi(-1, e)
    base = sqrt(r**2 - norm(sub(pr,c))**2)
    a = add(pr, multi(base, e))
    b = add(pr, multi(base, e2))
    ans = [a,b]
    ans = sorted(ans, key=lambda x: (x[0],x[1]))
    print(' '.join([str(round(i,8)) for o in ans for i in o]))