\datethis
@* Introduction. This program finds the length of the shortest strong
addition chain that leads to a given number~$n$.
A strong addition chain --- aka a Lucas chain or a Chebyshev chain ---
is a sequence of integers $1=a_0
int L[maxn]; /* the results */
int ub[maxm],lb[maxm]; /* upper and lower bounds on $a_k$ */
int choice[4*maxm]; /* current choices at each level */
int bound[4*maxm]; /* maximum possible choices at each level */
struct { int*ptr; int val; } assigned[8*maxm*maxm]; /* for undoing */
int undo_ptr; /* pointer to the |assigned| stack */
int save[4*maxm]; /* information for undoing at each level */
int hint[4*maxm]; /* additional information at each level */
int verbose=0; /* set nonzero when debugging */
int record=0; /* largest $L(n)$ seen so far */
@@;
main()
{
register int a,k,l,m,n,t,work;
int special=0, ap, app;
ub[0]=lb[0]=1;
n=1;@+ m=0;@+ work=0;
while (1) {
L[n]=m;
printf("L(%d)=%d:",n,m);
if (m>record) {
record=m; printf("*");
}
if (special) @@;
else for (k=0;k<=m;k++) printf(" %d",ub[k]);
printf(" [%d]\n",work);
n++;@+ work=0;
@;
}
}
@* Backtracking. The choices to be made on levels 0, 1, 2, \dots\ can
be grouped in fours, so that the actions at level $l$ depend on $l\bmod4$.
If $l\bmod4=0$, we're given a number $a$ and we want to decide how to
write it as a sum $a'+a''$. If $l\bmod4!=0$, we're given a number $b$ and
we want to place it in the chain if it isn't already present. The
cases $l\bmod4=1$,~2,~3 correspond respectively to $b=a'$, $b=a''$, and
$b=\vert a'-a''\vert$ in the previous choice $a=a'+a''$.
We keep the current state of the chain in arrays |lb| and |ub|.
If |lb[k]=ub[k]|,
their common value is $a_k$; otherwise $a_k$ has not yet been specified,
but its eventual value will satisfy $|lb|[k]\le a_k\le|ub|[k]$.
These bounds are maintained in such a way that
$$\hbox{|ub[k]=
@;
while (1) {
@;
@;
@;
not_found: m++;
}
found:
@ The obvious lower bound is $\lceil\lg n\rceil$.
@=
for (k=(n-1)>>1,m=1;k;m++) k>>=1;
@ A slightly subtle point arises here: We make |lb[k]| and |ub[k]| infinite for
$k>m$, because some of our routines below will look in positions $\ge L(a)$
when trying to insert an element~$a$.
@=
ub[m]=lb[m]=n;
for (k=m-1;k;k--) {
lb[k]=(lb[k+1]+1)>>1;
if (lb[k]<=k) lb[k]=k+1;
}
for (k=1;kn-(m-k)) ub[k]=n-(m-k);
}
l=0; t=m+1;
for (k=t;k<=record;k++) lb[k]=ub[k]=maxn;
undo_ptr=0;
@ At each level |l| we make all choices that lie between |choice[l]| and
|bound[l]|, inclusive. If successful, we advance |l| and go to |start_level|;
otherwise we go to |backup|.
@=
start_level: work++;
if (verbose) @;
if (l&3) @@;
else @;
@;
@ @=
{
printf("Entering level %d:\n",l);
for (k=1;k=
int lookup(int x) /* is $x$ already in the chain? */
{
register int k;
if (x<=2) return 1;
for (k=L[x];x>ub[k];k++) ;
return x==ub[k] && x==lb[k];
}
@ The values of $a_1$ and $a_2$ can never be a problem.
@=
save[l]=t; /* remember the current value of $t$, in case we fail */
decr_t: t--;
if (t<=2) goto found;
if (ub[t]>lb[t]) goto restore_t_and_backup;
a=ub[t];
for (k=t-1;;k--) if (ub[k]==lb[k]) {
ap=ub[k], app=a-ap;
if (app>ap) break;
if (lookup(app) && lookup(ap-app)) goto decr_t; /* yes, it's OK already */
}
choice[l]=(a+1)>>1; /* the minimum choice is $a'=\lceil a/2\rceil$ */
bound[l]=a-1; /* and the maximum choice is $a'=a-1$ */
vet_it: @;
advance: l++;@+goto start_level;
@ The |impossible| subroutine determines rapidly when there is no
``hole'' in which an element can be placed in the current chain.
@=
int impossible(int x) /* is there obviously no way to put $x$ in? */
{
register int k;
if (x<=2) return 0;
for (k=L[x];x>ub[k];k++) ;
return x=
ap=choice[l];@+ app=a-ap;
if (impossible(ap) || impossible(app) || impossible(ap-app))
goto next_choice;
hint[l+1]=ap;@+ hint[l+2]=app;@+ hint[l+3]=ap-app;
@* Placing the summands. Any change to the |ub| and |lb| table needs to
be recorded in the |assigned| array, because we may need to undo it.
@d assign(x,y) assigned[undo_ptr].ptr=x, assigned[undo_ptr++].val=*x, *x=y
@ The algorithm for adjusting upper and lower bounds is probably the most
interesting part of this whole program. I suppose I should prove it correct.
(Since this subroutine is called only in one place, I might want to try
experimenting to see how much faster this program runs when subroutine-call
overhead is avoided by converting to inline code. Subroutining might actually
turn out to be a win because of the limited number of registers
on x86-like computers.)
@=
place(int x,int k) /* set $a_k=x$ */
{
register int j=k, y=x;
if (ub[j]==y && lb[j]==y) return 0;
while (ub[j]>y) {
assign(&ub[j],y); /* the upper bound decreases */
j--, y--;
}
j=k+1, y=x+x;
while (ub[j]>y) {
assign(&ub[j],y); /* the upper bound decreases */
j++, y+=y;
}
j=k, y=x;
while (lb[j]>1;
}
j=k+1, y=x+1;
while (lb[j]=
int choice_lookup(int x) /* find the smallest viable place for $x$ */
{
register int k;
if (x<=2) return 0;
for (k=L[x];x>ub[k];k++) ;
return k;
}
@ In the special case that the entry to be placed is already present,
we avoid unnecessary computation by setting |bound[l]| to zero.
(Note: I thought this would be a good idea, but it didn't actually
decrease the observed running time.)
@=
{
a=hint[l];
save[l]=undo_ptr;
k=choice[l]=choice_lookup(a);
if (k==0 || (a==ub[k] && a==lb[k])) {
bound[l]=0; goto advance;
} else {
while (a>=lb[k]) k++;
bound[l]=k-1;
}
goto next_place;
}
@ @=
unplace: if (!bound[l]) goto backup;
while (undo_ptr>save[l]) {
--undo_ptr;
*assigned[undo_ptr].ptr=assigned[undo_ptr].val;
}
choice[l]++;
a=hint[l];
next_place:@+ if (choice[l]>bound[l]) goto backup;
place(a,choice[l]);
goto advance;
@ Finally, when we run out of steam on the current level, we reconsider
previous choices as follows.
@ @=
restore_t_and_backup: t=save[l];
backup: if (l==0) goto not_found;
--l;
if (l&3) goto unplace; /* $l\bmod4=1$, 2, or 3 */
a=ub[t]; /* $l\bmod4=0$ */
next_choice: choice[l]++;
if (choice[l]<=bound[l]) goto vet_it;
goto restore_t_and_backup;
@* Simple upper bounds. We can often save a lot of work by using the fact
that $L(mn)\le L(m)+L(n)$.
@=
{
for (k=2,a=n/k; k<=a; k++,a=n/k)
if (n%k==0 && m==L[k]+L[a]) {
special=k; goto found;
}
@;
}
@ Another simple upper bound,
$L(n)\le\lfloor\lg n\rfloor+\lfloor\lg{2\over3}n\rfloor$,
follows from the fact that a strong chain ending with $(a,a+1)$
can be extended by appending either $(2a,2a+1)$ or $(2a+1,2a+2)$.
I programmed it here just to see how often it helps, but I doubt if it will
be very effective. (Indeed, experience showed that it was the method of
choice only for $n=2$, 3, 5, 7, 11, and 23; probably not for any larger~$n$.)
Incidentally,
the somewhat plausible inequality $L(2n+1)\le L(n)+2$ is {\it not\/}
true, although the analogous inequality
$l(2n+1)\le l(n)+2$ obviously holds for ordinary addition chains.
Indeed, $L(17)=6$ and $L(8)=3$.
@=
if (m==lg(n)+lg((n+n)/3)) special=1;
@ @=
int lg(int n)
{
register int m, x;
for (x=n>>1, m=0; x; m++) x>>=1;
return m;
}
@ @=
{
if (special==1) printf(" Binary method");
else printf(" Factor method %d x %d", special, n/special);
special=0;
}
@ Experience showed that the factor method often gives an optimum result,
at least for small~$n$. Indeed, the factor method was optimum for all
nonprime $n<1219$. (The first exception, 1219, is $23\times53$, the
product of two primes that have worse-than-normal $L$~values.)
Yet the factoring shortcut reduced the total running time by only
about 4\%, because it didn't help with the hard cases --- the cases that
keep the computer working hardest. (These timing statistics are based
only on the calculations for $n\le1000$; larger values of~$n$ may well
be a different story. But I think most of the running time goes into
proving that shorter chains are impossible.)
@* Index.