有史以来最奇葩的事情发生了！！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 doubleusing 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-9struct 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_ptpn=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;}