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):


import sys

sys.setrecursionlimit(10 ** 9)
def input(): return sys.stdin.readline().strip()
def INT(): return int(input())
def MAP(): return map(int, input().split())
def LIST(): return list(map(int, input().split()))
INF=float('inf')

class Geometry:

    EPS = 10 ** -9

    def add(self, a, b):
        x1, y1 = a
        x2, y2 = b
        return (x1+x2, y1+y2)

    def sub(self, a, b):
        x1, y1 = a
        x2, y2 = b
        return (x1-x2, y1-y2)

    def mul(self, a, b):
        x1, y1 = a
        if not isinstance(b, tuple):
            return (x1*b, y1*b)
        x2, y2 = b 
        return (x1*x2, y1*y2)

    def div(self, a, b):
        x1, y1 = a
        if not isinstance(b, tuple):
            return (x1/b, y1/b)
        x2, y2 = b
        return (x1/x2, y1/y2)

    def abs(self, a):
        from math import hypot
        x1, y1 = a
        return hypot(x1, y1)

    def norm(self, a):
        x, y = a
        return x**2 + y**2

    def dot(self, a, b):
        x1, y1 = a
        x2, y2 = b
        return x1*x2 + y1*y2

    def cross(self, a, b):
        x1, y1 = a
        x2, y2 = b
        return x1*y2 - y1*x2

    def project(self, seg, p):
        """ 線分segに対する点pの射影 """
        p1, p2 = seg
        base = self.sub(p2, p1)
        r = self.dot(self.sub(p, p1), base) / self.norm(base)
        return self.add(p1, self.mul(base, r))

    def reflect(self, seg, p):
        """ 線分segを対称軸とした点pの線対称の点 """
        return self.add(p, self.mul(self.sub(self.project(seg, p), p), 2))

    def ccw(self, p0, p1, p2):
        """ 線分p0,p1から線分p0,p2への回転方向 """
        a = self.sub(p1, p0)
        b = self.sub(p2, p0)
        
        if self.cross(a, b) > self.EPS: return 1
        
        if self.cross(a, b) < -self.EPS: return -1
        
        if self.dot(a, b) < -self.EPS: return 2
        
        if self.norm(a) < self.norm(b): return -2
        
        return 0

    def intersect(self, seg1, seg2):
        """ 線分seg1と線分seg2の交差判定 """
        p1, p2 = seg1
        p3, p4 = seg2
        return (
            self.ccw(p1, p2, p3) * self.ccw(p1, p2, p4) <= 0
            and self.ccw(p3, p4, p1) * self.ccw(p3, p4, p2) <= 0
        )

    def get_distance_LP(self, line, p):
        """ 直線lineと点pの距離 """
        p1, p2 = line
        return abs(self.cross(self.sub(p2, p1), self.sub(p, p1)) / self.abs(self.sub(p2, p1)))

    def get_distance_SP(self, seg, p):
        """ 線分segと点pの距離 """
        p1, p2 = seg
        if self.dot(self.sub(p2, p1), self.sub(p, p1)) < 0: return self.abs(self.sub(p, p1))
        if self.dot(self.sub(p1, p2), self.sub(p, p2)) < 0: return self.abs(self.sub(p, p2))
        return self.get_distance_LP(seg, p)

    def get_distance_SS(self, seg1, seg2):
        """ 線分seg1と線分seg2の距離 """
        p1, p2 = seg1
        p3, p4 = seg2
        if self.intersect(seg1, seg2): return 0
        return min(
            self.get_distance_SP(seg1, p3), self.get_distance_SP(seg1, p4),
            self.get_distance_SP(seg2, p1), self.get_distance_SP(seg2, p2),
        )

    def get_cross_pointSS(self, seg1, seg2):
        """ 線分seg1と線分seg2の交点 """
        p1, p2 = seg1
        p3, p4 = seg2
        base = self.sub(p4, p3)
        dist1 = abs(self.cross(base, self.sub(p1, p3)))
        dist2 = abs(self.cross(base, self.sub(p2, p3)))
        t = dist1 / (dist1+dist2)
        return self.add(p1, self.mul(self.sub(p2, p1), t))

    def get_cross_pointCL(self, c, line):
        """ 円cと直線lineの交点 """
        from math import sqrt
        
        x, y, r = c
        p1, p2 = line
        pr = self.project(line, (x, y))
        e = self.div(self.sub(p2, p1), self.abs(self.sub(p2, p1)))
        base = sqrt(r*r - self.norm(self.sub(pr, (x, y))))
        return [self.add(pr, self.mul(e, base)), self.sub(pr, self.mul(e, base))]
    
    def arg(self, p):
        from math import atan2
        x, y = p
        return atan2(y, x)
    
    def polar(self, a, r):
        from math import sin, cos
        return (cos(r)*a, sin(r)*a)
    
    def get_cross_pointCC(self, c1, c2):
        """ 円c1と円c2の交点 """
        from math import acos
        
        x1, y1, r1 = c1
        x2, y2, r2 = c2
        d = self.abs(self.sub((x1, y1), (x2, y2)))
        a = acos((r1*r1+d*d-r2*r2) / (2*r1*d))
        t = self.arg(self.sub((x2, y2), (x1, y1)))
        return [self.add((x1, y1), self.polar(r1, t+a)), self.add((x1, y1), self.polar(r1, t-a))]

gm = Geometry()
x1, y1, r1 = MAP()
x2, y2, r2 = MAP()
c1 = (x1, y1, r1)
c2 = (x2, y2, r2)
res = gm.get_cross_pointCC(c1, c2)
res.sort()

res = sum(res, ())
for i in range(3):
    print('{:.10f}'.format(res[i]), end=' ')
print('{:.10f}'.format(res[-1]))

Solution (Python):
import math


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

    def __add__(self, v):
        return Vector2(self.x+v.x, self.y+v.y)

    def __sub__(self, v):
        return Vector2(self.x-v.x, self.y-v.y)

    def __mul__(self, v):
        return Vector2(self.x*v, self.y*v)

    def __truediv__(self, v):
        return Vector2(self.x/v, self.y/v)

    def __abs__(self):
        return math.sqrt(float(self.x*self.x+self.y*self.y))

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

    def __lt__(self, v):
        if math.isclose(self.x, v.x):
            return self.y < v.y
        return self.x < v.x

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

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

    def norm(self):
        d = abs(self)
        return Vector2(self.x/d, self.y/d)


def projection(v1, v2, p):
    v12 = (v2-v1).norm()
    v1p = p-v1
    d = v12.dot(v1p)
    return v1+v12*d


def distanceLP(v1, v2, p):
    '''
    v1 -> v2の直線とpとの距離
    '''
    return abs((v2-v1).cross(p-v1))/abs(v2-v1)


def distanceSP(v1, v2, p):
    '''
    v1 -> v2の線分とpとの距離
    '''
    if (v2-v1).dot(p-v1) < 0.0:
        return abs(p-v1)
    if (v1-v2).dot(p-v2) < 0.0:
        return abs(p-v2)
    return distanceLP(v1, v2, p)


def ccw(p0, p1, p2):
    c = (p1-p0).cross(p2-p0)
    if math.isclose(c, 0.0, abs_tol=1e-10):
        d = (p1-p0).dot(p2-p0)
        if d < 0.0:
            return -1
        else:
            d1 = abs(p1-p0)
            d2 = abs(p2-p0)
            if d1 < d2:
                return 1
            else:
                return 0
    elif c < 0.0:
        return -1
    else:
        return 1


def intersect(p1, p2, p3, p4):
    '''
    p1p2とp3p4の交差判定
    '''
    
    
    
    
    
    c1 = ccw(p1, p2, p3)*ccw(p1, p2, p4)
    c2 = ccw(p3, p4, p1)*ccw(p3, p4, p2)
    return c1 <= 0.0 and c2 <= 0.0


def distance(a1, a2, b1, b2):
    '''
    線分a1a2とb1b2の距離
    '''
    if intersect(a1, a2, b1, b2):
        return 0.0
    return min([
        min([distanceSP(a1, a2, b1), distanceSP(a1, a2, b2)]),
        min([distanceSP(b1, b2, a1), distanceSP(b1, b2, a2)])
    ])


cx1, cy1, r1 = map(int, input().split())
cx2, cy2, r2 = map(int, input().split())

c1 = Vector2(cx1, cy1)
c2 = Vector2(cx2, cy2)
c12 = c2-c1
d = abs(c12)

a = math.acos((r1*r1+d*d-r2*r2)/(2*r1*d))
t = math.atan2(c12.y, c12.x)
a1 = c1+Vector2(math.cos(t+a)*r1, math.sin(t+a)*r1)
a2 = c1+Vector2(math.cos(t-a)*r1, math.sin(t-a)*r1)
ans = ""
for a in sorted([a1, a2]):
    ans += f"{a.x:.8f} {a.y:.8f} "
print(ans.strip())

Solution (Python):


from typing import List, Tuple
from math import sqrt, sin, cos, acos, atan2
from enum import IntEnum

EPS = 1e-10


def float_equal(x: float, y: float) -> bool:
    return abs(x - y) < EPS


class PointLocation(IntEnum):
    COUNTER_CLOCKWISE = 1
    CCW = 1
    CLOCKWISE = -1
    CW = -1
    ONLINE_BACK = 2
    O_B = 2
    ONLINE_FRONT = -2
    O_F = -2
    ON_SEGMENT = 0
    O_S = 0


PL = PointLocation


class Point:

    def __init__(self, x: float=0.0, y: float=0.0) -> None:
        self.x = x
        self.y = y

    def __repr__(self) -> str:
        return "Point({}, {})".format(self.x, self.y)

    def __eq__(self, other: object) -> bool:
        if not isinstance(other, Point):
            
            return NotImplemented
        return float_equal(self.x, other.x) and \
            float_equal(self.y, other.y)

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

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

    def __mul__(self, k: float) -> 'Point':
        return Point(self.x * k, self.y * k)

    def __rmul__(self, k: float) -> 'Point':
        return self * k

    def __truediv__(self, k: float) -> 'Point':
        return Point(self.x / k, self.y / k)

    def __lt__(self, other: 'Point') -> bool:
        return self.y < other.y \
            if abs(self.x - other.x) < EPS \
            else self.x < other.x

    @staticmethod
    def polar(a: float, r: float) -> 'Point':
        return Point(a * cos(r), a * sin(r))

    def arg(self) -> float:
        return atan2(self.y, self.x)

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

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

    def dot(self, other: 'Point') -> float:
        return self.x * other.x + self.y * other.y

    def cross(self, other: 'Point') -> float:
        return self.x * other.y - self.y * other.x

    def is_orthogonal(self, other: 'Point') -> bool:
        return float_equal(self.dot(other), 0.0)

    def is_parallel(self, other: 'Point') -> bool:
        return float_equal(self.cross(other), 0.0)

    def distance(self, other: 'Point') -> float:
        return (self - other).abs()

    def in_side_of(self, seg: 'Segment') -> bool:
        return seg.vector().dot(
            Segment(seg.p1, self).vector()) >= 0

    def in_width_of(self, seg: 'Segment') -> bool:
        return \
            self.in_side_of(seg) and \
            self.in_side_of(seg.reverse())

    def distance_to_line(self, seg: 'Segment') -> float:
        return \
            abs((self - seg.p1).cross(seg.vector())) / \
            seg.length()

    def distance_to_segment(self, seg: 'Segment') -> float:
        if not self.in_side_of(seg):
            return self.distance(seg.p1)
        if not self.in_side_of(seg.reverse()):
            return self.distance(seg.p2)
        else:
            return self.distance_to_line(seg)

    def location(self, seg: 'Segment') -> PointLocation:
        p = self - seg.p1
        d = seg.vector().cross(p)
        if d > EPS:
            return PointLocation.COUNTER_CLOCKWISE
        if d < -EPS:
            return PointLocation.CLOCKWISE
        if seg.vector().dot(p) < 0.0:
            return PointLocation.ONLINE_BACK
        if seg.vector().norm() < p.norm():
            return PointLocation.ONLINE_FRONT
        return PointLocation.ON_SEGMENT


Vector = Point


class Segment:

    def __init__(self, p1: Point = None, p2: Point = None) -> None:
        self.p1: Point = Point() if p1 is None else p1
        self.p2: Point = Point() if p2 is None else p2

    def __repr__(self) -> str:
        return "Segment({}, {})".format(self.p1, self.p2)

    def __eq__(self, other: object) -> bool:
        if not isinstance(other, Segment):
            
            return NotImplemented
        return self.p1 == other.p1 and self.p2 == other.p2

    def vector(self) -> Vector:
        return self.p2 - self.p1

    def reverse(self) -> 'Segment':
        return Segment(self.p2, self.p1)

    def length(self) -> float:
        return self.p1.distance(self.p2)

    def is_orthogonal(self, other: 'Segment') -> bool:
        return self.vector().is_orthogonal(other.vector())

    def is_parallel(self, other: 'Segment') -> bool:
        return self.vector().is_parallel(other.vector())

    def projection(self, p: Point) -> Point:
        v = self.vector()
        vp = p - self.p1
        return v.dot(vp) / v.norm() * v + self.p1

    def reflection(self, p: Point) -> Point:
        x = self.projection(p)
        return p + 2 * (x - p)

    def intersects(self, other: 'Segment') -> bool:
        d0: PointLocation = self.p1.location(other)
        d1: PointLocation = self.p2.location(other)
        d2: PointLocation = other.p1.location(self)
        d3: PointLocation = other.p2.location(self)
        return d0 * d1 * d2 * d3 == 0 or \
            (d0 * d1 == -1 and d2 * d3 == -1)

    def intersection(self, other: 'Segment') -> Point:
        a = self.vector()
        b = other.vector()
        c = self.p1 - other.p1
        s = b.cross(c) / a.cross(b)
        return self.p1 + s * a

    def distance_with_segment(self, other: 'Segment') -> float:
        if not self.is_parallel(other) and \
                self.intersects(other):
            return 0
        else:
            return min(
                self.p1.distance_to_segment(other),
                self.p2.distance_to_segment(other),
                other.p1.distance_to_segment(self),
                other.p2.distance_to_segment(self))


Line = Segment


class Circle:

    def __init__(self, c: Point=None, r: float=0.0) -> None:
        self.c: Point = Point() if c is None else c
        self.r: float = r

    def __eq__(self, other: object) -> bool:
        if not isinstance(other, Circle):
            return NotImplemented
        return self.c == other.c and self.r == other.r

    def __repr__(self) -> str:
        return "Circle({}, {})".format(self.c, self.r)

    def cross_point_line(self, line: Line) -> List[Point]:
        proj = line.projection(self.c)
        dist = self.c.distance_to_line(line)
        tan = sqrt(self.r * self.r - dist * dist)
        u = line.vector() / line.vector().abs()
        return sorted([proj - tan * u, proj + tan * u])

    def cross_point_circle(self, other: 'Circle') -> List[Point]:
        vec = other.c - self.c
        d = vec.abs()
        r1 = self.r
        r2 = other.r
        a = acos((r1 * r1 + d * d - r2 * r2) / (2 * r1 * d))
        t = (other.c - self.c).arg()
        return sorted([self.c + Vector.polar(self.r, t + a),
                       self.c + Vector.polar(self.r, t - a)])


def main() -> None:
    c1x, c1y, c1r = [int(x) for x in input().split()]
    c1 = Circle(Point(c1x, c1y), c1r)
    c2x, c2y, c2r = [int(x) for x in input().split()]
    c2 = Circle(Point(c2x, c2y), c2r)
    left, right = c1.cross_point_circle(c2)
    print("{:.8f} {:.8f} {:.8f} {:.8f}".format(
        left.x, left.y, right.x, right.y))


if __name__ == "__main__":
    main()

Solution (Python):



from math import sqrt


def cross_point(circle1, circle2):
    x1, y1, r1 = circle1
    x2, y2, r2 = circle2

    if r1 < r2:
        x1, x2 = x2, x1
        y1, y2 = y2, y1
        r1, r2 = r2, r1

    p1p2 = (x2-x1, y2-y1)
    np = norm(p1p2)

    if np > (r1+r2)**2 or np < (r1-r2)**2:
        return []
    else:
        lp = sqrt(np)
        a = (p1p2[0]/lp, p1p2[1]/lp)
        b = ((y1-y2)/lp, (x2-x1)/lp)
        h = sqrt((r1+r2+lp)*(-r1+r2+lp)*(r1-r2+lp)*(r1+r2-lp)) / (2 * lp)
        t = sqrt(r1*r1 - h*h)
        if h < 1e-10:
            return [(x1 + t*a[0] + h*b[0], y1 + t*a[1] + h*b[1])]
        else:
            return [(x1 + t*a[0] + h*b[0], y1 + t*a[1] + h*b[1]),
                    (x1 + t*a[0] - h*b[0], y1 + t*a[1] - h*b[1])]


def norm(v):
    x, y = v
    return x**2 + y**2


def run():
    c1 = [int(i) for i in input().split()]
    c2 = [int(i) for i in input().split()]

    ps = cross_point(c1, c2)

    if len(ps) == 0:
        raise ValueError('two circles do not intersect')
    if len(ps) == 1:
        ps *= 2

    p1, p2 = ps
    if p1 > p2:
        p1, p2 = p2, p1

    print("{:.8f} {:.8f} {:.8f} {:.8f}".format(*p1, *p2))


if __name__ == '__main__':
    run()

Solution (Python):
import math
from typing import Union


class Point(object):
    __slots__ = ['x', 'y']

    def __init__(self, x, y):
        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: Union[int, float]):
        return Point(self.x * other, self.y * other)

    def norm(self):
        return pow(self.x, 2) + pow(self.y, 2)

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

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


class Vector(Point):
    __slots__ = ['x', 'y', 'pt1', 'pt2']

    def __init__(self, pt1: Point, pt2: Point):
        super().__init__(pt2.x - pt1.x, pt2.y - pt1.y)
        self.pt1 = pt1
        self.pt2 = pt2

    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 __repr__(self):
        return f"{self.pt1},{self.pt2}"


class Segment(Vector):
    __slots__ = ['x', 'y', 'pt1', 'pt2']

    def __init__(self, pt1: Point, pt2: Point):
        super().__init__(pt1, pt2)

    def projection(self, pt: Point)-> Point:
        t = self.dot(Vector(self.pt1, pt)) / self.norm()
        return self.pt1 + self * t

    def reflection(self, pt: Point) -> Point:
        return self.projection(pt) * 2 - pt

    def is_intersected_with(self, other) -> bool:
        if (self.point_geometry(other.pt1) * self.point_geometry(other.pt2)) <= 0\
                and other.point_geometry(self.pt1) * other.point_geometry(self.pt2) <= 0:
            return True
        else:
            return False

    def point_geometry(self, pt: Point) -> int:
        """
        [-2:"Online Back", -1:"Counter Clockwise", 0:"On Segment", 1:"Clockwise",  2:"Online Front"]
        """
        vec_pt1_to_pt = Vector(self.pt1, pt)
        cross = self.cross(vec_pt1_to_pt)
        if cross > 0:
            return -1   
        elif cross < 0:
            return 1    
        else:           
            dot = self.dot(vec_pt1_to_pt)
            if dot < 0:
                return -2    
            else:       
                if self.abs() < vec_pt1_to_pt.abs():
                    return 2    
                else:
                    return 0    

    def cross_point(self, other) -> Point:
        d1 = abs(self.cross(Vector(self.pt1, other.pt1)))       
        d2 = abs(self.cross(Vector(self.pt1, other.pt2)))       
        t = d1 / (d1 + d2)
        return other.pt1 + other * t

    def distance_to_point(self, pt: Point) -> Union[int, float]:
        vec_pt1_to_pt = Vector(self.pt1, pt)
        if self.dot(vec_pt1_to_pt) <= 0:
            return vec_pt1_to_pt.abs()
        vec_pt2_to_pt = Vector(self.pt2, pt)
        if Vector.dot(self * -1, vec_pt2_to_pt) <= 0:
            return vec_pt2_to_pt.abs()
        return (self.projection(pt) - pt).abs()

    def distance_to_segment(self, other) -> Union[int, float]:
        if self.is_intersected_with(other):
            return 0.0
        else:
            return min(
                self.distance_to_point(other.pt1),
                self.distance_to_point(other.pt2),
                other.distance_to_point(self.pt1),
                other.distance_to_point(self.pt2)
            )

    def __repr__(self):
        return f"{self.pt1},{self.pt2}"

class Circle(Point):
    __slots__ = ['x', 'y', 'r']
    def __init__(self, x, y, r):
        super().__init__(x, y)
        self.r = r

    def cross_point_with_circle(self, other):
        vec_self_to_other = Vector(self, other)
        vec_abs = vec_self_to_other.abs()
        t = ((pow(self.r, 2) - pow(other.r, 2)) / pow(vec_abs, 2) + 1) / 2
        pt = (other - self) * t
        abs_from_pt = math.sqrt(pow(self.r, 2) - pt.norm())
        inv = Point(vec_self_to_other.y / vec_abs,  - vec_self_to_other.x / vec_abs) * abs_from_pt
        pt_ = self + pt
        return (pt_ + inv), (pt_ - inv)

    def __repr__(self):
        return f"({self.x},{self.y}), {self.r}"


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

    pt1, pt2 = Circle(c1x, c1y, c1r).cross_point_with_circle(Circle(c2x, c2y, c2r))

    if pt1.x > pt2.x or (pt1.x == pt2.x and pt1.y > pt2.y):
        tmp = pt1
        pt1 = pt2
        pt2 = tmp

    print(" ".join(map('{:.10f}'.format, [pt1.x, pt1.y, pt2.x, pt2.y])))

    return


main()

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):
import math

EPS = 1e-10

def equals(a, b):
    return abs(a - b) < EPS


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

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

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

    def __mul__(self, a):
        return Point(self.x * a, self.y * a)

    def __rmul__(self, a):
        return self * a

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

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

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

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

    def __eq__(self, p):
        return equals(self.x, p.x) and equals(self.y, p.y)


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


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


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


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


def arg(p):
    return math.atan2(p.y, p.x)


def polar(r, a):
    return Point(math.cos(a) * r, math.sin(a) * r)


def getCrossPoints(c1, c2):
    d = (c1.c - c2.c).abs()
    a = math.acos((c1.r*c1.r + d*d - c2.r*c2.r) / (2 * c1.r * d))
    t = arg(c2.c - c1.c)
    return (c1.c + polar(c1.r, t + a), c1.c + polar(c1.r, t - a))


if __name__ == '__main__':
    c1x, c1y, c1r =  [int(v) for v in input().split()]
    c2x, c2y, c2r =  [int(v) for v in input().split()]
    c1 = Circle(Point(c1x, c1y), c1r)
    c2 = Circle(Point(c2x, c2y), c2r)
    v1, v2 = getCrossPoints(c1, c2)
    if v2 < v1:
        v1, v2 = v2, v1
    ans = [v1.x, v1.y, v2.x, v2.y]

    print('{0[0]:.8f} {0[1]:.8f} {0[2]:.8f} {0[3]:.8f}'.format(ans))

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

vx, vy = x2 - x1, y2 - y1
norm = vx * vx + vy * vy
add = math.sqrt(norm)
dec = (r1 * r1 - r2 * r2) / add
p = (add + dec) / 2 
cx = p * vx / add + x1
cy = p * vy / add + y1
q = math.sqrt(max(0, r1 * r1 - p * p))
dx = q * vy / add
dy = -q * vx / add
ans = sorted([(cx + dx, cy + dy), (cx - dx, cy - dy)])
print "%.10f %.10f %.10f %.10f" % (ans[0][0], ans[0][1], ans[1][0], ans[1][1])