reduce-historical/r35/lib/camal.tex

\documentstyle[fullpage]{article}
\title{REDUCE Meets CAMAL}
\author{J. P. Fitch \\
School of Mathematical Sciences\\
University of Bath\\
BATH, BA2 7AY, United Kingdom}
\def\today{}
\begin{document}\maketitle

\begin{abstract}
{\em It is generally accepted that special purpose algebraic systems
are more efficient than general purpose ones, but as machines get
faster this does not matter.  An experiment has been performed to see
if using the ideas of the special purpose algebra system CAMAL(F) it
is possible to make the general purpose system REDUCE perform
calculations in celestial mechanics as efficiently as CAMAL did twenty
years ago.  To this end a prototype Fourier module is created for
REDUCE, and it is tested on some small and medium-sized problems taken
from the CAMAL test suite. The largest calculation is the
determination of the Lunar Disturbing Function to the sixth order.  An
assessment is made as to the progress, or lack of it, which computer
algebra has made, and how efficiently we are using modern hardware.
}
\end{abstract}

\section{Introduction}

A number of years ago there emerged the divide between general-purpose
algebra systems and special purpose one.  Here we investigate how far
the improvements in software and more predominantly hardware have
enabled the general systems to perform as well as the earlier special
ones.  It is similar in some respects to the Possion program for
MACSYMA \cite{Fateman} which was written in response to a similar
challenge.

The particular subject for investigation is the Fourier series
manipulator which had its origins in the Cambridge University
Institute for Theoretical Astronomy, and later became the F subsystem
of CAMAL \cite{Barton67b,CAMALF}.  In the late 1960s this system was
used for both the Delaunay Lunar Theory \cite{Delaunay,Barton67a} and
the Hill Lunar Theory \cite{Bourne}, as well as other related
calculations.  Its particular area of application had a number of
peculiar operations on which the general speed depended.  These are
outlined below in the section describing how CAMAL worked.  There have
been a number of subsequent special systems for celestial mechanics,
but these tend to be restricted to the group of the originator.

The main body of the paper describes an experiment to create within
the REDUCE system a sub-system for the efficient manipulation of
Fourier series.  This prototype program is then assessed against both
the normal (general) REDUCE and the extant CAMAL results.  The tests
are run on a number of small problems typical of those for which CAMAL
was used, and one medium-sized problem, the calculation of the Lunar
Disturbing Function.  The mathematical background to this problem is
also presented for completeness.  It is important as a problem as it
is the first stage in the development of a Delaunay Lunar Theory.

The paper ends with an assessment of how close the performance of a
modern REDUCE on modern equipment is to the (almost) defunct CAMAL of
eighteen years ago.

\section{How CAMAL Worked}

The Cambridge Algebra System was initially written in assembler for
the Titan computer, but later was rewritten a number of times, and
matured in BCPL, a version which was ported to IBM mainframes and a
number of microcomputers.  In this section a brief review of the main
data structures and special algorithms is presented.

\subsection{CAMAL Data Structures}

CAMAL is a hierarchical system, with the representation of polynomials
being completely independent of the representations of the angular
parts.  

The angular part had to represent a polynomial coefficient, either a
sine or cosine function and a linear sum of angles.  In the problems
for which CAMAL was designed there are 6 angles only, and so the
design restricted the number, initially to six on the 24 bit-halfword
TITAN, and later to eight angles on the 32-bit IBM 370, each with
fixed names (usually u through z).  All that is needed is to remember
the coefficients of the linear sum.  As typical problems are
perturbations, it was reasonable to restrict the coefficients to small
integers, as could be represented in a byte with a guard bit.  This
allowed the representation to pack everything into four words.
\begin{verbatim}
    [ NextTerm, Coefficient, Angles0-3, Angles4-7 ]
\end{verbatim}
The function was coded by a single bit in the {\tt Coefficient} field.  This
gives a particularly compact representation.  For example the Fourier
term $\sin(u-2v+w-3x)$ would be represented as
\begin{verbatim}
    [ NULL, "1"|0x1, 0x017e017d, 0x00000000 ]
or
    [ NULL, "1"|0x1, 1:-2:1:-3, 0:0:0:0 ]
\end{verbatim}
where {\tt "1"} is a pointer to the representation of the polynomial
1.  In all this representation of the term took 48 bytes.  As the
complexity of a term increased the store requirements to no grow much;
the expression $(7/4) a e^3 f^5 \cos(u-2v+3w-4x+5y+6z)$ also takes 48
bytes.  There is a canonicalisation operation to ensure that the
leading angle is positive, and $\sin(0)$ gets removed.  It should be
noted that $\cos(0)$ is a valid and necessary representation.

The polynomial part was similarly represented, as a chain of terms
with packed exponents for a fixed number of variables.  There is no
particular significance in this except that the terms were held in
{\em increasing} total order, rather than the decreasing order which
is normal in general purpose systems.  This had a number of important
effects on the efficiency of polynomial multiplication in the presence
of a truncation to a certain order.  We will return to this point
later.  Full details of the representation can be found in
\cite{LectureNotes}.

The space administration system was based on explicit return rather
than garbage collection.  This meant that the system was sometimes
harder to write, but it did mean that much attention was focussed on
efficient reuse of space.  It was possible for the user to assist in
this by marking when an expression was needed no longer, and the
compiler then arranged to recycle the space as part of the actual
operation.  This degree of control was another assistance in running
of large problems on relatively small machines.

\subsection{Automatic Linearisation}

In order to maintain Fourier series in a canonical form it is
necessary to apply the transformations for linearising products of
sine and cosines.  These will be familiar to readers of the REDUCE
test program as
\begin{eqnarray}
\cos \theta \cos \phi & \Rightarrow & 
                (\cos(\theta+\phi)+\cos(\theta-\phi))/2, \\
\cos \theta \sin \phi & \Rightarrow & 
                (\sin(\theta+\phi)-\sin(\theta-\phi))/2, \\
\sin \theta \sin \phi & \Rightarrow & 
                (\cos(\theta-\phi)-\cos(\theta+\phi))/2, \\
\cos^2 \theta & \Rightarrow & (1+\cos(2\theta))/2,      \\
\sin^2 \theta & \Rightarrow & (1-\cos(2\theta))/2.
\end{eqnarray}
In CAMAL these transformations are coded directly into the
multiplication routines, and no action is necessary on the part of the
user to invoke them.  Of course they cannot be turned off either.

\subsection{Differentiation and Integration}

The differentiation of a Fourier series with respect to an angle is
particularly simple.  The integration of a Fourier series is a little
more interesting.  The terms like $\cos(n u + \ldots)$ are easily
integrated with respect to $u$, but the treatment of terms independent
of the angle would normally introduce a secular term.  By convention
in Fourier series these secular terms are ignored, and the constant of
integration is taken as just the terms independent of the angle in the
integrand.  This is equivalent to the substitution rules
\begin{eqnarray*}
\sin(n \theta) & \Rightarrow & -(1/n) \cos(n \theta) \\
\cos(n \theta) & \Rightarrow & (1/n) \sin(n \theta)
\end{eqnarray*}

In CAMAL these operations were coded directly, and independently of
the differentiation and integration of the polynomial coefficients.

\subsection{Harmonic Substitution}

An operation which is of great importance in Fourier operations is the
{\em harmonic substitution}.  This is the substitution of the sum of
some angles and a general expression for an angle.  In order to
preserve the format, the mechanism uses the translations
\begin{eqnarray*}
\sin(\theta + A) & \Rightarrow & \sin(\theta) \cos(A) + 
                                 \cos(\theta) \sin(A) \\
\cos(\theta + A) & \Rightarrow & \cos(\theta) \cos(A) -
                                 \sin(\theta) \sin(A) \\
\end{eqnarray*}
and then assuming that the value $A$ is small it can be replaced by
its expansion:
\begin{eqnarray*}
\sin(\theta + A) & \Rightarrow & \sin(\theta) \{1 - A^2/2! + A^4/4!\ldots\} +\\
                 &             & \cos(\theta) \{A - A^3/3! + A^5/5!\ldots\} \\
\cos(\theta + A) & \Rightarrow & \cos(\theta) \{1 - A^2/2! + A^4/4!\ldots\} -\\
                 &             & \sin(\theta) \{A - A^3/3! + A^5/5! \ldots\} \\
\end{eqnarray*}
If a truncation is set for large powers of the polynomial variables
then the series will terminate.  In CAMAL the {\tt HSUB} operation
took five arguments; the original expression, the angle for which
there is a substitution, the new angular part, the expression part
($A$ in the above), and the number of terms required.

The actual coding of the operation was not as expressed above, but by
the use of Taylor's theorem.  As has been noted above the
differentiation of a harmonic series is particularly easy.

\subsection{Truncation of Series}

The main use of Fourier series systems is in generating perturbation
expansions, and this implies that the calculations are performed to
some degree of the small quantities.  In the original CAMAL all
variables were assumed to be equally small (a restriction removed in
later versions).  By maintaining polynomials in increasing maximum
order it is possible to truncate the multiplication of two
polynomials.  Assume that we are multiplying the two polynomials
\begin{eqnarray*}
        A = a_0 + a_1 + a_2 + \ldots \\
        B = b_0 + b_1 + b_2 + \ldots
\end{eqnarray*}
If we are generating the partial answer
\[
        a_i (b_0 + b_1 + b_2 + \ldots)
\]
then if for some $j$ the product $a_i b_j$ vanishes, then so will all
products $a_i b_k$ for $k>j$.  This means that the later terms need
not be generated.  In the product of $1+x+x^2+x^3+\ldots+x^{10}$ and
$1+y+y^2+y^3+\ldots+y^10$ to a total order of 10 instead of generating
100 term products only 55 are needed.  The ordering can also make the
merging of the new terms into the answer easier.

\section{Towards a CAMAL Module}

For the purposes of this work it was necessary to reproduce as many of
the ideas of CAMAL as feasible within the REDUCE framework and
philosophy.  It was not intended at this stage to produce a complete
product, and so for simplicity a number of compromises were made with
the ``no restrictions'' principle in REDUCE and the space and time
efficiency of CAMAL.  This section describes the basic design
decisions. 

\subsection{Data Structures}

In a fashion similar to CAMAL a two level data representation is used.
The coefficients are the standard quotients of REDUCE, and their
representation need not concern us further.  The angular part is
similar to that of CAMAL, but the ability to pack angle multipliers
and use a single bit for the function are not readily available in
Standard LISP, so instead a longer vector is used.  Two versions were
written.  One used a balanced tree rather than a linear list for the
Fourier terms, this being a feature of CAMAL which was considered but
never coded.  The other uses a simple linear representation for sums.
The angle multipliers are held in a separate vector in order to allow
for future flexibility.  This leads to a representation as a vector of
length 6 or 4;
\begin{verbatim}
Version1:  [ BalanceBits, Coeff, Function, Angles, LeftTree, RightTree ]
Version2:  [ Coeff, Function, Angles, Next ]
\end{verbatim}
where the {\tt Angles} field is a vector of length 8, for the
multipliers.  It was decided to forego packing as for portability we
do not know how many to pack into a small integer.  The tree system
used is AVL, which needs 2 bits to maintain balance information, but
these are coded as a complete integer field in the vector.  We can
expect the improvements implicit in a binary tree to be advantageous
for large expressions, but the additional overhead may reduce its
utility for smaller expressions.

A separate vector is kept relating the position of an angle to its
print name, and on the property list of each angle the allocation of
its position is kept.  So long as the user declares which variables
are to be treated as angles this mechanism gives flexibility which was
lacking in CAMAL.

\subsection{Linearisation}

As in the CAMAL system the linearisation of products of sines and
cosines is done not by pattern matching but by direct calculation at
the heart of the product function, where the transformations (1)
through (3) are made in the product of terms function.  A side effect
of this is that there are no simple relations which can be used from
within the Fourier multiplication, and so a full addition of partial
products is required.  There is no need to apply linearisations
elsewhere as a special case.  Addition, differentiation and
integration cannot generate such products, and where they can occur in
substitution the natural algorithm uses the internal multiplication
function anyway.

\subsection{Substitution}

Substitution is the main operation of Fourier series.  It is useful to
consider three different cases of substitutions.
\begin{enumerate}
\item Angle Expression for Angle:
\item Angle Expression + Fourier Expression for Angle:
\item Fourier Expression for Polynomial Variable.
\end{enumerate}

The first of these is straightforward, and does not require any
further comment.  The second substitution requires a little more care,
but is not significantly difficult to implement.  The method follows
the algorithm used in CAMAL, using TAYLOR series.  Indeed this is the
main special case for substitution.

The problem is the last case.  Typically many variables used in a
Fourier series program have had a WEIGHT assigned to them.  This means
that substitution must take account of any possible WEIGHTs for
variables.  The standard code in REDUCE does this in effect by
translating the expression to prefix form, and recalculating the value.
A Fourier series has a large number of coefficients, and so this
operations are repeated rather too often.  At present this is the
largest problem area with the internal code, as will be seen in the
discussion of the Disturbing Function calculation.

\section{Integration with REDUCE}

The Fourier module needs to be seen as part of REDUCE rather than as a
separate language.  This can be seen as having internal and external
parts.

\subsection{Internal Interface}

The Fourier expressions need to co-exist with the normal REDUCE syntax
and semantics.  The prototype version does this by (ab)using the
module method, based in part on the TPS code \cite{Barnes}.  Of course
Fourier series are not constant, and so are not really domain
elements.  However by asserting that Fourier series form a ring of
constants REDUCE can arrange to direct basic operations to the Fourier
code for addition, subtraction, multiplication and the like.

The main interface which needs to be provided is a simplification
function for Fourier expressions.  This needs to provide compilation
for linear sums of angles, as well as constructing sine and cosine
functions, and creating canonical forms.

\subsection{User Interface}

The creation of {\tt HDIFF} and {\tt HINT} functions for
differentiation disguises this.  An unsatisfactory aspect of the
interface is that the tokens {\tt SIN} and {\tt COS} are already in
use.  The prototype uses the operator form
\begin{verbatim}
        fourier sin(u)
\end{verbatim}
to introduce harmonically represented sine functions.  An alternative of
using the tokens {\tt F\_SIN} and {\tt F\_COS} is also available.  

It is necessary to declare the names of the angles, which is achieved
with the declaration 
\begin{verbatim}
        harmonic theta, phi;
\end{verbatim}

At present there is no protection against using a variable as both an
angle and a polynomial varaible.  This will nooed to be done in a
user-oriented version.

\section{The Simple Experiments}

The REDUCE test file contains a simple example of a Fourier
calculation, determining the value of $(a_1 \cos({wt}) + a_3
\cos(3{wt}) + b_1 \sin({wt}) + b_3 \sin(3{wt}))^3$.  For the purposes
of this system this is too trivial to do more than confirm the correct
answers. 

The simplest non-trivial calculation for a Fourier series manipulator
is to solve Kepler's equation for the eccentric anomoly E in terms of
the mean anomoly u, and the eccentricity of an orbit e, considered as a
small quantity
\[
        E = u + e \sin E
\]
The solution procedes by repeated approximation.  Clearly the initial
approximation is $E_0 = u$.  The $n^{th}$ approximation can be written
as $u + A_n$, and so $A_n$ can be calculated by
\[
        A_k = e \sin (u + A_{k-1})
\]
This is of course precisely the case for which the HSUB operation is
designed, and so in order to calculate $E_n - u$ all one requires is
the code
\begin{verbatim}
        bige := fourier 0;
        for k:=1:n do <<
          wtlevel k;
          bige:=fourier e * hsub(fourier(sin u), u, u, bige, k);
        >>;
        write "Kepler Eqn solution:", bige$
\end{verbatim}

It is possible to create a regular REDUCE program to simulate this (as
is done for example in Barton and Fitch\cite{Barton72}, page 254).
Comparing these two programs indicates substantial advantages to the
Fourier module, as could be expected.
\medskip
\begin{center}
\begin{tabular}{ | c | l l |}
\multicolumn{3}{c}{\bf Solving Kepler's Equation} \\
\hline
Order   &       REDUCE  &       Fourier Module \\
5       &       9.16    &       2.48    \\
6       &       17.40   &       4.56    \\
7       &       33.48   &       8.06    \\
8       &       62.76   &       13.54   \\
9       &       116.06  &       21.84   \\
10      &       212.12  &       34.54   \\
11      &       381.78  &       53.94   \\
12      &       692.56  &       82.96   \\
13      &       1247.54 &       125.86  \\
14      &       2298.08 &       187.20  \\
15      &       4176.04 &       275.60  \\
16      &       7504.80 &       398.62  \\
17      &       13459.80        &       569.26  \\
18      &       ***     &       800.00  \\
19      &       ***     &       1116.92 \\
20      &       ***     &       1536.40 \\
\hline
\end{tabular}
\end{center}
\medskip
These results were with the linear representation of Fourier series.
The tree representation was slightly slower.  The ten-fold speed-up
for the 13th order is most satisfactory.

\section{A Medium-Sized Problem}

Fourier series manipulators are primarily designed for large-scale
calculations, but for the demonstration purposes of this project a
medium problem is considered.  The first stage in calculating the
orbit of the Moon using the Delaunay theory (of perturbed elliptic
motion for the restricted 3-body problem) is to calculate the energy
of the Moon's motion about the Earth --- the Hamiltonian of the
system.   This is the calculation we use for comparisons.

\subsection{Mathematical Background}

The full calculation is described in detail in \cite{Brown}, but a
brief description is given here for completeness, and to grasp the
extent of the calculation.

Referring to the figure 1 which gives the cordinate system, the basic
equations are 
\begin{eqnarray}
S  & = & (1-\gamma ^2)\cos(f + g +h -f' -g' -h')
+ \gamma ^2 cos(f + g -h +f' +g' +h') \\
r & = & a (1 - e \cos E) \\
l & = & E - e \sin E \\
a & = & r {{\bf d} E} \over {{\bf d} l} \\
r ^2 {{\bf d} f} \over {{\bf d} l} & = & a^2 (1 - e^2)^{1 \over 2}\\
R & = & m' {a^2 \over {a'} ^3} {{a'}\over {r
'}} \left \{ \left ({r \over a}\right )^2
\left ({{a'} \over {r'}}\right )^2 P_2(S) +
\left ({a \over {a'}}\right )\left
({r \over a}\right )^3 \left ({{a'} \over {r'}}\right )^3 P_3(S)
+ \ldots \right \}
\end{eqnarray}

There are similar equations to (7) to (10) for the quantities $r'$,
$a'$, $e'$, $l'$, $E'$ and $f'$ which refer to the position of the Sun
rather than the Moon.  The problem is to calculate the expression $R$
as an expansion in terms of the quantities $e$, $e'$, $\gamma$,
$a/a'$, $l$, $g$, $h$, $l'$, $g'$ and $h'$.  The first three
quantities are small quantities of the first order, and $a/a'$ is of
second order.

The steps required are
\begin{enumerate}
\item Solve the Kepler equation (8)
\item Substiture into (7) to give $r/a$ in terms of $e$ and $l$.
\item Calculate $a/r$ from (9) and $f$ from (10)
\item Substitute for $f$ and $f'$ into $S$ using (6)
\item Calculate $R$ from $S$, $a'/r'$ and $r/a$
\end{enumerate}

The program is given in the Appendix.

\subsection{Results}

The Lunar Disturbing function was calculated by a direct coding of the
previous sections' mathematics.  The program was taken from Barton
and Fitch \cite{Barton72} with just small changes to generalise it for
any order, and to make it acceptable for Reduce3.4.  The Fourier
program followed the same pattern, but obviously used the {\tt HSUB}
operation as appropriate and the harmonic integration.  It is very
similar to the CAMAL program in \cite{Barton72}.   

The disturbing function was calculated to orders 2, 4 and 6 using
Cambridge LISP on an HLH Orion 1/05 (Intergraph Clipper), with the
three programs $\alpha$) Reduce3.4, $\beta$) Reduce3.4 + Camal Linear
Module and $\gamma$) Reduce3.4 + Camal AVL Module.  The timings for
CPU seconds (excluding garbage collection time) are summarised the
following table:
\medskip
\begin{center}
\begin{tabular}{ | c || l | l | l |}
\hline
Order of DDF    & Reduce        & Camal Linear  & Camal Tree \\
\hline
2       &       23.68   &       11.22   &       12.9    \\
4       &       429.44  &       213.56  &       260.64  \\
6       &       $>$7500 &       3084.62 &       3445.54 \\
\hline
%%% Linear n=4 138.72 (4Mb + unsafe vector access + recurrance)
%%% Linear n=6 1870.10 (4Mb + unsafe vector access + recurrance)
\end{tabular}
\end{center}
\medskip

If these numbers are normalised so REDUCE calculating the DDF is 100
units for each order the table becomes
\medskip
\begin{center}
\begin{tabular}{ | c || l | l | l |}
\hline
Order of DDF    & Reduce        & Camal Linear  & Camal Tree \\ \hline
2       &       100     &       47.38   &       54.48   \\
4       &       100     &       49.73   &       60.69   \\
6       &       100     &       $<$41.13        &       $<$45.94 \\
\hline
\end{tabular}
\end{center}
\medskip

From this we conclude that a doubling of speed is about correct, and
although the balanced tree system is slower as the problem size
increases the gap between it and the simpler linear system is
narrowing. 

It is disappointing that the ratio is not better, nor the absolute
time less.  It is worth noting in this context that Jefferys claimed
that the sixth order DDF took 30s on a CDC6600 with TRIGMAN in 1970
\cite{Jefferys}, and Barton and Fitch took about 1s for the second
order DDF on TITAN with CAMAL \cite{Barton72}.  A closer look at the
relative times for individual sections of the program shows that the
substitution case of replacing a polynomial variable by a Fourier
series is only marginally faster than the simple REDUCE program.  In
the DDF program this operation is only used once in a major form,
substituting into the Legendre polynomials, which have been previously
calculated by Rodrigues formula.  This suggests that we replace this
with the recurrence relationship.

Making this change actually slows down the normal REDUCE by a small
amount but makes a significant change to the Fourier module; it
reduces the run time for the 6th order DDF from 3084.62s to 2002.02s.
This gives some indication of the problems with benchmarks.  What is
clear is that the current implementation of substitution of a Fourier
series for a polynomial variable is inadequate.

\section{Conclusion}

The Fourier module is far from complete.  The operations necessary for
the solution of Duffing's and Hill's equations are not yet written,
although they should not cause much problem.  The main defficiency is
the treatment of series truncation; at present it relies on the REDUCE
WTLEVEL mechanism, and this seems too coarse for efficient truncation.
It would be possible to re-write the polynomial manipulator as well,
while retaining the REDUCE syntax, but that seems rather more than one
would hope.

The real failure so far is the large time lag between the REDUCE-based
system on a modern workstation against a mainframe of 25 years ago
running a special system.  The CAMAL Disturbing function program could
calculate the tenth order with a maximum of 32K words (about
192Kbytes) whereas this system failed to calculate the eigth order in
4Mbytes (taking 2000s before failing).  I have in my archives the
output from the standard CAMAL test suite, which includes a sixth
order DDF on an IBM 370/165 run on 2 June 1978, taking 22.50s and
using a maximum of 15459 words of memory for heap --- or about
62Kbytes.  A rough estimate is that the Orion 1/05 is comparable in
speed to the 360/165, but with more real memory and virtual memory.

However, a simple Fourier manipulator has been created for REDUCE which
performs between twice and three times the speed of REDUCE using
pattern matching.  It has been shown that this system is capable of
performing the calculations of celestial mechanics, but it still
seriously lags behind the efficiency of the specialist systems of
twenty years before.  It is perhaps fortunate that it was not been
possible to compare it with a modern specialist system.

There is still work to do to provide a convenient user interface, but
it is intended to develop the system in this direction.  It would be
pleasant to have again a system of the efficiency of CAMAL(F).

I would like to thank Codemist Ltd for the provision of computing
resources for this project, and David Barton who taught be so much
about Fourier series and celstial mechanics.  Thank are also due to
the National Health Service, without whom this work and paper could not
have been produced.

\section*{Appendix: The DDF Function}
\begin{verbatim}
array p(n/2+2);
harmonic u,v,w,x,y,z;
weight e=1, b=1, d=1, a=1;

%% Generate Legendre Polynomials to sufficient order
for i:=2:n/2+2 do <<
  p(i):=(h*h-1)^i;
  for j:=1:i do p(i):=df(p(i),h)/(2j)
>>;

%%%%%%%%%%%%%%%% Step1: Solve Kepler equation
bige := fourier 0;
for k:=1:n do <<
  wtlevel k;
  bige:=fourier e * hsub(fourier(sin u), u, u, bige, k);
>>;

%% Ensure we do not calculate things of too high an order
wtlevel n;

%%%%%%%%%%%%%%%% Step 2: Calculate r/a in terms of e and l
dd:=-e*e; hh:=3/2; j:=1; cc := 1;
for i:=1:n/2 do <<
  j:=i*j; hh:=hh-1; cc:=cc+hh*(dd^i)/j
>>;
bb:=hsub(fourier(1-e*cos u), u, u, bige, n);
aa:=fourier 1+hdiff(bige,u); ff:=hint(aa*aa*fourier cc,u);

%%%%%%%%%%%%%%%% Step 3: a/r and f
uu := hsub(bb,u,v); uu:=hsub(uu,e,b);
vv := hsub(aa,u,v); vv:=hsub(vv,e,b);
ww := hsub(ff,u,v); ww:=hsub(ww,e,b);

%%%%%%%%%%%%%%%% Step 4: Substitute f and f' into S
yy:=ff-ww; zz:=ff+ww;
xx:=hsub(fourier((1-d*d)*cos(u)),u,u-v+w-x-y+z,yy,n)+
    hsub(fourier(d*d*cos(v)),v,u+v+w+x+y-z,zz,n);

%%%%%%%%%%%%%%%% Step 5: Calculate R
zz:=bb*vv; yy:=zz*zz*vv;

on fourier;
for i := 2:n/2+2 do <<
  wtlevel n+4-2i; p(i) := hsub(p(i), h, xx) >>;

wtlevel n;
for i:=n/2+2 step -1 until 3 do
    p(n/2+2):=fourier(a*a)*zz*p(n/2+2)+p(i-1);
yy*p(n/2+2);

\end{verbatim}
\newpage
\bibliographystyle{plain}
\bibliography{camal}

\end{document}