// Calculate pi in decimal by Dave Coffin  1/25/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[3];

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

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

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

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

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

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

  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 (3, size);
  if (!reg[0]) {
    fprintf (stderr, "Not enough memory");
    exit(1);
  }
  reg[1] = reg[0] + size;
  reg[2] = reg[1] + size;
  puts("This program calculates π to arbitrarily many decimal places.\n"
	"It works best with three abaci, but if you have only one, calculate\n"
	"and write down all the A values, then the B values, then add them up.\n"
	"\"Undivide by 4\" is actually an in-place division by 25 (× 0.04).\n");
  puts ("Initialize to calculate 4 * atan(1/5):");
  maxden = floor(digits/log10(5)/2)*2-3;
  dwide = floor(log10(maxden)) + 1;
  reg[0][3+dwide] = reg[2][2] = 8;
  show (0, digits, dwide+1);
  for (denom=3; denom <= maxden; denom+=2) {
    printf ("Undivide A with \"Q 4\", copy to B, divide by %0*d, %s C:\n",
	dwide, denom, (denom & 2) ? "subtract from":"add to");
    for (i=size-3; i > 0; i--) {
      reg[0][i+2] = reg[0][i]*4;
      reg[0][i] = 0;
    }
    carry (reg[0], size);
    memcpy (reg[1], reg[0], size);
    divide (reg[1], digits, dwide, denom);
    sign = (denom & 2) ? -1:1;
    for (i=1; i <= digits; i++)
      reg[2][i] += sign * reg[1][i];
    carry (reg[2], size);
    show (0, digits, dwide+1);
  }
  puts ("Now calculate and subtract atan(1/239):");
  maxden = floor(digits/log10(239)/2)*2-1;
  dwide = floor(log10(maxden)) + 1;
  memset (reg[0], 0, 2*size);
  reg[0][5] = 1;
  show (0, digits, dwide+1);
  puts ("Divide A by 239 and subtract from C:");
  divide (reg[0], digits, 3, 239);
  for (i=1; i < digits; i++)
    reg[2][i] -= reg[0][i];
  carry (reg[2], digits);
  show (0, digits, 4);
  shift = 7+dwide;
  for (denom=3; denom <= maxden; denom+=2) {
    printf ("Divide A by 57121 (239²), shift right %d places, copy to B, divide by %0*d, %s C:\n",
	shift, dwide, denom, (denom & 2) ? "add to":"subtract from");
    divide (reg[0], digits, 5, 57121);
    for (i = size; i--; )
      reg[0][i] = i < shift ? 0 : reg[0][i-shift];
    shift = 6;
    memcpy (reg[1], reg[0], size);
    divide (reg[1], digits, dwide, denom);
    sign = (denom & 2) ? 1:-1;
    for (i=1; i <= digits; i++)
      reg[2][i] += sign * reg[1][i];
    carry (reg[2], size);
    show (0, digits, dwide+1);
  }
  puts ("That's π/4. Now add C to itself to get π/2, π, and 2π:");
  for (i=0; i < 3; i++) {
    for (j=1; j <= digits; j++)
      reg[2][j] *= 2;
    carry (reg[2], size);
    show (2, digits, 0);
  }
  free(reg[0]);
}