ZOJ3825(2014年牡丹江现场赛G题)

比赛时读这题一开始感觉并不难，思路好像蛮容易想的，但是现场没有一支队过了这题……最后一小时上机打了这题，在终场前十秒提交，可惜结果自然是WA……

比完赛这几天好好想了下这题，WA了好多次，期间还尝试用JAVA来写，结果都没有成功……

最后特么终于用C++把这题AC了……虽然是用Mathematica打了个对拍查了很久错误，不算是完全自己做出来的……

这题题意大致可以转化为：给定一个圆x^2+y^2<=R^2和一条直线A x+B y+C=0，其中A,B,R是数量级10^8的整数，C是数量级10^16的整数，求这条直线在这个圆内及边界上经过的整点个数。这里我们可以假设gcd(A,B)=1，因为若gcd(A,B)不能整除C，这个方程无解；若能整除，A,B,C同时除以gcd(A,B)后可使x,y的系数互质。

思路其实很直观。首先可以用拓展欧几里得算法求出A x+B y+C=0的一个特解(-CX,-CY)，这里的X,Y满足AX+BY=1（当然因为精度问题程序中我们求出(X,Y)后不能直接求出特解，这里(-CX,-CY)只作为记号）。显然求出特解后，A x+B y+C=0上所有整点都必有形式(x,y)=k(B,-A)+(-CX,-CY)，k是整数。那么，就应该有(kB-CX)^2+(-kA-CY)^2=0，化简即有k要在[(C(BX-AY)-S)/(A^2+B^2),(C(BX-AY)+S)/(A^2+B^2)]:=[T1,T2]中，这里S=sqrt(R^2(A^2+B^2)-C^2)，接下来只要求这个闭区间里整数点的个数即可。

但是，因为S^2已经到了10^32的数量级，C(BX-AY)的数量级也达到了10^24，正常方法直接用double 甚至是long double精度都是不够的，而因为这个问题我也WA了很多次。要做出这题，后来我发现必须完全避开浮点数的运算，才能保证精度。注意到，首先如果将整个区间整体左移或者右移整数个单位，整数点的个数是不变的。这样，首先我们可以给区间端点把分子中的C(BX-AY)用C(BX-AY)%(A^2+B^2)替代，而这个取余虽然被除数超出long long精度但仍是可以用程序实现的。其次，注意到S用floor(S)替代对结果是没有影响的，因为[T1,T2]中整数点个数等于floor(T2)-ceil(T1)+1。其中若k=floor(T2),一定有k(A^2+B^2)<=S+...，由左式是整数可知右式向下取整不影响不等式的正确性。对于ceil(T1)的讨论类似。

那么现在，我们需要做的事情就是，写一个高精度整数类表示至少达到10^32整数，并且对于一个整数可以快速求出不超过它平方根的最大整数（这个整数可直接用long long表示）。最后，因为T1，T2都是有理数，那么求floor和ceil时我们不需要转化为浮点数运算，直接相除即可，当然要注意T1，T2的符号这些细节。具体实现可见代码。

最后感叹一句啊，为什么我比赛要开这么道神题啊……

#include<stdio.h>
#include<math.h>
#include<string.h>
#include<algorithm>
#include<iostream>
using namespace std;

typedef long long ll;
typedef long double lf;
const ll mod=10000000;
const int maxn=6;

struct bigint
{
	ll a[maxn];
	bigint(ll tet=0)
	{
		memset(a,0,sizeof(a));
		for(int i=0;i<maxn;i++)
		{
			a[i]=tet%mod;
			tet/=mod;
		}
	}
};
bool operator <(bigint a,bigint b)
{
	int i,j;
	for(i=maxn-1;i>=0;i--)
		if(a.a[i]!=b.a[i])
			return a.a[i]<b.a[i];
	return false;
}
bool operator ==(bigint a,bigint b)
{
	return !(a<b)&&!(b<a);
}
bigint operator +(bigint a,bigint b)
{
	bigint res;
	for(int i=0;i<maxn;i++) res.a[i]=a.a[i]+b.a[i];
	for(int i=0;i<maxn-1;i++) res.a[i+1]+=(res.a[i]/mod),res.a[i]%=mod;
	return res;
}
bigint operator -(bigint a,bigint b)
{
	bigint res;
	if(a==b||a<b) return res;
	for(int i=0;i<maxn;i++) res.a[i]=a.a[i]-b.a[i];
	for(int i=0;i<maxn-1;i++)
	{
		if(res.a[i]<0) {res.a[i]+=mod;res.a[i+1]--;}
	}
	return res;
}
bigint operator *(bigint a,bigint b)
{
	bigint res;
	int i,j;
	for(i=0;i<maxn/2;i++) for(j=0;j<maxn/2;j++) res.a[i+j]+=(a.a[i]*b.a[j]);
	for(i=0;i<maxn-1;i++) res.a[i+1]+=(res.a[i]/mod),res.a[i]%=mod;
	return res;
}
bigint c2(bigint a)
{
	for(int i=maxn-1;i>=0;i--) {if(i&&a.a[i]%2) {a.a[i-1]+=mod;}a.a[i]/=2;}
	return a;
}
ll zh(bigint a)
{
	ll te=0;
	for(int i=maxn-1;i>=0;i--) te=te*mod+a.a[i];
	return te;
}
void print(bigint a)
{
	int i,j;
	for(i=maxn-1;i>=0;i--) if(a.a[i]||!i) {printf("%lld",a.a[i]);break;}
	for(--i;i>=0;i--) printf("%07lld",a.a[i]);
	printf("\n");
}

void Sqrt(bigint a,ll &r1,ll &r2)
{
	if(a==bigint(0)) {r1=r2=0;return;}
	bigint at(0),bt(100000000),tt,ss;
	bt=bt*bt*bigint(100);
	while(bigint(1)<bt-at)
	{
		tt=c2(at+bt);ss=tt*tt;
		if(ss==a) {r1=zh(tt);r2=r1;return;}
		if(ss<a) at=tt;
		else bt=tt;
	}
	r1=zh(at);r2=r1+1;
}

void gcd(ll a,ll b,ll &d,ll &x,ll &y)
{
	if(!b) {d=a;x=1;y=0;}
	else {gcd(b,a%b,d,y,x);y-=(a/b)*x;}
}

ll tabs(ll a) {return a>0?a:-a;}

ll mo(ll x, ll y, ll mo){
    ll sgn=(x<0?-1:(x>0))*(y<0?-1:(y>0));
    x=tabs(x);y=tabs(y);
    if(!sgn) return 0;
    ll n=mo;
    ll T=floorl(sqrtl(n)+0.5);
    ll t=T*T-n;
    ll a=x/T, b=x%T;
    ll c=y/T, d=y%T;
    ll e=a*c/T, f=a*c%T;
    ll v=((a*d+b*c)%n+e*t)%n;
    ll g=v/T, h=v%T;
    ll ret=(((f+g)*t%n + b*d)%n+h*T)%n;
    ret=(ret%n+n)%n;
    return ret*sgn;
}

ll xt[5][2],r,s;
ll x,y;

ll sqz(ll a,ll b)
{
	if(a>0) return (a+b-1)/b;
	if(a==0) return 0;
	return -((-a)/b);
}
ll xqz(ll a,ll b)
{
	if(a>0) return a/b;
	if(a==0) return 0;
	return -(-a+b-1)/b;
}

ll cal(ll A,ll B,ll C)
{
	//cout<<C<<" "<<-(A*y-B*x)<<endl;
	ll t1=A*A+B*B;
	ll fm=mo(C,-A*y+B*x,t1);
	C=tabs(C);
	ll tt1,tt2;
	if(bigint(t1)*bigint(r*r)<bigint(C)*bigint(C)) return 0;
	//print(bigint(t1)*bigint(r*r)-bigint(C)*bigint(C));
	Sqrt(bigint(t1)*bigint(r*r)-bigint(C)*bigint(C),tt1,tt2);
	ll kst=sqz(fm-tt1,t1),ken=xqz(fm+tt1,t1);
	//cout<<fm<<" "<<t1<<" "<<kst<<" "<<ken<<endl;
	return kst>ken?0:ken-kst+1;
}



int main()
{
	ll T;
	scanf("%lld",&T); 
	while(T--)
	{
		int i,j;
		scanf("%lld",&s);
		for(i=0;i<3;i++)
		{
			scanf("%lld %lld",&xt[i][0],&xt[i][1]);
			if(!i) scanf("%lld",&r);
		}
		for(i=1;i<3;i++) for(j=0;j<2;j++) xt[i][j]=xt[i][j]-xt[0][j];
		ll A=xt[2][1]-xt[1][1],B=xt[1][0]-xt[2][0],C=ll(xt[2][0])*xt[1][1]-ll(xt[1][0])*xt[2][1];
		if(A==0&&B==0){printf("0\n");continue;}
		ll d;
		gcd(tabs(A),tabs(B),d,x,y);
		if(s%d) {printf("0\n");continue;}
		if(A<0) x=-x;
		if(B<0) y=-y;
		A/=d;B/=d;C/=d;s/=d;
		ll res=cal(A,B,C+s)+cal(A,B,C-s);
		printf("%lld\n",res);
	}
	return 0;
}