// Abacus teaching program by Dave Coffin  9/11/2018
// Explains and demonstrates various algorithms I've learned over the years.

// $Revision: 1.11 $
// $Date: 2020/07/17 01:57:07 $

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

// This program is GPL if Readline is used, otherwise do as you please with it.
#ifdef READLINE
#include <readline/readline.h>
#include <readline/history.h>
#endif

// This should be big enough not to matter
#define SIZE 512

char desc[SIZE*2+16], explain[SIZE+16];			// null-terminated ASCII
char board[SIZE], arg[3][SIZE], multab[129][SIZE];	// raw arrays of digits
int size, len[3], base=10, colons=0, pitman=0, dwide=1, verbose=1, noob=1;
char prompt[16]="abacus-10> ";
jmp_buf error;

#define SWAP(a,b) { a=a+b; b=a-b; a=a-b; }

// All put_*() functions return the number of bytes written

int put_udec (char *str, unsigned val)
{
  char buf[16];
  int n=0, bytes;

  do {
    buf[n++] = val % 10;
    val /= 10;
  } while (val);
  for (bytes=n; n--; )
    *str++ = buf[n] + '0';
  return bytes;
}

// 128 means "not a digit". Every abacus should have a bead pattern for this.
int put_digit (char *str, unsigned char dig)
{
  if (colons) {
    if (dig == 128)
      str[0] = str[1] = str[2] = '-';
    else {
      str[0]       = '0' + dig / 100;
      str[dwide-3] = '0' + dig / 10 % 10;
      str[dwide-2] = '0' + dig % 10;
    }
    str[dwide-1] = ':';
    return dwide;
  }
  if (dig == 128) {
    str[0] = '|';
    return 1;
  }
  if (dig > 9) {
    if (pitman && dig < 12) {
      str[0] = 0xe2;
      str[1] = 0x86;
      str[2] = 0x80 + dig;
      return 3;
    }
    str[0] = dig-10 + 'A';
  } else str[0] = dig + '0';
  return 1;
}

int put_int (char *str, int val)
{
  char buf[40];
  int n, bytes;

  n = 0; do {
    buf[n++] = val % base;
    val /= base;
  } while (val);
  for (bytes=0; n--; )
    bytes += put_digit (str+bytes, buf[n]);
  return bytes-colons;		// drop the last colon
}

int bprintf (char *str, const char *fmt, ...)
{
  const char *fp;
  va_list argp;
  int i, len;
  char buf[SIZE*4], *op, *cp;

  op = str ? str : buf;
  va_start (argp, fmt);
  for (fp = fmt; *fp; fp++) {
    if (*fp != '%')
      *op++ = *fp;
    else switch (*++fp) {
      case '%':
	*op++ = '%';
	break;
      case '*':			// just print spaces
	i = va_arg (argp, int);
	while (i--) *op++ = ' ';
	break;
      case 'q':			// %q to mark quotient digit
      case 'Q':
	for (i=0; i < dwide-colons; i++)
	  *op++ = 'Q';		// could use × instead of Q
	for (; i < dwide << (*fp == 'q'); i++)
	  *op++ = ' ';
	break;
      case 'c':			// %c for single char
	i = va_arg (argp, int);
	*op++ = i;
	break;
      case 'u':			// %u for unsigned decimal
	i = va_arg (argp, int);
	op += put_udec (op, i);
	break;
      case 'i':			// %i for integer in abacus base
	i = va_arg (argp, int);
	op += put_int (op, i);
	break;
      case 'n':			// %n for array of digits
	len = va_arg (argp, int);
	cp = va_arg (argp, char *);
	for (i=0; i < len; i++)
	  op += put_digit (op, cp[i]);
	op -= colons;
	break;
      case 's':			// %s for null-terminated ASCII
	cp = va_arg (argp, char *);
	while (*cp) *op++ = *cp++;
	break;
    }
  }
  *op = 0;
  if (!str) fputs (buf, stdout);
  va_end (argp);
  return op - (str ? str : buf);
}

int input_num (unsigned char *ip, char *op, int base, int colons, int shift)
{
  char *sp;
  int dval, i, o, zero=1;

  for (dval=o=0; *ip && o < SIZE-shift-1; ip++) {
    if (colons) {
      if (isdigit(*ip)) dval = dval*10 + *ip - '0';
      if (!isdigit(ip[1]))
	if (dval < base) {
	  op[o++] = dval;
	  dval = 0;
	} else {
	  bprintf (0, "Digit \'%u\' is too big!\n", dval);
	  longjmp (error, 2);
	}
      goto repch;
    } else {
      if ((unsigned) (*ip-'0') < (base < 10 ? base:10))
	op[o++] = *ip-'0';
      else if (toupper(*ip) >= 'A' && toupper(*ip) < 'A'+base-10)
	op[o++] = toupper(*ip)-'A'+10;
      else if (*ip == 0xe2 && *++ip == 0x86 && (*++ip | 1) == 0x8b)	// Pitman digits
	op[o++] = *ip - 0x80;
      else repch:
      if (isalpha(*ip) && (sp = strchr ("h2t3m6sH", tolower(*ip))))
	for (i = sp[1]-'0'; i-- && o < SIZE-shift-1; )
	  op[o++] = 0;
      else if (!colons && *ip != ',') {
	bprintf (0, "Illegal digit \"%c\".\n", *ip);
	longjmp (error, 2);
      }
    }
    if (o && op[o-1]) zero = 0;
    if (zero && shift && ip[1]) o = 0;
  }
  if (shift && (shift = (720720 - o) % shift))
    for (i = o += shift; i >= 0; i--)
      op[i] = (i < shift) ? 0: op[i-shift];
  return o;
}

void user_help()
{
  bprintf (0,
"This is the recommended order to learn (<N> is always base ten):"
"\n   add <num> <num>"
"\n   sub <num> <num>"
"\n   neg <num>"
"\n   mul <num> <num>"
"\n   div <num> <num>"
"\n   undiv <num> <num> [<rem>]"
"\n   square <num>"
"\n   cube <num>"
"\n   sqrt <num>"
"\n   cbrt <num> [<N>] - cube root, no more than N digits per cycle"
"\n   perm - test your knowledge of the %i×%i multiplication table"
"\n   perd, pers - same thing, but with division and square root"
"\n   rand <N> [<num>] - show <num> times a random N-digit number"
"\n   rand2 <N> - show the square of a random N-digit number"
"\n   rand3 <N> - show the cube a random N-digit number"
"\n   base <N> - change the base of computation (BoC)"
"\n   tbase <N> <num> - convert <num> to base <N> (from the BoC)"
"\n   fbase <N> <num> - convert <num> from base <N> (to the BoC)"
"\n   quit (or \"q\")"
"\nNumbers may include \"h\" for \"00\", \"t\" for \"000\","
"\n\"m\" for \"000000\", and \"s\" for two dozen zeroes.\n", base-1, base-1);
}

void set_size (int req)
{
  size = req;
  if (size >= SIZE) {
    bprintf (0, "%u digits needed, only %u available!\n", size, SIZE-1);
    longjmp (error, 1);
  }
}

// board[0] is never displayed; it's only there to absorb carries.
void show (int underline, int normal)
{
  int i, j;
  char line[SIZE*4+16];

  if (verbose) {
    if (desc[0]) puts (desc);
    line[0] = ' ';
    for (j=i=1; i <= size; i++) {
      if (i == underline) memcpy (line+j, "\e[4m", 4), j+=4;
      j += put_digit (line+j, board[i]) - colons;
      if (i == normal) memcpy (line+j, "\e[0m", 4), j+=4;
      if (colons && i < size) line[j++] = ':';
    }
    line[j] = 0;
    puts (line);
  }
  desc[0] = 0;
}

void carry (char *board, int size, int sign)
{
  for ( ; size > 0; size--)
    if ((unsigned) board[size] >= base) {
      board[size-1] += sign;
      board[size] -= sign*base;
    }
}

void add_num (char *board, int pos, int sign, int len, char *num)
{
  int i;

  for (i=0; i < len; i++)
    board[pos+i] += sign * (num[i]);
  carry (board, pos+i-1, sign);
}

// val is in units of one-half and uses one's-complement
// for negative numbers so we can display "-0".
void add_int (int ones, int val)
{
  int sign=1, half, disp, i;

  if (val < 0) {
    val = ~val;
    sign = -1;
  }
  if (half = val & 1)
    board[ones+1] += base/2 * sign;	// Don't try this in odd-numbered bases!
  disp = val >>= 1;
  i = ones; do
    board[i] += val % base * sign;
  while (i--, val /= base);
  bprintf (desc, "%*%c%i%s %s", i*dwide, ','-sign, disp, half ? "½":"", explain);
  carry (board, ones+half, sign);
  explain[0] = 0;
}

void do_add (int sign)
{
  int pos, i;

  set_size ((len[0] > len[1] ? len[0]:len[1]) + 1);
  pos = size - len[0] + 1;
  memcpy (board+pos, arg[0], len[0]);
  show (0, 0);
  pos = size - len[1] + 1;
  for (i=0; i < len[1]; i++) {
    add_int (pos+i, (arg[1][i] << 1) ^ sign);
    show (0, 0);
  }
}

void user_add()
{
  static int hint=1;

  if (hint)
    puts ("\nAs you push beads with one hand, keep a finger"
	"\nof the other hand firmly on the digit being added"
	"\nso as not to lose your place while carrying.\n");
  hint = 0;
  do_add (0);
}

void user_sub()
{
  static int hint=1;

  if (hint)
    bprintf (0,"\nIf the answer starts with a \"%i\", it is negative."
	"\nSee the \"neg\" instruction for more info.\n", base-1);
  hint = 0;
  do_add (-1);
}

void user_neg()
{
  int i;
  static int hint=1;

  if (hint)
    bprintf (0,"\nNegation is very easy:"
	"\n%u\'s-complement the rightmost non-zero digit, then"
	"\n%u\'s-complement all digits to its left."
	"\n\nNegation is useful when you have to calculate e.g."
	"\nA - B × C. Do B × C first, then negate it and add A.\n", base, base-1);
  hint = 0;
  memcpy (board+1, arg[0], size = len[0]);
  show (0, 0);
  desc[size+1] = 0;
  for (i = size; i > 0 && !board[i]; i--)
    desc[i] = 'u';
  desc[i] = 't';
  board[i] = base - board[i];
  while (--i > 0) {
    desc[i] = 'n';
    board[i] = base-1 - board[i];
  }
  desc[0] = ' ';
  show (0, 0);
}

// multab[] may be used when noob == 0
void noob_table (int verbose)
{
  char board[SIZE];
  int i;

  memset (board, 0, SIZE);
  if (verbose) bprintf (0, "\nRepeatedly add %n on the abacus\n"
	"and write down your answers like so:\n", len[1], arg[1]);
  for (i=0; i <= base; i++) {
    memcpy (multab[i], board, len[1]+1);
    add_num (board, 1, 1, len[1], arg[1]);
    if (verbose && i)
      bprintf (0, "\n%*%i    %n", dwide*(i < base), i, len[1]+1, multab[i]);
  }
  if (verbose) puts (" (check for mistakes here)");
}

// Multiply whatever's on the board by arg[1]
void do_mul (char *board, int size, int noob)
{
  int i, j, mul;

  if (verbose) putchar('\n');
  for (i=1; i <= size; i++) {
    if (!(mul = board[i])) continue;
    show (i, size);
    if (noob) {
      bprintf (desc, "%*+%n (%n × %i)", i-len[1]-1, len[1]+1, multab[mul], len[1], arg[1], mul);
      board[i] = 0;
      add_num (board, i-len[1], 1, len[1]+1, multab[mul]);
    } else {
      if (verbose) bprintf (0, "%*%n (paper)\n", (i-len[1])*dwide+1, len[1], arg[1]);
      for (j=0; j < len[1]; j++) {
	bprintf (explain, "(%i × %i)", mul, arg[1][j]);
	if (j == len[1]-1) board[i] = 0;
	add_int (i-len[1]+1+j, mul * arg[1][j] * 2);
	if (j < len[1]-1) show (i, size);
      }
    }
  }
  show (0, 0);
}

void user_mul()
{
  static int hint=1;

  set_size (len[0] + len[1]);
  noob_table (verbose && noob);
  if (hint && verbose) {
    bprintf (0, "\nPlace all beads of the multiplicand \e[4m%n\e[0m"
	"\nto the right side of the abacus and slightly off"
	"\nthe bar as a reminder that they are to be replaced,"
	"\nnot added to.\n", len[0], arg[0]);

    if (!noob) {
      bprintf (0, "\nWrite the multiplier %n on a slip of"
	"\npaper with the last digit under the first digit"
	"\nof the multiplicand. Move it to the right as the"
	"\nmultiplicand is consumed.\n", len[1], arg[1]);
      puts ("\nYou can add simple products one digit at a time,\n"
	"but try to gain speed by learning their aliases,\n"
	"e.g. 9 × 3 = 30-3, 7 × 4 = 30-2, 7 × 8 = 100-44");
    }
  }
  hint = 0;
  memcpy (board+len[1]+1, arg[0], len[0]);
  do_mul (board, size, noob);
}

void do_div (int qpos0, int steps, int noob)
{
  int look, rup, qpos, guess, i;

  look = len[1]+1;
  if (!noob) {
    look = 2;
    rup = arg[1][0];
    for (i = 1; i < len[1]; i++)
      if (arg[1][i] && rup++) break;
    bprintf (0, "   %n rounded up to a single digit is %i.\n", len[1], arg[1], rup);
  }
  for (qpos = qpos0; qpos < qpos0+steps; qpos++) {
    bprintf (0, "%*%q%n\n", (qpos-1)*dwide+1, len[1], arg[1]);
    if (noob) {
      for (guess = base+1; --guess; )
	if (memcmp (multab[guess], board+qpos+1, len[1]+1) <= 0) break;
    } else {
      guess = (board[qpos+1]*base + board[qpos+2])/rup;
      if (guess < 2) guess = 0;
    }
    board[qpos] = guess;
    if (noob < 2) {
      if (guess) {
	i = bprintf (desc, "%*%i = floor(%n ÷ ",
		(qpos-1)*dwide+1, guess, look, board+qpos+1);
	if (noob) bprintf (desc+i, "%n)", len[1], arg[1]);
	else      bprintf (desc+i, "%i)", rup);
      }
      show (qpos0, qpos);
    }
    if (guess) {
      if (noob) {
	if (noob < 2)
	  bprintf (desc, "%*-%n (%n × %i)", qpos,
		len[1]+1, multab[guess], len[1], arg[1], guess);
	add_num (board+qpos, 1, -1, len[1]+1, multab[guess]);
	show (qpos0, qpos);
      } else
	for (i=0; i < len[1]; i++) {
	  bprintf (explain, "(%i × %i)", guess, arg[1][i]);
	  add_int (qpos+2+i, ~(arg[1][i] * guess * 2));
	  show (qpos0, qpos);
	}
    } else if (noob == 2) show (qpos0, qpos);
    while (memcmp (multab[1], board+qpos+1, len[1]+1) <= 0) {
      board[qpos]++;
      add_num (board, qpos+2, -1, len[1], arg[1]);
      bprintf (desc, "%*+1%*-%n", qpos*dwide-1-colons, dwide-1+colons, len[1], arg[1]);
      show (qpos0, qpos);
    }
  }
}

void user_div()
{
  static int hint=1;

  set_size (len[0] + 1 + (memcmp (arg[0], arg[1], len[1]) >= 0));
  if (!arg[1][0]) {
    puts ("Division by zero is not allowed.");
    return;
  }
  if (len[0] < len[1] || (len[0] == len[1] && memcmp(arg[0],arg[1],len[0]) < 0)) {
    bprintf (0, "%n / %n is zero with remainder %n. For more precision,\n"
	"append 0\'s, h\'s, t\'s, or m\'s to the dividend.\n",
	len[0], arg[0], len[1], arg[1], len[1], arg[1]);
    return;
  }
  noob_table (noob);
  if (hint) {
    if (!noob)
      bprintf (0, "\nNotice how this program guesses a low quotient digit, then"
	"\nsteps it up until the remainder is less than the divisor."
	"\nIn most cases, it's faster to guess high and then step"
	"\ndown until the remainder is non-negative.\n");
  }
  hint = 0;
  bprintf (0, "\nWrite \"%q%n\" on a slip of paper.\n\n", len[1], arg[1]);
  memcpy (board+size-len[0]+1, arg[0], len[0]);
  show (0, 0);
  do_div (1, size-len[1]-1, noob);
}

void do_undiv (int qpos0, int steps, int noob)
{
  int qpos, qval, i;

  for (qpos=qpos0; qpos > qpos0-steps; qpos--) {
    if ((qval = board[qpos]) == 0) continue;
    show (qpos0-steps+1, qpos);
    bprintf (0, "%*%q%n\n", (qpos-1)*dwide+1, len[1], arg[1]);
    if (noob) {
      if (noob < 2)
	bprintf (desc, "%*+%n (%n × %i)", qpos, len[1]+1, multab[qval], len[1], arg[1], qval);
      add_num (board, qpos+1, 1, len[1]+1, multab[qval]);
    } else
      for (i=0; i < len[1]; i++) {
	bprintf (explain, "(%i × %i)", qval, arg[1][i]);
	add_int (qpos+i+2, qval * arg[1][i] * 2);
	if (i != len[1]-1) show (qpos0-steps+1, qpos);
      }
    board[qpos] = 0;
  }
}

void user_undiv()
{
  set_size (len[0] + len[1] + 1);
  if (len[2] > len[1]) {
    puts ("Remainder is too big!");
    return;
  }
  noob_table (noob);
  bprintf (0, "\nWrite \"%q%n\" on a slip of paper.\n\n", len[1], arg[1]);
  memcpy (board+1, arg[0], len[0]);
  memcpy (board+size-len[2]+1, arg[2], len[2]);
  do_undiv (len[0], len[0], noob);
  show (0, 0);
}

void do_square (char *board, int size, int pos)
{
  int i, j;

  memcpy (arg[1], arg[0], len[0]);
  len[1] = len[0];
  noob_table (verbose && noob);
  if (noob) {
    memcpy (board+pos+len[1]+1, arg[0], len[0]);
    do_mul (board, size, noob);
    return;
  }
  if (verbose) {
    bprintf (0, "\nMultiplication requires that every digit of one number be"
	"\nmultiplied by every digit of the other, but there's a shortcut"
	"\nwhen multiplying a number by itself.\n"
	"\nZero the abacus and write the number on a slip of paper with"
	"\nblank spaces between the digits. This paper does not move.\n"
	"\nThen multiply each paper digit by all digits to the right,"
	"\nand add each product exactly half-way between them:\n %*",pos*dwide);
    for (i=0; i < len[1]; i++)
      bprintf (0, "%*%n", (i && colons)+dwide, 1, arg[1]+i);
    putchar ('\n');
  }
  show (0, 0);
  for (i=0; i < len[1]; i++) {
    if (!arg[1][i]) continue;
    for (j = i+1; j < len[0]; j++) {
      bprintf (explain, "(%i × %i)", arg[1][i], arg[1][j]);
      add_int (pos+i+j+2, arg[1][i] * arg[1][j] * 2);
      show (0, 0);
    }
  }
  memcpy (arg[0], board+pos+2, size-pos-2);
  len[0] = size-pos-2;
  if (verbose) {
    bprintf (0, "\nNext, multiply each paper digit by all digits to the left."
	"\nThis yields the same result, so just double what you have:\n");
    show (0, 0);
    bprintf (0, "%*+%n\n", (pos+1)*dwide, len[0], arg[0]);
  }
  add_num (board, pos+2, 1, len[0], arg[0]);
  show (0, 0);
  if (verbose) {
    bprintf (0, "\nLastly, multiply each paper digit by itself:\n %*", pos*dwide);
    for (i=0; i < len[1]; i++)
      bprintf (0, "%*%n", (i && colons)+dwide, 1, arg[1]+i);
    putchar ('\n');
  }
  show (0, 0);
  for (i=0; i < len[1]; i++) {
    bprintf (explain, "(%i²)", arg[1][i]);
    add_int (pos+i*2+2, arg[1][i] * arg[1][i] * 2);
    show (0, 0);
  }
}

void user_square()
{
  set_size (len[0] * 2);
  do_square (board, size, 0);
}

void user_cube()
{
  set_size (len[0] * 3);
  do_square (board, size, len[0]);
  do_mul (board, size, noob);
}

void user_sqrt()
{
  int val, guess, half, root, i, j;
  long long num, den;
  char sroot[SIZE];
  static int hint=1;

  if (base & 1) {
    puts ("Sorry, square root in odd bases is not supported.");
    return;
  }
  set_size (len[0] + 2);
  if (hint)
    bprintf (0, "\nAppend h's, m's, or s's if you want more precision."
	"\nThis is Kato Fukutaro's method, based on the formula:\n"
	"\n(a + b)² = a² + 2ab + b² = a² + 2(ab + ½b²)\n"
	"\nWe'll subtract a², halve the remainder, and subtract"
	"\nsimple products and half-squares from then on.\n"
	"\nAs with long division, it's usually faster to guess"
	"\nhigh and step down.\n");
  hint = 0;
  bprintf (0, "\nCalculating the square root of");
  for (i=0; i < len[0]; i += 2)
    bprintf (0, "%c%n", i ? ',':' ', 2, arg[0]+i);
  putchar ('\n');
  memcpy (board+2, arg[0], len[0]);
  show (0, 0);
  for (val=0, i=2; i < 6; i++)
    val = val*base + board[i];
  for (guess = base-1; guess*guess*base*base > val; guess--);
  board[1] = guess;
  bprintf (desc," %i = floor(sqrt(%i))", guess, val/base/base);
  show (1, 1);
  bprintf (explain, "(%i²)", guess);
  add_int (3, ~(guess * guess * 2));
  show (1, 1);
  for (i = size-1; i > 1; i--)
    if (board[i]) {
      if (half = board[i] & 1)
	board[i+1] += base/2;
      board[i] = board[i] >> 1;
      bprintf (desc, "%*-%i%s", i-1, board[i]+half, half ? "+½":"");
      show (1, 1);
    }
  for (i=guess*base, guess=base-1; (i+guess)*(i+guess) > val; guess--);
  for (root=2; root <= len[0]/2; root++) {
    if ((j = root) > 6) j = 6;
    for (num=i=0; i < j; i++)
      num = num*base + board[root+1+i];
    for (den=i=0; i < j-1; i++)
      den = den*base + board[i+1];
    if (root > 2) guess = (num*8 / den - 3) >> 3;
    if (guess < 2) guess = 0;
    bprintf (desc, "%*%i < %n / %n", (root-1)*dwide+1, guess, j, board+1+root, j-1, board+1);
    board[root] = guess;
    show (1, root);
    if (guess) {
      if (noob) {
	memcpy (arg[1], board+1, root-1);
	arg[1][len[1] = root-1] = 0;
	noob_table(0);
	bprintf (desc, "%*-%n (%n × %i)", root, len[1]+1, multab[guess], len[1], arg[1], guess);
	add_num (board, root+1, -1, len[1]+1, multab[guess]);
	show (1, root);
      } else
	for (i=1; i < root; i++) {
	  bprintf (explain, "(%i × %i)", guess, board[i]);
	  add_int (i+1+root, ~(guess * board[i] * 2));
	  show (1, root);
	}
      bprintf (explain, "(½%i²)", guess);
      add_int (root*2+1, ~(guess*guess));
      show (1, root);
    }
    while (1) {
      memcpy (sroot, board, root+1);
      sroot[root+1] = base/2;
      bprintf (desc, "%*+1%*-%n½", root*dwide-1-colons, dwide-1+colons, root, sroot+1);
      if (memcmp (board+root+1, sroot, root+2) < 0) break;
      add_num (board, root+2, -1, root+1, sroot+1);
      board[root] = ++guess;
      show (1, root);
    }
  }
}

void user_cbrt()
{
  int kato, third, grow, guess, val, root=1, maxbite, bite, b, i;
  const char vfrac[3][4] = { "", "⅓", "⅔" };
  char save[SIZE];
  static int hint[2] = { 1, 1 };

  if (base == 2) {
    puts ("Sorry, cube root in base two is not supported.");
    return;
  }
  kato = !(base % 3);
  if (strchr (arg[1], 'k')) kato = 0;
  if (hint[kato]) {
    bprintf (0, "\nAppend t's, m's, or s's if you want more precision."
	"\nThis algorithm is based on the formula:\n");
    if (kato)
      bprintf (0, "\n(a + b)³ = a³ + 3a²b + 3ab² + b³ = a³ + 3(a(ab + b²) + ⅓b³)\n"
	"\nWe'll subtract a³, multiply the remainder by 1/3 = 0.%i, divide"
	"\nby a, subtract ab and b², undivide by a, and subtract ⅓b³.", base/3);
    else
      bprintf (0, "\n(a + b)³ = a³ + 3a²b + 3ab² + b³ = a³ + 3a(ab + b²) + b³\n"
	"\nWe'll subtract a³, divide by 3a, subtract ab and b², then"
	"\nundivide by 3a and subtract b³.");
    bprintf (0, "\nAppend b's digits to a and divide by the new %sa to extract"
	"\nmore digits. b must never have more digits than a.\n"
	"\nFor a \"fence\", I move all beads to the center of the wire.\n",
	kato ? "":"3");
  }
  hint[kato] = 0;
  bprintf (0, "\nCalculating the cube root of");
  for (i=0; i < len[0]; i += 3)
    bprintf (0, "%c%n", i ? ',':' ', 3, arg[0]+i);
  putchar ('\n');
  if (kato) {
    bprintf (0, "\nBecause the root will not be multiplied by three,\n");
    grow = 0;
  } else {
    for (b=1, i=0; i < 30 || i < len[0]; i++, b %= 27)
      save[i] = (b *= base) / 27;
    grow = (i = memcmp (arg[0], save, len[0])) > 0;
    if ((i | base % 3) == 0) grow++;
    bprintf (0, "\nBecause this number is %s 1/3³\n(",
	grow ? "greater than or equal to":"less than");
    for (i=0; i < 30; i += 3)
      bprintf (0, "%n,", 3, save+i);
    bprintf (0, grow ? "...),\n3a must always be one digit longer than a, and\n" :
		 "...),\n3a will always be the same length as a, and\n");
  }
  bprintf (0, "there must be %u digits to the left of the first comma.\n\n", 5+grow);
  set_size (len[0] + 2 + grow + kato);
  if (!(maxbite = atoi(arg[1])))
    maxbite = SIZE;
  memcpy (board+3+grow, arg[0], len[0]);
  show (0, 0);
  val = (arg[0][0]*base + arg[0][1])*base + arg[0][2];
  for (guess = base-1; guess*guess*guess > val; guess--);
  board[1] = guess;
  bprintf (desc," %i = floor(cbrt(%i))", guess, val);
  show (1, 1);
  bprintf (explain, "(%i³)", guess);
  add_int (5+grow, ~(guess*guess*guess*2));
  show (0, 0);
  for (third=0, root=1; root < len[0]/3; root += bite) {
    if (kato && !third && (third = 1)) {
      puts ("Third the remainder in place. This will not be undone.");
      show (0, 0);
      for (i = size-1; i > root+1; i--) {
	board[i+1] += b = board[i] % 3 * base / 3;
	bprintf (desc, "%*%i%s%i = %i × %i", (i-1)*dwide+1,
		board[i]/3, colons ? ":":"", b, base/3, board[i]);
	board[i] /= 3;
	show (0, 0);
      }
    }
    if ((bite = root) > maxbite) bite = maxbite;
    if (bite >= len[0]/3 - root)
	bite  = len[0]/3 - root;
    else if (root + bite*2 > len[0]/3)
      bite = (len[0]/3 - root) >> 1;
    if (kato)
      memcpy (arg[1], board+1, len[1] = root);
    else {
      for (arg[1][0] = i=0; i < root; i++) {
	arg[1][i+grow-1] += 3*board[i+1] / base;
	arg[1][i+grow  ]  = 3*board[i+1] % base;
      }
      len[1] = root+grow;
      carry (arg[1] - 1, len[1], 1);
    }
    noob_table (0);
    do_div (root+2, root+bite*2, 2);
    bprintf (0, "Set a fence between quotient and remainder\n");
    board[(root+bite+1)*2] = 128;
    show (0, 0);
    for (b=1; b <= bite; b++) {
      if (board[root+b]) {
	bprintf (0, "We need a zero for the next root digit!\n");
	board[root+b]--;
	board[root+b+1] += base;
	show (0, 0);
      }
      memcpy (save, board, SIZE);
      for (verbose=0, guess = base-1; guess >= 0; guess--) {
	memcpy (board, save, SIZE);
	board[root+b] = guess;
	bprintf (desc, "%*%i < %n / %n",
	  (root+b-1)*dwide+1, guess, root+1, board+root+b+1, root, board+1);
	show (root+1, root+b);
	if (guess) {
	  for (i=1; i <= root; i++) {
	    bprintf (explain, "(%i × %i)", guess, board[i]);
	    add_int (root+b+1+i, ~(guess * board[i] * 2));
	    show (root+1, root+b);
	  }
	  for (; i < root+b; i++) {
	    bprintf (explain, "(%i × %i × 2)", guess, board[i]);
	    add_int (root+b+1+i, ~(guess * board[i] * 4));
	    show (root+1, root+b);
	  }
	  bprintf (explain, "(%i²)", guess);
	  add_int (root+b+1+i, ~(guess * guess * 2));
	  show (root+1, root+b);
	}
	if (board[root+b] != guess) continue;	// remainder went negative
	if (verbose) break;
	guess += verbose = 1;
      }
    }
    bprintf (0, "Remove the fence and undo the divison\n");
    board[(root+bite+1)*2] = 0;
    do_undiv ((root+bite)*2+1, root+bite+1, 2);
    show (root+1, root+bite);
    memcpy (arg[1], board+root+1, bite);
recover:
    len[1] = bite;
    memset (save, 0, SIZE);
    if (third) save[bite*3+1] = base/3;
    else save[bite*3] = 1;
    noob_table (verbose = 0);
    for (i=0; i < 3; i++)
      do_mul (save, bite*3+1, 1);
    bprintf (desc, "%*-%n%s (%s%n³)", (root*3+grow+2)*dwide, bite*3, save+1,
		vfrac[save[bite*3+1]/(base/3)], vfrac[third], bite, arg[1]);
    add_num (board, root*3+3+grow, -1, bite*3+1, save+1);
    verbose = 1;
    if (!board[1])
      memcpy (board+root+1, arg[1], bite);
    else if (board[root+bite] != guess) {
      show (0, 0);
      bprintf (0, "Oops, the remainder went negative! Time for a reset:\n");
      bite += root;
      third = root = 0;
      memcpy (arg[1], board+1, bite);
      memset (board, 0, SIZE);
      memcpy (board+3+grow, arg[0], len[0]);
      show (0, 0);
      goto recover;
    }
    show (0, 0);
  }
}

void solution (char *op, int sarg)
{
#ifdef READLINE
  char solve[SIZE*2];
  int i;

  i = bprintf (solve, "%s %n", op, size, board+1);
  if (sarg) bprintf (solve+i, " %n", len[1], arg[1]);
  add_history (solve);
  puts ("Press ↑ or Ctrl-P for the solution.");
#endif
}

void do_perm()
{
  int a, i, j, t;

  len[0] = len[1] = base;
  for (a=0; a < 2; ) {
    for (i=0; i < base; i++)
      arg[a][i] = i;
    for (i=0; i < base; i++) {
      j = random() % base;
      t = arg[a][i];
      arg[a][i] = arg[a][j];
      arg[a][j] = t;
    }
    if (arg[a][0] && arg[a][base-1]) a++;
  }
}

void user_perm()
{
  do_perm();
  bprintf (0, "Multiply these numbers: %n %n\n", len[0], arg[0], len[1], arg[1]);
  memcpy (board+1, arg[0], size=base);
  solution ("undiv", 1);
}

void user_perd()
{
  do_perm();
  verbose = 0;
  user_mul();
  verbose = 1;
  bprintf (0, "Divide these numbers: %n %n\n", base*2, board+1, base, arg[1]);
  solution ("div", 1);
}

void user_pers()
{
  do_perm();
  memcpy (arg[0], arg[1], base);
  verbose = 0;
  user_mul();
  verbose = 1;
  bprintf (0, "Square root this number: %n\n", base*2, board+1);
  solution ("sqrt", 0);
}

void do_rand()
{
  int i;

  len[0] = atoi(arg[0]);
  for (i=0; i < len[0] && i < SIZE-1; i++)
    if (i && i < len[0]-1)
      arg[0][i] = random() % base;
    else
      arg[0][i] = random() % (base-1) + 1;
}

void user_rand()
{
  if (!arg[1][0]) arg[1][0] = 1;
  do_rand();
  verbose = 0;
  user_mul();
  verbose = 1;
  show (0, 0);
  solution ("div", 1);
}

void user_rand2()
{
  do_rand();
  verbose = 0;
  user_square();
  verbose = 1;
  show (0, 0);
  solution ("sqrt", 0);
}

void user_rand3()
{
  do_rand();
  verbose = 0;
  user_cube();
  verbose = 1;
  show (0, 0);
  solution ("cbrt", 0);
}

void parse_base (int *base, int *colons)
{
  int new = atoi(arg[0]);
  if (new < 2 || new > 128) {
    bprintf (0, "Valid bases are 2 to 128.\n");
    longjmp (error, 3);
  }
  *base = new;
  *colons = (strchr (arg[0], ':') && new > 10) || new > 36;
}

void user_base()
{
  parse_base (&base, &colons);
  pitman = (base == 12) && strchr (arg[0], 'p');
  if (base > 10 && base <= 36 && !colons)
    bprintf (0, "Type \"base %u:\" if you want colons.\n", base);
  if (base == 12 && !pitman) 
    bprintf (0, "Type \"base 12p\" if you want Pitman digits.\n");
  bprintf (prompt, "abacus-%u%c ", base, colons ? ':':'>');
  dwide = colons ? (base > 100 ? 4:3) : 1;
}

void user_tbase()
{
  int tbase, tcol, dw, qstop, qpos, u0, dig, i, j, k;
  char slider[8];

  parse_base (&tbase, &tcol);
  for (i = tbase-1, dw=0; i; dw++, i/= base);	// dw = # of digits in tbase-1
  for (j = tbase, i=dw; i >= 0; i--, j /= base)
    slider[i] = j % base;
  if (dw > 1)
    bprintf (0, "\nEach digit in base %u occupies %u digits in base %u.\n", tbase, dw, base);
  set_size (ceil (len[1] * log(base)/log(tbase)) * dw);
  bprintf (0, "\nTo convert out of the base of computation (%u), repeatedly divide by the"
	"\ntarget base (%u), piling up remainders on the right side of the board.\n\n",
	base, tbase);
  memcpy (board+1+size-len[1], arg[1], len[1]);
  show (0, 0);
  if (slider[0])
    bprintf (0, "\n%u is an exact power of %u, so we're done.\n", tbase, base);
  else {
    memcpy (arg[1], slider+1, len[1] = dw);
    noob_table (0);
    for (qstop = size-dw; qstop >= dw; qstop -= dw) {
      for (u0 = 0, qpos = 1; qpos <= qstop; qpos++) {
	for (dig = base-1; dig; dig--)
	  if (memcmp (board+qpos, multab[dig], dw+1) >= 0) break;
	if (!u0)
	  if (dig) u0 = qpos;
	  else continue;
	add_num (board, qpos, -1, len[1]+1, multab[dig]);
	board[qpos] = dig;
	bprintf (0, "%*%Q%n\n", (qpos-1)*dwide+1, dw, arg[1]);
	if (qpos < qstop) show (u0, qpos);
      }
      if (u0) show (0, 0);
    }
  }
  if (tbase > base || tcol != colons) {
    for (k=0, i=1; i <= size; i += dw) {
      for (dig=j=0; j < dw; j++)
	dig = dig*base + board[i+j];
      board[k++] = dig;
    }
    SWAP (base, tbase);
    SWAP (colons, tcol);
    dwide = colons ? (base > 100 ? 4:3) : 1;
    bprintf (0, "\nIn \"base %u%s\" syntax, this number is %n .\n\n",
		base, colons && base <= 36 ? ":":"", k, board);
    SWAP (base, tbase);
    SWAP (colons, tcol);
    dwide = colons ? (base > 100 ? 4:3) : 1;
  }
}

void user_fbase()
{
  int fbase, fcol, dw, qstop, qpos, dig, i, j;
  char slider[8];

  parse_base (&fbase, &fcol);
  len[0] = input_num (arg[1], arg[0], fbase, fcol, 1);
  for (i = fbase-1, dw=0; i; dw++, i/= base);	// dw = # of digits in fbase-1
  for (j = fbase, i=dw; i >= 0; i--, j /= base)
    slider[i] = j % base;
  if (dw > 1)
    bprintf (0, "\nEach digit in base %u occupies %u digits in base %u.\n", fbase, dw, base);
  set_size (len[0] * dw);
  bprintf (0, "\nTo convert into the base of computation (%u), repeatedly undivide by the"
	"\noriginal base (%u), taking one more base-%u digit as a remainder each time.\n\n",
	base, fbase, fbase);
  for (i=0; i < len[0]; i++)
    for (dig=arg[0][i], j=dw; j--; dig /= base)
      board[1+i*dw+j] = dig % base;
  show (0, 0);
  if (slider[0])
    bprintf (0, "\n%u is an exact power of %u, so we're done.\n\n", fbase, base);
  else {
    memcpy (arg[1], slider+1, len[1] = dw);
    noob_table (0);
    for (qstop = dw; qstop < size; qstop += dw)
      for (qpos = qstop; qpos; qpos--) {
	if (!(dig = board[qpos])) continue;
	board[qpos] = 0;
	add_num (board, qpos, 1, len[1]+1, multab[dig]);
	bprintf (0, "%*%Q%n\n", (qpos-1)*dwide+1, dw, arg[1]);
	if (qpos > 1) show (1, qpos-1);
	else show (0, 0);
      }
  }
}

#define NOPS (sizeof opcode / sizeof *opcode)

int main()
{
  char yn[8], *cp;
// # = no arg, n = numeric arg, a = ASCII arg, capitalized if mandatory
// 0 = leave as is, 1-3 = strip zeroes and pad to a multiple of that many digits
  char opcode[][12] = {
    "###0q", "###0quit", "###0help",
    "NN#0add", "NN#0sub", "N##0neg", "NN#1mul", "NN#1div", "NNn1undiv",
    "N##1square", "N##1cube", "N##2sqrt", "Na#3cbrt",
    "###1perm", "###1perd", "###1pers", "An#1rand", "A##1rand2", "A##1rand3",
    "A##0base", "AN#1tbase", "AA#1fbase",
  };
  void (*func[])() = {
    0, 0, user_help, user_add, user_sub, user_neg, user_mul,
    user_div, user_undiv, user_square, user_cube, user_sqrt, user_cbrt,
    user_perm, user_perd, user_pers, user_rand, user_rand2, user_rand3,
    user_base, user_tbase, user_fbase,
  };
  int op, a;

  srandom(time(0));
  bprintf (0, "Have you memorized the multiplication table up to %i × %i? ", base-1, base-1);
  cp = fgets (yn, 6, stdin);
  if (tolower(yn[0]) == 'y') noob = 0;
  for (;;) {
    if (setjmp (error)) goto next;
#ifdef READLINE
    char *cmd;
    if (cmd = readline (prompt))
      add_history (cmd);
    else {
#else
    char cmd[SIZE*4];
    bprintf (0, prompt);
    fgets (cmd, SIZE*4-1, stdin);
    if (feof(stdin)) {
#endif
      puts ("^D");
      break;
    }
    if (!cmd[0]) continue;
    cp = strtok (cmd, " \t\n");
    for (op = cp ? 0:NOPS; op < NOPS; op++)
      if (!strcmp (cp, opcode[op]+4)) break;
    if (op == NOPS) {
      puts ("Unknown command, try \"help\"");
      continue;
    }
    if (!func[op]) break;
    for (a=0; isalpha(opcode[op][a]); a++) {
      if (!(cp = strtok (0, " \t\n")))
	if (isupper(opcode[op][a])) {
	  bprintf (0, "Missing argument to %s!\n", opcode[op]+4);
	  goto next;
	} else if (opcode[op][a] == 'n') {
	  arg[a][0] = 0;
	  len[a] = 1;
	  break;
	}
      if (toupper(opcode[op][a]) == 'A') {
	strcpy (arg[a], cp ? cp:"0");
	continue;
      }
      len[a] = input_num (cp, arg[a], base, colons, opcode[op][3]-'0');
    }
    explain[0] = 0;
    memset (desc, 0, SIZE);
    memset (board, 0, SIZE);
    (*func[op])();
    next: ;
#ifdef READLINE
    free (cmd);
#endif
  }
}