Problem: Find the tangent lines between a point $p$ and a circle $c$.

Constraints

* $-1,000 \leq px, py, cx, cy \leq 1,000$
* $1 \leq r \leq 1,000$
* Distance between $p$ and the center of $c$ is greater than the radius of $c$.

Input

The input is given in the following format.

$px \; py$
$cx \; cy \; r$

$px$ and $py$ represents the coordinate of the point $p$. $cx$, $cy$ and $r$ represents the center coordinate and radius of the circle $c$ respectively. All input values are given in integers.

Output

Print coordinates of the tangent points on the circle $c$ based on the following rules.

* 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.00001.

Examples

Input

0 0
2 2 2


Output

0.0000000000 2.0000000000
2.0000000000 0.0000000000


Input

-3 0
2 2 2


Output

0.6206896552 3.4482758621
2.0000000000 0.0000000000

Solution (in CPP):

Solution:
#include <iostream>
#include <vector>
#include <cstdio>
#include <cmath>
#include <algorithm>
using namespace std;
const int N = 200000 + 10;
const double EPS = 1e-7;
const double INF = 1e12;
int sgn(double x) {
    if (x < -EPS) return -1;
    if (x > +EPS) return +1;
    return 0;
}
struct Vec {
    double x,y;
    Vec() {}
    Vec(double x, double y) { 
        this->x = x, this->y = y; 
    }
    void read() { 
        scanf("%lf %lf", &x, &y); 
    }
    void prt() {
        printf("%.8f %.8f\n", x, y); 
    }
    double len() { 
        return sqrt(x * x + y * y); 
    }
    double len2() { 
        return x * x + y * y; 
    }
    bool operator <(const Vec& o)const {
        if(x == o.x)
            return y < o.y;
        return x < o.x;
    }
    Vec operator +(const Vec& o)const { 
        return Vec(x + o.x, y + o.y); 
    }
    Vec operator -(const Vec& o)const { 
        return Vec(x - o.x, y - o.y); 
    }
    Vec operator *(const double& k)const { 
        return Vec(k * x, k * y); 
    }
    double operator *(const Vec& o)const {
        return x * o.x + y * o.y; 
    }
    double operator ^(const Vec& o)const {
        return x * o.y - y * o.x; 
    }
    Vec rotate(double ang){ 
        return Vec(x * cos(ang) - y * sin(ang), x * sin(ang) + y * cos(ang)); 
    }
    Vec change(double l) { 
        if(len() < EPS) return *this; 
        return (*this) * (l/len()); 
    }
};
bool cmpang(Vec A, Vec B) {
    if ((A ^ B) == 0) {
        return A < B;
    }
    return (A ^ B) > 0;
}

struct Line {
    Vec A1, A2;
    Line() {}
    Line(Vec A1, Vec A2) {
        this->A1 = A1, this->A2 = A2;
    }
    double len() {
        return (A2-A1).len();
    }
    double Len2() {
        return (A2-A1).len2();
    }
    double getAng() const {
        return atan2(A2.y-A1.y, A2.x-A1.x);
    }
    void prt() {
        printf("(%.5f %.5f) - (%.5f %.5f) %.7f\n", A1.x,A1.y,A2.x,A2.y,(*this).getAng());
    }
    bool operator <(const Line& o)const {
        Line l = o;
        double A = (*this).getAng();
        double B = l.getAng();
        if (abs(A-B) < EPS)
            return ((A2-A1)^(l.A1-A1)) > EPS;
        return A < B;
    }   
};
// P to L
double disToLine(Vec P, Line L) {
    return abs((P-L.A1)^(L.A2-L.A1)) / L.len();
}
Vec projection(Vec P, Line L) {
    double shadowLen = (P - L.A1) * (L.A2 - L.A1);
    return L.A1 + (L.A2-L.A1).change(shadowLen);
}
Vec rejection(Vec P, Line L) {
    return projection(P,L) * 2 - P;
}
bool onseg(Vec P, Line L) {
    if(disToLine(P, L) > EPS) return 0;
    return (P-L.A1) * (P-L.A2) <= 0;
}
Vec Lineintersect(Line L1, Line L2) {
    // 0 ~ (L2.A1 - L1.A1) ^ (L2.A2 - L1.A1);
    // 1 ~ (L2.A1 - L1.A2) ^ (L2.A2 - L1.A2);
    double F0 = (L2.A1 - L1.A1) ^ (L2.A2 - L1.A1);
    double F1 = (L2.A1 - L1.A2) ^ (L2.A2 - L1.A2);
    if (abs(F1 - F0) < EPS) return Vec(INF, INF);
    return L1.A1 + (L1.A2 - L1.A1) * (- F0 / (F1 - F0));
}
Vec Segintersect(Line L1, Line L2) {
    Vec P = Lineintersect(L1, L2);
    if (P.x == INF) return P;
    if (onseg(P, L1) && onseg(P, L2)) return P;
    return Vec(INF, INF);
}
bool disToSeg(Vec P, Line L) {
    double ans = min((P-L.A1).len(), (P-L.A2).len());
    ans = min(ans, (P - projection(P,L)).len());
    return ans;
}

// 
vector<Vec> shiftConvexHull(vector<Vec> A) {
    double my = INF, mx = INF; int st = -1;
    for (int i = 0; i < A.size(); i ++) {
        if (A[i].y < my || (A[i].y == my && A[i].x < mx) ) {
            my = A[i].y; mx = A[i].x; st = i;
        }
    }
    vector<Vec> ans;
    for (int i = 0; i < A.size(); i ++) {
        ans.push_back(A[(i+st)%A.size()]);
    }
    return ans;
}

vector<Vec> convexHull(vector<Vec> v) {
    sort(v.begin(), v.end());
    vector<Vec> ans; ans.resize(2*v.size());
    int top = 0, n = v.size();
    if(n <= 1) return v;

    for(int i = 0; i < n; ans[top++] = v[i++]) {
        while (top > 1 && ((ans[top-1] - ans[top-2]) ^ (v[i] - ans[top-1])) < EPS) -- top;
    }
    for(int i = n-2, t = top; i >= 0; ans[top++] = v[i--]) {
        while (top > t && ((ans[top-1] - ans[top-2]) ^ (v[i] - ans[top-1])) < EPS) -- top;
    } 
    while(ans.size() >= top) ans.pop_back();
    return ans;
}
double PolyArea(vector<Vec> v) {
    double ans = 0;
    int n = v.size();
    for(int i = 0; i < v.size(); i ++) {
        ans = ans + (v[i] ^ v[(i+1)%n]);
    }
    return abs(ans) / 2;
}
vector<Vec> convexCut(vector<Vec> v, Vec P1, Vec P2) {
    vector<Vec> ans;
    for(int i=0;i<v.size();i++) {
        int d1 = sgn((P2-P1)^(v[i]-P1));
        int d2 = sgn((P2-P1)^(v[(i+1)%v.size()]-P1));
        if (d1 >= 0) ans.push_back(v[i]);
        if (d1 * d2 < 0) ans.push_back(Lineintersect(Line(P1,P2), Line(v[i],v[(i+1)%v.size()])));
    }
    return ans;
}
double cloestPair(vector<Vec> v) {
    if (v.size() <= 1) return INF;

    sort(v.begin(), v.end());
    vector<Vec> vL, vR;
    int mid = v.size() / 2;
    for (int i = 0; i < v.size(); i ++) {
        if (i < mid)
            vL.push_back(v[i]);
        else
            vR.push_back(v[i]);
    }
    double d = min(cloestPair(vL), cloestPair(vR));
    
    vector<Vec> strip;
    for (int i = 0; i < v.size(); i ++) {
        if (abs(v[mid].x - v[i].x) <= d + EPS) 
            strip.push_back(v[i]);
    }

    for (int i = 0; i < strip.size(); i ++) swap(strip[i].x, strip[i].y);
    sort(strip.begin(), strip.end());
    for (int i = 0; i < strip.size(); i ++) swap(strip[i].x, strip[i].y);

    for (int i = 0; i < strip.size(); i ++) {
        for (int j = i + 1; j < strip.size(); j ++) {
            if (strip[j].y - strip[i].y > d + EPS) break;
            d = min(d, (strip[j]-strip[i]).len());
        }
    }
    return d;
}
bool onRight(Line L1, Line L2, Line L3) {
    Vec P = Lineintersect(L2, L3);
    return ((L1.A2 - L1.A1) ^ (P - L1.A1)) < 0;
}
vector<Vec> halfPlaneIntersection(vector<Line> L) {
    sort(L.begin(), L.end());

    int head = 0, tail = 0, cnt = 0; //[head, tail)
    Line que[L.size()+1];
    for (int i = 0; i < (int)L.size() - 1; i ++) {
        if (abs(L[i].getAng() - L[i+1].getAng()) < EPS) continue;
        L[cnt ++] = L[i];
    }

    L[cnt ++] = L[(int)L.size() - 1];

    for (int i = 0; i < cnt; i ++) {
        while (tail - head > 1 && onRight(L[i], que[tail-1], que[tail-2])) -- tail;
        while (tail - head > 1 && onRight(L[i], que[head], que[head+1])) ++ head;
        que[tail++] = L[i];
    }
    while (tail - head > 1 && onRight(que[tail-1], que[head], que[head+1])) ++ head;
    while (tail - head > 1 && onRight(que[head], que[tail-1], que[tail-2])) -- tail;
    
    vector<Vec> ans;
    for (int i = head; i < tail - 1; i ++) {
        ans.push_back(Lineintersect(que[i], que[i+1]));
    }
    ans.push_back(Lineintersect(que[tail - 1], que[head]));
    
    return ans;
}
vector<Vec> Minkowski(vector<Vec> a, vector<Vec> b) {
    vector<Line> L;
    for (int i = 0; i < a.size(); i ++) {
        L.push_back(Line(a[i], a[(i+1)%(int)a.size()]));
    }
    for (int i = 0; i < b.size(); i ++) {
        L.push_back(Line(b[i], b[(i+1)%(int)b.size()]));
    }
    sort(L.begin(), L.end());
    vector<Vec> ans;

    Vec now = Vec(0, 0); 
    ans.push_back(now);
    for (int i = 0; i < L.size(); i ++) {
        ans.push_back(now + L[i].A2 - L[i].A1);
        now = now + L[i].A2 - L[i].A1;
    }

    Vec LD = Vec(INF,INF);
    for (int i = 0; i < ans.size(); i ++) {
        LD = min(LD, ans[i]);
    }
    for (int i = 0; i < ans.size(); i ++) {
        ans[i] = ans[i] + (a[0] + b[0] - LD);        
    }
    return ans;
}
/*
Circle:
*/

struct Circle {
    Vec O; double r;
    void read() {
        O.read(); scanf("%lf", &r);
    }
};

vector<Vec> isCC(Circle C1, Circle C2) {
    vector<Vec> ans;
    double d = (C1.O - C2.O).len();
    if (C1.r + C2.r > d + EPS) return ans;
    if (abs(C1.r - C2.r) < d - EPS) return ans;
    double cos1 = (C1.r * C1.r + d * d - C2.r * C2.r) / (2 * C1.r * d);
    Vec v1 = (C2.O - C1.O).change(C1.r * cos1);
    Vec v2 = v1.rotate(acos(-1)/2).change(C1.r * sqrt(1 - cos1*cos1));
    ans.push_back(C1.O + v1 + v2);
    ans.push_back(C1.O + v1 - v2);
    return ans;
} 
vector<Vec> tanCP(Circle C, Vec P) {
    double x = (C.O - P).len();
    double d = x * x - C.r * C.r;
    vector<Vec> ans;
    if (d < -EPS) return ans;
    if (d < 0) d = 0;
    Vec q1 = (P - C.O) * ( (C.r * C.r) / (x * x) );
    Vec q2 = q1.rotate(acos(-1) / 2).change(C.r / x * (sqrt(x * x - C.r * C.r)));
    ans.push_back(C.O + q1 + q2);
    ans.push_back(C.O + q1 - q2);
    return ans;
}

vector<Line> tanCC(Circle C1, Circle C2) {

}
double AreaCC(Circle C1, Circle C2) {
    double d = (C1.O - C2.O).len();
    if (C1.r + C2.r > d + EPS) {
        return C1.r * C1.r * acos(-1) + C2.r * C2.r * acos(-1);
    }
    if (abs(C1.r - C2.r) < d - EPS) {
        return max(C1.r, C2.r) * max(C1.r, C2.r) * acos(-1);
    }
    double ans = C1.r * C1.r * acos(-1) + C2.r * C2.r * acos(-1);
    double cos1 = (C1.r * C1.r + d * d - C2.r * C2.r) / (2 * C1.r * d);
    double the1 = 2 * acos(cos1);

    double cos2 = (C2.r * C2.r + d * d - C1.r * C1.r) / (2 * C2.r * d);
    double the2 = 2 * acos(cos2);

    ans = ans - ( (2 * acos(-1) - the1) * C1.r * C1.r - sin(the1) * C1.r * C1.r / 2 );
    ans = ans - ( (2 * acos(-1) - the2) * C2.r * C2.r - sin(the2) * C2.r * C2.r / 2 );
    return ans;
}   
double AreaCT(Circle C1, Vec P1, Vec P2) {

}

int main(){
    Vec P; Circle C;
    P.read(); C.read();
    vector<Vec> ans = tanCP(C, P);
    sort(ans.begin(), ans.end());
    for (auto p: ans)
        p.prt();
}