\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; 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.