uva 11139 格点计数问题

#include<stdio.h>
#include<string.h>
#include<math.h>
#include<algorithm>
#include<map>
using namespace std;
typedef long long LL;
int n;
#define N 125
int gc

,sum1

,sum2

;
LL dp

;

LL fff(int x,int y)
{   if (x>y) swap(x,y);
if (dp[x][y]!=-1) return dp[x][y];
LL s=x*y+2,ss=x*y*2+2,tmp,sss,nowt;
LL re1=0,re2=0,re3=0;

tmp=(sum2[y-1][x]-sum2[0][x])*2+(y-1)*y;
re1+=(s*(y-1)-tmp)/2;

tmp=(sum2[x-1][y]-sum2[0][y])*2+(x-1)*x;
re1+=(s*(x-1)-tmp)/2;

LL t=(x-1)*(y-1);
tmp=(sum2[x-1][y]-sum2[0][y])*(y-1)+(sum2[y-1][x]-sum2[0][x])*(x-1)+sum1[x-1][y-1]-sum1[0][y-1]-sum1[x-1][0];
sss=(x-1)*(y-1)*x*y+x*(x-1)/2*y*(y-1)/2;
re2+=(ss*t-sss-tmp)/2;

t=0;ss=0,nowt=0;tmp=0;sss=0;
for (int j=1;j<=y;j++)
{
int be=(y*x-j*x)/y+1;
nowt=x-be+1;
t+=nowt;
tmp+=gc[x][y]*nowt;
tmp+=sum2[x][y-j];
if (be>0) tmp-=sum2[be-1][y-j];
tmp+=sum2[x-be][j];
sss+=j*x*nowt+y*(x+be)*(x-be+1)/2;
}
re3+=(-x*y+2)*t;
re3-=tmp;
re3+=sss;
re3/=2;
dp[x][y]=re1*2+re2*4+re3*4;

return dp[x][y];
}
int gcd(int a,int b)
{   if (b==0) return a;
return gcd(b,a%b);
}
LL C4(LL a)
{   if (a<4) return 0;

return a*(a-1)*(a-2)/6*(a-3)/4;
}
LL C3(LL a)
{   if (a<3) return 0;
return a*(a-1)*(a-2)/6;
}

void doit()
{   LL ans1=0,ans2=0,x,tmp,q;
LL tot=(n+1)*(n+1);
for (int i=1;i<=n;i++)
for (int j=1;j<=n;j++)
{tmp=gc[i][j];
tmp--;
q=(n-i+1)*(n-j+1);
ans1+=q*tmp;
ans2+=q*tmp*(tmp-1)/2;
}

LL ans=C4(tot)-C4(n+1)*(n+1)*2-C3(n+1)*n*(n+1)*(n+1)*2;
ans-=(ans1*(tot-3)-3*ans2)*2;
for (int i=2;i<=n;i++)
for (int j=2;j<=n;j++)
{
ans+=fff(i,j)*2*(n+1-i)*(n+1-j);
}
printf("%lld\n",ans);
}

int main()
{   for (int i=0;i<N;i++)
for (int j=0;j<N;j++)
dp[i][j]=-1;
gc[0][0]=0;
for (int i=0;i<N;i++)
for (int j=i;j<N;j++) if (i||j)
gc[j][i]=gc[i][j]=gcd(i,j);

for (int j=0;j<N;j++)
{
sum2[0][j]=gc[0][j];
for (int i=1;i<N;i++)
sum2[i][j]=sum2[i-1][j]+gc[i][j];
}
for (int i=0;i<N;i++)
sum1[i][0]=sum2[i][0];
for (int j=1;j<N;j++)
for (int i=0;i<N;i++)
sum1[i][j]=sum1[i][j-1]+sum2[i][j];
while (scanf("%d",&n),n) {printf("%d ",n),doit();}
return 0;
}

LL fff(int x,int y)
{   if (x>y) swap(x,y);
if (dp[x][y]!=-1) return dp[x][y];
int s=x*y+2,ss=x*y*2+2,tmp,sss;
LL re1=0,re2=0,re3=0;
for (int i=1;i<y;i++)
{
tmp=gc[i][x]+gc[y-i][x]+y;
re1+=(s-tmp)/2;
}

for (int i=1;i<x;i++)
{
tmp=gc[i][y]+gc[x-i][y]+x;
re1+=(s-tmp)/2;
}

for (int i=1;i<x;i++)
for (int j=1;j<y;j++)
{
tmp=gc[i][y]+gc[j][x]+gc[x-i][y-j];
sss=ss-(i*y+j*x+(x-i)*(y-j));
re2+=(sss-tmp)/2;
}
s-=y*x*2;
for (int i=1;i<=x;i++)
for (int j=1;j<=y;j++)if ((y-j)*x<y*i)
{  tmp=gc[x][y]+gc[i][y-j]+gc[x-i][j];
sss=s+j*x+i*y;
re3+=(sss-tmp)/2;
}
dp[x][y]=re1*2+re2*4+re3*4;

return dp[x][y];
}