// Calculate trig functions in decimal by Dave Coffin  1/28/2020

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

signed char *reg[5];
int ones;

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

  size += pad;
  for (; i < 5; i++) {
    if (i==2) size -= pad;
    printf ("%c ", "AASCE"[i]);
    for (j=1; j <= size; j++) {
      putchar (reg[i][j]+'0');
      if (j == ones) 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, q, rem;

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

int main (int argc, char **argv)
{
  double fd=10, x=1, ln=0;
  int digits=0, size, maxden=1, dwide, denom, i, j, sign;
  char *arg = "10";

  if (argc > 1) digits = atoi(argv[1]);
  if (!digits) {
    fprintf (stderr, "Usage: %s <number of digits> <argument>\n", argv[0]);
    exit(1);
  }
  if (argc > 2) {
    arg = argv[2];
    x=i=j=0;  do {
      arg[i] = arg[j];
      if (isdigit(arg[i]))
	x += (arg[i++]-'0') * (fd /= 10);
    } while (arg[j++]);
  }
  ones = floor(log10(exp(x))) + 1;
  size = (digits + ones) + 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;
  while ((ln += log(maxden) - log(x)) < digits * log(10)) maxden++;
  dwide = floor(log10(maxden)) + 1;
  printf ("Taylor series for x=%c.%s to %d decimal digits\n\n", arg[0], arg+1, digits);
  puts ("In each step, you may change the number of zeros after the Q's,");
  puts ("provided you do the same to both Q numbers. If there isn't at");
  puts ("least one zero, you'll have to hold quotient digits in your head");
  puts ("while dividing and undividing.\n");
  puts ("If you only have one abacus, write down the final A value for");
  puts ("each step and sum them up later.");
  puts ("\nInitial setting independent of x:");
  reg[0][ones] = reg[1][ones] = reg[3][ones] = reg[4][ones] = 1;
  show (0, digits += ones, dwide);
  for (denom=1; denom <= maxden; denom++) {
    for (i=size; i--; ) {
      for (j=0; arg[j] && i+j+dwide < size; j++)
	reg[0][i+j+dwide] += reg[0][i] * (arg[j]-'0');
      reg[0][i] = 0;
      carry (reg[0], size);
    }
    printf ("Undivide A with \"Q");
    for (i=dwide; --i; ) putchar ('0');
    printf ("%s\", divide with \"Q%0*d\", add to E, %s %c:\n",
	arg, dwide, denom, (denom & 2) ? "subtract from":"add to", "CS"[denom & 1]);
    memcpy (reg[1], reg[0], size);
    divide (reg[1], digits, dwide, denom);
    sign = (denom & 2) ? -1:1;
    for (i=1; i <= digits; i++) {
      reg[3-(denom & 1)][i] += sign * reg[1][i];
      reg[4][i] += reg[1][i];
    }
    carry (reg[3-(denom & 1)], size);
    carry (reg[4], size);
    show (0, digits, dwide);
    memcpy (reg[0], reg[1], size);
  }
  printf ("\nThere you have it, the sine, cosine, and exponential of %c.%s\n", arg[0], arg+1);
  free(reg[0]);
}