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;}