camal.tex 28 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632
  1. \documentstyle[11pt]{article}
  2. \title{REDUCE Meets CAMAL}
  3. \author{J. P. Fitch \\
  4. School of Mathematical Sciences\\
  5. University of Bath\\
  6. BATH, BA2 7AY, United Kingdom}
  7. \def\today{}
  8. \begin{document}\maketitle
  9. \begin{abstract}
  10. {\em It is generally accepted that special purpose algebraic systems
  11. are more efficient than general purpose ones, but as machines get
  12. faster this does not matter. An experiment has been performed to see
  13. if using the ideas of the special purpose algebra system CAMAL(F) it
  14. is possible to make the general purpose system REDUCE perform
  15. calculations in celestial mechanics as efficiently as CAMAL did twenty
  16. years ago. To this end a prototype Fourier module is created for
  17. REDUCE, and it is tested on some small and medium-sized problems taken
  18. from the CAMAL test suite. The largest calculation is the
  19. determination of the Lunar Disturbing Function to the sixth order. An
  20. assessment is made as to the progress, or lack of it, which computer
  21. algebra has made, and how efficiently we are using modern hardware.
  22. }
  23. \end{abstract}
  24. \section{Introduction}
  25. A number of years ago there emerged the divide between general-purpose
  26. algebra systems and special purpose one. Here we investigate how far
  27. the improvements in software and more predominantly hardware have
  28. enabled the general systems to perform as well as the earlier special
  29. ones. It is similar in some respects to the Possion program for
  30. MACSYMA \cite{Fateman} which was written in response to a similar
  31. challenge.
  32. The particular subject for investigation is the Fourier series
  33. manipulator which had its origins in the Cambridge University
  34. Institute for Theoretical Astronomy, and later became the F subsystem
  35. of CAMAL \cite{Barton67b,CAMALF}. In the late 1960s this system was
  36. used for both the Delaunay Lunar Theory \cite{Delaunay,Barton67a} and
  37. the Hill Lunar Theory \cite{Bourne}, as well as other related
  38. calculations. Its particular area of application had a number of
  39. peculiar operations on which the general speed depended. These are
  40. outlined below in the section describing how CAMAL worked. There have
  41. been a number of subsequent special systems for celestial mechanics,
  42. but these tend to be restricted to the group of the originator.
  43. The main body of the paper describes an experiment to create within
  44. the REDUCE system a sub-system for the efficient manipulation of
  45. Fourier series. This prototype program is then assessed against both
  46. the normal (general) REDUCE and the extant CAMAL results. The tests
  47. are run on a number of small problems typical of those for which CAMAL
  48. was used, and one medium-sized problem, the calculation of the Lunar
  49. Disturbing Function. The mathematical background to this problem is
  50. also presented for completeness. It is important as a problem as it
  51. is the first stage in the development of a Delaunay Lunar Theory.
  52. The paper ends with an assessment of how close the performance of a
  53. modern REDUCE on modern equipment is to the (almost) defunct CAMAL of
  54. eighteen years ago.
  55. \section{How CAMAL Worked}
  56. The Cambridge Algebra System was initially written in assembler for
  57. the Titan computer, but later was rewritten a number of times, and
  58. matured in BCPL, a version which was ported to IBM mainframes and a
  59. number of microcomputers. In this section a brief review of the main
  60. data structures and special algorithms is presented.
  61. \subsection{CAMAL Data Structures}
  62. CAMAL is a hierarchical system, with the representation of polynomials
  63. being completely independent of the representations of the angular
  64. parts.
  65. The angular part had to represent a polynomial coefficient, either a
  66. sine or cosine function and a linear sum of angles. In the problems
  67. for which CAMAL was designed there are 6 angles only, and so the
  68. design restricted the number, initially to six on the 24 bit-halfword
  69. TITAN, and later to eight angles on the 32-bit IBM 370, each with
  70. fixed names (usually u through z). All that is needed is to remember
  71. the coefficients of the linear sum. As typical problems are
  72. perturbations, it was reasonable to restrict the coefficients to small
  73. integers, as could be represented in a byte with a guard bit. This
  74. allowed the representation to pack everything into four words.
  75. \begin{verbatim}
  76. [ NextTerm, Coefficient, Angles0-3, Angles4-7 ]
  77. \end{verbatim}
  78. The function was coded by a single bit in the {\tt Coefficient} field. This
  79. gives a particularly compact representation. For example the Fourier
  80. term $\sin(u-2v+w-3x)$ would be represented as
  81. \begin{verbatim}
  82. [ NULL, "1"|0x1, 0x017e017d, 0x00000000 ]
  83. or
  84. [ NULL, "1"|0x1, 1:-2:1:-3, 0:0:0:0 ]
  85. \end{verbatim}
  86. where {\tt "1"} is a pointer to the representation of the polynomial
  87. 1. In all this representation of the term took 48 bytes. As the
  88. complexity of a term increased the store requirements to no grow much;
  89. the expression $(7/4) a e^3 f^5 \cos(u-2v+3w-4x+5y+6z)$ also takes 48
  90. bytes. There is a canonicalisation operation to ensure that the
  91. leading angle is positive, and $\sin(0)$ gets removed. It should be
  92. noted that $\cos(0)$ is a valid and necessary representation.
  93. The polynomial part was similarly represented, as a chain of terms
  94. with packed exponents for a fixed number of variables. There is no
  95. particular significance in this except that the terms were held in
  96. {\em increasing} total order, rather than the decreasing order which
  97. is normal in general purpose systems. This had a number of important
  98. effects on the efficiency of polynomial multiplication in the presence
  99. of a truncation to a certain order. We will return to this point
  100. later. Full details of the representation can be found in
  101. \cite{LectureNotes}.
  102. The space administration system was based on explicit return rather
  103. than garbage collection. This meant that the system was sometimes
  104. harder to write, but it did mean that much attention was focussed on
  105. efficient reuse of space. It was possible for the user to assist in
  106. this by marking when an expression was needed no longer, and the
  107. compiler then arranged to recycle the space as part of the actual
  108. operation. This degree of control was another assistance in running
  109. of large problems on relatively small machines.
  110. \subsection{Automatic Linearisation}
  111. In order to maintain Fourier series in a canonical form it is
  112. necessary to apply the transformations for linearising products of
  113. sine and cosines. These will be familiar to readers of the REDUCE
  114. test program as
  115. \begin{eqnarray}
  116. \cos \theta \cos \phi & \Rightarrow &
  117. (\cos(\theta+\phi)+\cos(\theta-\phi))/2, \\
  118. \cos \theta \sin \phi & \Rightarrow &
  119. (\sin(\theta+\phi)-\sin(\theta-\phi))/2, \\
  120. \sin \theta \sin \phi & \Rightarrow &
  121. (\cos(\theta-\phi)-\cos(\theta+\phi))/2, \\
  122. \cos^2 \theta & \Rightarrow & (1+\cos(2\theta))/2, \\
  123. \sin^2 \theta & \Rightarrow & (1-\cos(2\theta))/2.
  124. \end{eqnarray}
  125. In CAMAL these transformations are coded directly into the
  126. multiplication routines, and no action is necessary on the part of the
  127. user to invoke them. Of course they cannot be turned off either.
  128. \subsection{Differentiation and Integration}
  129. The differentiation of a Fourier series with respect to an angle is
  130. particularly simple. The integration of a Fourier series is a little
  131. more interesting. The terms like $\cos(n u + \ldots)$ are easily
  132. integrated with respect to $u$, but the treatment of terms independent
  133. of the angle would normally introduce a secular term. By convention
  134. in Fourier series these secular terms are ignored, and the constant of
  135. integration is taken as just the terms independent of the angle in the
  136. integrand. This is equivalent to the substitution rules
  137. \begin{eqnarray*}
  138. \sin(n \theta) & \Rightarrow & -(1/n) \cos(n \theta) \\
  139. \cos(n \theta) & \Rightarrow & (1/n) \sin(n \theta)
  140. \end{eqnarray*}
  141. In CAMAL these operations were coded directly, and independently of
  142. the differentiation and integration of the polynomial coefficients.
  143. \subsection{Harmonic Substitution}
  144. An operation which is of great importance in Fourier operations is the
  145. {\em harmonic substitution}. This is the substitution of the sum of
  146. some angles and a general expression for an angle. In order to
  147. preserve the format, the mechanism uses the translations
  148. \begin{eqnarray*}
  149. \sin(\theta + A) & \Rightarrow & \sin(\theta) \cos(A) +
  150. \cos(\theta) \sin(A) \\
  151. \cos(\theta + A) & \Rightarrow & \cos(\theta) \cos(A) -
  152. \sin(\theta) \sin(A) \\
  153. \end{eqnarray*}
  154. and then assuming that the value $A$ is small it can be replaced by
  155. its expansion:
  156. \begin{eqnarray*}
  157. \sin(\theta + A) & \Rightarrow & \sin(\theta) \{1 - A^2/2! + A^4/4!\ldots\} +\\
  158. & & \cos(\theta) \{A - A^3/3! + A^5/5!\ldots\} \\
  159. \cos(\theta + A) & \Rightarrow & \cos(\theta) \{1 - A^2/2! + A^4/4!\ldots\} -\\
  160. & & \sin(\theta) \{A - A^3/3! + A^5/5! \ldots\} \\
  161. \end{eqnarray*}
  162. If a truncation is set for large powers of the polynomial variables
  163. then the series will terminate. In CAMAL the {\tt HSUB} operation
  164. took five arguments; the original expression, the angle for which
  165. there is a substitution, the new angular part, the expression part
  166. ($A$ in the above), and the number of terms required.
  167. The actual coding of the operation was not as expressed above, but by
  168. the use of Taylor's theorem. As has been noted above the
  169. differentiation of a harmonic series is particularly easy.
  170. \subsection{Truncation of Series}
  171. The main use of Fourier series systems is in generating perturbation
  172. expansions, and this implies that the calculations are performed to
  173. some degree of the small quantities. In the original CAMAL all
  174. variables were assumed to be equally small (a restriction removed in
  175. later versions). By maintaining polynomials in increasing maximum
  176. order it is possible to truncate the multiplication of two
  177. polynomials. Assume that we are multiplying the two polynomials
  178. \begin{eqnarray*}
  179. A = a_0 + a_1 + a_2 + \ldots \\
  180. B = b_0 + b_1 + b_2 + \ldots
  181. \end{eqnarray*}
  182. If we are generating the partial answer
  183. \[
  184. a_i (b_0 + b_1 + b_2 + \ldots)
  185. \]
  186. then if for some $j$ the product $a_i b_j$ vanishes, then so will all
  187. products $a_i b_k$ for $k>j$. This means that the later terms need
  188. not be generated. In the product of $1+x+x^2+x^3+\ldots+x^{10}$ and
  189. $1+y+y^2+y^3+\ldots+y^10$ to a total order of 10 instead of generating
  190. 100 term products only 55 are needed. The ordering can also make the
  191. merging of the new terms into the answer easier.
  192. \section{Towards a CAMAL Module}
  193. For the purposes of this work it was necessary to reproduce as many of
  194. the ideas of CAMAL as feasible within the REDUCE framework and
  195. philosophy. It was not intended at this stage to produce a complete
  196. product, and so for simplicity a number of compromises were made with
  197. the ``no restrictions'' principle in REDUCE and the space and time
  198. efficiency of CAMAL. This section describes the basic design
  199. decisions.
  200. \subsection{Data Structures}
  201. In a fashion similar to CAMAL a two level data representation is used.
  202. The coefficients are the standard quotients of REDUCE, and their
  203. representation need not concern us further. The angular part is
  204. similar to that of CAMAL, but the ability to pack angle multipliers
  205. and use a single bit for the function are not readily available in
  206. Standard LISP, so instead a longer vector is used. Two versions were
  207. written. One used a balanced tree rather than a linear list for the
  208. Fourier terms, this being a feature of CAMAL which was considered but
  209. never coded. The other uses a simple linear representation for sums.
  210. The angle multipliers are held in a separate vector in order to allow
  211. for future flexibility. This leads to a representation as a vector of
  212. length 6 or 4;
  213. \begin{verbatim}
  214. Version1: [ BalanceBits, Coeff, Function, Angles, LeftTree, RightTree ]
  215. Version2: [ Coeff, Function, Angles, Next ]
  216. \end{verbatim}
  217. where the {\tt Angles} field is a vector of length 8, for the
  218. multipliers. It was decided to forego packing as for portability we
  219. do not know how many to pack into a small integer. The tree system
  220. used is AVL, which needs 2 bits to maintain balance information, but
  221. these are coded as a complete integer field in the vector. We can
  222. expect the improvements implicit in a binary tree to be advantageous
  223. for large expressions, but the additional overhead may reduce its
  224. utility for smaller expressions.
  225. A separate vector is kept relating the position of an angle to its
  226. print name, and on the property list of each angle the allocation of
  227. its position is kept. So long as the user declares which variables
  228. are to be treated as angles this mechanism gives flexibility which was
  229. lacking in CAMAL.
  230. \subsection{Linearisation}
  231. As in the CAMAL system the linearisation of products of sines and
  232. cosines is done not by pattern matching but by direct calculation at
  233. the heart of the product function, where the transformations (1)
  234. through (3) are made in the product of terms function. A side effect
  235. of this is that there are no simple relations which can be used from
  236. within the Fourier multiplication, and so a full addition of partial
  237. products is required. There is no need to apply linearisations
  238. elsewhere as a special case. Addition, differentiation and
  239. integration cannot generate such products, and where they can occur in
  240. substitution the natural algorithm uses the internal multiplication
  241. function anyway.
  242. \subsection{Substitution}
  243. Substitution is the main operation of Fourier series. It is useful to
  244. consider three different cases of substitutions.
  245. \begin{enumerate}
  246. \item Angle Expression for Angle:
  247. \item Angle Expression + Fourier Expression for Angle:
  248. \item Fourier Expression for Polynomial Variable.
  249. \end{enumerate}
  250. The first of these is straightforward, and does not require any
  251. further comment. The second substitution requires a little more care,
  252. but is not significantly difficult to implement. The method follows
  253. the algorithm used in CAMAL, using TAYLOR series. Indeed this is the
  254. main special case for substitution.
  255. The problem is the last case. Typically many variables used in a
  256. Fourier series program have had a WEIGHT assigned to them. This means
  257. that substitution must take account of any possible WEIGHTs for
  258. variables. The standard code in REDUCE does this in effect by
  259. translating the expression to prefix form, and recalculating the value.
  260. A Fourier series has a large number of coefficients, and so this
  261. operations are repeated rather too often. At present this is the
  262. largest problem area with the internal code, as will be seen in the
  263. discussion of the Disturbing Function calculation.
  264. \section{Integration with REDUCE}
  265. The Fourier module needs to be seen as part of REDUCE rather than as a
  266. separate language. This can be seen as having internal and external
  267. parts.
  268. \subsection{Internal Interface}
  269. The Fourier expressions need to co-exist with the normal REDUCE syntax
  270. and semantics. The prototype version does this by (ab)using the
  271. module method, based in part on the TPS code \cite{Barnes}. Of course
  272. Fourier series are not constant, and so are not really domain
  273. elements. However by asserting that Fourier series form a ring of
  274. constants REDUCE can arrange to direct basic operations to the Fourier
  275. code for addition, subtraction, multiplication and the like.
  276. The main interface which needs to be provided is a simplification
  277. function for Fourier expressions. This needs to provide compilation
  278. for linear sums of angles, as well as constructing sine and cosine
  279. functions, and creating canonical forms.
  280. \subsection{User Interface}
  281. The creation of {\tt HDIFF} and {\tt HINT} functions for
  282. differentiation disguises this. An unsatisfactory aspect of the
  283. interface is that the tokens {\tt SIN} and {\tt COS} are already in
  284. use. The prototype uses the operator form
  285. \begin{verbatim}
  286. fourier sin(u)
  287. \end{verbatim}
  288. to introduce harmonically represented sine functions. An alternative of
  289. using the tokens {\tt F\_SIN} and {\tt F\_COS} is also available.
  290. It is necessary to declare the names of the angles, which is achieved
  291. with the declaration
  292. \begin{verbatim}
  293. harmonic theta, phi;
  294. \end{verbatim}
  295. At present there is no protection against using a variable as both an
  296. angle and a polynomial varaible. This will nooed to be done in a
  297. user-oriented version.
  298. \section{The Simple Experiments}
  299. The REDUCE test file contains a simple example of a Fourier
  300. calculation, determining the value of $(a_1 \cos({wt}) + a_3
  301. \cos(3{wt}) + b_1 \sin({wt}) + b_3 \sin(3{wt}))^3$. For the purposes
  302. of this system this is too trivial to do more than confirm the correct
  303. answers.
  304. The simplest non-trivial calculation for a Fourier series manipulator
  305. is to solve Kepler's equation for the eccentric anomoly E in terms of
  306. the mean anomoly u, and the eccentricity of an orbit e, considered as a
  307. small quantity
  308. \[
  309. E = u + e \sin E
  310. \]
  311. The solution procedes by repeated approximation. Clearly the initial
  312. approximation is $E_0 = u$. The $n^{th}$ approximation can be written
  313. as $u + A_n$, and so $A_n$ can be calculated by
  314. \[
  315. A_k = e \sin (u + A_{k-1})
  316. \]
  317. This is of course precisely the case for which the HSUB operation is
  318. designed, and so in order to calculate $E_n - u$ all one requires is
  319. the code
  320. \begin{verbatim}
  321. bige := fourier 0;
  322. for k:=1:n do <<
  323. wtlevel k;
  324. bige:=fourier e * hsub(fourier(sin u), u, u, bige, k);
  325. >>;
  326. write "Kepler Eqn solution:", bige$
  327. \end{verbatim}
  328. It is possible to create a regular REDUCE program to simulate this (as
  329. is done for example in Barton and Fitch\cite{Barton72}, page 254).
  330. Comparing these two programs indicates substantial advantages to the
  331. Fourier module, as could be expected.
  332. \medskip
  333. \begin{center}
  334. \begin{tabular}{ | c | l l |}
  335. \multicolumn{3}{c}{\bf Solving Kepler's Equation} \\
  336. \hline
  337. Order & REDUCE & Fourier Module \\
  338. 5 & 9.16 & 2.48 \\
  339. 6 & 17.40 & 4.56 \\
  340. 7 & 33.48 & 8.06 \\
  341. 8 & 62.76 & 13.54 \\
  342. 9 & 116.06 & 21.84 \\
  343. 10 & 212.12 & 34.54 \\
  344. 11 & 381.78 & 53.94 \\
  345. 12 & 692.56 & 82.96 \\
  346. 13 & 1247.54 & 125.86 \\
  347. 14 & 2298.08 & 187.20 \\
  348. 15 & 4176.04 & 275.60 \\
  349. 16 & 7504.80 & 398.62 \\
  350. 17 & 13459.80 & 569.26 \\
  351. 18 & *** & 800.00 \\
  352. 19 & *** & 1116.92 \\
  353. 20 & *** & 1536.40 \\
  354. \hline
  355. \end{tabular}
  356. \end{center}
  357. \medskip
  358. These results were with the linear representation of Fourier series.
  359. The tree representation was slightly slower. The ten-fold speed-up
  360. for the 13th order is most satisfactory.
  361. \section{A Medium-Sized Problem}
  362. Fourier series manipulators are primarily designed for large-scale
  363. calculations, but for the demonstration purposes of this project a
  364. medium problem is considered. The first stage in calculating the
  365. orbit of the Moon using the Delaunay theory (of perturbed elliptic
  366. motion for the restricted 3-body problem) is to calculate the energy
  367. of the Moon's motion about the Earth --- the Hamiltonian of the
  368. system. This is the calculation we use for comparisons.
  369. \subsection{Mathematical Background}
  370. The full calculation is described in detail in \cite{Brown}, but a
  371. brief description is given here for completeness, and to grasp the
  372. extent of the calculation.
  373. Referring to the figure 1 which gives the cordinate system, the basic
  374. equations are
  375. \begin{eqnarray}
  376. S & = & (1-\gamma ^2)\cos(f + g +h -f' -g' -h')
  377. + \gamma ^2 cos(f + g -h +f' +g' +h') \\
  378. r & = & a (1 - e \cos E) \\
  379. l & = & E - e \sin E \\
  380. a & = & r {{\bf d} E} \over {{\bf d} l} \\
  381. r ^2 {{\bf d} f} \over {{\bf d} l} & = & a^2 (1 - e^2)^{1 \over 2}\\
  382. R & = & m' {a^2 \over {a'} ^3} {{a'}\over {r
  383. '}} \left \{ \left ({r \over a}\right )^2
  384. \left ({{a'} \over {r'}}\right )^2 P_2(S) +
  385. \left ({a \over {a'}}\right )\left
  386. ({r \over a}\right )^3 \left ({{a'} \over {r'}}\right )^3 P_3(S)
  387. + \ldots \right \}
  388. \end{eqnarray}
  389. There are similar equations to (7) to (10) for the quantities $r'$,
  390. $a'$, $e'$, $l'$, $E'$ and $f'$ which refer to the position of the Sun
  391. rather than the Moon. The problem is to calculate the expression $R$
  392. as an expansion in terms of the quantities $e$, $e'$, $\gamma$,
  393. $a/a'$, $l$, $g$, $h$, $l'$, $g'$ and $h'$. The first three
  394. quantities are small quantities of the first order, and $a/a'$ is of
  395. second order.
  396. The steps required are
  397. \begin{enumerate}
  398. \item Solve the Kepler equation (8)
  399. \item Substiture into (7) to give $r/a$ in terms of $e$ and $l$.
  400. \item Calculate $a/r$ from (9) and $f$ from (10)
  401. \item Substitute for $f$ and $f'$ into $S$ using (6)
  402. \item Calculate $R$ from $S$, $a'/r'$ and $r/a$
  403. \end{enumerate}
  404. The program is given in the Appendix.
  405. \subsection{Results}
  406. The Lunar Disturbing function was calculated by a direct coding of the
  407. previous sections' mathematics. The program was taken from Barton
  408. and Fitch \cite{Barton72} with just small changes to generalise it for
  409. any order, and to make it acceptable for Reduce3.4. The Fourier
  410. program followed the same pattern, but obviously used the {\tt HSUB}
  411. operation as appropriate and the harmonic integration. It is very
  412. similar to the CAMAL program in \cite{Barton72}.
  413. The disturbing function was calculated to orders 2, 4 and 6 using
  414. Cambridge LISP on an HLH Orion 1/05 (Intergraph Clipper), with the
  415. three programs $\alpha$) Reduce3.4, $\beta$) Reduce3.4 + Camal Linear
  416. Module and $\gamma$) Reduce3.4 + Camal AVL Module. The timings for
  417. CPU seconds (excluding garbage collection time) are summarised the
  418. following table:
  419. \medskip
  420. \begin{center}
  421. \begin{tabular}{ | c || l | l | l |}
  422. \hline
  423. Order of DDF & Reduce & Camal Linear & Camal Tree \\
  424. \hline
  425. 2 & 23.68 & 11.22 & 12.9 \\
  426. 4 & 429.44 & 213.56 & 260.64 \\
  427. 6 & $>$7500 & 3084.62 & 3445.54 \\
  428. \hline
  429. %%% Linear n=4 138.72 (4Mb + unsafe vector access + recurrance)
  430. %%% Linear n=6 1870.10 (4Mb + unsafe vector access + recurrance)
  431. \end{tabular}
  432. \end{center}
  433. \medskip
  434. If these numbers are normalised so REDUCE calculating the DDF is 100
  435. units for each order the table becomes
  436. \medskip
  437. \begin{center}
  438. \begin{tabular}{ | c || l | l | l |}
  439. \hline
  440. Order of DDF & Reduce & Camal Linear & Camal Tree \\ \hline
  441. 2 & 100 & 47.38 & 54.48 \\
  442. 4 & 100 & 49.73 & 60.69 \\
  443. 6 & 100 & $<$41.13 & $<$45.94 \\
  444. \hline
  445. \end{tabular}
  446. \end{center}
  447. \medskip
  448. From this we conclude that a doubling of speed is about correct, and
  449. although the balanced tree system is slower as the problem size
  450. increases the gap between it and the simpler linear system is
  451. narrowing.
  452. It is disappointing that the ratio is not better, nor the absolute
  453. time less. It is worth noting in this context that Jefferys claimed
  454. that the sixth order DDF took 30s on a CDC6600 with TRIGMAN in 1970
  455. \cite{Jefferys}, and Barton and Fitch took about 1s for the second
  456. order DDF on TITAN with CAMAL \cite{Barton72}. A closer look at the
  457. relative times for individual sections of the program shows that the
  458. substitution case of replacing a polynomial variable by a Fourier
  459. series is only marginally faster than the simple REDUCE program. In
  460. the DDF program this operation is only used once in a major form,
  461. substituting into the Legendre polynomials, which have been previously
  462. calculated by Rodrigues formula. This suggests that we replace this
  463. with the recurrence relationship.
  464. Making this change actually slows down the normal REDUCE by a small
  465. amount but makes a significant change to the Fourier module; it
  466. reduces the run time for the 6th order DDF from 3084.62s to 2002.02s.
  467. This gives some indication of the problems with benchmarks. What is
  468. clear is that the current implementation of substitution of a Fourier
  469. series for a polynomial variable is inadequate.
  470. \section{Conclusion}
  471. The Fourier module is far from complete. The operations necessary for
  472. the solution of Duffing's and Hill's equations are not yet written,
  473. although they should not cause much problem. The main defficiency is
  474. the treatment of series truncation; at present it relies on the REDUCE
  475. WTLEVEL mechanism, and this seems too coarse for efficient truncation.
  476. It would be possible to re-write the polynomial manipulator as well,
  477. while retaining the REDUCE syntax, but that seems rather more than one
  478. would hope.
  479. The real failure so far is the large time lag between the REDUCE-based
  480. system on a modern workstation against a mainframe of 25 years ago
  481. running a special system. The CAMAL Disturbing function program could
  482. calculate the tenth order with a maximum of 32K words (about
  483. 192Kbytes) whereas this system failed to calculate the eigth order in
  484. 4Mbytes (taking 2000s before failing). I have in my archives the
  485. output from the standard CAMAL test suite, which includes a sixth
  486. order DDF on an IBM 370/165 run on 2 June 1978, taking 22.50s and
  487. using a maximum of 15459 words of memory for heap --- or about
  488. 62Kbytes. A rough estimate is that the Orion 1/05 is comparable in
  489. speed to the 360/165, but with more real memory and virtual memory.
  490. However, a simple Fourier manipulator has been created for REDUCE which
  491. performs between twice and three times the speed of REDUCE using
  492. pattern matching. It has been shown that this system is capable of
  493. performing the calculations of celestial mechanics, but it still
  494. seriously lags behind the efficiency of the specialist systems of
  495. twenty years before. It is perhaps fortunate that it was not been
  496. possible to compare it with a modern specialist system.
  497. There is still work to do to provide a convenient user interface, but
  498. it is intended to develop the system in this direction. It would be
  499. pleasant to have again a system of the efficiency of CAMAL(F).
  500. I would like to thank Codemist Ltd for the provision of computing
  501. resources for this project, and David Barton who taught be so much
  502. about Fourier series and celstial mechanics. Thank are also due to
  503. the National Health Service, without whom this work and paper could not
  504. have been produced.
  505. \section*{Appendix: The DDF Function}
  506. \begin{verbatim}
  507. array p(n/2+2);
  508. harmonic u,v,w,x,y,z;
  509. weight e=1, b=1, d=1, a=1;
  510. %% Generate Legendre Polynomials to sufficient order
  511. for i:=2:n/2+2 do <<
  512. p(i):=(h*h-1)^i;
  513. for j:=1:i do p(i):=df(p(i),h)/(2j)
  514. >>;
  515. %%%%%%%%%%%%%%%% Step1: Solve Kepler equation
  516. bige := fourier 0;
  517. for k:=1:n do <<
  518. wtlevel k;
  519. bige:=fourier e * hsub(fourier(sin u), u, u, bige, k);
  520. >>;
  521. %% Ensure we do not calculate things of too high an order
  522. wtlevel n;
  523. %%%%%%%%%%%%%%%% Step 2: Calculate r/a in terms of e and l
  524. dd:=-e*e; hh:=3/2; j:=1; cc := 1;
  525. for i:=1:n/2 do <<
  526. j:=i*j; hh:=hh-1; cc:=cc+hh*(dd^i)/j
  527. >>;
  528. bb:=hsub(fourier(1-e*cos u), u, u, bige, n);
  529. aa:=fourier 1+hdiff(bige,u); ff:=hint(aa*aa*fourier cc,u);
  530. %%%%%%%%%%%%%%%% Step 3: a/r and f
  531. uu := hsub(bb,u,v); uu:=hsub(uu,e,b);
  532. vv := hsub(aa,u,v); vv:=hsub(vv,e,b);
  533. ww := hsub(ff,u,v); ww:=hsub(ww,e,b);
  534. %%%%%%%%%%%%%%%% Step 4: Substitute f and f' into S
  535. yy:=ff-ww; zz:=ff+ww;
  536. xx:=hsub(fourier((1-d*d)*cos(u)),u,u-v+w-x-y+z,yy,n)+
  537. hsub(fourier(d*d*cos(v)),v,u+v+w+x+y-z,zz,n);
  538. %%%%%%%%%%%%%%%% Step 5: Calculate R
  539. zz:=bb*vv; yy:=zz*zz*vv;
  540. on fourier;
  541. for i := 2:n/2+2 do <<
  542. wtlevel n+4-2i; p(i) := hsub(p(i), h, xx) >>;
  543. wtlevel n;
  544. for i:=n/2+2 step -1 until 3 do
  545. p(n/2+2):=fourier(a*a)*zz*p(n/2+2)+p(i-1);
  546. yy*p(n/2+2);
  547. \end{verbatim}
  548. \newpage
  549. \bibliographystyle{plain}
  550. \bibliography{camal}
  551. \end{document}