Problem:
Cross Points of Circles For given two circles $c1$ and $c2$, print the coordinates of the cross points of them. Input The input is given in the following format. $c1x\; c1y\; c1r$ $c2x\; c2y\; c2r$  $c1x$, $c1y$ and $c1r$ represent the coordinate and radius of the first circle.
  $c2x$, $c2y$ and $c2r$ represent the coordinate and radius of the second circle.
  All input values are given in integers. Output Print the coordinates ($x1$, $y1$) and ($x2$, $y2$) of the cross points $p1$ and $p2$ respectively 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  The given circle have at least one cross point and have different center coordinates.  $-10,000 \leq c1x, c1y, c2x, c2y \leq 10,000$  $1 \leq c1r, c2r \leq 10,000$  Sample Input and Output Sample Input 1 0 0 2
2 0 2 Sample Output 1 1.00000000 -1.73205080 1.00000000 1.73205080 Sample Input 2 0 0 2
0 3 1 Sample Output 2 0.00000000 2.00000000 0.00000000 2.00000000

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

c1x, c1y, c1r = map(int, input().split())
c2x, c2y, c2r = map(int, input().split())
a = complex(c1x,c1y)
b = complex(c2x,c2y)
ab=abs(b-a)
ac=(c1r**2 + ab**2-c2r**2)/(2*ab)
c=a+(b-a)*(ac/ab)
pc=(c1r**2-ac**2)**0.5
p=c+pc*unit(b-a)*1j
q=c+pc*unit(b-a)*-1j

if p.real>q.real:
  p,q=q,p
if abs(p.real - q.real) < 1e-6 and p.imag > q.imag:
  p, q = q, p
print('{:.10f} {:.10f} {:.10f} {:.10f}'.format(p.real, p.imag, q.real, q.imag))

Solution (Python):
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())

    def rotate(self, theta):
        return Vector(self.x * math.cos(theta)
                               - self.y * math.sin(theta),
                               self.x * math.sin(theta)
                               + self.y * math.cos(theta))

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 get_cross_point(c1, c2):
    v = c2.c - c1.c
    d = v.abs()
    assert d <= (c1.r + c2.r)
    a = math.acos((c1.r * c1.r + d * d - c2.r * c2.r) / (2 * c1.r  * d))
    return c1.c + v.rotate(a) * (c1.r/d), c1.c + v.rotate(-a) * (c1.r/d)

x1, y1, r1 = [ int(i) for i in input().split() ]
c1 = Circle(Point(x1, y1), r1)
x2, y2, r2 = [ int(i) for i in input().split() ]
c2= Circle(Point(x2, y2), r2)
p1, p2 = get_cross_point(c1, c2)
x1, y1 = min((p1.x, p1.y), (p2.x, p2.y))
x2, y2 = max((p1.x, p1.y), (p2.x, p2.y))
print("{0:0.8f} {1:0.8f} {2:0.8f} {3:0.8f}".format(x1, y1, x2, y2))

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


def cosine_formula(a, b, c):
    return (b ** 2 + c ** 2 - a ** 2) / (2 * b * c)


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

import math

def main():
    ax, ay, ar = map(int, readline().split())
    bx, by, br = map(int, readline().split())
    a, b = complex(ax, ay), complex(bx, by)

    ab = abs(b - a)
    ac = (ar ** 2 + ab ** 2 - br ** 2) / (2 * ab)
    c = a + (b - a) * (ac / ab) 
    pc = math.sqrt(ar ** 2 - ac ** 2)
    p, q = c + pc * unit(b - a) * 1j, c + pc * unit(b - a) * -1j
    
    if p.real > q.real:
        p, q = q, p
    if abs(p.real - q.real) < 1e-6 and p.imag > q.imag:
        p, q = q, p

    print('{:.10f} {:.10f} {:.10f} {:.10f}'.format(p.real, p.imag, q.real, q.imag))
main()

Solution (Python):



"""
input:
0 0 2
2 0 2

output:
1.00000000 -1.73205080 1.00000000 1.73205080
"""

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


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

    def __init__(self, centre=float('inf'), radius=-1):
        self.centre = complex(centre)
        self.radius = float(radius)


def arg(p_vector):
    return atan2(p_vector.imag, p_vector.real)


def polar(a, r):
    return cos(r) * a + sin(r) * a * 1j


def get_circles_cross_points(c1, c2):
    centre_distance = abs(c1.centre - c2.centre)
    a = acos((pow(c1.radius, 2) + pow(centre_distance, 2) - pow(c2.radius, 2)) /
             (2 * c1.radius * centre_distance))
    t = arg(c2.centre - c1.centre)
    p1, p2 = c1.centre + polar(c1.radius, t + a), c1.centre + polar(c1.radius, t - a)

    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(_circle_info):
    for i, circle in enumerate(_circle_info):
        cx, cy, r = map(int, circle)
        circles[i].centre = cx + cy * 1j
        circles[i].radius = r

    cp1, cp2 = get_circles_cross_points(*circles)
    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()
    circle_info = map(lambda x: x.split(), _input)

    circles = [Circle() for _ in range(2)]
    solve(circle_info)

Solution (Python):



"""
input:
0 0 2
2 0 2

output:
1.00000000 -1.73205080 1.00000000 1.73205080
"""

import sys
from math import sin, cos, atan2, acos


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

    def __init__(self, centre=float('inf'), radius=-1):
        self.centre = complex(centre)
        self.radius = float(radius)


def arg(p_vector):
    return atan2(p_vector.imag, p_vector.real)


def polar(a, r):
    return cos(r) * a + sin(r) * a * 1j


def get_circles_cross_points(c1, c2):
    centre_distance = abs(c1.centre - c2.centre)
    a = acos((pow(c1.radius, 2) + pow(centre_distance, 2) - pow(c2.radius, 2)) /
             (2 * c1.radius * centre_distance))
    t = arg(c2.centre - c1.centre)
    p1, p2 = c1.centre + polar(c1.radius, t + a), c1.centre + polar(c1.radius, t - a)

    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(_circle_info):
    for i, circle in enumerate(_circle_info):
        cx, cy, r = map(int, circle)
        circles[i].centre = cx + cy * 1j
        circles[i].radius = r

    cp1, cp2 = get_circles_cross_points(*circles)
    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()
    circle_info = map(lambda x: x.split(), _input)

    circles = [Circle() for _ in range(2)]
    solve(circle_info)

Solution (Python):



"""
input:
0 0 2
2 0 2

output:
1.00000000 -1.73205080 1.00000000 1.73205080
"""

import sys
from math import sin, cos, atan2, acos


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

    def __init__(self, centre=float('inf'), radius=-1):
        self.centre = complex(centre)
        self.radius = float(radius)


def arg(p_vector):
    return atan2(p_vector.imag, p_vector.real)


def polar(a, r):
    return cos(r) * a + sin(r) * a * 1j


def get_circles_cross_points(c1, c2):
    centre_distance = abs(c1.centre - c2.centre)
    a = acos((pow(c1.radius, 2) + pow(centre_distance, 2) - pow(c2.radius, 2)) /
             (2 * c1.radius * centre_distance))
    t = arg(c2.centre - c1.centre)
    p1, p2 = c1.centre + polar(c1.radius, t + a), c1.centre + polar(c1.radius, t - a)

    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(_circle_info):
    for i, circle in enumerate(_circle_info):
        cx, cy, r = map(int, circle)
        circles[i].centre = cx + cy * 1j
        circles[i].radius = r

    cp1, cp2 = get_circles_cross_points(*circles)
    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()
    circle_info = map(lambda x: x.split(), _input)

    circles = [Circle() for _ in range(2)]
    solve(circle_info)

Solution (Python):



"""
input:
0 0 2
2 0 2

output:
1.00000000 -1.73205080 1.00000000 1.73205080
"""

import sys
from math import sin, cos, atan2, acos


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

    def __init__(self, centre=float('inf'), radius=-1):
        self.centre = complex(centre)
        self.radius = float(radius)


def arg(p_vector):
    return atan2(p_vector.imag, p_vector.real)


def polar(a, r):
    return a * cos(r) + a * sin(r) * 1j


def get_circles_cross_points(c1, c2):
    centre_distance = abs(c1.centre - c2.centre)
    a = acos((pow(c1.radius, 2) + pow(centre_distance, 2) - pow(c2.radius, 2)) /
             (2 * c1.radius * centre_distance))
    t = arg(c2.centre - c1.centre)
    p1, p2 = c1.centre + polar(c1.radius, t + a), c1.centre + polar(c1.radius, t - a)

    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(_circle_info):
    for i, circle in enumerate(_circle_info):
        cx, cy, r = map(int, circle)
        circles[i].centre = cx + cy * 1j
        circles[i].radius = r

    cp1, cp2 = get_circles_cross_points(*circles)
    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()
    circle_info = map(lambda x: x.split(), _input)

    circles = [Circle() for _ in range(2)]
    solve(circle_info)

Solution (Python):



"""
input:
0 0 2
2 0 2

output:
1.00000000 -1.73205080 1.00000000 1.73205080
"""

import sys
from math import sin, cos, atan2, acos


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

    def __init__(self, centre=float('inf'), radius=-1):
        self.centre = complex(centre)
        self.radius = float(radius)


def arg(p_vector):
    return atan2(p_vector.imag, p_vector.real)


def polar(a, r):
    return a * cos(r) + a * sin(r) * 1j


def get_circles_cross_points(c1, c2):
    centre_distance = abs(c1.centre - c2.centre)
    a = acos((pow(c1.radius, 2) + pow(centre_distance, 2) - pow(c2.radius, 2)) /
             (2 * c1.radius * centre_distance))
    t = arg(c2.centre - c1.centre)
    p1, p2 = c1.centre + polar(c1.radius, t + a), c1.centre + polar(c1.radius, t - a)

    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(_circle_info):
    for i, circle in enumerate(_circle_info):
        cx, cy, r = map(int, circle)
        circles[i].centre = cx + cy * 1j
        circles[i].radius = r

    cp1, cp2 = get_circles_cross_points(*circles)
    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()
    circle_info = map(lambda x: x.split(), _input)

    circles = [Circle() for _ in range(2)]
    solve(circle_info)

Solution (Python):



"""
input:
0 0 2
2 0 2

output:
1.00000000 -1.73205080 1.00000000 1.73205080
"""

import sys
from math import sin, cos, atan2, acos
from operator import attrgetter


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

    def __init__(self, centre=float('inf'), radius=-1):
        self.centre = complex(centre)
        self.radius = float(radius)


def arg(p_vector):
    return atan2(p_vector.imag, p_vector.real)


def polar(a, r):
    return a * cos(r) + a * sin(r) * 1j


def get_circles_cross_points(c1, c2):
    centre_distance = abs(c1.centre - c2.centre)
    a = acos((pow(c1.radius, 2) + pow(centre_distance, 2) - pow(c2.radius, 2)) /
             (2 * c1.radius * centre_distance))
    t = arg(c2.centre - c1.centre)
    ans = [c1.centre + polar(c1.radius, t + a), c1.centre + polar(c1.radius, t - a)]
    ans.sort(key=attrgetter('real', 'imag'))

    return ans


def solve(_circle_info):
    for i, circle in enumerate(_circle_info):
        cx, cy, r = map(int, circle)
        circles[i].centre = cx + cy * 1j
        circles[i].radius = r

    cp1, cp2 = get_circles_cross_points(*circles)
    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()
    circle_info = map(lambda x: x.split(), _input)

    circles = [Circle() for _ in range(2)]
    solve(circle_info)

Solution (Python):



"""
input:
0 0 2
2 0 2

output:
1.00000000 -1.73205080 1.00000000 1.73205080
"""

import sys
from math import sin, cos, atan2, acos
from operator import attrgetter
from collections import namedtuple


def arg(p_vector):
    return atan2(p_vector.imag, p_vector.real)


def polar(a, r):
    return a * cos(r) + a * sin(r) * 1j


def get_circles_cross_points(c1, c2):
    centre_distance = abs(c1.centre - c2.centre)
    
    a = acos((pow(c1.radius, 2) + pow(centre_distance, 2) - pow(c2.radius, 2)) /
             (2 * c1.radius * centre_distance))
    t = arg(c2.centre - c1.centre)
    ans = [c1.centre + polar(c1.radius, t + a), c1.centre + polar(c1.radius, t - a)]
    ans.sort(key=attrgetter('real', 'imag'))

    return ans


def solve(_circles):
    cp1, cp2 = get_circles_cross_points(*_circles)
    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()
    circle_info = map(lambda x: x.split(), _input)

    Circle = namedtuple('Circle', 'centre radius')
    circles = [Circle(centre=int(cx) + int(cy) * 1j, radius=int(r)) for cx, cy, r in circle_info]
    solve(circles)

Solution (Python):
x1,y1,r1 = map(float,input().split())
x2,y2,r2 = map(float,input().split())
a,b = x2-x1,y2-y1
d = (a**2+b**2+r1**2-r2**2)/2
x = (d*a+b*(((a**2+b**2)*r1**2 - d**2)**0.5))/(a**2+b**2)
y = (d*b-a*(((a**2+b**2)*r1**2 - d**2)**0.5))/(a**2+b**2)

z = (d*a-b*(((a**2+b**2)*r1**2 - d**2)**0.5))/(a**2+b**2)
w = (d*b+a*(((a**2+b**2)*r1**2 - d**2)**0.5))/(a**2+b**2)
if x < z:
    print(x1+x,y1+y,end=' ')
    print(x1+z,y1+w)

elif x > z:
    print(x1+z,y1+w,end=' ')
    print(x1+x,y1+y)
elif y < w:
    print(x1+x,y1+y,end=' ')
    print(x1+z,y1+w)
else:
    print(x1+z,y1+w,end=' ')
    print(x1+x,y1+y)

Solution (Python):
import math, cmath

def print_cross_points(c1, r1, c2, r2):
    c1c2 = c2 - c1
    d = abs(c1c2)
    a = math.acos((r1 ** 2 + d ** 2 - r2 ** 2) / (2 * r1 * d))
    t = cmath.phase(c1c2)
    cp1 = c1 + cmath.rect(r1, t + a)
    cp2 = c1 + cmath.rect(r1, t - a)
    print("{0:.8f} {1:.8f} {2:.8f} {3:.8f}".format(cp1.real, cp1.imag,
                                                       cp2.real, cp2.imag))

c1x, c1y, c1r = map(int, input().split())
c2x, c2y, c2r = map(int, input().split())

if c1y < c2y:
    print_cross_points(c1x + c1y * 1j, c1r, c2x + c2y * 1j, c2r)
elif c1y > c2y or c1x < c2x:
    print_cross_points(c2x + c2y * 1j, c2r, c1x + c1y * 1j, c1r)
else:
    print_cross_points(c1x + c1y * 1j, c1r, c2x + c2y * 1j, c2r)

Solution (Python):
c1x, c1y, c1r = map(int, input().split())
c2x, c2y, c2r = map(int, input().split())
dx, dy = c2x - c1x, c2y - c1y
s2 = dx ** 2 + dy ** 2
a = (s2 + c1r ** 2 - c2r ** 2) / 2
df = (s2 * c1r ** 2 - a ** 2) ** 0.5
ix1, ix2 = (a * dx + dy * df) / s2, (a * dx - dy * df) / s2
iy1, iy2 = (a * dy - dx * df) / s2, (a * dy + dx * df) / s2
if (ix1, iy1) > (ix2, iy2):
    ix1, ix2 = ix2, ix1
    iy1, iy2 = iy2, iy1
print(ix1 + c1x, iy1 + c1y, ix2 + c1x, iy2 + c1y)

Solution (Python):
from math import sqrt
x1, y1, r1 = map(int, raw_input().split())
x2, y2, r2 = map(int, raw_input().split())

def solve(x1, y1, r1, x2, y2, r2):
    rr0 = (x2 - x1)**2 + (y2 - y1)**2
    xd = x2 - x1
    yd = y2 - y1
    rr1 = r1**2; rr2 = r2**2
    cv = (rr0 + rr1 - rr2)
    sv = sqrt(4*rr0*rr1 - cv**2)
    return (
        (x1 + (cv*xd - sv*yd)/(2.*rr0), y1 + (cv*yd + sv*xd)/(2.*rr0)),
        (x1 + (cv*xd + sv*yd)/(2.*rr0), y1 + (cv*yd - sv*xd)/(2.*rr0)),
    )

p0, p1 = sorted(solve(x1, y1, r1, x2, y2, r2))
print "%.08f %.08f %.08f %.08f" % (p0 + p1)

Solution (Python):
from math import sqrt
x1, y1, r1 = map(int, input().split())
x2, y2, r2 = map(int, input().split())

def intersection_cc(C0, r0, C1, r1):
    x0, y0 = C0; x1, y1 = C1
    xd = x1 - x0; yd = y1 - y0
    rr0 = xd**2 + yd**2
    rr1 = r0**2; rr2 = r1**2
    v_cos = (rr0 + rr1 - rr2)
    if 4*rr0*rr1 < v_cos**2:
        
        return tuple()
    """
    if 4*rr0*rr1 == v_cos**2:
        
        return tuple((x0 + v_cos*xd/(2.*rr0), y0 + v_cos*yd/(2.*rr0)))
    """
    
    v_sin = sqrt(4*rr0*rr1 - v_cos**2)
    return (
        (x0 + (v_cos*xd - v_sin*yd)/(2.*rr0), y0 + (v_cos*yd + v_sin*xd)/(2.*rr0)),
        (x0 + (v_cos*xd + v_sin*yd)/(2.*rr0), y0 + (v_cos*yd - v_sin*xd)/(2.*rr0))
    )

p0, p1 = sorted(intersection_cc((x1, y1), r1, (x2, y2), r2))
print("%.08f %.08f %.08f %.08f" % (p0 + p1))

Solution (Python):
import sys
from operator import itemgetter, attrgetter
from itertools import starmap
import cmath
from math import isinf, sqrt, acos, atan2
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 = sqrt(c.r * c.r - norm(pr - c.c))
    return Segment(pr + e * base, pr - e * base)
def getCrossPointsCL(c, l):
    pr = project(l, c.c)
    e = (l.se - l.fi) / abs(l.se - l.fi)
    base = sqrt(c.r * c.r - norm(pr - c.c))
    return Segment(pr + e * base, pr - e * base)
def getCrossPointsCC(c1, c2):
    d = abs(c1.c - c2.c)
    a = acos((c1.r * c1.r + d * d - c2.r * c2.r) / (2.0 * c1.r * d))
    t = cmath.phase(c2.c - c1.c)
    return Segment(c1.c + cmath.rect(c1.r, t + a), c1.c + cmath.rect(c1.r, t - a))
x, y, r = map(int, readline().split())
c1 = Circle(complex(x, y), r)
x, y, r = map(int, readline().split())
c2 = Circle(complex(x, y), r)
cpc = getCrossPointsCC(c1, c2)
if cpc.se.real < cpc.fi.real or (cpc.se.real == cpc.fi.real and cpc.se.imag < cpc.fi.imag):
    cpc.fi, cpc.se = cpc.se, cpc.fi
print(*map('{:.10f}'.format, (cpc.fi.real, cpc.fi.imag, cpc.se.real, cpc.se.imag)))

Solution (Python):
import sys
from operator import itemgetter, attrgetter
from itertools import starmap
import cmath
from math import isinf, sqrt, acos, atan2
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 = sqrt(c.r * c.r - norm(pr - c.c))
    return Segment(pr + e * base, pr - e * base)
def getCrossPointsCL(c, l):
    pr = project(l, c.c)
    e = (l.se - l.fi) / abs(l.se - l.fi)
    base = sqrt(c.r * c.r - norm(pr - c.c))
    return Segment(pr + e * base, pr - e * base)
def getCrossPointsCC(c1, c2):
    d = abs(c1.c - c2.c)
    a = acos((c1.r * c1.r + d * d - c2.r * c2.r) / (2.0 * c1.r * d))
    t = cmath.phase(c2.c - c1.c)
    return Segment(*sorted((c1.c + cmath.rect(c1.r, t + a), c1.c + cmath.rect(c1.r, t - a)), key=attrgetter('real', 'imag')))
x, y, r = map(int, readline().split())
c1 = Circle(complex(x, y), r)
x, y, r = map(int, readline().split())
c2 = Circle(complex(x, y), r)
cpc = getCrossPointsCC(c1, c2)
print(*map('{:.10f}'.format, (cpc.fi.real, cpc.fi.imag, cpc.se.real, cpc.se.imag)))

Solution (Python):
import sys
from operator import itemgetter, attrgetter
from itertools import starmap
import cmath
from math import isinf, sqrt, acos, atan2
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 getCrossPointsCL(c, l):
    pr = project(l, c.c)
    e = (l.se - l.fi) / abs(l.se - l.fi)
    base = sqrt(c.r * c.r - norm(pr - c.c))
    return Segment(*sorted((pr + e * base, pr - e * base)), key=attrgetter('real', 'imag'))
def getCrossPointsCC(c1, c2):
    d = abs(c1.c - c2.c)
    a = acos((c1.r * c1.r + d * d - c2.r * c2.r) / (2.0 * c1.r * d))
    t = cmath.phase(c2.c - c1.c)
    return Segment(*sorted((c1.c + cmath.rect(c1.r, t + a), c1.c + cmath.rect(c1.r, t - a)), key=attrgetter('real', 'imag')))
x, y, r = map(int, readline().split())
c1 = Circle(complex(x, y), r)
x, y, r = map(int, readline().split())
c2 = Circle(complex(x, y), r)
cpc = getCrossPointsCC(c1, c2)
print(*map('{:.10f}'.format, (cpc.fi.real, cpc.fi.imag, cpc.se.real, cpc.se.imag)))