@*Intro. This program finds all ``perfect digital invariants'' of order~$m$
in the decimal system, namely all integers that satisfy $\pi_m x=x$,
where $\pi_m$ takes an integer into the sum of the $m$th powers of its digits.
It can be shown without difficulty that such integers have at most $m+1$ digits.
Indeed, if $10^p\le x<10^{p+1}$ we have $\pi_m x<10^p$ whenever $p>m$.
(The proof follows from the fact that $(m+1)9^m<10^{m+1}$.)
It's an interesting backtrack program, in which I successively choose the
digits $9\ge x_1\ge x_2\ge\cdots\ge x_{m+1}\ge0$ that will be the
digits of~$x$ (in some order). Lower bounds and upper bounds on $x$
are sufficiently sharp to rule out lots of cases before very many
of those digits have been specified. (And if $m$ is small,
I could even run through {\it all\/} such sequences of digits, because there
are only $m+10\choose 9$ of them. That's about 2.5 billion when $m=40$.)
The only high-precision arithmetic needed here is addition. I implement it
with binary-coded decimal representation (15 digits per octabyte),
using bitwise techniques as suggested in exercise 7.1.3--100.
Memory references (mems) are counted as if an optimizing compiler were
doing things like inlining subroutines, and as if the distribution arrays were
packed into a single octabyte. I actually keep the elements unpacked, to keep
debugging simple.
@d maxm 1000
@d maxdigs (1+(maxm/15)) /* octabytes per binary-coded decimal number */
@d o mems++
@d oo mems+=2
@c
#include
#include
int m; /* command-line parameter */
typedef unsigned long long ull;
ull mems;
ull nodes;
ull thresh=10000000000; /* reporting time */
ull profile[maxm+3];
int count;
int vbose; /* level of verbosity */
@;
@;
main(int argc,char*argv[]) {
register int j,k,l,p,r,t,pd,alt,blt,xl,change;
@;
@;
@;
fprintf(stderr,"Altogether %d solutions for m=%d (%llu nodes, %llu mems).\n",
count,m,nodes,mems);
if (vbose) @;
}
@ @=
{
fprintf(stderr,"Profile: 1\n");
for (k=2;k<=m+2;k++) fprintf(stderr,"%19lld\n",profile[k]);
}
@ @=
if (argc<2 || sscanf(argv[1],"%d",&m)!=1) {
fprintf(stderr,"Usage: %s m [profile] [verbose] [extraverbose]\n",argv[0]);
exit(-1);
}
vbose=argc-2;
if (m<2 || m>maxm) {
fprintf(stderr,"Sorry, m should be between 2 and %d, not %d!\n",maxm,m);
exit(-2);
}
mdigs=1+(m/15);
@*Tricky arithmetic. I've got to deal with biggish numbers and inspect their
decimal digits. But I'm using a binary computer and I don't want to
be repeatedly dividing by powers of~10. So I have an addition routine
that computes (say) the sum of hexadecimal-coded numbers
|0x344159959| and |0x271828043|, giving |0x615988002| as if
the numbers were decimal instead.
@=
void add(ull*p,ull*q,ull*r) { /* add |p| to |q|, giving |r| */
register int k,c;
register ull t,w,x,y;
for (k=c=0;k;
if (c) { /* this shouldn't happen */
fprintf(stderr,"Overflow!\n");
exit(-999);
}
}
@ It's interesting that I must add |c| to |x| here, {\it not\/} to~|y|.
Otherwise the nondecimal digit \.{a} might appear in the result.
@=
{
o,x=*(p+k)+c; /* |x| might have a nondecimal digit now */
o,y=*(q+k)+0x666666666666666; /* no cross-digit carries occur */
t=x+y;
w=(t^x^y)&0x1111111111111110; /* this is where cross-digit carries happen */
w=(w^0x1111111111111110)>>3;
t-=w+(w<<1); /* subtract 6 where there were no carries */
o,*(r+k)=t&0xfffffffffffffff;
c=t>>60;
}
@ At the beginning of this program, I need a table of $0^m$, $1^m$, $2^m$,
\dots,~$9^m$. So why not compute it via addition?
@=
void kmult(int k,ull*a) { /* multiply |a| by $k$ */
switch(k) {
case 8: add(a,a,a);
case 4: add(a,a,a);
case 2: add(a,a,a);@+break;
case 6: add(a,a,a);
case 3: add(a,a,z);@+add(a,z,a);@+break;
case 5: add(a,a,z);@+add(z,z,z);@+add(a,z,a);@+break;
case 9: add(a,a,z);@+add(z,z,z);@+add(z,z,z);@+add(a,z,a);@+break;
case 7: add(a,a,z);@+add(a,z,z);@+add(z,z,z);@+add(a,z,a);@+break;
case 0: case 1: break;
}
}
@ @=
for (k=1;k<10;k++) {
table[1][k][0]=k;
for (j=2;j<=m;j++) kmult(k,table[1][k]); /* compute $k^m$ */
for (j=2;j<=m+1;j++)
add(table[1][k],table[j-1][k],table[j][k]); /* compute $j\cdot k^m$ */
}
@ @=
int mdigs; /* our multiprecision arithmetic routine uses this many octabytes */
ull table[maxm+2][10][maxdigs]; /* precomputed tables of $j\cdot k^m$ */
ull z[maxdigs]; /* temporary buffer for bignums */
@ Here's a macro that delivers a given digit (nybble) of a mutlibyte number.
@d nybb(a,p) (int)((a[p/15]>>(4*(p%15)))&0xf)
@ When debugging, or operating verbosely, I want to see all digits of a
multiprecise number, with a vertical bar just before digit number~|t|.
@=
void printnum(ull*a,int t) {
register int k;
for (k=m;k>=0;k--) {
if (t==k) fprintf(stderr,"|");
fprintf(stderr,"%d",nybb(a,k));
}
}
@*The algorithm. This program has the overall structure of a typical
backtrack program, with a few twists. One of those twists is the
state parameter~|pd|, which is nonzero when the move at level~|l-1| was forced.
(Such cases are rare, but important.)
@=
b1: @;
b2: profile[l]++,nodes++;
@=thresh|@>;
for (k=0;k<10;k++) {
pdist[l][k]=pdist[l-1][k];
dist[l][k]=dist[l-1][k]+(k==xl?1:0);
}
oo,oo; /* two mems to copy |pdist| and |dist|, which could have been packed */
if (pd) @@;
else {
if (r==0) goto b5; /* we haven't room to accept a new digit |xl| */
r--,add(sig[l-1],table[1][xl],sig[l]);
}
if (l>m+1) @;
b3:@+if (vbose>1) fprintf(stderr,"Level %d, trying %d (%lld mems)\n",
l,xl,mems);
@;
move: @;
b4:@+if (xl) {
xl--;
o,pd=pdist[l][xl]; /* |dist[l][xl]| was zero */
goto b3;
}
b5:@+if (--l) {
o,pd=pdsave[l];
if (pd) goto b5;
@;
goto b4;
}
@ @=thresh|@>=
if (mems>=thresh) {
thresh+=10000000000;
fprintf(stderr,"After %lld mems:",mems);
for (k=2;k<=l;k++) fprintf(stderr," %lld",profile[k]);
fprintf(stderr,"\n");
}
@ The purpose of backtrack level |l| is to compute the $l$th largest digit,
$x_l$, of a solution~$x$, assuming that $x_1$, \dots,~$x_{l-1}$ have already
been specified.
The main idea is to compute bounds $a_l$ and $b_l$ such that $a_l\le x\le b_l$
must be valid, whenever $x_1$, \dots,~$x_{l-1}$ have the given values and
$x_l$ is at most a given threshold value~|xl|. Those bounds, like all of the
multiprecise numbers in this computation, are $(m+1)$-digit numbers whose
individual digits are $a_{lm}\ldots a_{l0}$ and $b_{lm}\ldots b_{l0}$.
They share a common prefix $p_m\ldots p_{t+1}$ of length $m+1-t$; thus
if $a_lxl| we must have |pdist[l][d]<=dist[l][d]|.
If |d=xl| we set $|pd|=\max(0,|pdist|[l][d]-|dist|[l][d])$.
Thus, if |xl| occurs thrice in $D$ but only
once in~$P$, we have $pd=0$; but if |xl| occurs thrice in $P$ but only
once in~$D$, we have $pd=2$. In the latter case we must choose $x_l=|xl|$
and also $x_{l+1}=|xl|$.
Let |r| be the number of unknown digits of |x|. (When |pd=0|, this is
$m+1$ minus $\vert S\vert$, the number of known digits.) If
$a_{lt}0$ and that one of the unknown digits lies between
$a_{lt}$ and~$b_{lt}$, inclusive.
When |xl| decreases, the bounds get tighter, hence the prefix can
become longer. And that's good.
These are the key facts governing our bounds $a_l$ and $b_l$. In order
to do the computations conveniently we maintain the sum of known digits,
$|sig|[l]=\sum_{k=0}^{xl-1}|pdist|[l][k]\cdot k^m+
\sum_{k=xl}^9|dist|[l][k]\cdot k^m+|pd|\cdot|xl|^m$.
@ @=
int dist[maxm+1][16],pdist[maxm+1][16];
ull a[maxm+1][maxdigs],b[maxm+1][maxdigs],sig[maxm+1][maxdigs];
int x[maxm+1],rsave[maxm+1],tsave[maxm+1],pdsave[maxm+1];
@ @=
l=1;
pd=pdsave[1]=0;
alt=0,blt=9;
t=m,r=m+1;
xl=9;
profile[1]=1;
goto b3; /* I really don't want to do step |b2| at root level! */
@ @=
oo,tsave[l]=t,rsave[l]=r;
o,pdsave[l]=pd;
o,x[l++]=xl;
goto b2;
@ @=
{
count++;
printf("%d: ",count);
for (k=1;k<=m+1;k++) printf("%d",x[k]);
printf("->");
for (k=m;k>=0;k--) printf("%d",nybb(sig[l],k));
printf("\n");
goto b5;
}
@ When this code is performed, |sig[l]| and |dist[l]| and |pdist[l]|
are supposed to be up to date, as well as |xl|, |t|, |r|, |alt|, and |blt|.
@=
loop:@+if (t>=0) {
change=0;
@;
if (vbose>2) {
fprintf(stderr," a=");
printnum(a[l],t);
fprintf(stderr,",b=");
printnum(b[l],t);
fprintf(stderr,"\n");
}
if (change) goto loop; /* either $a_l$ or $b_l$ or both can be improved */
while (alt==blt) @;
if (change) goto loop;
}
@ The numbers |alt| and |blt| just past the prefix give important constraints
on what the future can bring. If we can improve them, we can often
improve them further yet, and possibly even extend the prefix.
@=
if (bltnybb(a[l],t)) {
fprintf(stderr,"Confusion (a decreased)!\n");
exit(-13);
}
alt=nybb(a[l],t);
if (blt=
{
o,p=pdist[l][blt];
if (blt>=xl) {
if (o,pxl) goto b5; /* oops, we've already saturated that digit */
pd=p+1-dist[l][blt]; /* |pd| becomes positive, if it wasn't already */
}
if (--r<0) goto b5;
add(sig[l],table[1][blt],sig[l]); /* newly known digit less than |xl| */
okay: o,pdist[l][blt]=p+1;
t--,change=1;
if (t<0) break;
oo,alt=nybb(a[l],t),blt=nybb(b[l],t);
}
@ @=
oo,t=tsave[l],r=rsave[l];
if (t>=0) oo,alt=nybb(a[l],t),blt=nybb(b[l],t);
else alt=blt=9;
o,xl=x[l];
@ When |dist| is ``catching up'' with |pdist|, we don't change |sig|,
because a digit that occurred in the prefix was already accounted for;
we knew that an |xl| would be coming, and it has finally arrived.
(Also |t| and |r| remain unchanged.)
@=
{
if (vbose>1) fprintf(stderr,"Level %d, that %d was forced\n",l,xl);
for (k=0;k