/*
* countshapes.c
* Calculates the number of circles, diamond, and hearts present in the nth
* generation of the Q-toothpick cellular automaton described by Sloane's A187210
*
* Created by Nathaniel Johnston (nathaniel@nathanieljohnston.com)
* http://www.nathanieljohnston.com/2011/03/the-q-toothpick-cellular-automaton/
* March 27, 2011
*/

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


unsigned int prevPts = 1, numPts = 1, 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 char 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 getOrien(pt[j]);
    }
    return 4;
}

unsigned short countCircles(void){
    unsigned int j;
    unsigned short px, py, po, numCirc=0;
    for(j=0;j<prevPts;j++)
    {
        px = getCoord(pt[j],1);
        py = getCoord(pt[j],2);
        po = getOrien(pt[j]);

        if(po == 2 && isOccupied(px+1,py)==3 && isOccupied(px,py+1)==1 && isOccupied(px+1,py+1)==0)numCirc++;
    }
    return numCirc;
}

unsigned short countDiamonds(void){
    unsigned int j;
    unsigned short px, py, po, numDiam=0;
    for(j=0;j<prevPts;j++)
    {
        px = getCoord(pt[j],1);
        py = getCoord(pt[j],2);
        po = getOrien(pt[j]);

        if(po == 0 && isOccupied(px+1,py)==1 && isOccupied(px,py+1)==3 && isOccupied(px+1,py+1)==2)numDiam++;
    }
    return numDiam;
}

unsigned short countHearts(void){
    unsigned int j;
    unsigned short px, py, po, numHeart=0;
    for(j=0;j<prevPts;j++)
    {
        px = getCoord(pt[j],1);
        py = getCoord(pt[j],2);
        po = getOrien(pt[j]);

        if(po == 0 && isOccupied(px,py+1)==4 && isOccupied(px+1,py+1)==4 && isOccupied(px-1,py+1)==2 && isOccupied(px-1,py+2)==1 && isOccupied(px,py+2)==0 && isOccupied(px+1,py+2)==1 && isOccupied(px+2,py+2)==0 && isOccupied(px+2,py+1)==3 && isOccupied(px+1,py)==1)numHeart++;
        if(po == 3 && isOccupied(px+1,py)==4 && isOccupied(px+1,py-1)==4 && isOccupied(px+1,py+1)==1 && isOccupied(px+2,py+1)==0 && isOccupied(px+2,py)==3 && isOccupied(px+2,py-1)==0 && isOccupied(px+2,py-2)==3 && isOccupied(px+1,py-2)==2 && isOccupied(px,py-1)==0)numHeart++;
        if(po == 2 && isOccupied(px,py-1)==4 && isOccupied(px-1,py-1)==4 && isOccupied(px+1,py-1)==0 && isOccupied(px+1,py-2)==3 && isOccupied(px,py-2)==2 && isOccupied(px-1,py-2)==3 && isOccupied(px-2,py-2)==2 && isOccupied(px-2,py-1)==1 && isOccupied(px-1,py)==3)numHeart++;
        if(po == 1 && isOccupied(px-1,py)==4 && isOccupied(px-1,py+1)==4 && isOccupied(px-1,py-1)==3 && isOccupied(px-2,py-1)==2 && isOccupied(px-2,py)==1 && isOccupied(px-2,py+1)==2 && isOccupied(px-2,py+2)==1 && isOccupied(px-1,py+2)==0 && isOccupied(px,py+1)==2)numHeart++;
    }
    return numHeart;
}

int main ()
{
    unsigned short xi, yi, oi, numCircs, numDiams, numHearts, tOr;
    unsigned char gen, maxGens;
    pt[0] = getPt(8192,8192,0);

    printf("This tool will calculate the number of circles, diamonds and hearts in the Q-toothpick sequence (Sloane's A187210).\nPlease enter the number of terms to compute (an integer from 1 to 8191): ");
    scanf("%d",&maxGens);

    file = fopen("numCircles.txt","w");
    fprintf(file,"0 0\n");
    fclose(file);
    file = fopen("numDiamonds.txt","w");
    fprintf(file,"0 0\n");
    fclose(file);
    file = fopen("numHearts.txt","w");
    fprintf(file,"0 0\n");
    fclose(file);
    printf("0 Circles: 0, Diamonds: 0, Hearts: 0\n");

    for(gen=1;gen<maxGens;gen++){
        numCircs = countCircles();
        numDiams = countDiamonds();
        numHearts = countHearts();
        file = fopen("numCircles.txt","a");
        fprintf(file,"%d %d \n",gen,numCircs);
        fclose(file);
        file = fopen("numDiamonds.txt","a");
        fprintf(file,"%d %d \n",gen,numDiams);
        fclose(file);
        file = fopen("numHearts.txt","a");
        fprintf(file,"%d %d \n",gen,numHearts);
        fclose(file);
        printf("%d Circles: %d, Diamonds: %d, Hearts: %d\n",gen,numCircs,numDiams,numHearts);

        for(xi=8192-gen;xi<=8192+gen;xi++){
            for(yi=8192-gen;yi<=8192+gen;yi++){
                tOr = isOccupied(xi,yi);
                for(oi=0;oi<4;oi++){
                    if(tOr==4 && addOrien(xi,yi,oi))pt[numPts++]=getPt(xi,yi,oi);
                }
            }
        }
        prevPts = numPts;
    }

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