next previous
Up: Expansions for nearly Gaussian


Appendix C: Algorithm for computing indices km

 To find all the solutions of Eq. (31) it is desirable to order all sets of non-negative integers $\{k_m\}$ satisfying it. We first rewrite (31) as  
 \begin{displaymath}
nk_n + ... +2k_2+ k_1=n \; .
 \end{displaymath} (C1)
It is natural to establish the ordering of the sets $S_i\equiv\{k_m\}$ satisfying Eq. (31) according to the order of numbers in their decimal representation. For n=3, say, we have 3 non-negative solutions

\begin{displaymath}
S_1=\{k_3=0,k_2=0,k_1=3\}\end{displaymath}

\begin{displaymath}
S_2=\{k_3=0,k_2=1,k_1=1\}\end{displaymath}

\begin{displaymath}
S_3=\{k_3=1,k_2=0,k_1=0\}\end{displaymath}

which can simply be written as

\begin{displaymath}
S_1=003, \; S_2=011, \; S_3=100 \; ,\end{displaymath}

and we say that S1<S2<S3 since 3<11<100. For $n\ge 10$,when the base 10 is no longer convenient, the sets of solutions can be ordered according to the order of the integer numbers

 
kn rn-1 + ... +k2 r+ k1 (C2)

for any natural base r. Those numbers are not important in themselves, what matters is thinking about the sets of $\{k_m\}$ as numbers in an abstract representation to base r.

On entry to our algorithm we have an integer ${\tt n} \gt 0$, and on exit we wish to have the number nsol of all solutions of Eq. (31) and the solutions themselves. For a set Si we introduce an abstract variable S to describe the array k(1:n) of n integer elements ordered as ``quasi-decimal digits''.

So, in an abstract form the algorithm looks like:

  --     Statement QS: 
  --             all S <= S(CURRENT) 
  --             are out if and only if 
  --             k(1)+2*k(2)+...+n*k(n)=n 
  <*initiate: make QS true for S=S(INITIAL) *>; 
  _while <*condition:  S ^= S(FINAL) *> _do 
    <*nextset: find next  S  
                        keeping QS invariant *> 
  _od; -- Here S=S(FINAL) and QS=.TRUE., i.e.
       -- all needed solutions are found
It is easy to initiate QS:

       k(1)=n; 
       mold=1; -- keeps largest m 
               -- for non-zero k(m) 
       _do m=2,n; 
         k(m)=0
       _od; 
       nsol=1;
The integer mold keeps track of the progress of the algorithm and allows us to establish the condition for the continuation of the main loop,

\begin{displaymath}
{\tt mold}\ne{\tt n} \; .\end{displaymath}

Finally, the core of the algorithm is the node nextset where we work with our sets as with digits, which is like doing the addition of decimal numbers by hand, column by column from right to left.
 
Table C1: Table of solutions to Eq. (31) for n=1 to 8. Zeros for higher orders are left blank

\begin{tabular}
{rrrrr}
\hline
n & 1 & 2 & 3 & 4 \\ \hline
&1 & 2 & 3 & 4 \\ & &...
 ...2 \\ & & & & 100010 \\ & & & & 1000001 \\ & & & &10000000 \\ \hline\end{tabular}

  --   advance S(CURRENT), i.e. consider adding
  --   n-( 2*k(2) +3*k(3) + ...+mold*k(mold) )
  --   to k(1) trying to add 1 to the lowest of
  --   k(2),...,k(mold) possible.
  --   If adding 1 to the lowest k(m) 
  --   makes the sum > n then
  --   put this k(m)=0 and try k(m+1) 
  -- 
   m=1;                                            
   sumcur=n; -- integer sum of current set S       
   _repeat                                         
       sumcur=sumcur-k(m)*m+m+1;                   
       k(m)=0;                                     
       k(m+1)=k(m+1)+1;                            
       m=m+1;                                      
   _until sumcur <= n _or  m>mold;                 
   _if m>mold _then mold=m _fi;                    
   k(1)=n-sumcur;                                  
   nsol=nsol+1;

Here the node nextset ends, and the whole algorithm is finished. We have written it here in a pseudocode which can be translated mechanically into any machine language. We actually use a special preprocessor Trefor which automatically transforms the text above to standard Fortran (see Weinstein & Blinnikov 1984, and Bartunov et al. 1997).

For reference, we present the first 8 sets of solutions found by this algorithm. In practice it is easier not to use the Table C1 even for low n, but generate all coefficients in the code.


next previous
Up: Expansions for nearly Gaussian

Copyright The European Southern Observatory (ESO)