#include <math.h>
#include <stdio.h>

/*
    Created by Nathaniel Johnston (nathaniel@nathanieljohnston.com), Mar 25 2011
*/

unsigned int prevPts = 1, numPts = 1, numXYaxis = 1, pt[268435455];

unsigned int getCoord (unsigned int pt, unsigned short coord)
{
    unsigned short shift = 10*(coord-1);
    return (pt & (1023*(1 << shift))) >> shift;
}

unsigned int getPt (unsigned short x, unsigned short y, unsigned short z)
{
    return x | (y << 10) | (z << 20);
}

unsigned short numNeigh (unsigned short x, unsigned short y, unsigned short z)
{
    unsigned int j, ndist, nn = 0;
    short px, py, pz;
    for(j=0;j<prevPts;j++)
    {
        px = (getCoord(pt[j],1)-x);
        py = (getCoord(pt[j],2)-y);
        pz = (getCoord(pt[j],3)-z);
        ndist = px*px + py*py + pz*pz;
        if(ndist == 2)nn++;
        else if(ndist == 0)return 0;
    }
    return nn;
}

int main ()
{
    unsigned short xi, yi, zi;
    unsigned char gen, maxGens;
    pt[0] = getPt(512,512,512);

    printf("This tool will calculate terms of Sloane's A147552.\nPlease enter the number of terms to compute (an integer from 1 to 255): ");
    scanf("%d",&maxGens);

    printf("0: 0\n1: 1\n");
    for(gen=1;gen<maxGens;gen++){
        for(xi=512-gen;xi<=512+gen;xi++){
            for(yi=512-gen;yi<=512+gen;yi++){
                for(zi=512-gen + ((xi+yi+gen) % 2);zi<=512+gen;zi+=2){
                    if(numNeigh(xi,yi,zi) == 1){
                        pt[numPts++]=getPt(xi,yi,zi);
                        if(zi==512)numXYaxis++;
                    }
                }
            }
        }
        prevPts = numPts;
        printf("%d: %d\n",gen+1,numXYaxis);
    }

    printf("Done!");
    getchar();
    return 0;
}