/* 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 multiplication.) 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 #include #include #include #include #include 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; isize; i++) h->table[i] = (digit *) NULL; } void digithash_free (digithash *h) { unsigned int i; if (!h->table) return; for (i=0; isize; 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; itable[i] = (digit *) NULL; } h->size = newsize; for (i=0; ifull = (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 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; klen && 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-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=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=0; i--) { for (j=0; j nmatch[j+1]) { SWAP (h.numorder[j], h.numorder[j+1], t); SWAP (nmatch[j], nmatch[j+1], t); } } } for (i=0; inum = 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; idismal_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; }