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

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

"""
import sys
from math import sqrt, atan2, acos, sin, cos
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_points(self, s):
        if isinstance(s, Segment):
            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
        elif isinstance(s, Cirle):
            c2 = s
            d = abs(self.c - c2.c)
            a = acos((self.r * self.r + d * d - c2.r * c2.r) / (2 * self.r * d))
            t = atan2(c2.c.y - self.c.y, c2.c.x - self.c.x)
            temp1 = Point(cos(t+a)*self.r, sin(t+a)*self.r)
            temp2 = Point(cos(t-a)*self.r, sin(t-a)*self.r)
            return self.c + temp1, self.c + temp2


def main(args):
    x1, y1, r1 = map(int, input().split())
    c1 = Cirle(Point(x1, y1), r1)
    x2, y2, r2 = map(int, input().split())
    c2 = Cirle(Point(x2, y2), r2)

    ans1, ans2 = c1.cross_points(c2)
    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, sin, cos, acos, atan2
c1x,c1y,r1 = map(int,input().split())
c2x,c2y,r2 = map(int,input().split())
c1 = [c1x, c1y]
c2 = [c2x, c2y]
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))
def polar(a,r):
    return [cos(r) * a, sin(r) * a]

d = norm(sub(c2,c1))
a = acos((r1**2 + d **2 - r2**2)/(2 * r1 * d))
c2c1 = sub(c2,c1)
t = atan2(c2c1[1],c2c1[0])
ans1 = add(c1, polar(r1, t+a))
ans2 = add(c1, polar(r1, t-a))

ans = sorted([ans1] + [ans2], key=lambda x: (x[0],x[1]))
print(' '.join([str(round(i, 9)) for o in ans for i in o]))

Solution (Python):


import collections
import math


class Vector2(collections.namedtuple("Vector2", ["x",  "y"])):

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

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

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

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

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

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

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

    def abs2(self):
        return float(self.x * self.x + self.y * self.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 getDistanceSP(segment, point):
    p = point
    p1, p2 = segment
    if (p2 - p1).dot(p - p1) < 0:
        return abs(p - p1)
    if (p1 - p2).dot(p - p2) < 0:
        return abs(p - p2)
    return abs((p2 - p1).cross(p - p1)) / abs(p2 - p1)


def getDistance(s1, s2):
    a, b = s1
    c, d = s2
    if intersect(s1, s2):  
        return 0
    return min(getDistanceSP(s1, c), getDistanceSP(s1, d), getDistanceSP(s2, a), getDistanceSP(s2, b))


def ccw(p0, p1, p2):
    a = p1 - p0
    b = p2 - p0
    if a.cross(b) > 0:
        return 1  
    elif a.cross(b) < 0:
        return -1  
    elif a.dot(b) < 0:
        return 2  
    elif abs(a) < abs(b):
        return -2  
    else:
        return 0  


def intersect(s1, s2):
    a, b = s1
    c, d = s2
    return ccw(a, b, c) * ccw(a, b, d) <= 0 and ccw(c, d, a) * ccw(c, d, b) <= 0


def project(l, p):
    p1, p2 = l
    base = p2 - p1
    hypo = p - p1
    return p1 + base * (hypo.dot(base) / abs(base)**2)


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


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


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


def getCrossPointsCircle(c1, c2):
    base = c2.c - c1.c
    d = abs(base)
    a = math.acos((c1.r**2 + d**2 - c2.r**2) / (2 * c1.r * d))
    t = math.atan2(base.y, base.x)
    return [c1.c + polar(c1.r, t + a), c1.c + polar(c1.r, t - a)]


if __name__ == '__main__':

    a, b, r = map(int, input().split())
    c1 = Circle(Vector2(a, b), r)
    a, b, r = map(int, input().split())
    c2 = Circle(Vector2(a, b), r)

    ans = getCrossPointsCircle(c1, c2)
    ans = sorted(ans, key=lambda x: (x.x, x.y))

    print("{:.8f} {:.8f} {:.8f} {:.8f}".format(
        ans[0].x, ans[0].y, ans[1].x, ans[1].y))

Solution (Python):
import math

class Vector:
    def __init__(self,x,y):
        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,scalar):
        return Vector(self.x*scalar,self.y*scalar)
    def __rmul__(self, scalar):
        return Vector(self.x*scalar,self.y*scalar)
    
   
    def __repr__(self):
        return str([self.x,self.y])
        
    def norm_2(self):
        return dot(self,self)
    
    def norm(self):
        return math.sqrt(self.norm_2())
    
def v_sum(v1,v2):
    return Vector(v1.x+v2.x,v1.y+v2.y) 
def scalar_multi(k,v):
    return Vector(k*v.x,k*v.y)
def v_diff(v1,v2):
    return v_sum(v1,scalar_multi(-1,v2))
    
def dot(vector1,vector2):
    return vector1.x*vector2.x+vector1.y*vector2.y

def cross(vector1,vector2):
    return vector1.x*vector2.y-vector1.y*vector2.x

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

d = c2-c1

cos = (r1**2+d.norm_2()-r2**2)/(2*r1*d.norm())
h = r1*cos/d.norm()*d + c1
l = r1*math.sqrt(1-cos**2)
n = Vector(d.y,-d.x)

a = h + l/n.norm()*n
b = h - l/n.norm()*n

if abs(a.x)<10**(-9):
    a.x = 0
if abs(a.y)<10**(-9):
    a.y = 0
if abs(b.x)<10**(-9):
    b.x = 0
if abs(b.y)<10**(-9):
    b.y = 0

if a.x<b.x:
    print('{} {} {} {}'.format(a.x,a.y,b.x,b.y))
elif a.x>b.x:
    print('{} {} {} {}'.format(b.x,b.y,a.x,a.y))
else:
    if a.y<b.y:
        print('{} {} {} {}'.format(a.x,a.y,b.x,b.y))
    else:
        print('{} {} {} {}'.format(b.x,b.y,a.x,a.y))

Solution (Python):
import math
EPS = 10**(-9)
def is_equal(a,b):
    return abs(a-b) < EPS

def norm(v,i=2):
    ret = 0
    n = len(v)
    for j in range(n):
        ret += abs(v[j])**i
    return math.pow(ret,1/i)

class Vector(list):
    """
    ベクトルクラス
    対応演算子
    +  : ベクトル和 
    -  : ベクトル差
    *  : スカラー倍、または内積
    /  : スカラー除法
    ** : 外積
    += : ベクトル和
    -= : ベクトル差
    *= : スカラー倍
    /= : スカラー除法

    メソッド
    self.norm(i) : L{i}ノルムを計算
    """
    def __add__(self,other):
        n = len(self)
        ret = [0]*n
        for i in range(n):
            ret[i] = super().__getitem__(i) + other.__getitem__(i)
        return self.__class__(ret)
    
    def __radd__(self,other):
        n = len(self)
        ret = [0]*n
        for i in range(n):
            ret[i] = other.__getitem__(i) + super().__getitem__(i)
        return self.__class__(ret)
    
    def __iadd__(self, other):
        n = len(self)
        for i in range(n):
            self[i] += other.__getitem__(i)
        return self

    def __sub__(self,others):
        n = len(self) 
        ret = [0]*n
        for i in range(n):
            ret[i] = super().__getitem__(i) - others.__getitem__(i)
        return self.__class__(ret)

    def __isub__(self, other):
        n = len(self)
        for i in range(n):
            self[i] -= other.__getitem__(i)
        return self

    def __rsub__(self,others):
        n = len(self) 
        ret = [0]*n
        for i in range(n):
            ret[i] = others.__getitem__(i) - super().__getitem__(i)
        return self.__class__(ret)
    
    def __mul__(self,other):
        n = len(self)
        if isinstance(other,list):
            ret = 0
            for i in range(n):
                ret += super().__getitem__(i)*other.__getitem__(i)
            return ret
        else:
            ret = [0]*n
            for i in range(n):
                ret[i] = super().__getitem__(i)*other
            return self.__class__(ret)

    def __rmul__(self,other):
        n = len(self)
        if isinstance(other,list):
            ret = 0
            for i in range(n):
                ret += super().__getitem__(i)*other.__getitem__(i)
            return ret
        else:
            ret = [0]*n
            for i in range(n):
                ret[i] = super().__getitem__(i)*other
            return self.__class__(ret)
    
    
    def __truediv__(self,other):
        """
        ベクトルのスカラー除法
        Vector/scalar
        """
        n = len(self)
        ret = [0]*n
        for i in range(n):
            ret[i] = super().__getitem__(i)/other
        return self.__class__(ret)
    
    def norm(self,i):
        """
        L{i}ノルム
        self.norm(i)
        """
        return norm(self,i)
    
    def __pow__(self,other):
        """
        外積
        self**other
        """
        n = len(self)
        ret = [0]*3
        x = self[:]
        y = other[:]
        if n == 2:
            x.append(0)
            y.append(0)
        if n == 2 or n == 3:
            for i in range(3):
                ret[0],ret[1],ret[2] = x[1]*y[2]-x[2]*y[1],x[2]*y[0]-x[0]*y[2],x[0]*y[1]-x[1]*y[0]
            ret = Vector(ret)
            if n == 2:
                return ret
            else:
                return ret
    def angle(self):
        if len(self) != 2:
            return False
        else:
            if is_equal(super().__getitem__(0), 0) and is_equal(super().__getitem__(1), 0):
                return False
            elif is_equal(super().__getitem__(0), 0):
                if super().__getitem__(1) > 0:
                    return math.pi/2
                elif super().__getitem__(1) < 0:
                    return -math.pi/2
            else:
                if super().__getitem__(0) > 0:
                    return math.atan(super().__getitem__(1)/super().__getitem__(0))
                elif super().__getitem__(0) < 0:
                    return math.atan(super().__getitem__(1)/super().__getitem__(0)) + math.pi


class Segment:
    """
    線分クラス

    methods
    length()                 :
    get_unit_vec()           : 
    projection(vector)       :
    is_vertical(segment)     :
    is_horizontal(segment)   :
    reflection()             :
    include()                :
    distance()               :
    ccw()                    :
    intersect()              :
    intersect_point          :
    """
    def __init__(self,v1,v2):
        if v2[0] < v1[0]:
            v1,v2 = v2,v1
        self.v1 = v1
        self.v2 = v2
    
    def length(self):
        return norm(self.v1-self.v2)

    def get_unit_vec(self):
        
        dist = norm(self.v2-self.v1)
        if dist != 0:
            return (self.v2-self.v1)/dist
        else:
            return False
    
    def projection(self,vector):
        
        unit_vec = self.get_unit_vec()
        t = unit_vec*(vector-self.v1)
        return self.v1 + t*unit_vec
    
    def is_vertical(self,other):
        
        return is_equal(0,self.get_unit_vec()*other.get_unit_vec())
    
    def is_horizontal(self,other):
        
        return is_equal(0,self.get_unit_vec()**other.get_unit_vec())
    
    def reflection(self,vector):
        
        projection = self.projection(vector)
        v = projection - vector
        return projection + vector
    
    def include(self,vector):
        
        proj = self.projection(vector)
        if not is_equal(norm(proj-vector),0):
            return False
        else:
            n = len(self.v1)
            f = True
            for i in range(n):
                f &= ((self.v1[i] <= vector[i] <= self.v2[i]) or (self.v2[i] <= vector[i] <=self.v1[i]))
            return f
    
    def distance(self,other):
        
        if isinstance(other,Vector):
            proj = self.projection(other)
            if self.include(proj):
                return norm(proj-other)
            else:
                ret = []
                ret.append(norm(self.v1-other))
                ret.append(norm(self.v2-other))
                return min(ret)

    def ccw(self,vector):
        """
        線分に対して点が反時計回りの位置にある(1)か時計回りの位置にある(-1)か線分上にある(0)か
        ただし、直線上にはあるが線分上にはない場合は反時計回りの位置にあると判定して1を返す。
        """
        direction = self.v2 - self.v1
        v = vector - self.v1
        if self.include(vector):
            return 0
        else:
            cross = direction**v
            if cross[2] <= 0:
                return 1
            else:
                return -1
    
    def intersect(self,segment):
        """
        線分の交差判定
        """
        ccw12 = self.ccw(segment.v1)
        ccw13 = self.ccw(segment.v2)
        ccw20 = segment.ccw(self.v1)
        ccw21 = segment.ccw(self.v2)

        if ccw12*ccw13*ccw20*ccw21 == 0:
            return True
        else:
            if ccw12*ccw13 < 0 and ccw20*ccw21 < 0:
                return True
            else:
                return False
    
    def intersect_point(self,segment):
        """
        線分の交点
        """
        d1 = norm(self.projection(segment.v1) - segment.v1)
        d2 = norm(self.projection(segment.v2) - segment.v2)
        return segment.v1 + (d1/(d1 + d2)) * (segment.v2 - segment.v1)
    
    def angle(self):
        if is_equal(self.v2[0] - self.v1[0], 0):
            return math.pi/2
        else:
            return math.atan((self.v2[1] - self.v1[1])/(self.v2[0] - self.v1[0]))



class Line(Segment):
    """
    直線クラス
    """
    
    def include(self,vector):
        proj = self.projection(vector)
        return is_equal(norm(proj-vector),0)
    
    def distance(self,vector):
        return norm(self.projection(vector) - vector)


class Circle():
    def __init__(self,vector,r):
        self.c = vector
        self.r = r
    
    def intersect_point(self,other):
        if isinstance(other, Line):
            if other.distance(self.c) > self.r:
                return False
            else:
                d = other.distance(self.c)
                proj = other.projection(self.c)
                unit_v = other.get_unit_vec()
                t = (self.r**2 - d**2)**0.5
                return proj + t*unit_v, proj - t*unit_v

        if isinstance(other,Circle):
            if norm(self.c - other.c) > self.r + other.r:
                return False
            else:
                d = norm(self.c - other.c)
                alpha = math.acos((self.r**2 + d**2 - other.r**2)/(2*d*self.r))
                L = other.c - self.c
                theta = L.angle()
                to1 = self.r*Vector([math.cos(theta + alpha), math.sin(theta + alpha)])
                to2 = self.r*Vector([math.cos(theta - alpha), math.sin(theta - alpha)])
                return self.c+ to1, self.c + to2

x1,y1,r1 = list(map(int,input().split()))
x2,y2,r2 = list(map(int,input().split()))
circle1 = Circle(Vector([x1,y1]), r1)
circle2 = Circle(Vector([x2,y2]), r2)
p1,p2 = circle1.intersect_point(circle2)
if p2[0] < p1[0]:
    p1,p2 = p2,p1
elif p2[0] == p1[0]:
    p1,p2 = p2,p1
print("{:.7f} {:.7f} {:.7f} {:.7f}".format(*p1,*p2))

Solution (Python):
import math

class Point:
    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):
        return Point(self.x * other, self.y * other)
    def __truediv__(self, other):
        return Point(self.x / other, self.y / other)
    def __repr__(self):
        return str(round(self.x, 8)) + ' ' + str(round(self.y, 8))
    def __lt__(self, other):
        if self.x == other.x:
            return self.y < other.y
        else:
            return self.x < other.x

class Vector(Point):
    pass

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

class Segment(Line):
    pass

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

def points_to_vector(p1, p2):
    x = p1.x - p2.x
    y = p1.y - p2.y
    return Vector(x, y)

def vector(p):
    return Vector(p.x, p.y)

def dot(v1, v2):
    return v1.x * v2.x + v1.y * v2.y

def cross(v1, v2):
    return v1.x * v2.y - v1.y * v2.x

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

def distance(v):
    return math.sqrt(norm(v))

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

def reflect(s, p):
    return p + (project(s, p) -p) * 2

def get_distance(s1, s2):
    if intersect_s(s1, s2):
        return 0
    d1 = get_distance_sp(s1, s2.p1)
    d2 = get_distance_sp(s1, s2.p2)
    d3 = get_distance_sp(s2, s1.p1)
    d4 = get_distance_sp(s2, s1.p2)
    return min(d1, min(d2, min(d3, d4)))

def get_distance_pp(p1, p2):
    return distance(p1 - p2)

def get_distance_lp(l, p):
    return abs(cross(l.p2 - l.p1, p - l.p1) / distance(l.p2 - l.p1))

def get_distance_sp(s, p):
    if dot(s.p2 - s.p1, p - s.p1) < 0:
        return distance(p - s.p1)
    elif dot(s.p1 - s.p2, p - s.p2) < 0:
        return distance(p - s.p2)
    else:
        return get_distance_lp(s, p)

def ccw(p0, p1, p2):
    EPS = 1e-10
    COUNTER_CLOCKWISE = 1
    CLOCKWISE = -1
    ONLINE_BACK = 2
    ONLINE_FRONT = -2
    ON_SEGMENT = 0
    v1 = p1 - p0
    v2 = p2 - p0
    if cross(v1, v2) > EPS:
        return COUNTER_CLOCKWISE
    elif cross(v1, v2) < -EPS:
        return CLOCKWISE
    elif dot(v1, v2) < -EPS:
        return ONLINE_BACK
    elif norm(v1) < norm(v2):
        return ONLINE_FRONT
    else:
        return ON_SEGMENT

def intersect_p(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 intersect_s(s1, s2):
    return intersect_p(s1.p1, s1.p2, s2.p1, s2.p2)

def get_cross_point(s1, s2):
    base = s2.p2 - s2.p1
    d1 = abs(cross(base, s1.p1 - s2.p1))
    d2 = abs(cross(base, s1.p2 - s2.p1))
    t = d1 / (d1 + d2)
    return s1.p1 + (s1.p2 - s1.p1) * t

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

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

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

def get_cross_points_cirlces(c1, c2):
    d = distance(c1.c - c2.c)
    a = math.acos((c1.r**2 + d * d - c2.r**2) / (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)

import sys


temp = list(map(int, input().split()))
c1 = Circle(Point(temp[0], temp[1]), temp[2])
temp = list(map(int, input().split()))
c2 = Circle(Point(temp[0], temp[1]), temp[2])
ans = get_cross_points_cirlces(c1, c2)
print(min(ans), max(ans))

Solution (Python):
from functools import singledispatch
import math

EPS = 1e-10

COUNTER_CLOCKWISE = 1
CLOCKWISE = -1
ONLINE_BACK = 2
ONLINE_FRONT = -2
ON_SEGMENT = 0

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

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

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

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

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

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

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

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

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

    def __eq__(self, point):
        
        return abs(self.x-point.x) < EPS and abs(self.y-point.y) <EPS

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 isOrthogonal(a, b):
    return dot(a, b) == 0

def isParallel(a, b):
    return cross(a, b) == 0

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

def reflect(s, p):
    return p+(project(s, p)-p)*2

@singledispatch
def get_distance(a, b):
    return abs(a-b)

def get_distance_lp(l, p):
    return abs(cross(l.p2-l.p1, p-l.p1)/abs(l.p2-l.p1))

def get_distance_sp(s, p):
    if dot(s.p2-s.p1, p-s.p1) < 0:
        return abs(p-s.p1)
    if dot(s.p1-s.p2, p-s.p2) < 0:
        return abs(p-s.p2)
    return get_distance_lp(s, p)
    
@get_distance.register(Segment)
def _(s1, s2):
    if intersect(s1, s2):
        return 0
    return min([get_distance_sp(s1, s2.p1), get_distance_sp(s1, s2.p2),
                get_distance_sp(s2, s1.p1), get_distance_sp(s2, s1.p2)])

def ccw(p0, p1, p2):
    a = p1 - p0
    b = p2 - p0
    if cross(a, b) > EPS:
        return COUNTER_CLOCKWISE
    if cross(a, b) < -EPS:
        return CLOCKWISE
    if dot(a, b) < -EPS:
        return ONLINE_BACK
    if a.norm() < b.norm():
        return ONLINE_FRONT
    return ON_SEGMENT

@singledispatch
def intersect(p1, p2, p3, p4):
    return (ccw(p1, p2, p3)*ccw(p1, p2, p4) <= 0
            and ccw(p3, p4, p1)*ccw(p3, p4, p2) <= 0)

@intersect.register(Segment)
def _(s1, s2):
    return intersect(s1.p1, s1.p2, s2.p1, s2.p2)

def get_cross_point(s1, s2):
    base = s2.p2 - s2.p1
    d1 = abs(cross(base, s1.p1-s2.p1))
    d2 = abs(cross(base, s1.p2-s2.p1))
    t = d1/(d1+d2)
    return s1.p1 + (s1.p2-s1.p1)*t

def get_cross_points(c, l):
    
    
    pr = project(l, c.c)
    e = (l.p2 - l.p1)/abs(l.p2-l.p1)
    base = (c.r**2-(pr-c.c).norm())**0.5
    return (pr+e*base, pr-e*base)

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

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

def get_cross_points_(c1, c2):
    
    
    
    
    d = abs(c1.c-c2.c)
    a = math.acos((c1.r**2+d**2-c2.r**2)/(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__':
    from sys import stdin
    input = stdin.readline

    cx, cy, r = map(int, input().split())
    center = Point(cx, cy)
    c1 = Circle(center, r)

    cx, cy, r = map(int, input().split())
    center = Point(cx, cy)
    c2 = Circle(center, r)

    ans1, ans2 = get_cross_points_(c1, c2)
    
    
    for ans in [ans1, ans2]:
        if round(ans.x, 10) == 0:
            ans.x = 0
        if round(ans.y, 10) == 0:
            ans.y = 0
    if ans1.x < ans2.x or (ans1.x == ans2.x and ans1.y < ans2.y):
        
        print(ans1.x, ans1.y, ans2.x, ans2.y)
    else:
        print(ans2.x, ans2.y, ans1.x, ans1.y)

Solution (Python):

from collections import defaultdict, deque
from heapq import heappush, heappop
from itertools import permutations, accumulate
import sys
import math
import bisect
def LI(): return [int(x) for x in sys.stdin.readline().split()]
def I(): return int(sys.stdin.readline())
def LS():return [list(x) for x in sys.stdin.readline().split()]
def S():
    res = list(sys.stdin.readline())
    if res[-1] == "\n":
        return res[:-1]
    return res
def IR(n):
    return [I() for i in range(n)]
def LIR(n):
    return [LI() for i in range(n)]
def SR(n):
    return [S() for i in range(n)]
def LSR(n):
    return [LS() for i in range(n)]

sys.setrecursionlimit(1000000)
mod = 1000000007

def solve():
    def D2(x1,y1,x2,y2):
        return (x1-x2)**2+(y1-y2)**2

    
    def check(R1,R2,eps = 1e-6):
        x1,y1,r1 = R1
        x2,y2,r2 = R2
        d2 = D2(x1,y1,x2,y2)
        if (r1-r2)**2 <= d2 <= (r1+r2)**2:
            return 1
        else:
            return 0

    
    def point(R1,R2):
        if not check(R1,R2):
            return []
        else:
            x1,y1,r1 = R1
            x2,y2,r2 = R2
            x2 -= x1
            y2 -= y1
            r = x2**2+y2**2
            a = (r+r1**2-r2**2)/2
            d = (r*r1**2-a**2)**0.5
            X,Y = a*x2/r, a*y2/r
            k = d/r
            X1,Y1 = X+y2*k+x1, Y-x2*k+y1
            X2,Y2 = X-y2*k+x1, Y+x2*k+y1
            return [(X1,Y1),(X2,Y2)]

    r1 = LI()
    r2 = LI()
    p1,p2 = point(r1,r2)
    if p2 < p1:
        p1,p2 = p2,p1
    print(*p1,*p2)
    return


if __name__ == "__main__":
    solve()