Language: C++

#include <bits/stdc++.h>

#define _overload(_1,_2,_3,name,...) name
#define _rep(i,n) _range(i,0,n)
#define _range(i,a,b) for(int i=int(a);i<int(b);++i)
#define rep(...) _overload(__VA_ARGS__,_range,_rep,)(__VA_ARGS__)

#define _rrep(i,n) _rrange(i,n,0)
#define _rrange(i,a,b) for(int i=int(a)-1;i>=int(b);--i)
#define rrep(...) _overload(__VA_ARGS__,_rrange,_rrep,)(__VA_ARGS__)

#define _all(arg) begin(arg),end(arg)
#define uniq(arg) sort(_all(arg)),(arg).erase(unique(_all(arg)),end(arg))
#define getidx(ary,key) lower_bound(_all(ary),key)-begin(ary)
#define clr(a,b) memset((a),(b),sizeof(a))
#define bit(n) (1LL<<(n))
#define popcount(n) (__builtin_popcountll(n))

using namespace std;

template<class T>bool chmax(T &a, const T &b) { return (a<b)?(a=b,1):0;}
template<class T>bool chmin(T &a, const T &b) { return (b<a)?(a=b,1):0;}

using ll=long long;
using R=long double;
const R EPS=1e-9L; // [-1000,1000]->EPS=1e-8 [-10000,10000]->EPS=1e-7
inline int sgn(const R& r){return(r > EPS)-(r < -EPS);}
inline R sq(R x){return sqrt(max(x,0.0L));}

const int dx[8]={1,0,-1,0,1,-1,-1,1};
const int dy[8]={0,1,0,-1,1,1,-1,-1};

// Description: ????????????
// Verifyed: various problem 
const R INF = 1E40L;
const R PI = acos(-1.0L);
using P=complex<R>;
const P O=0.0L;
using L=struct{P s,t;};
using VP=vector<P>;
using C=struct{P c;R r;};

auto& operator >> (istream& is,P& p){ R x,y;is >> x >> y; p=P(x,y); return is;}
auto& operator << (ostream& os,P& p){ os << real(p) << " " << imag(p); return os;}

inline R dot(P o,P a,P b){return real(conj(a-o)*(b-o));}
inline R det(P o,P a,P b){return imag(conj(a-o)*(b-o));}
inline int sdot(P o,P a,P b){return sgn(dot(o,a,b));}
inline int sdet(P o,P a,P b){return sgn(det(o,a,b));}

//projection verify AOJ CGL_1_A
P proj(L l,P p){ R u=real((p-l.s)/(l.t-l.s)); return (1-u)*l.s+u*l.t;} 
 
namespace std{
    bool operator <  (const P& a,const P& b){ return sgn(real(a-b))?real(a-b)<0:sgn(imag(a-b))<0;}
    bool operator == (const P& a,const P& b){ return sgn(real(a-b))==0 && sgn(imag(a-b))==0;}
}

bool operator <  (const C& a,const C& b){ return a.c==b.c?sgn(a.r-b.r)<0:a.c<b.c;}
bool operator == (const C& a,const C& b){ return a.c==b.c&&sgn(a.r-b.r)==0;}

L tan(C c,P p){P d=1il*(p-c.c);return L{p-d,p+d};}
VP contact(C c,P p){
    VP ret;
    P d=p-c.c;
    for(P s:{-1il,1il}) ret.emplace_back(c.c+c.r*(c.r+sq(norm(d)-norm(c.r))*s)*d/norm(d));
    return uniq(ret);
}


int main(void){
    P p;
    C c;
    cin >> p >> c.c >> c.r;
    auto res=contact(c,p);
    cout.precision(20);
    for(auto &it:res) cout << fixed << real(it) << " " << imag(it) << endl;
    return 0;
}

Compiler output:
a.cc: In function 'VP contact(C, P)':
a.cc:13:34: error: could not convert '(std::sort<__gnu_cxx::__normal_iterator<complex<long double>*, vector<complex<long double> > > >(std::begin<vector<complex<long double> > >(ret), std::end<vector<complex<long double> > >(ret)), ret.std::vector<std::complex<long double> >::erase(__gnu_cxx::__normal_iterator<const std::complex<long double>*, std::vector<std::complex<long double> > >(std::unique<__gnu_cxx::__normal_iterator<complex<long double>*, vector<complex<long double> > > >(std::begin<vector<complex<long double> > >(ret), std::end<vector<complex<long double> > >(ret))), __gnu_cxx::__normal_iterator<const std::complex<long double>*, std::vector<std::complex<long double> > >(std::end<vector<complex<long double> > >(ret))))' from 'std::vector<std::complex<long double> >::iterator' to 'VP' {aka 'std::vector<std::complex<long double> >'}
   13 | #define uniq(arg) sort(_all(arg)),(arg).erase(unique(_all(arg)),end(arg))
      |                   ~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      |                                  |
      |                                  std::vector<std::complex<long double> >::iterator
a.cc:67:12: note: in expansion of macro 'uniq'
   67 |     return uniq(ret);
      |            ^~~~