計算幾何，狗都不寫

題意簡述：在平面上有一個座標 \((x_c,y_c)\) 和半徑 \(r\) 都是整數的圓 \((1\leq r_c\leq \sqrt{x_c^2+y_c^2}-1)\)，你可以詢問不超過 \(60\) 次一個從原點出發的射線 \((0,0)\to (x,y)\) 與圓的距離（若相交，則返回 \(0\)）。確定圓的座標和半徑。

\(|x_c|,|y_c|,r\le 10^5\)，你需要保證 \(|x|,|y|\leq 10^6\)。

一道有趣的計算幾何，思維含量並不是很高（*3300 有點離譜），就是寫起來麻煩了些。比賽的時候花了四個小時寫這題，結果發現前面兩個半小時寫的東西完全沒有用（狂笑）。

---

Step 1：旋轉座標軸規避分類討論

首先我們得確定這個圓的大致位置，也就是在哪個象限。因為如果對於四個象限分別分類討論很麻煩，不如旋轉座標軸使得這個圓落在某些固定的象限，每次旋轉 \(90^{\circ}\)。不妨將圓心旋轉到 \(x\) 軸上方，且確保這個圓不經過第三象限，這是可以同時做到的：

我們可以通過詢問 \((0,\pm i)\) 與 \((\pm i,0)\) 來確定圓和座標軸是否有交。

• 若沒有交點，將座標軸旋轉使得圓落在第一象限。
• 若有一個交點，將座標軸旋轉使得圓落在第一、二象限。
• 若有兩個交點，將座標軸旋轉使得圓落在第一、二、四象限。

不要忘了旋轉後詢問要轉回來。比如我們將座標軸順時針旋轉 \(90^{\circ}\)

後，詢問時需要將座標逆時針旋轉 \(90^{\circ}\)。

---

Step 2：二分找到與圓相切的射線

由於圓心在 \(x\) 軸上方，詢問互動庫 \((0,x)\ (x<0)\) 的返回值 \(d_0\) 顯然是 \(\sqrt{x_c^2+y_c^2}-r\)：因為過圓心且垂直於射線的直線與射線沒有交點，故最短距離為原點與圓心距離減去半徑。

注意到根據三個方程我們可以直接解出 \((x,y,r)\)，\(d_0=\sqrt{x_c^2+y_c^2}-r\) 已經是其中一個了。只需要再找到兩個和圓距離 \(\neq d_0\) 且 \(\neq 0\) 的射線即可，稱其為合法射線。

但是我們不能隨便詢問！因為可以構造一個與 \(y=kx\)

非常相近且半徑非常大的圓，例如與 \(y=-x\) 幾乎相切的 \((50000,50000,70709)\)，使得合法的射線範圍非常小。因為一旦稍有偏移，我們詢問到的是 \(d_0\)，要麼就是 \(0\)（上述例子中，只有 \(y=-x\) 和過原點的圓的兩條切線之間所夾住的區域合法）。而我們不知道出題人會構造什麼毒瘤資料，這也意味著我們不能詢問固定的射線，這樣是沒有前途的。

但注意到距離切線非常近（指與切線夾角非常小且與圓相離的射線範圍）的那一塊區域是合法的，這啟發我們二分出幾乎與圓相切（但與圓相離）的兩條射線 \(a:(0,0)\to (x_a,y_a)\) 與 \(b:(0,0)\to (x_b,y_b)\)。具體怎麼二分呢，想象一個 \((2\times10^6)\times (2\times 10^6)\) 的盒子包住了整個座標軸。我們在從左逼近時，若返回值 \(>\rm eps\) 則認為與圓相離，向順時針方向二分，否則認為與圓相交或相切，向逆時針方向二分。從右逼近則相反。

程式碼中尋找射線 \(a\) 只在第二象限內（原因在接下來會提到），因為圓心在 \(x\) 軸上方導致就算二分不到與切線很相近的位置，我們也可以找到不同於 \(d_0\) 和 \(0\) 的其它合法射線（想一想，為什麼？Hint：\(10^6>10^5\)）。

Step 3：列舉 \(r\) 並根據 \(d_0,d_a\) 求出 \((x,y)\)

根據圓與 \(a\) 的距離 \(d_a\) 以及 \(d_0\) 仍不足以確定這個圓：兩個方程，三個變數。如果再加上 \(d_b\) 那麼算起來太麻煩，精度爆炸怎麼辦？沒有關係，注意到 \(r\) 是整數並且 \(\leq 10^5\)，我們直接列舉 \(r\) 並根據 \(d_0,d_a\) 算出 \(x,y\)，再檢查當前的圓 \((x,y,r)\) 與 \(b\) 的距離是否等於 \(d_b\)。若相等則說明 \((x,y,r)\) 即為答案。這樣精度就很優秀了。

怎麼算 \(x,y\)：注意到射線 \(a\) 在第二象限內，設其與 \(x\) 軸負半軸的夾角為 \(\alpha\)，先將座標軸逆時針旋轉 \(\alpha\) 使得 \(a\) 落在 \(x\) 軸負半軸上，我們設原來的 \(x,y\) 經過這樣的變換後新的座標為 \(x_T,y_T\)。我們很容易求出 \(x_T,y_T\)：\(y_T\) 就是 \(r+d_a\)，而 \(d_0\) 在這樣的變換下是不會變的（旋轉變換不改變兩點之間的距離），因此 \(x_T^2+y_T^2=(r+d_0)^2\)。解得 \(|x_T|\) 後還需要確定符號：究竟是在新座標系的第一象限還是第二象限。這個可以通過求得 \(a\) 後詢問互動庫 \((-x_a,-y_a)\) 得到 \(d_a'\)，若 \(d_a<d_a'\) 說明 \(T\) 在第二象限，令 \(x_T\gets -|x_T|\)，否則 \(T\) 在第一象限，\(x_T\gets |x_T|\)。

注意到我們將座標軸逆時針旋轉了 \(\alpha\) 角度，因此還要順時針轉回來。將其乘以

\[\begin{bmatrix}\cos (-\alpha)&-\sin(-\alpha)\\\sin (-\alpha)&\cos(-\alpha)\end{bmatrix} \]

即

\[\begin{bmatrix}\cos \alpha&\sin\alpha\\-\sin \alpha&\cos\alpha\end{bmatrix} \]

即可。根據 \(\alpha\in\left[0,\dfrac \pi 2\right]\) 進行了化簡。至於怎麼算圓到射線的距離，點積搞搞就好了吧，實在不行可以像我一樣上網現學。

至此，本題被我們完美解決。時間複雜度 \(\mathcal{O}(r)\)，詢問次數為 \(2\log V\)，其中 \(V\) 是詢問最大範圍。程式碼中使用 \(V=2\times 10^5\)，一個測試點大約詢問 \(50\) 次。

const ld eps=1e-5;
const ll N=2e5;

ld max(ld x,ld y){return x>y?x:y;}
bool Eq(ld x,ld y){return fabs(x-y)<eps;}
bool Isz(ld x){return fabs(x)<eps;}

ld up,down,left,right;
ld dis(ld x1,ld y1,ld x2,ld y2){
	return sqrt((x1-x2)*(x1-x2)+((y1-y2)*(y1-y2)));
}
ld dis(ld x1,ld y1,ld x2,ld y2,ld r){
	ld dot=x1*x2+y1*y2,d=dis(0,0,x2,y2);
	if(dot<eps)return max((ld)0,d-r);
	ld need=dot/dis(0,0,x1,y1);
	return max((ld)0,sqrt(d*d-need*need)-r);
}

int swp,type;
void ref(int &x,int &y){int t=y; y=x,x=-t;}
ld query(int x,int y){
	for(int j=1;j<=swp;j++)ref(x,y);
	std::cout<<"? "<<x<<" "<<y<<std::endl;
	ld res; std::cin>>res;
	return res;
}
void answ(int x,int y,int r){
	for(int j=1;j<=swp;j++)ref(x,y);
	std::cout<<"! "<<x<<" "<<y<<" "<<r<<std::endl,exit(0);
}
void answer(ld x,ld y,ld r){
	int xx=x<0?x-eps:x+eps;
	int yy=y<0?y-eps:y+eps;
	return answ(xx,yy,r+eps);
}

int main(){
	up=query(0,10),down=query(0,-10);
	left=query(-10,0),right=query(10,0);
	type=Isz(up)+Isz(down)+Isz(left)+Isz(right);
	if(type==0){
		ld mx=max(max(up,down),max(left,right));
		if(down==mx&&right==mx)swp=1;
		if(right==mx&&up==mx)swp=2;
		if(up==mx&&left==mx)swp=3;
	}
	if(type==1){
		if(Isz(left))swp=1;
		if(Isz(down))swp=2;
		if(Isz(right))swp=3;
	}
	if(type==2){
		if(Isz(left)&&Isz(up))swp=1;
		if(Isz(down)&&Isz(left))swp=2;
		if(Isz(right)&&Isz(down))swp=3;
	}
	up=query(0,10),down=query(0,-10);
	left=query(-10,0),right=query(10,0);
	ll l=1,r=N*2-1,xx,yy; ld res;
	while(l<r){
		ll m=l+r>>1;
		xx=-std::min(N,m),yy=N-max(0ll,m-N);
		(res=query(xx,yy))>1e-6?r=m:l=m+1;
	}
	xx=-std::min(N,l),yy=N-max(0ll,l-N);
	ld upl=query(xx,yy),downr=query(-xx,-yy);
	ld mdis=std::min(upl,downr);
	
	ll xx2,yy2; l=1,r=N*3-1;
	while(l<r){
		ll m=(l+r>>1)+1;
		xx2=std::min(N,m),yy2=-N+max(0ll,m-N);
		(res=query(xx2,yy2)>1e-6)?l=m:r=m-1;
	}
	xx2=std::min(N,l),yy2=-N+max(0ll,l-N);
	res=query(xx2,yy2);
	for(int r=1;r<=1e5;r++){
		ld y=mdis+r,x=sqrt((down+r)*(down+r)-y*y);
		ld sn=(ld)fabs(yy)/(ld)(sqrt(xx*xx+yy*yy));
		ld cs=(ld)fabs(xx)/(ld)(sqrt(xx*xx+yy*yy));
		if(Eq(mdis,upl))x=-x;
		ld x0=x*cs+y*sn,y0=-x*sn+y*cs; x=x0,y=y0;
		if(Eq(res,dis(xx2,yy2,x,y,r)))answer(x,y,r);
	}
	return 0;
}
// qwq Kelly qwq