有史以来最奇葩的事情发生了！！lightOJ 1130 Intersection between Circle and Rectangle

http://www.lightoj.com/volume_showproblem.php?problem=1130

题目是求一个矩形和一个圆的相交面积，这不是重点，接下来把非重点都略过。

有一个地方要用到sqrt函数，由于精度问题改了好久才过，于是就想起有个QUAKE III里用的高速开根号函数，不知道会不会比sqrt更好用，就拿来用了，代码如下：

float Q_rsqrt( float number ){
   int i;
   float x2, y;
   const float threehalfs = 1.5F;

   x2 = number * 0.5F;
   y   = number;
   i   = * ( int * ) &y;   // evil doubleing point bit level hacking
   i   = 0x5f3759df - ( i >> 1 ); // what the fuck?
   y   = * ( float * ) &i;
   y   = y * ( threehalfs - ( x2 * y * y ) ); // 1st iteration
   y   = y * ( threehalfs - ( x2 * y * y ) ); // 2nd iteration, this can be removed

   return y;
}

一换上果然也能A，于是奇葩的事情发生了！我发现这个函数不是用来求sqrt()的！！！！而是求1/sqrt()的！！！这也能A啊？！！！我把它取了个倒数再用！！尼玛！！居然还就WA了！！！！

更更奇葩的是！！为了精度需要我将它拿过来时把float都改成了double，而其实这个函数只能用于float！！！如果改成double它只会输出0啊！！！！！！！只会输出0他也能A啊！！！！！！！

把那句关键的试了几种情况如下：

ds=cj/sqrt(da*db);AC

ds=cj/sqrt(da)/sqrt(db);WA

ds=cj/Q_rsqrt(da*db);AC

ds=cj*Q_rsqrt(da*db);WA

估计应该有精度问题啊，数据也不会弱啊，不然1、2就不会有区别了。

3、4情况就完全不理解了……

#include <set>
#include <map>
#include <list>
#include <cmath>
#include <ctime>
#include <deque>
#include <queue>
#include <stack>
#include <cctype>
#include <cstdio>
#include <string>
#include <vector>
#include <cassert>
#include <cstdlib>
#include <cstring>
#include <sstream>
#include <iostream>
#include <algorithm>
#define ll __int64
//#define ll long long
#define PI acos(-1.0)
#define PRN 105
#define INF 2000000000
#define typ double
using namespace std;
int f_min(int x,int y){return x<y?x:y;}
//int f_max(int x,int y){return x>y?x:y;}
typ f_max(typ x,typ y){return x>y?x:y;}
int f_abs(int x){return x<0?-x:x;}
typ f_abs(typ x){return x<0?-x:x;}
int prime[PRN],Pn;
bool unprime[PRN];
void get_prime(){
	int i,j,lim=sqrt(PRN*1.0);
	memset(unprime,0,sizeof(unprime));Pn=0;
	for(i=2;i<PRN;i++){
		if(unprime[i])continue;
		prime[Pn++]=i;
		if(i>lim)continue;
		for(j=i*i;j<PRN;j+=i)unprime[j]=1;
	}
}
struct Edge{
	int u,v,w,next;
};
Edge edge[210000];
int head[100005],en;
void insert(int u,int v,int w){
	edge[en].u=u;edge[en].v=v;edge[en].w=w;
	edge[en].next=head[u];head[u]=en++;
}
#define eps 1e-9
struct Point{
	double x,y;
	void disp(){
		printf("x=%lf y=%lf\n",x,y);
	}
	Point friend mi(Point a,Point b){
		Point res;
		res.x=a.x-b.x;
		res.y=a.y-b.y;
		return res;
	}
};
struct Circle{
	Point c;
	double r;
};
struct Line{
	Point s,e;
	void set(double sx,double sy,double ex,double ey){
		s.x=sx;s.y=sy;e.x=ex;e.y=ey;
	}
};

Circle cir;
Line line[4];
Point pt[20];
int pn;
void put_pt(Line l){
	Point temp;
	double t,d;
	if(f_abs(l.s.x-l.e.x)<eps){
		if(f_abs(cir.c.x-l.s.x)>cir.r)return;
		d=cir.c.x-l.s.x;
		t=sqrt(f_abs(cir.r*cir.r-d*d));
		temp.x=l.s.x;
		if(l.s.y<l.e.y){
			temp.y=cir.c.y-t;
			if(temp.y>l.s.y&&temp.y<l.e.y)pt[pn++]=temp;
			temp.y=cir.c.y+t;
			if(temp.y>l.s.y&&temp.y<l.e.y)pt[pn++]=temp;
		}else{
			temp.y=cir.c.y+t;
			if(temp.y<l.s.y&&temp.y>l.e.y)pt[pn++]=temp;
			temp.y=cir.c.y-t;
			if(temp.y<l.s.y&&temp.y>l.e.y)pt[pn++]=temp;
		}
	}else{
		if(f_abs(cir.c.y-l.s.y)>cir.r)return;
		d=cir.c.y-l.s.y;
		t=sqrt(f_abs(cir.r*cir.r-d*d));
		temp.y=l.s.y;
		if(l.s.x<l.e.x){
			temp.x=cir.c.x-t;
			if(temp.x>l.s.x&&temp.x<l.e.x)pt[pn++]=temp;
			temp.x=cir.c.x+t;
			if(temp.x>l.s.x&&temp.x<l.e.x)pt[pn++]=temp;
		}else{
			temp.x=cir.c.x+t;
			if(temp.x<l.s.x&&temp.x>l.e.x)pt[pn++]=temp;
			temp.x=cir.c.x-t;
			if(temp.x<l.s.x&&temp.x>l.e.x)pt[pn++]=temp;
		}
	}
}
void get_data(){
	scanf("%lf%lf%lf",&cir.c.x,&cir.c.y,&cir.r);
	double lx,ly,rx,ry;
	scanf("%lf%lf%lf%lf",&lx,&ly,&rx,&ry);
	if(f_abs(lx-rx)<eps||f_abs(ly-ry)<eps)while(1);
	line[0].set(lx,ry,rx,ry);
	line[1].set(rx,ry,rx,ly);
	line[2].set(rx,ly,lx,ly);
	line[3].set(lx,ly,lx,ry);
	
	//get_pt
	pn=0;
	pt[pn++]=line[0].s;
	put_pt(line[0]);
	pt[pn++]=line[1].s;
	put_pt(line[1]);
	pt[pn++]=line[2].s;
	put_pt(line[2]);
	pt[pn++]=line[3].s;
	put_pt(line[3]);
	//	for(int i=0;i<pn;i++)pt[i].disp();
}
bool out[20];
double f_sqrt(double x){
	double h=f_max(2,x),l=0,mid;
	while(h>l+(1e-12)){
		mid=(h+l)/2;
		if(mid*mid<x)l=mid;
		else h=mid;
	}
	return h;
}
double Q_rsqrt( double number ){
   int i;
   double x2, y;
   const double threehalfs = 1.5F;

   x2 = number * 0.5F;
   y   = number;
   i   = * ( int * ) &y;   // evil floating point bit level hacking
   i   = 0x5f3759df - ( i >> 1 ); // what the fuck?
   y   = * ( double * ) &i;
   y   = y * ( threehalfs - ( x2 * y * y ) ); // 1st iteration
   y   = y * ( threehalfs - ( x2 * y * y ) ); // 2nd iteration, this can be removed

   return y;
} 
double get_dis(Point a,Point b){
	return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
double get_dis2(Point a,Point b){
	return (a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y);
}
void get_out(){
	int i;
	double t;
	for(i=0;i<pn;i++){
		t=get_dis(pt[i],cir.c);
		if(t>cir.r+eps)out[i]=1;
		else out[i]=0;
	}
}

double get_cj(Point a,Point b){
	return a.x*b.y-a.y*b.x;
}
double get_dj(Point a,Point b){
	return a.x*b.x+a.y*b.y;
}
double get_tri(Point a,Point b){
	return get_cj(mi(a,cir.c),mi(b,cir.c))/2;
}
double get_shan(Point a,Point b){
	double deg,da=get_dis2(a,cir.c),db=get_dis2(b,cir.c);
	double cj=get_cj(mi(a,cir.c),mi(b,cir.c)),dj=get_dj(mi(a,cir.c),mi(b,cir.c)),dc,ds;
	ds=cj/Q_rsqrt(da*db);
	dc=dj;
	deg=asin(ds);
	if(dc<-eps){
		if(deg>0)deg=PI-deg;
		else deg=-PI-deg;
	}
	return deg/2*cir.r*cir.r;
}
void run(){
	if(f_abs(cir.r)<eps){
		printf("0.00000000\n");
		return;
	}
	get_out();
	int i,tn=0;
	double res=0,deg=0,td[20];
	for(i=0;i<pn;i++){
		if(!out[i]&&!out[(i+1)%pn])res+=get_tri(pt[i],pt[(i+1)%pn]);
		else res+=get_shan(pt[i],pt[(i+1)%pn]);
	}
	printf("%.12lf\n",f_abs(res));
}
int main(){
	int i,t;
	scanf("%d",&t);
	for(i=1;i<=t;i++){
		get_data();
		printf("Case %d: ",i);
		run();
	}
	return 0;
}