/*
* A187216.c
* Calculates the number of Q-toothpicks in the nth generation of the Q-toothpick cellular automaton, starting from two opposite Q-toothpicks
*
* Created by Nathaniel Johnston (nathaniel@nathanieljohnston.com)
* April 15, 2011
*/

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


unsigned int prevPts = 2, numPts = 2, pt[268435455];
FILE *file;

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

unsigned int getOrien (unsigned int pt)
{
    return (pt & (3 << 28)) >> 28;
}

unsigned int getPt (unsigned short x, unsigned short y, unsigned short orien)
{
    return x | (y << 14) | (orien << 28);
}

unsigned char addOrien (unsigned short x, unsigned short y, unsigned short orien)
{
    unsigned int j, nh = 0, nv = 0;
    unsigned short px, py, po;
    for(j=0;j<prevPts;j++)
    {
        px = getCoord(pt[j],1);
        py = getCoord(pt[j],2);
        po = getOrien(pt[j]);

        if(orien == 0 && ((po == 2 && px == x-1 && py == y+1) || (po == 1 && px == x-1 && py == y)))nh++;
        if(orien == 0 && ((po == 3 && px == x && py == y-1) || (po == 2 && px == x+1 && py == y-1)))nv++;

        if(orien == 1 && ((po == 3 && px == x+1 && py == y+1) || (po == 0 && px == x+1 && py == y)))nh++;
        if(orien == 1 && ((po == 3 && px == x-1 && py == y-1) || (po == 2 && px == x && py == y-1)))nv++;

        if(orien == 2 && ((po == 3 && px == x+1 && py == y) || (po == 0 && px == x+1 && py == y-1)))nh++;
        if(orien == 2 && ((po == 0 && px == x-1 && py == y+1) || (po == 1 && px == x && py == y+1)))nv++;

        if(orien == 3 && ((po == 2 && px == x-1 && py == y) || (po == 1 && px == x-1 && py == y-1)))nh++;
        if(orien == 3 && ((po == 0 && px == x && py == y+1) || (po == 1 && px == x+1 && py == y+1)))nv++;
    }
    return (nh == 1 || nv == 1);
}

unsigned char isOccupied (unsigned short x, unsigned short y)
{
    unsigned int j;
    for(j=0;j<prevPts;j++)
    {
        if(x == getCoord(pt[j],1) && y == getCoord(pt[j],2))return 1;
    }
    return 0;
}

int main ()
{
    unsigned short xi, yi, oi;
    unsigned short gen, maxGens;
    pt[0] = getPt(8192,8192,0);
    pt[1] = getPt(8191,8191,2);

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

    file = fopen("b187216.txt","w");
    fprintf(file,"0 0\n1 2\n");
    fclose(file);
    printf("0 0\n1 2\n");

    for(gen=1;gen<maxGens;gen++){
        for(xi=8191-gen;xi<=8192+gen;xi++){
            for(yi=8191-gen;yi<=8192+gen;yi++){
                for(oi=0;oi<4;oi++){
                    if(!isOccupied(xi,yi) && addOrien(xi,yi,oi))pt[numPts++]=getPt(xi,yi,oi);
                }
            }
        }
        file = fopen("b187216.txt","a");
        fprintf(file,"%d %d\n",gen+1,numPts);
        fclose(file);
        printf("%d %d\n",gen+1,numPts);
        prevPts = numPts;
    }

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