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[N][N],sum1[N][N],sum2[N][N];
LL dp[N][N];



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