// Calculate pi in dozenal by Dave Coffin  1/27/2020

// $Revision$
// $Date$

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <time.h>
#include <setjmp.h>

signed char *reg[5];

void put_digit (int dig)
{
//char sym[][8] = { "A","B" };
//char sym[][8] = { "T","E" };
//char sym[][8] = { "Ⅹ","Ɛ" };  // Dwiggins
  char sym[][8] = { "↊","↋" };	// Pitman

  if (dig > 9)
    fputs (sym[dig-10], stdout);
  else putchar (dig+'0');
}

void put_int (int wide, int val)
{
  int pow[] = { 1,12,144,1728,20736,248832,2985984,35831808,429981696 };

  while (--wide >= 0)
    put_digit (val / pow[wide] % 12);
}

void show (int i, int size, int pad)
{
  int j;

  size += pad;
  for (; i < 5; i++) {
    if (i==4) size -= pad;
    printf ("%c ", 'A'+i-(i>2));
    for (j=1; j <= size; j++) {
      put_digit (reg[i][j]);
      if (j == 1) putchar ('.');
    }
    putchar ('\n');
  }
}

void carry (signed char *board, int size)
{
  while (--size) {
    while (board[size] < 0) {
      board[size-1]--;
      board[size] += 12;
    }
    while (board[size] > 11) {
      board[size-1]++;
      board[size] -= 12;
    }
  }
}

void divide (signed char *reg, int digits, int dwide, int denom)
{
  int i, j, k, p, rem;

  for (i=1; i <= digits; i++) {
    for (rem=0, j=i+1; j < i+dwide+2; j++)
      rem = rem*12 + reg[j];
    reg[i] = rem / denom;
    for (k=denom, j=i+1+dwide; j > i; k /= 12, j--) {
      p = reg[i] * (k % 12);
      reg[j]   -= p % 12;
      reg[j-1] -= p / 12;
    }
    carry (reg+i+1, dwide+1);
  }
}

int main (int argc, char **argv)
{
  int digits=0, size, maxden, dwide, denom, i, j, sign;

  if (argc > 1) digits = atoi(argv[1]);
  if (!digits) {
    fprintf (stderr, "Usage: %s <number of digits>\n", argv[0]);
    exit(1);
  }
  size = ++digits+10;
  reg[0] = calloc (5, size);
  if (!reg[0]) {
    fprintf (stderr, "Not enough memory");
    exit(1);
  }
  for (i=0; i < 4; i++)
    reg[i+1] = reg[i] + size;
  puts("This program calculates π to arbitrarily many dozenal places.\n"
	"It works best with four abaci, but if you have only one, calculate\n"
	"and write down all the A values, then all the B values, etc.\n"
	"\"Undiv Q3\" and \"Undiv Q14\" do in-place division by four and nine\n"
	"respectively, i.e. multiplying by 0.3 and 0.14.\n");
  puts ("Initialize to calculate atan(1/2) + atan(1/3):");
  maxden = floor(digits*log(12)/log(2)/2)*2-5;
  dwide = floor(log(maxden)/log(12)) + 1;
  reg[0][3+dwide] = 6;
  reg[1][3+dwide] = 4;
  reg[4][2] = 10;
  show (0, digits, dwide+1);
  for (denom=3; denom <= maxden; denom+=2) {
    printf ("Undivide A with Q3, B with Q14, sum into C, divide with \"Q ");
    put_int (dwide, denom);
    printf ("\", %s D:\n", (denom & 2) ? "subtract from":"add to");
    for (i=size-2; i > 0; i--) {
      reg[0][i+1] = reg[0][i]*3;
      reg[0][i] = 0;
    }
    carry (reg[0], size);
    reg[1][size-1] = 0;
    for (i=size-3; i > 0; i--) {
      reg[1][i+2] += reg[1][i]*4;
      reg[1][i+1] = reg[1][i];
      reg[1][i] = 0;
    }
    carry (reg[1], size);
    for (i=1; i < size; i++)
      reg[2][i] = reg[1][i] + reg[0][i];
    carry (reg[2], size);
    memcpy (reg[3], reg[2], size);
    divide (reg[3], digits, dwide, denom);
    sign = (denom & 2) ? -1:1;
    for (i=1; i <= digits; i++)
      reg[4][i] += sign * reg[3][i];
    carry (reg[4], size);
    show (0, digits, dwide+1);
  }
  puts ("That's π/4. Now add D to itself to get π/2, π, and 2π:");
  for (i=0; i < 3; i++) {
    for (j=1; j <= digits; j++)
      reg[4][j] *= 2;
    carry (reg[4], size);
    show (4, digits, 0);
  }
  free(reg[0]);
}