/*

              PROGRAM FOR DISMAL ARITHMETIC

  David Applegate (david(AT)research.att.com), Nov 11 2003

 *************************************************************************

DESCRIPTION

Dismal arithmetic was invented by Marc LeBrun in 2003.
The following is a C program for doing dismal arithmetic
and dismal number theory.

The operations of dismal addition and dismal multiplication are
described in sequences A087061 (addition) and A087062 (multiplication)
in the On-Line Encyclopedia of Integer Sequences.

(A087062 also has Maple code for addition and multiplcation.)

Once this C program is compiled, you can type

dismal help

and it will reply

Usage: dismal cmd [args...]
available cmds are:
   div a b           Outputs c such that b*c = a
   divs a            Outputs b such that b*c = a for some c
   add a b           Outputs a+b
   mul a b           Outputs a*b
   pcount a b        Counts the primes in [a,b]
   dinfo a           Outputs info about divisors of a
   dinfo a b         Outputs info about divisors of n, a<=n<=b
   pdinfo a          Outputs info about prime divisors of a
   pdinfo a b        Outputs info about prime divisors of n, a<=n<=b
   help_dinfo        Describes the output of dinfo and pdinfo
   help              Outputs this text

For example,

dismal divs 10

will produce all dismal divisors of 10 in some random order:

1
10
20
30
40
50
60
70
80
90
2
3
4
5
6
7
8
9

(Note that 10 has just one dismal prime divisor, 90 !)

**************************************************************************


INSTALLATION:

call this file  dismal_eis.c  and compile it, for example using

   cc -o dismal_eis dismal_eis.c


**************************************************************************

                                           */


#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <malloc.h>
#include <string.h>
#include <strings.h>

typedef char digit;

#define SWAP(a,b,t) ((t)=(a),(a)=(b),(b)=(t))

char *my_strdup (char *s)
{
  char *p = (char *) malloc ((strlen(s)+1) * sizeof (char));
  
  if (p == (char *) NULL) {
    fprintf (stderr, "Out of memory\n");
    exit (1);
  }
  strcpy (p, s);
  return p;
}

int isprime (unsigned int p)
{
  unsigned int i;
  
  if ((p&1) == 0) return 0;
  for (i=3; i*i<=p; i+=2) {
    if (p%i == 0) return 0;
  }
  return 1;
}

unsigned int nextprime (unsigned int x)
{
  if (x < 3) return 3;
  x |= 1;
  while (!isprime (x)) x += 2;
  return x;
}

typedef struct digithash {
  unsigned int size;
  unsigned int nelems;
  unsigned int full;
  digit **table;
} digithash;

#define DIGITHASH_INITSIZE 1019
#define DIGITHASH_HIDENS 0.75
#define DIGITHASH_LODENS 0.35

void digithash_init (digithash *h)
{
  unsigned int i;
  
  h->size = nextprime (DIGITHASH_INITSIZE);
  h->nelems = 0;
  h->full = (int) (DIGITHASH_HIDENS * h->size);
  h->table = (digit **) malloc (h->size * sizeof (digit *));
  if (!h->table) {
    fprintf (stderr, "Out of memory\n");
    exit (1);
  }
  for (i=0; i<h->size; i++) h->table[i] = (digit *) NULL;
}

void digithash_free (digithash *h)
{
  unsigned int i;
  
  if (!h->table) return;
  
  for (i=0; i<h->size; i++) {
    if (h->table[i]) free (h->table[i]);
  }
  free (h->table);
}

unsigned int digithash_hash (digit *n)
{
  unsigned int v = 0;
  while (*n) v = v * 11 + *n++;
  return v;
}

int digithash_add_work (digithash *h, digit *n, int add)
{
  unsigned int loc = digithash_hash (n) % h->size;

  while (h->table[loc]) {
    if (!strcmp (h->table[loc], n)) return 0;
    loc++;
    if (loc >= h->size) loc = 0;
  }
  if (add) {
    h->table[loc] = n;
    h->nelems++;
  }
  return 1;
}

void digithash_resize (digithash *h)
{
  unsigned int newsize = nextprime ((unsigned int) (h->nelems /
						    DIGITHASH_LODENS) + 1);
  unsigned int oldsize = h->size;
  digit **oldtable = h->table;
  unsigned int i;

  h->table = (digit **) malloc (newsize * sizeof (digit *));
  if (!h->table) {
    fprintf (stderr, "Out of memory\n");
    exit (1);
  }

  for (i=0; i<newsize; i++) {
    h->table[i] = (digit *) NULL;
  }
  h->size = newsize;

  for (i=0; i<oldsize; i++) {
    if (oldtable[i]) {
      digithash_add_work (h, oldtable[i], 1);
    }
  }

  free (oldtable);

  h->full = (int) (DIGITHASH_HIDENS * h->size);
}
  
int digithash_add (digithash *h, digit *n)
{
  digit *ncopy = my_strdup (n);
  
  if (h->nelems >= h->full) {
    digithash_resize (h);
  }

  if (digithash_add_work (h, ncopy, 1)) {
    return 1;
  } else {
    free (ncopy);
    return 0;
  }
}

int digithash_find (digithash *h, digit *n)
{
  return !digithash_add_work (h, n, 0);
}

int nfound = 0;
int ndeadend = 0;
int ndup = 0;

typedef struct dig_range {
  digit min;
  digit max;
} dig_range;

typedef struct div_num {
  int len;
  dig_range *digs;
} div_num;

typedef struct divide_help {
  div_num a;
  int numlen;
  int *numorder;
  digit *res;
  digithash reshash;
} divide_help;

typedef struct stack_elem {
  int idx;
  digit d;
} stack_elem;

typedef struct stack {
  int top;
  int maxtop;
  stack_elem *ents;
} stack;

#define PUSH(s,id,dig) { \
    if ((s)->top >= (s)->maxtop) { \
      fprintf (stderr, "Out of stack space\n"); \
      exit (1); \
    } \
    (s)->ents[(s)->top].idx = (id); \
    (s)->ents[(s)->top].d = (dig); \
    (s)->top++; \
}

#define POP(s,id,dig) { \
    if ((s)->top <= 0) { \
      fprintf (stderr, "stack underflow\n"); \
      exit (1); \
    } \
    (s)->top--; \
    (id) = (s)->ents[(s)->top].idx; \
    (dig) = (s)->ents[(s)->top].d; \
}
     
typedef struct divisor_help {
  div_num a;
  div_num b;
  int numlen;
  int *numorder;
  digit *res;
  stack amaxstack;
  stack bmaxstack;
  digithash reshash;
} divisor_help;

typedef int (div_callback) (digit *, void *);

void create_stack (stack *s, int size)
{
  s->ents = (stack_elem *) malloc (size * sizeof (stack_elem));
  if (!s->ents) {
    fprintf (stderr, "Out of memory\n");
    exit (1);
  }
  s->maxtop = size;
  s->top = 0;
}

void free_stack (stack *s)
{
  free (s->ents);
  s->maxtop = 0;
  s->top = 0;
}

void dismal_mul (digit *a, digit *b, digit *res)
{
  int reslen = strlen(a) + strlen(b) - 1;
  int i;
  int j;

  for (i=0; i<reslen; i++) res[i] = '0';
  res[reslen] = '\0';

  for (i=0; a[i]; i++) {
    for (j=0; b[j]; j++) {
      if (a[i] <= b[j] && a[i] > res[i+j]) res[i+j] = a[i];
      if (b[j] < a[i] && b[j] > res[i+j]) res[i+j] = b[j];
    }
  }
}

void dismal_add (digit *a, digit *b, digit *res)
{
  int alen = strlen(a);
  int blen = strlen(b);

  while (*a && alen > blen) {
    *res++ = *a++;
    alen--;
  }
  while (*b && blen > alen) {
    *res++ = *b++;
    blen--;
  }

  while (*a) {
    if (*a > *b) *res = *a;
    else         *res = *b;
    a++;
    b++;
    res++;
  }
  *res = '\0';
}

int spin (digit *res, div_num *n, digithash *reshash, int depth,
		 div_callback *callback, void *u_data)
{
  digit d;
  int nsol = 0;
  
  if (depth >= n->len) {
    res[depth] = '\0';
    while (res[0] == '0') res++;
    if (digithash_add (reshash, res)) {
      nfound++;
      nsol = (*callback)(res, u_data);
      return nsol;
    } else {
      ndup++;
      return 0;
    }
  }

  for (d = n->digs[depth].min; d <= n->digs[depth].max; d++) {
    res[depth] = d;
    nsol += spin (res, n, reshash, depth+1, callback, u_data);
  }
  return nsol;
}

int divide_work (digit *num, digit *den, divide_help *h,
			int depth, div_callback *callback, void *u_data)
{
  int i;
  int j;
  int alen = h->a.len;
  int omin;
  int omax;
  int nsol = 0;
  int rval;
  int ncall = 0;

  if (depth >= h->numlen) {
    nsol = spin (h->res, &h->a, &h->reshash, 0, callback, u_data);
    return nsol;
  }

  i = h->numorder[depth];

  for (j=(i>=alen?i-alen+1:0); j<=i && den[j]; j++) {
    if (den[j] > num[i] && h->a.digs[i-j].min <= num[i] &&
	h->a.digs[i-j].max >= num[i]) {
      omin = h->a.digs[i-j].min;
      omax = h->a.digs[i-j].max;
      h->a.digs[i-j].min = num[i];
      h->a.digs[i-j].max = num[i];
      ncall++;
      rval = divide_work (num, den, h, depth+1, callback, u_data);
      if (rval < 0) return rval;
      nsol += rval;
      h->a.digs[i-j].min = omin;
      h->a.digs[i-j].max = omax;
    } else if (den[j] == num[i] && h->a.digs[i-j].max >= num[i]) {
      omin = h->a.digs[i-j].min;
      if (omin < num[i]) h->a.digs[i-j].min = num[i];
      ncall++;
      rval = divide_work (num, den, h, depth+1, callback, u_data);
      if (rval < 0) return rval;
      nsol += rval;
      h->a.digs[i-j].min = omin;
    }
  }

  if (ncall == 0) ndeadend++;
  return nsol;
}

void min_change (digit *num, int numlen, div_num *a, int j,
			div_num *b, stack *s)
{
  int k;
  digit m = a->digs[j].min;

  for (k=0; k<b->len && j+k < numlen; k++) {
    if (m > num[j+k] && b->digs[k].max > num[j+k]) {
      PUSH (s, k, b->digs[k].max);
      b->digs[k].max = num[j+k];
    }
  }
}


void unroll_stack (stack *s, int loc, div_num *a)
{
  int i;
  digit d;
  
  while (s->top > loc) {
    POP (s, i, d);
    a->digs[i].max = d;
  }
}

int divisor_work (digit *num, divisor_help *h, int depth,
			 div_callback *callback, void *u_data)
{
  int rval;
  int rval2;
  int i;
  int j;
  int nsol = 0;
  int ncall = 0;
  int omina;
  int omaxa;
  int ominb;
  int omaxb;
  int astackloc;
  int bstackloc;
  
  if (depth >= h->numlen) {
    rval = spin (h->res, &h->a, &h->reshash, 0, callback, u_data);
    if (rval < 0) return rval;
    rval2 = spin (h->res, &h->b, &h->reshash, 0, callback, u_data);
    if (rval2 < 0) return rval2;
    return rval + rval2;
  }

  i = h->numorder[depth];
  
  for (j=(i>=h->a.len?i-h->a.len+1:0); j<=i && j < h->b.len; j++) {
    if (h->a.digs[i-j].max < num[i] || h->b.digs[j].max < num[i]) continue;
    if (h->a.digs[i-j].min <= num[i]) {
      omina = h->a.digs[i-j].min;
      omaxa = h->a.digs[i-j].max;
      ominb = h->b.digs[j].min;
      h->a.digs[i-j].min = num[i];
      h->a.digs[i-j].max = num[i];
      astackloc = h->amaxstack.top;
      bstackloc = h->bmaxstack.top;
      min_change (num, h->numlen, &h->a, i-j, &h->b, &h->bmaxstack);
      if (h->b.digs[j].min < num[i]) {
	h->b.digs[j].min = num[i];
	min_change (num, h->numlen, &h->b, j, &h->a, &h->amaxstack);
      }
      ncall++;
      rval = divisor_work (num, h, depth+1, callback, u_data);
      if (rval < 0) return rval;
      nsol += rval;
      unroll_stack (&h->amaxstack, astackloc, &h->a);
      h->b.digs[j].min = ominb;
      unroll_stack (&h->bmaxstack, bstackloc, &h->b);
      h->a.digs[i-j].min = omina;
      h->a.digs[i-j].max = omaxa;
    }
    if (h->b.digs[j].min <= num[i]) {
      ominb = h->b.digs[j].min;
      omaxb = h->b.digs[j].max;
      omina = h->a.digs[i-j].min;
      h->b.digs[j].min = num[i];
      h->b.digs[j].max = num[i];
      astackloc = h->amaxstack.top;
      bstackloc = h->bmaxstack.top;
      min_change (num, h->numlen, &h->b, j, &h->a, &h->amaxstack);
      if (h->a.digs[i-j].min < num[i]) {
	h->a.digs[i-j].min = num[i];
	min_change (num, h->numlen, &h->a, i-j, &h->b, &h->bmaxstack);
      }
      ncall++;
      rval = divisor_work (num, h, depth+1, callback, u_data);
      if (rval < 0) return rval;
      nsol += rval;
      unroll_stack (&h->bmaxstack, bstackloc, &h->b);
      h->a.digs[i-j].min = omina;
      unroll_stack (&h->amaxstack, astackloc, &h->a);
      h->b.digs[j].min = ominb;
      h->b.digs[j].max = omaxb;
    }
  }

  if (ncall == 0) ndeadend++;

  return nsol;
}

int dismal_divide (digit *num, digit *den, div_callback *callback,
		   void *u_data)
{
  int numlen = strlen(num);
  int reslen = numlen + 1 - strlen(den);
  divide_help h;
  int i;
  int j;
  int *nmatch;
  int *nloc;
  int t;
  int rval;

  if (reslen <= 0) {
    fprintf (stderr, "result length %d -> no solutions\n", reslen);
    return 0;
  }

  h.res = (digit *) malloc ((reslen + 1) * sizeof (digit));
  if (!h.res) {
    fprintf (stderr, "Out of memory\n");
    exit (1);
  }
  for (i=0; i<=reslen; i++) {
    h.res[i] = '\0';
  }
  
  h.a.len = reslen;
  h.numlen = numlen;
  h.a.digs = (dig_range *) malloc (reslen * sizeof (dig_range));
  h.numorder = (int *) malloc (numlen * sizeof (int));
  nmatch = (int *) malloc (numlen * sizeof (int));
  nloc = (int *) malloc (numlen * sizeof (int));
  if (!h.a.digs || !h.numorder || !nmatch || !nloc) {
    fprintf (stderr, "out of memory\n");
    exit (1);
  }

  for (i=0; i<reslen; i++) {
    h.a.digs[i].min = '0';
    h.a.digs[i].max = '9';
  }
  for (i=0; i<numlen; i++) {
    h.numorder[i] = i;
    nmatch[i] = 0;
    nloc[i] = i;
  }

  for (i=0; i<numlen; i++) {
    for (j=(i>=reslen?i-reslen+1:0); j<=i && den[j]; j++) {
      if (den[j] > num[i] && h.a.digs[i-j].max > num[i]) {
	h.a.digs[i-j].max = num[i];
      }
    }
  }

  for (i=0; i<numlen; i++) {
    for (j=(i>=reslen?i-reslen+1:0); j<=i && den[j]; j++) {
      if (den[j] >= num[i] && h.a.digs[i-j].max >= num[i]) {
	nmatch[i]++;
      }
    }
  }

#if 0
  printf ("result digit max: ");
  for (i=0; i<reslen; i++) {
    printf ("%c", h.a.digs[i].max);
  }
  printf ("\n");

  printf ("num nmatch:");
  for (i=0; i<numlen; i++) {
    printf (" %d", nmatch[i]);
  }
  printf ("\n");
#endif
  
  for (i=0; i<numlen; i++) {
    if (nmatch[i] == 0) {
      free (nloc);
      free (nmatch);
      free (h.res);
      free (h.a.digs);
      free (h.numorder);
      return 0;
    }
  }

  for (i=numlen-1; i>=0; i--) {
    for (j=0; j<i; j++) {
      if (nmatch[j] > nmatch[j+1]) {
	SWAP (h.numorder[j], h.numorder[j+1], t);
	SWAP (nmatch[j], nmatch[j+1], t);
      }
    }
  }

  for (i=0; i<numlen; i++) {
    h.numorder[nloc[i]] = i;
  }
  
#if 0
  printf ("num assign order:");
  for (i=0; i<numlen; i++) {
    printf (" %d", h.numorder[i]);
  }
  printf ("\n");
#endif

  free (nloc);
  free (nmatch);

  digithash_init (&h.reshash);
  
  rval = divide_work (num, den, &h, 0, callback, u_data);

  digithash_free (&h.reshash);

  free (h.res);
  free (h.a.digs);
  free (h.numorder);

  return rval;
}

int dismal_divisors (digit *num, div_callback *callback, void *u_data)
{
  int numlen = strlen(num);
  divisor_help h;
  int i;
  int j;
  int rval;
  int nsol = 0;

  h.numlen = numlen;
  h.a.digs = (dig_range *) malloc (numlen * sizeof (dig_range));
  h.b.digs = (dig_range *) malloc (numlen * sizeof (dig_range));
  h.numorder = (int *) malloc ((numlen+1) * sizeof (int));
  h.res = (digit *) malloc ((numlen+1) * sizeof (digit));

  if (!h.a.digs || !h.numorder || !h.b.digs || !h.res) {
    fprintf (stderr, "out of memory\n");
    exit (1);
  }

  create_stack (&h.amaxstack, 10*numlen);
  create_stack (&h.bmaxstack, 10*numlen);
  
  for (i=0; i<numlen; i++) {
    h.a.digs[i].min = '0';
    h.a.digs[i].max = '9';
    h.b.digs[i].min = '0';
    h.b.digs[i].max = '9';
  }
  j=0;
  for (i=0; 2*i<numlen; i++) {
    h.numorder[j++] = i;
    h.numorder[j++] = numlen - 1 - i;
  }    

  digithash_init (&h.reshash);

  for (i=0; 2*i<numlen; i++) {
    h.a.len = i+1;
    h.b.len = numlen - i;
    rval = divisor_work (num, &h, 0, callback, u_data);
    if (rval < 0) return rval;
    nsol += rval;
  }

  digithash_free (&h.reshash);

  free_stack (&h.amaxstack);
  free_stack (&h.bmaxstack);
  
  free (h.res);
  free (h.b.digs);
  free (h.a.digs);
  free (h.numorder);

  return nsol;
}

div_callback div_print;

#define DIV_PROD_MUL 16

int div_print (digit *res, void *u_data)
{
  printf ("%s\n", res);
  return 0;
}

void usage (char *fname)
{
  fprintf (stderr, "Usage: %s cmd [args...]\n", fname);
  fprintf (stderr, "available cmds are:\n");
  fprintf (stderr, "   div a b           Outputs c such that b*c = a\n");
  fprintf (stderr, "   divs a            Outputs b such that b*c = a for some c\n");
  fprintf (stderr, "   add a b           Outputs a+b\n");
  fprintf (stderr, "   mul a b           Outputs a*b\n");
  fprintf (stderr, "   pcount a b        Counts the primes in [a,b]\n");
  fprintf (stderr, "   dinfo a           Outputs info about divisors of a\n");
  fprintf (stderr, "   dinfo a b         Outputs info about divisors of n, a<=n<=b\n");
  fprintf (stderr, "   pdinfo a          Outputs info about prime divisors of a\n");
  fprintf (stderr, "   pdinfo a b        Outputs info about prime divisors of n, a<=n<=b\n");
  fprintf (stderr, "   help_dinfo        Describes the output of dinfo and pdinfo\n");
  fprintf (stderr, "   help              Outputs this text\n");
  exit (1);
}

digit *result_space (digit *a, digit *b)
{
  digit *s = (digit *) malloc ((strlen (a) + strlen (b) + 1) * sizeof (digit));
  
  if (!s) {
      fprintf (stderr, "Out of memory\n");
      exit (1);
  }
  return s;
}

typedef struct div_inf {
  digit *num;
  int len;
  int is_prime;
  int is_pseudoprime;
  int ndivisors;
  int nbounded_divisors;
  digit *dismal_sum_divisors;
  digit *dismal_sum_bounded_divisors;
  digit *dismal_sum_bounded_divisors2;
  digit *dismal_sum_ne_divisors;
  digit *sum_work;
  double sum_divisors;
  double sum_bounded_divisors;
  double sum_bounded_divisors2;
  double sum_ne_divisors;
  int nprime_divisors;
  digit *dismal_sum_prime_divisors;
  digit *dismal_prod_prime_divisors;
  digithash *primehash;
} div_inf;

void init_div_inf (div_inf *di, digit *a)
{
  di->num = a;
  di->len = strlen (a);
  if (!strcmp (a, "8") || !strcmp (a, "9")) {
    di->is_prime = 0;
  } else {
    di->is_prime = 1;
  }
  di->is_pseudoprime = 1;
  di->ndivisors = 0;
  di->nbounded_divisors = 0;
  di->sum_divisors = 0.0;
  di->sum_bounded_divisors = 0.0;
  di->sum_bounded_divisors2 = 0.0;
  di->sum_ne_divisors = 0.0;
  di->dismal_sum_divisors = (digit *) malloc ((di->len+1) * sizeof (digit));
  di->dismal_sum_bounded_divisors = (digit *) malloc ((di->len+1) * sizeof (digit));
  di->dismal_sum_bounded_divisors2 = (digit *) malloc ((di->len+1) * sizeof (digit));
  di->dismal_sum_ne_divisors = (digit *) malloc ((di->len+1) * sizeof (digit));
  di->nprime_divisors = 0;
  di->dismal_sum_prime_divisors = (digit *) malloc ((di->len+1) * sizeof (digit));
  di->dismal_prod_prime_divisors = (digit *) malloc ((DIV_PROD_MUL*di->len+1) * sizeof (digit));
  di->sum_work = (digit *) malloc ((DIV_PROD_MUL*di->len+1) * sizeof (digit));
  if (!di->dismal_sum_divisors || !di->dismal_sum_bounded_divisors ||
      !di->dismal_sum_bounded_divisors2 || !di->dismal_sum_ne_divisors ||
      !di->dismal_sum_prime_divisors || !di->dismal_prod_prime_divisors ||
      !di->sum_work) {
    fprintf (stderr, "Out of memory\n");
    exit (1);
  }

  di->dismal_sum_divisors[0] = '0';
  di->dismal_sum_bounded_divisors[0] = '0';
  di->dismal_sum_bounded_divisors2[0] = '0';
  di->dismal_sum_ne_divisors[0] = '0';
  di->dismal_sum_prime_divisors[0] = '0';
  di->dismal_prod_prime_divisors[0] = '9';
  di->sum_work[0] = '0';
  di->dismal_sum_divisors[1] = '\0';
  di->dismal_sum_bounded_divisors[1] = '\0';
  di->dismal_sum_bounded_divisors2[1] = '\0';
  di->dismal_sum_ne_divisors[1] = '\0';
  di->dismal_sum_prime_divisors[1] = '\0';
  di->dismal_prod_prime_divisors[1] = '\0';
  di->sum_work[1] = '\0';
}

void free_div_inf (div_inf *di)
{
  free (di->dismal_sum_divisors);
  free (di->dismal_sum_bounded_divisors);
  free (di->dismal_sum_bounded_divisors2);
  free (di->dismal_sum_ne_divisors);
  free (di->dismal_sum_prime_divisors);
  free (di->dismal_prod_prime_divisors);
  free (di->sum_work);
}

div_callback gather_inf;

int gather_inf (digit *res, void *u_data)
{
  div_inf *di = (div_inf *) u_data;
  unsigned int reslen = strlen (res);
  unsigned int numlen = strlen (di->num);
  int numcmp = strcmp (res, di->num);
  double resf = atof (res);
  unsigned int i;

  if (numcmp && strcmp (res, "9")) {
    di->is_prime = 0;
  }
  if (reslen > 1 && reslen < numlen) di->is_pseudoprime = 0;
  di->ndivisors++;
  dismal_add (res, di->dismal_sum_divisors, di->sum_work);
  strcpy (di->dismal_sum_divisors, di->sum_work);
  di->sum_divisors += resf;
  if (reslen < numlen || numcmp <= 0) {
    di->nbounded_divisors++;
    dismal_add (res, di->dismal_sum_bounded_divisors, di->sum_work);
    strcpy (di->dismal_sum_bounded_divisors, di->sum_work);
    di->sum_bounded_divisors += resf;
    if (reslen < numlen || numcmp < 0) {
      dismal_add (res, di->dismal_sum_bounded_divisors2, di->sum_work);
      strcpy (di->dismal_sum_bounded_divisors2, di->sum_work);
      di->sum_bounded_divisors2 += resf;
    }
  }
  if (numcmp != 0) {
    dismal_add (res, di->dismal_sum_ne_divisors, di->sum_work);
    strcpy (di->dismal_sum_ne_divisors, di->sum_work);
    di->sum_ne_divisors += resf;
  }

  if (digithash_find (di->primehash, res)) {
    di->nprime_divisors++;
    dismal_add (res, di->dismal_sum_prime_divisors, di->sum_work);
    strcpy (di->dismal_sum_prime_divisors, di->sum_work);
    if (reslen + strlen (di->dismal_prod_prime_divisors) <= DIV_PROD_MUL*numlen) {
      dismal_mul (res, di->dismal_prod_prime_divisors, di->sum_work);
      strcpy (di->dismal_prod_prime_divisors, di->sum_work);
    } else {
      for (i=0; i<DIV_PROD_MUL*numlen; i++) di->dismal_prod_prime_divisors[i] = '9';
      di->dismal_prod_prime_divisors[DIV_PROD_MUL*numlen] = '\0';
    }
  }

  return 0;
}

div_callback proper_divisor;

int proper_divisor (digit *res, void *u_data)
{
  digit *num = (digit *) u_data;

  if (!strcmp (res, num) || !strcmp (res, "9")) {
    return 0;
  } else {
    return 1;
  }
}

void incr_digit_num (digit *n)
{
  int i = strlen (n) - 1;

  while (i >= 0) {
    if (n[i] != '9') {
      n[i]++;
      return;
    } else {
      n[i] = '0';
      i--;
    }
  }

  n[0] = '1';
  i = strlen (n);
  n[i] = '0';
  n[i+1] = '\0';
  return;
}

void help_dinfo (void)
{
  printf ("dinfo and pdinfo output a header row, and then a row for each\n");
  printf ("number requested.  Each row consists of '|' separated fields.\n");
  printf ("\n");
  printf ("The fields are:\n");
  printf ("  n                   the number\n");
  printf ("  9ish                1 if n contains a 9, 0 otherwise\n");
  printf ("  prime               1 if n is prime, 0 otherwise\n");
  printf ("                         n is prime if n != 9, and if n=a*b, \n");
  printf ("                         then either a==9 or b==9\n");
  printf ("  pseudoprime         1 if n has a divisor with length in [2,n)\n");
  printf ("  divisors            the number of divisors of n\n");
  printf ("  bd_divisors         the number of divisors of n in [1,n]\n");
  printf ("  dsum_divisors       the dismal sum of the divisors of n\n");
  printf ("  dsum_bd_divisors    the dismal sum of the divisors of n in [1,n]\n");
  printf ("  dsum_bd_divisors2   the dismal sum of the divisors of n in [1,n)\n");
  printf ("  dsum_ne_divisors    the dismal sum of the divisors of n that are != n\n");
  printf ("  2*n                 the dismal product 2*n\n");
  printf ("  sum_divisors        the normal sum of the divisors of n\n");
  printf ("  sum_bd_divisors     the normal sum of the divisors of n in [1,n]\n");
  printf ("  sum_bd_divisors2    the normal sum of the divisors of n in [1,n)\n");
  printf ("  sum_ne_divisors     the normal sum of the divisors of n that are != n\n");
  printf ("\n");
  printf ("In addition, pdinfo includes the following fields about prime divisors\n");
  printf ("  prime_divisors      the number of prime divisors of n\n");
  printf ("  sum_prime_divisors  the dismal sum of the prime divisors of n\n");
  printf ("  prod_prime_divisors the dismal product of the prime divisors of n\n");
  printf ("                        because numbers like 111...111 have very many\n");
  printf ("                        prime divisors, this product is capped at\n");
  printf ("                        length %d*length(n)\n", DIV_PROD_MUL);
}

void collect_primes (digithash *primehash, digit *hi, digit *x)
{
  unsigned int hilen = strlen (hi);

  if (!x) {
    fprintf (stderr, "Out of memory\n");
    exit (1);
  }
  x[0] = '1';
  x[1] = '\0';

  while (strlen (x) <= hilen) {
    if (strcmp (x, "8") && strcmp (x, "9") &&
	dismal_divisors (x, proper_divisor, x) == 0) {
      digithash_add (primehash, x);
    }
    incr_digit_num (x);
  }
}

void dinfo_range (digit *lo, digit *hi, int primefields)
{
  digit *x = (digit *) malloc ((strlen (hi)+2)*sizeof (digit));
  digit *x2 = (digit *) malloc ((strlen (hi)+2)*sizeof (digit));
  div_inf di;
  digithash primehash;

  if (!x || !x2) {
    fprintf (stderr, "Out of memory\n");
    exit (1);
  }
  digithash_init (&primehash);

  if (primefields) {
    collect_primes (&primehash, hi, x);
  }

  strcpy (x, lo);

  printf ("n|9ish|prime|pseudoprime|divisors|bd_divisors|dsum_divisors|dsum_bd_divisors|dsum_bd_divisors2|dsum_ne_divisors|2*n|sum_divisors|sum_bd_divisors|sum_bd_divisors2|sum_ne_divisors");
  if (primefields) printf ("|prime_divisors|sum_prime_divisors|prod_prime_divisors");
  printf ("\n");

  while (strlen (x) < strlen (hi) ||
	 (strlen (x) == strlen (hi) && strcmp (x, hi) <= 0)) {
    init_div_inf (&di, x);
    di.primehash = &primehash;
    
    dismal_divisors (x, gather_inf, &di);
    dismal_mul (x, "2", x2);

    printf ("%s|%d|%d|%d|%d|%d|%s|%s|%s|%s|%s|%.0f|%.0f|%.0f|%.0f",
	    x, index(x, '9') != NULL, di.is_prime, di.is_pseudoprime,
	    di.ndivisors, di.nbounded_divisors, di.dismal_sum_divisors,
	    di.dismal_sum_bounded_divisors, di.dismal_sum_bounded_divisors2,
	    di.dismal_sum_ne_divisors, x2, di.sum_divisors,
	    di.sum_bounded_divisors, di.sum_bounded_divisors2,
	    di.sum_ne_divisors);
    if (primefields) {
      printf ("|%d|%s|%s", di.nprime_divisors, di.dismal_sum_prime_divisors,
	      di.dismal_prod_prime_divisors);
    }
    printf ("\n");

    free_div_inf (&di);
    incr_digit_num (x);
  }

  digithash_free (&primehash);
  free (x2);
  free (x);
}

void pcount_range (digit *lo, digit *hi)
{
  digit *x = (digit *) malloc ((strlen (hi)+2)*sizeof (digit));
  int pcount = 0;
  unsigned int output_length = strlen (lo);

  if (!x) {
    fprintf (stderr, "Out of memory\n");
    exit (1);
  }
  
  strcpy (x, lo);

  while (strlen (x) < strlen (hi) ||
	 (strlen (x) == strlen (hi) && strcmp (x, hi) <= 0)) {
    if (strcmp (x, "8") && strcmp (x, "9") &&
	dismal_divisors (x, proper_divisor, x) == 0) {
      pcount++;
    }
    incr_digit_num (x);
    if (strlen (x) > output_length) {
      printf ("%d primes with <= %d digits\n", pcount, output_length);
      fflush (stdout);
      output_length++;
    }
  }
  printf ("%d primes in [%s,%s]\n", pcount, lo, hi);

  free (x);
}

int main (int ac, char **av)
{
  digit *res;
  
  if (ac < 2) usage (av[0]);
  if (!strcmp (av[1], "div") && ac == 4) {
    dismal_divide (av[2], av[3], div_print, (void *) NULL);
  } else if (!strcmp (av[1], "divs") && ac == 3) {
    dismal_divisors (av[2], div_print, (void *) NULL);
  } else if (!strcmp (av[1], "add") && ac == 4) {
    res = result_space (av[2], av[3]);
    dismal_add (av[2], av[3], res);
    div_print (res, (void *) NULL);
    free (res);
  } else if (!strcmp (av[1], "mul") && ac == 4) {
    res = result_space (av[2], av[3]);
    dismal_mul (av[2], av[3], res);
    div_print (res, (void *) NULL);
    free (res);
  } else if (!strcmp (av[1], "dinfo") && ac == 3) {
    dinfo_range (av[2], av[2], 0);
  } else if (!strcmp (av[1], "dinfo") && ac == 4) {
    dinfo_range (av[2], av[3], 0);
  } else if (!strcmp (av[1], "pdinfo") && ac == 3) {
    dinfo_range (av[2], av[2], 1);
  } else if (!strcmp (av[1], "pdinfo") && ac == 4) {
    dinfo_range (av[2], av[3], 1);
  } else if (!strcmp (av[1], "pcount") && ac == 4) {
    pcount_range (av[2], av[3]);
  } else if (!strcmp (av[1], "help_dinfo")) {
    help_dinfo ();
  } else {
    usage (av[0]);
  }

  return 0;
}