Merge Eigenmath modules

All eigenmath modules were imported using:

	git cvsimport -vik -C /tmp/eigenmath/modules/$MODULE -A /tmp/cvsuser.txt \
	-d :pserver:anonymous@eigenmath.cvs.sourceforge.net:/cvsroot/eigenmath $MODULE

/tmp/cvsuser.txt contained:

	gweigt=George Weigt <gweigt@users.sourceforge.net> Etc/GMT+7

The files are copied into subdirectories and the commits
manually added as parents.

George Weigt's last commits of the individual modules:

7fea6f26a2 c
44e026d257 eigenmath
29403fc655 lisp
f053a059fe man
74a1f5ec30 ross
804bd98a9c tex

doc/manual/16.tex

@@ -0,0 +1,37 @@
+
+\newpage
+
+\subsection{}
+Consider the graph of
+$$y={1\over\sqrt{1+x^2}}$$
+over the domain from 0 to infinity.
+The following graph demonstrates the general idea.
+
+\medskip
+\verb$xrange=(0,10)$
+
+\verb$yrange=(-1,1)$
+
+\verb$draw(1/sqrt(1+x^2))$
+
+\begin{center}
+\includegraphics[scale=0.4]{16.png}
+\end{center}
+
+\medskip
+\noindent
+Now imagine that the graph is rotated about the $x$ axis.
+What is the volume subtended by the area under the curve?
+
+\medskip
+\noindent
+To solve, use cylindrical coordinates.
+The $x$ axis becomes the $z$ axis.
+We want $x$ to become a function of $y$ so invert the function.
+$$y={1\over\sqrt{1+x^2}}
+\quad\rightarrow\quad
+{1\over y^2}=1+x^2
+\quad\rightarrow\quad
+x=\sqrt{{1\over y^2}-1}
+$$
+

doc/manual/Eigenmath.tex

@@ -0,0 +1,74 @@
+\documentclass[11pt]{article}
+\pagestyle{myheadings}
+
+\usepackage{graphicx}
+\usepackage{makeidx}
+
+
+\title{Eigenmath Manual}
+%\author{George Weigt}
+\date{\today}
+
+%\makeindex
+
+\begin{document}
+
+\maketitle
+
+\newpage
+
+\tableofcontents
+
+\newpage
+
+\include{nabokov}
+
+\include{symbols}
+
+%\include{geometric-series}
+
+%\include{units-of-measure}
+
+\include{draw}
+
+\include{scripting}
+
+\include{complex}
+
+\include{linear-algebra}
+
+
+
+
+\include{derivative}
+\include{integrals}
+\include{integral-trick}
+\include{fund-thm-of-calculus}
+\include{line-integral}
+\include{surface-area}
+\include{surface-integral}
+\include{greens-theorem}
+\include{stokes-theorem}
+
+
+
+\include{francois-viete}
+\include{curl-in-tensor-form}
+\include{qho}
+\include{hydrogen-wavefunctions}
+\include{space-shuttle-and-corvette}
+\include{avogadro}
+\include{zerozero}
+%\include{16}
+
+\newpage
+
+\include{list-of-functions}
+
+\newpage
+
+\include{syntax}
+
+%\printindex
+
+\end{document}

doc/manual/Makefile

@@ -0,0 +1,6 @@
+all :
+	pdflatex Eigenmath
+	pdflatex Eigenmath
+	pdflatex Eigenmath
+	makeindex Eigenmath
+	pdflatex Eigenmath

doc/manual/avogadro.tex

@@ -0,0 +1,37 @@
+\subsection{Avogadro's constant}
+
+There is a proposal to define Avogadro's constant as exactly
+84446886 to the third power.\footnote{Fox, Ronald and Theodore Hill.
+``An Exact Value for Avogadro's Number.''
+{\it American Scientist} 95 (2007): 104--107.
+The proposed number in the article is actually $84446888^3$.
+In a subsequent addendum the authors reduced it to $84446886^3$ to make the
+number divisible by 12. See {\tt www.physorg.com/news109595312.html}}
+This number corresponds to an ideal cube of atoms with 84,446,886
+atoms along each edge.
+Let us check the difference between the proposed value and the measured value
+of $(6.0221415\pm0.0000010)\times10^{23}$ atoms.
+
+\medskip
+\verb$A=84446886^3$
+
+\verb$B=6.0221415*10^23$
+
+\verb$A-B$
+
+$$-5.17173\times10^{16}$$
+
+\verb$0.0000010*10^23$
+
+$$1\times10^{17}$$
+
+\medskip
+\noindent
+We see that the proposed value is within the experimental error.
+Just for the fun of it, let us factor the proposed value.
+
+\medskip
+\verb$factor(A)$
+
+$$2^3\times3^3\times1667^3\times8443^3$$
+

doc/manual/complex-numbers.tex

@@ -0,0 +1,5 @@
+\documentclass{article}
+\pagestyle{empty}
+\begin{document}
+\include{complex}
+\end{document}

doc/manual/complex.tex

@@ -0,0 +1,95 @@
+\section{Complex numbers}
+When Eigenmath starts up, it defines the symbol $i$ as $i=\sqrt{-1}$.
+Other than that, there is nothing special about $i$.
+It is just a regular symbol that can be redefined and used for some other purpose if need be.
+
+\medskip
+\noindent
+Complex quantities can be entered in rectangular or polar form.
+
+\medskip
+\verb$a+i*b$
+$$a+ib$$
+
+\verb$exp(i*pi/3)$
+$$\exp({1\over3}i\pi)$$
+
+\medskip
+\noindent
+Converting to rectangular or polar coordinates simplifies mixed forms.
+
+\medskip
+\verb$A=1+i$
+
+\verb$B=sqrt(2)*exp(i*pi/4)$
+
+\verb$A-B$
+$$1+i-2^{1/2}\exp({1\over4}i\pi)$$
+
+\verb$rect$
+$$0$$
+
+\medskip
+\noindent
+Rectangular complex quantities, when raised to a power, are multiplied out.
+
+\medskip
+\verb$(a+i*b)^2$
+$$a^2-b^2+2iab$$
+
+\medskip
+\noindent
+When $a$ and $b$ are numerical and the power is negative, the evaluation is done as follows.
+$$(a+ib)^{-n}=\left[{a-ib\over(a+ib)(a-ib)}\right]^n=\left[{a-ib\over a^2+b^2}\right]^n$$
+Of course, this causes $i$ to be removed from the denominator.
+%For $n=1$ we have
+%$${1\over a+ib}={a-ib\over a^2+b^2}$$
+Here are a few examples.
+
+\medskip
+\verb$1/(2-i)$
+$${2\over5}+{1\over5}i$$
+
+\verb$(-1+3i)/(2-i)$
+$$-1+i$$
+
+\newpage
+
+\noindent
+The absolute value of a complex number returns its magnitude.
+
+\medskip
+\verb$abs(3+4*i)$
+$$5$$
+
+\medskip
+\noindent
+In light of this, the following result might be unexpected.
+
+\medskip
+\verb$abs(a+b*i)$
+$$\mathop{\rm abs}(a+ib)$$
+
+\medskip
+\noindent
+The result is not $\sqrt{a^2+b^2}$ because that would assume that
+$a$ and $b$ are real.
+For example, suppose that $a=0$ and $b=i$.
+Then
+$$|a+ib|=|-1|=1$$
+and
+$$\sqrt{a^2+b^2}=\sqrt{-1}=i$$
+Hence
+$$|a+ib|\ne\sqrt{a^2+b^2}\quad\hbox{for some $a,b\in\cal C$}$$
+
+\medskip
+\noindent
+The {\it mag} function is an alternative.
+It treats symbols like $a$ and $b$ as real.
+
+\medskip
+\verb$mag(a+b*i)$
+
+$$(a^2+b^2)^{1/2}$$
+
+

doc/manual/curl-in-tensor-form.tex

@@ -0,0 +1,94 @@
+\subsection{Curl in tensor form}
+The curl of a vector function can be expressed in tensor form as
+$$\mathop{\rm curl}{\bf F}=\epsilon_{ijk}\,{\partial F_k\over\partial x_j}$$
+where $\epsilon_{ijk}$ is the Levi-Civita tensor.
+The following script demonstrates that this formula is equivalent
+to computing curl the old fashioned way.
+
+\medskip
+\verb$-- Define epsilon.$
+
+\verb$epsilon=zero(3,3,3)$
+
+\verb$epsilon[1,2,3]=1$
+
+\verb$epsilon[2,3,1]=1$
+
+\verb$epsilon[3,1,2]=1$
+
+\verb$epsilon[3,2,1]=-1$
+
+\verb$epsilon[1,3,2]=-1$
+
+\verb$epsilon[2,1,3]=-1$
+
+\verb$-- F is a generic vector function.$
+
+\verb$F=(FX(),FY(),FZ())$
+
+\verb$-- A is the curl of F.$
+
+\verb$A=outer(epsilon,d(F,(x,y,z)))$
+
+\verb$A=contract(A,3,4) --sum across k$
+
+\verb$A=contract(A,2,3) --sum across j$
+
+\verb$-- B is the curl of F computed the old fashioned way.$
+
+\verb$BX=d(F[3],y)-d(F[2],z)$
+
+\verb$BY=d(F[1],z)-d(F[3],x)$
+
+\verb$BZ=d(F[2],x)-d(F[1],y)$
+
+\verb$B=(BX,BY,BZ)$
+
+\verb$-- Are A and B equal? Subtract to find out.$
+
+\verb$A-B$
+
+$$\left(\matrix{0\cr0\cr0}\right)$$
+
+\newpage
+
+\noindent
+The following is a variation on the previous script.
+The product $\epsilon_{ijk}\,\partial F_k/\partial x_j$
+is computed in just one line of code.
+In addition, the outer product and the contraction across $k$
+are now computed with a dot product.
+
+\medskip
+\verb$F=(FX(),FY(),FZ())$
+
+\verb$epsilon=zero(3,3,3)$
+
+\verb$epsilon[1,2,3]=1$
+
+\verb$epsilon[2,3,1]=1$
+
+\verb$epsilon[3,1,2]=1$
+
+\verb$epsilon[3,2,1]=-1$
+
+\verb$epsilon[1,3,2]=-1$
+
+\verb$epsilon[2,1,3]=-1$
+
+\verb$A=contract(dot(epsilon,d(F,(x,y,z))),2,3)$
+
+\verb$BX=d(F[3],y)-d(F[2],z)$
+
+\verb$BY=d(F[1],z)-d(F[3],x)$
+
+\verb$BZ=d(F[2],x)-d(F[1],y)$
+
+\verb$B=(BX,BY,BZ)$
+
+\verb$-- Are A and B equal? Subtract to find out.$
+
+\verb$A-B$
+
+$$\left(\matrix{0\cr0\cr0}\right)$$
+

doc/manual/derivative.tex

@@ -0,0 +1,104 @@
+\section{Calculus}
+\subsection{Derivative}
+\index{derivative}
+$d(f,x)$ returns the derivative of $f$ with respect to $x$.
+The $x$ can be omitted for expressions in $x$.
+
+\medskip
+\verb$d(x^2)$
+$$2x$$
+
+\bigskip
+\noindent
+The following table summarizes the various ways to obtain multiderivatives.
+
+\begin{center}
+\begin{tabular}{cllllll}
+%& & & & {\it alternate form} \\
+%\\
+$\displaystyle{\partial^2f\over\partial x^2}$ & & \verb$d(f,x,x)$ & & \verb$d(f,x,2)$ \\
+\\
+$\displaystyle{\partial^2f\over\partial x\,\partial y}$ & & \verb$d(f,x,y)$ \\
+\\
+$\displaystyle{\partial^{m+n+\cdot\cdot\cdot} f\over\partial x^m\,\partial y^n\cdots}$ & &
+\verb$d(f,x,...,y,...)$ & & \verb$d(f,x,m,y,n,...)$ \\
+\end{tabular}
+\end{center}
+
+%\medskip
+%\verb$r=sqrt(x^2+y^2)$
+
+%\verb$d(r,x,y)$
+%$${-{xy\over(x^2+y^2)^{3/2}}}$$
+
+\subsection{Gradient}
+\index{gradient}
+
+\noindent
+The gradient of $f$ is obtained by using a vector for $x$ in $d(f,x)$.
+
+\medskip
+\verb$r=sqrt(x^2+y^2)$
+
+\verb$d(r,(x,y))$
+$$\left(\matrix{
+\displaystyle{{x\over(x^2+y^2)^{1/2}}}\cr
+\cr
+\displaystyle{{y\over(x^2+y^2)^{1/2}}}\cr
+}\right)$$
+
+\medskip
+\noindent
+The $f$ in $d(f,x)$ can be a tensor function.
+Gradient raises the rank by one.
+
+\medskip
+\verb$F=(x+2y,3x+4y)$
+
+\verb$X=(x,y)$
+
+\verb$d(F,X)$
+$$\left(\matrix{1&2\cr3&4}\right)$$
+
+\subsection{Template functions}
+The function $f$ in $d(f)$ does not have to be defined.
+It can be a template function with just a name and an argument list.
+Eigenmath checks the argument list to figure out what to do.
+For example, $d(f(x),x)$ evaluates to itself because $f$ depends on $x$.
+However, $d(f(x),y)$ evaluates to zero because $f$ does not depend on $y$.
+
+\medskip
+\verb$d(f(x),x)$
+$$\partial(f(x),x)$$
+
+\verb$d(f(x),y)$
+$$0$$
+
+\verb$d(f(x,y),y)$
+$$\partial(f(x,y),y)$$
+
+\verb$d(f(),t)$
+$$\partial(f(),t)$$
+
+\medskip
+\noindent
+As the final example shows, an empty argument list causes
+$d(f)$ to always evaluate to itself, regardless
+of the second argument.
+
+\medskip
+\noindent
+Template functions are useful for experimenting with differential forms.
+For example, let us check the identity
+$$\mathop{\rm div}(\mathop{\rm curl}{\bf F})=0$$
+for an arbitrary vector function $\bf F$.
+
+\medskip
+\verb$F=(F1(x,y,z),F2(x,y,z),F3(x,y,z))$
+
+\verb$curl(U)=(d(U[3],y)-d(U[2],z),d(U[1],z)-d(U[3],x),d(U[2],x)-d(U[1],y))$
+
+\verb$div(U)=d(U[1],x)+d(U[2],y)+d(U[3],z)$
+
+\verb$div(curl(F))$
+$$0$$

doc/manual/draw.tex

@@ -0,0 +1,100 @@
+\section{Draw}
+$draw(f,x)$ draws a graph of the function $f$ of $x$.
+The second argument can be omitted when the dependent variable
+is literally $x$ or $t$.
+The vectors $xrange$ and $yrange$ control the scale of the graph.
+
+\medskip
+\verb$draw(x^2)$
+
+\medskip
+\begin{center}
+\includegraphics[scale=0.4]{parabola.png}
+\end{center}
+
+\verb$xrange=(-1,1)$
+
+\verb$yrange=(0,2)$
+
+\verb$draw(x^2)$
+
+\medskip
+\begin{center}
+\includegraphics[scale=0.4]{parabola2.png}
+\end{center}
+
+\verb$clear$
+
+\medskip
+\noindent
+The clear command (or a click of the Clear button)
+resets $xrange$ and $yrange$.
+This needs to be done so that the next graph
+appears as shown.
+
+\newpage
+
+\noindent
+Parametric drawing occurs when a function returns a vector.
+The vector $trange$ controls the parameter range.
+The default range is $(-\pi,\pi)$.
+
+\medskip
+\verb$f=(cos(t),sin(t))$
+
+\verb$draw(5*f)$
+
+\medskip
+\begin{center}
+\includegraphics[scale=0.4]{circle.png}
+\end{center}
+
+\verb$trange=(0,pi/2)$
+
+\verb$draw(5*f)$
+
+\medskip
+\begin{center}
+\includegraphics[scale=0.4]{circle2.png}
+\end{center}
+
+\newpage
+
+\noindent
+Here are a couple of interesting curves and the code for drawing them.
+First is a lemniscate.
+
+\medskip
+\verb$clear$
+
+\verb$X=cos(t)/(1+sin(t)^2)$
+
+\verb$Y=sin(t)*cos(t)/(1+sin(t)^2)$
+
+\verb$draw(5*(X,Y))$
+
+\medskip
+\begin{center}
+\includegraphics[scale=0.4]{lemniscate.png}
+\end{center}
+
+\medskip
+\noindent
+Next is a cardioid.
+
+\medskip
+\verb$r=(1+cos(t))/2$
+
+\verb$u=(cos(t),sin(t))$
+
+\verb$xrange=(-1,1)$
+
+\verb$yrange=(-1,1)$
+
+\verb$draw(r*u)$
+
+\medskip
+\begin{center}
+\includegraphics[scale=0.4]{cardioid.png}
+\end{center}
+

doc/manual/francois-viete.tex

@@ -0,0 +1,55 @@
+\section{More examples}
+\subsection{Fran\c cois Vi\`ete}
+Fran\c cois Vi\`ete was the first to discover an exact formula for $\pi$.
+Here is his formula.
+\begin{displaymath}
+{2\over\pi}={\sqrt2\over2}\times{\sqrt{2+\sqrt2}\over2}\times
+{\sqrt{2+\sqrt{2+\sqrt2}}\over2}\times\cdots
+\end{displaymath}
+%We can flip it around and write the formula like this.
+%\begin{displaymath}
+%\pi=2\times{2\over\sqrt2}\times{2\over\sqrt{2+\sqrt2}}\times
+%{2\over\sqrt{2+\sqrt{2+\sqrt2}}}\times\cdots
+%\end{displaymath}
+Let $a_0=0$ and $a_{n}=\sqrt{2+a_{n-1}}$.
+Then we can write
+\begin{displaymath}
+{2\over\pi}={a_1\over2}\times{a_2\over2}\times
+{a_3\over2}\times\cdots
+\end{displaymath}
+%
+Solving for $\pi$ we have
+\begin{displaymath}
+\pi=2\times{2\over a_1}\times{2\over a_2}\times{2\over a_3}\times\cdots=2\prod_{k=1}^\infty
+{2\over a_k}
+\end{displaymath}
+%
+Let us now use Eigenmath to compute $\pi$ according to Vi\`ete's formula.
+Of course, we cannot calculate all the way out to infinity, we have to stop somewhere.
+It turns out that nine factors are just enough to get six digits of accuracy.
+
+\medskip
+\verb$a(n)=test(n=0,0,sqrt(2+a(n-1)))$
+
+\verb$float(2*product(k,1,9,2/a(k)))$
+
+$$3.14159$$
+
+\medskip
+\noindent
+The function $a(n)$ calls itself $n$ times so overall there are
+54 calls to $a(n)$.
+By using a different algorithm with temporary variables, we can get the
+answer in just nine steps.
+
+\medskip
+\verb$a=0$
+
+\verb$b=2$
+
+\verb$for(k,1,9,a=sqrt(2+a),b=b*2/a)$
+
+\verb$float(b)$
+
+$$3.14159$$
+

doc/manual/fund-thm-of-calculus.tex

@@ -0,0 +1,39 @@
+
+\newpage
+
+\index{fundamental theorem of calculus}
+
+\noindent
+The fundamental theorem of calculus was established by James Gregory,
+a contemporary of Newton.
+The theorem is a formal expression of the inverse relation between
+integrals and derivatives.
+$$\int_a^b f'(x)\,dx=f(b)-f(a)$$
+Here is an Eigenmath demonstration of the fundamental theorem of calculus.
+
+\medskip
+\verb$f=x^2/2$
+
+\verb$xrange=(-1,1)$
+
+\verb$yrange=xrange$
+
+\verb$draw(d(f))$
+% no blank line, otherwise goes to two pages
+\begin{center}
+\includegraphics[scale=0.4]{funda1.png}
+\end{center}
+
+\verb$draw(integral(d(f)))$
+
+\begin{center}
+\includegraphics[scale=0.4]{funda2.png}
+\end{center}
+
+\noindent
+The first graph shows that $f'(x)$ is antisymmetric, therefore the total
+area under the curve from $-1$ to $1$ sums to zero.
+The second graph shows that $f(1)=f(-1)$.
+Hence for $f(x)={1\over2}x^2$ we have
+$$\int_{-1}^1f'(x)\,dx=f(1)-f(-1)=0$$
+

doc/manual/geometric-series.tex

@@ -0,0 +1,54 @@
+\section{Details}
+\subsection{User-defined symbols}
+
+A geometric series converges according to the formula
+$$\sum_{k=0}^\infty a^k={1\over1-a},\qquad|a|<1$$
+If we use $a=-1/2$ and for practical purposes only count up to nine instead of infinity,
+we should have
+$$\sum_{k=0}^9\left(-{1\over2}\right)^k\approx{2\over3}$$
+The above calculation can be done in one line of code using the $sum$ function.
+
+\medskip
+\verb$sum(k,0,9,(-0.5)^k)$
+$$0.666016$$
+
+\medskip
+\noindent
+The following example uses an intermediate variable.
+
+\medskip
+\verb$f=sum(k,0,9,a^k)$
+
+\verb$f$
+$$f=1+a+a^2+a^3+a^4+a^5+a^6+a^7+a^8+a^9$$
+
+\verb$eval(f,a,-1/2)$
+$$341\over512$$
+
+\verb$float(last)$
+$$0.666016$$
+
+\medskip
+\noindent
+As seen on the first line, no result is printed when a symbol is defined.
+When you do in fact want to see the value of a symbol,
+just enter it as shown on the second line.
+
+\medskip
+\noindent
+When a result is displayed, it is also stored in the symbol $last$.
+
+\subsection{User-defined functions}
+
+\noindent
+The following example shows how to define a function.
+
+\medskip
+\verb$f(a)=sum(k,0,9,a^k)$
+
+\verb$f(-1/2)$
+$$341\over512$$
+
+\verb$f(-0.5)$
+$$0.666016$$
+

doc/manual/greens-theorem.tex

@@ -0,0 +1,123 @@
+\subsection{Green's theorem}
+\index{Green's theorem}
+Green's theorem tells us that
+$$\oint P\,dx+Q\,dy=\int\!\!\!\int
+\left({\partial Q\over\partial x}-{\partial P\over\partial y}\right)
+dx\,dy$$
+
+\medskip
+\noindent
+Example 1.
+Evaluate $\oint (2x^3-y^3)\,dx+(x^3+y^3)\,dy$ around the circle
+$x^2+y^2=1$ using Green's theorem.\footnote{
+Wilfred Kaplan, {\it Advanced Calculus, 5th Edition,} 287.}
+
+\medskip
+\noindent
+It turns out that Eigenmath cannot solve the double integral over
+$x$ and $y$ directly.
+Polar coordinates are used instead.
+
+\medskip
+\verb$P=2x^3-y^3$
+
+\verb$Q=x^3+y^3$
+
+\verb$f=d(Q,x)-d(P,y)$
+
+\verb$x=r*cos(theta)$
+
+\verb$y=r*sin(theta)$
+
+\verb$defint(f*r,r,0,1,theta,0,2pi)$
+
+$${3\over2}\pi$$
+
+\medskip
+\noindent
+The $defint$ integrand is $f{*}r$ because $r\,dr\,d\theta=dx\,dy$.
+
+\medskip
+\noindent
+Now let us try computing the line integral side of Green's theorem
+and see if we get the same result.
+We need to use the trick of converting sine and cosine to exponentials
+so that Eigenmath can find a solution.
+
+\medskip
+\verb$x=cos(t)$
+
+\verb$y=sin(t)$
+
+\verb$P=2x^3-y^3$
+
+\verb$Q=x^3+y^3$
+
+\verb$f=P*d(x,t)+Q*d(y,t)$
+
+\verb$f=circexp(f)$
+
+\verb$defint(f,t,0,2pi)$
+
+$${3\over2}\pi$$
+
+\newpage
+
+\noindent
+Example 2.
+Compute both sides of Green's theorem for
+$F=(1-y,x)$ over the disk $x^2+y^2\le4$.
+
+\medskip
+\noindent
+First compute the line integral along the boundary of the disk.
+Note that the radius of the disk is 2.
+
+\medskip
+\verb$--Line integral$
+
+\verb$P=1-y$
+
+\verb$Q=x$
+
+\verb$x=2*cos(t)$
+
+\verb$y=2*sin(t)$
+
+\verb$defint(P*d(x,t)+Q*d(y,t),t,0,2pi)$
+
+$$8\pi$$
+
+\verb$--Surface integral$
+
+\verb$x=quote(x) --remove parametrization of x$
+
+\verb$y=quote(y) --remove parametrization of y$
+
+\verb$h=sqrt(4-x^2)$
+
+\verb$defint(d(Q,x)-d(P,y),y,-h,h,x,-2,2)$
+
+$$8\pi$$
+
+\verb$--Bonus point: Compute the surface integral using polar coordinates.$
+
+\verb$f=d(Q,x)-d(P,y) --do before change of coordinates$
+
+\verb$x=r*cos(theta)$
+
+\verb$y=r*sin(theta)$
+
+\verb$defint(f*r,r,0,2,theta,0,2pi)$
+
+$$8\pi$$
+
+\verb$defint(f*r,theta,0,2pi,r,0,2) --try integrating over theta first$
+
+$$8\pi$$
+
+\medskip
+\noindent
+In this case, Eigenmath solved both forms of the polar integral.
+However, in cases where Eigenmath fails to solve a double integral, try
+changing the order of integration.

doc/manual/how-it-works.tex

@@ -0,0 +1,177 @@
+\chapter{How it works}
+
+Eigenmath is written in the C programming language.
+Its central data structure is
+an array of 100 thousand blocks of memory.
+The size of each block is 12 bytes for Mac OS X
+and 16 bytes for Windows.
+\begin{verbatim}
+         _______
+      0 |_______|  \
+      1 |_______|   |
+      2 |_______|   |
+        |       |   |
+        .       .   | 1,200,000 bytes (Mac OS X)
+        .       .   |
+        .       .   |
+        |_______|   |
+ 99,999 |_______|  /
+\end{verbatim}
+
+\bigskip
+\noindent
+If necessary, Eigenmath will allocate additional memory in increments of
+100 thousand blocks, up to a maximum of 10 million blocks
+(one hundred times what is shown above).
+
+%\begin{verbatim}
+%100000 blocks (12 bytes/block)
+%99706 free
+%294 used
+%\end{verbatim}
+
+\newpage
+
+\noindent
+Blocks are C structs defined as follows.
+
+\medskip
+\begin{verbatim}
+typedef struct U {
+     union {
+          struct {
+               struct U *car;          // pointing down
+               struct U *cdr;          // pointing right
+          } cons;
+          char *printname;
+          char *str;
+          struct tensor *tensor;
+          struct {
+               unsigned int *a, *b;    // rational number a over b
+          } q;
+          double d;
+     } u;
+     unsigned char k, tag;
+} U;
+\end{verbatim}
+
+\medskip
+\noindent
+The member $k$ identifies the union that is stored in the block.
+The value of $k$ is one of the following.
+
+\medskip
+\begin{verbatim}
+enum {
+        CONS,
+        NUM,
+        DOUBLE,
+        STR,
+        TENSOR,
+        SYM,
+};
+\end{verbatim}
+
+\newpage
+
+\noindent
+Blocks are linked together to store mathematical expressions.
+For example, the following shows how the expression
+$A\times B+C$ is stored.
+\begin{verbatim}
+ _______      _______                                _______
+|CONS   |--->|CONS   |----------------------------->|CONS   |
+|       |    |       |                              |       |
+|_______|    |_______|                              |_______|
+    |            |                                      |
+ ___V___      ___V___      _______      _______      ___V___
+|SYM    |    |CONS   |--->|CONS   |--->|CONS   |    |SYM    |
+|add    |    |       |    |       |    |       |    |C      |
+|_______|    |_______|    |_______|    |_______|    |_______|
+                 |            |            |
+              ___V___      ___V___      ___V___
+             |SYM    |    |SYM    |    |SYM    |
+             |times  |    |A      |    |B      |
+             |_______|    |_______|    |_______|
+\end{verbatim}
+
+\bigskip
+\noindent
+Only CONS blocks contain pointers to other blocks.
+Every other kind of block is a terminal node.
+Fundamentally, this approach is a subset of the S-expression data structure
+invented by John McCarthy for the LISP programming language.
+
+\newpage
+
+\noindent
+Confining pointers to CONS blocks differs from the more traditional linked
+list data structure in a significant way.
+In a linked list, blocks contain both data and pointers simultaneously.
+For example, this is how one might store $A+B$ using a linked list.
+\begin{verbatim}
+ _______       _______       _______
+|_______|---->|_______|---->|_______|
+|SYM    |     |SYM    |     |SYM    |
+|add    |     |A      |     |B      |
+|_______|     |_______|     |_______|
+\end{verbatim}
+
+\bigskip
+\noindent
+Now suppose we want to store an additional expression $A+C$.
+Using a linked list, we have to make copies of the SYM-add and
+SYM-A blocks
+because of the pointers.
+If we were to rewrite the pointers we would destroy the original
+expression.
+However, with S-expressions all we have to do is allocate new CONS
+blocks to create new expressions.
+We never have to copy data.
+We can store thousands of expressions with CONS
+blocks all pointing to a common SYM-A block.
+
+\newpage
+
+\noindent
+As it runs, Eigenmath allocates memory blocks to build new expressions.
+For example, $A+A$ is changed to $2A$ by $evel\_add$.
+Here is the input expression to $eval\_add$.
+\begin{verbatim}
+ _______      _______      _______
+|CONS   |--->|CONS   |--->|CONS   |
+|       |    |       |    |       |
+|_______|    |_______|    |_______|
+    |            |            |
+ ___V___         |         ___V___
+|SYM    |        |------->|SYM    |
+|add    |                 |A      |
+|_______|                 |_______|
+\end{verbatim}
+
+\bigskip
+\noindent
+And here is the output expression from $eval\_add$.
+\begin{verbatim}
+ _______      _______      _______
+|CONS   |--->|CONS   |--->|CONS   |
+|       |    |       |    |       |
+|_______|    |_______|    |_______|
+    |            |            |
+ ___V___      ___V___      ___V___
+|SYM    |    |NUM    |    |SYM    |
+|times  |    |2      |    |A      |
+|_______|    |_______|    |_______|
+\end{verbatim}
+
+\bigskip
+\noindent
+The three CONS blocks in the output expression are new.
+They are not the same as the CONS blocks in the input expression.
+In this particular example, Eigenmath allocates a total of four blocks
+to build the new expression,
+three CONS blocks and one NUM block.
+Eventually all of the available blocks get used up, at which point Eigenmath
+does garbage collection to recover unused blocks.
+Returning to the example, garbage collection would recover the original
+three CONS blocks that formed $A+A$, if nothing points to that expression.

doc/manual/hydrogen-wavefunctions.tex

@@ -0,0 +1,53 @@
+\subsection{Hydrogen wavefunctions}
+\index{hydrogen wavefunctions}
+Hydrogen wavefunctions $\psi$ are solutions to the differential equation
+$${\psi\over n^2}=\nabla^2\psi+{2\psi\over r}$$
+where $n$ is an integer representing the quantization of total energy and
+$r$ is the radial distance of the electron.
+The Laplacian operator in spherical coordinates is
+$$\nabla^2={1\over r^2}{\partial\over\partial r}
+\left(r^2{\partial\over\partial r}\right)
++{1\over r^2\sin\theta}{\partial\over\partial\theta}
+\left(\sin\theta{\partial\over\partial\theta}\right)
++{1\over r^2\sin^2\theta}{\partial^2\over\partial\phi^2}$$
+The general form of $\psi$ is
+$$\psi=r^le^{-r/n}L_{n-l-1}^{2l+1}(2r/n)
+P_l^{|m|}(\cos\theta)e^{im\phi}$$
+where $L$ is a Laguerre polynomial, $P$ is a Legendre polynomial and
+$l$ and $m$ are integers such that
+$$1\le l\le n-1,\qquad -l\le m\le l$$
+The general form can be expressed as the product of a radial
+wavefunction $R$ and a spherical harmonic $Y$.
+$$\psi=RY,\qquad R=r^le^{-r/n}L_{n-l-1}^{2l+1}(2r/n),\qquad
+Y=P_l^{|m|}(\cos\theta)e^{im\phi}$$
+The following code checks $E=K+V$ for $n,l,m=7,3,1$.
+
+\medskip
+\verb$laplacian(f)=1/r^2*d(r^2*d(f,r),r)+$
+
+\verb$1/(r^2*sin(theta))*d(sin(theta)*d(f,theta),theta)+$
+
+\verb$1/(r*sin(theta))^2*d(f,phi,phi)$
+
+\verb$n=7$
+
+\verb$l=3$
+
+\verb$m=1$
+
+\verb$R=r^l*exp(-r/n)*laguerre(2*r/n,n-l-1,2*l+1)$
+
+\verb$Y=legendre(cos(theta),l,abs(m))*exp(i*m*phi)$
+
+\verb$psi=R*Y$
+
+\verb$E=psi/n^2$
+
+\verb$K=laplacian(psi)$
+
+\verb$V=2*psi/r$
+
+\verb$simplify(E-K-V)$
+
+$$0$$
+

doc/manual/integral-trick.tex

@@ -0,0 +1,34 @@
+
+\newpage
+
+\noindent
+Here is a useful trick.
+Difficult integrals involving sine and cosine
+can often be solved by using exponentials.
+Trigonometric simplifications involving powers
+and multiple angles turn into simple algebra in the
+exponential domain.
+For example, the definite integral
+$$\int_0^{2\pi}\left(\sin^4t-2\cos^3(t/2)\sin t\right)dt$$
+can be solved as follows.
+
+\medskip
+\verb$f=sin(t)^4-2*cos(t/2)^3*sin(t)$
+
+\verb$f=circexp(f)$
+
+\verb$defint(f,t,0,2*pi)$
+
+$$-{16\over5}+{3\over4}\pi$$
+
+\medskip
+\noindent
+Here is a check.
+
+\medskip
+\verb$g=integral(f,t)$
+
+\verb$f-d(g,t)$
+
+$$0$$
+

doc/manual/integrals.tex

@@ -0,0 +1,44 @@
+\subsection{Integral}
+\index{integral}
+\noindent
+$integral(f,x)$ returns the integral of $f$ with respect to $x$.
+The $x$ can be omitted for expressions in $x$.
+A multi-integral can be obtained by extending the argument list.
+
+\medskip
+\verb$integral(x^2)$
+$${1\over3}x^3$$
+
+\verb$integral(x*y,x,y)$
+$${1\over4}x^2y^2$$
+
+\medskip
+\noindent
+$defint(f,x,a,b,\ldots)$
+computes the definite integral of $f$ with respect to $x$ evaluated from $a$ to $b$.
+The argument list can be extended for multiple integrals.
+
+\medskip
+\noindent
+The following example computes the integral of $f=x^2$
+over the domain of a semicircle.
+For each $x$ along the abscissa, $y$ ranges from 0 to $\sqrt{1-x^2}$.
+
+\medskip
+\verb$defint(x^2,y,0,sqrt(1-x^2),x,-1,1)$
+
+$${1\over8}\pi$$
+
+\medskip
+\noindent
+As an alternative, the $eval$ function can be used to compute a definite integral step by step.
+
+\medskip
+\verb$I=integral(x^2,y)$
+
+\verb$I=eval(I,y,sqrt(1-x^2))-eval(I,y,0)$
+
+\verb$I=integral(I,x)$
+
+\verb$eval(I,x,1)-eval(I,x,-1)$
+$${1\over8}\pi$$

doc/manual/line-integral.tex

@@ -0,0 +1,257 @@
+\subsection{Arc length}
+
+Let $g(t)$ be a function that draws a curve.
+The arc length from $g(a)$ to $g(b)$ is given by
+$$\int_a^b|g'(t)|\,dt$$
+where $|g'(t)|$ is the length of the tangent vector at $g(t)$.
+The integral sums over all of the tangent lengths to arrive at the total length
+from $a$ to $b$.
+For example, let us measure the length of
+
+\medskip
+\verb$xrange=(0,1)$
+
+\verb$yrange=(0,1)$
+
+\verb$draw(x^2)$
+
+\begin{center}
+\includegraphics[scale=0.4]{arc.png}
+\end{center}
+
+\medskip
+\noindent
+A suitable $g(t)$ for the arc is
+$$g(t)=(t,t^2),\quad0\le t\le1$$
+The Eigenmath solution is
+
+\medskip
+\verb$x=t$
+
+\verb$y=t^2$
+
+\verb$g=(x,y)$
+
+\verb$defint(abs(d(g,t)),t,0,1)$
+
+$$\hbox{$1\over4$}\log(2+5^{1/2})+\hbox{$1\over2$}5^{1/2}$$
+
+\verb$float$
+
+$$1.47894$$
+
+\medskip
+\noindent
+As we would expect, the result is greater than $\sqrt2$, the length of the
+diagonal.
+
+\medskip
+\noindent
+The result seems rather complicated given that we
+started with a simple parabola.
+Let us inspect $|g'(t)|$ to see why.
+
+\medskip
+\verb$g$
+
+$$g=\left(\matrix{t\cr t^2}\right)$$
+
+\medskip
+\verb$d(g,t)$
+
+$$\left(\matrix{1\cr2t}\right)$$
+
+\medskip
+\verb$abs(d(g,t))$
+
+$$(4t^2+1)^{1/2}$$
+
+\medskip
+\noindent
+The following script does a discrete computation of the arc length by dividing
+the curve into 100 pieces.
+
+\medskip
+\verb$g(t)=(t,t^2)$
+
+\verb$h(t)=abs(g(t)-g(t-0.01))$
+
+\verb$L=0$
+
+\verb$for(k,1,100,L=L+h(k/100.0))$
+
+\verb$L$
+
+$$L=1.47894$$
+
+\newpage
+
+\noindent
+Find the length of the curve $y=x^{3/2}$ from the origin to
+$x={4\over3}$.
+
+\medskip
+\verb$x=t$
+
+\verb$y=x^(3/2)$
+
+\verb$g=(x,y)$
+
+\verb$defint(abs(d(g,x)),x,0,4/3)$
+
+$$\hbox{$56\over27$}$$
+
+\medskip
+\noindent
+Because of the way $t$ is substituted for $x$, the previous solution is
+really no different from the following.
+
+\medskip
+\verb$g=(t,t^(3/2))$
+
+\verb$defint(abs(d(g,t)),t,0,4/3)$
+
+$$\hbox{$56\over27$}$$
+
+
+\newpage
+
+\subsection{Line integrals}
+
+There are two different kinds of line integrals,
+one for scalar fields and one
+for vector fields.
+The following table shows how both are based on the calculation of
+arc length.
+
+\bigskip
+
+\begin{center}
+\begin{tabular}{|lll|}
+\hline
+ & & \\
+& Abstract form
+& Computable form
+\\
+ & & \\
+Arc length
+& $\displaystyle{\int_C ds}$
+& $\displaystyle{\int_a^b |g'(t)|\,dt}$
+\\
+ & & \\
+Line integral, scalar field
+& $\displaystyle{\int_C f\,ds}$
+& $\displaystyle{\int_a^b f(g(t))\,|g'(t)|\,dt}$
+\\
+ & & \\
+Line integral, vector field
+& $\displaystyle{\int_C(F\cdot u)\,ds}$
+& $\displaystyle{\int_a^b F(g(t))\cdot g'(t)\,dt}$
+\\
+ & & \\
+\hline
+\end{tabular}
+\end{center}
+
+\bigskip
+\noindent
+For the vector field form, the symbol $u$ is the unit tangent vector
+$$u={g'(t)\over|g'(t)|}$$
+The length of the tangent vector cancels with $ds$
+as follows.
+$$\int_C(F\cdot u)\,ds
+=\int_a^b\bigg(F(g(t))\cdot{g'(t)\over|g'(t)|}\bigg)\,\bigg(|g'(t)|\,dt\bigg)
+=\int_a^b F(g(t))\cdot g'(t)\,dt
+$$
+
+\newpage
+
+\noindent
+Evaluate
+$$\int_Cx\,ds\quad\hbox{and}\quad\int_Cx\,dx$$
+where $C$ is a straight line from $(0,0)$ to $(1,1)$.
+
+\medskip
+\noindent
+What a difference the measure makes.
+The first integral is over a scalar field and the second is over a vector field.
+This can be understood when we recall that
+$$ds=|g'(t)|\,dt
+%\quad\hbox{and}\quad
+%\int_Cx\,dx=\int_Cx\,dx+0\,dy
+$$
+Hence for $\int_Cx\,ds$ we have
+
+\medskip
+\verb$x=t$
+
+\verb$y=t$
+
+\verb$g=(x,y)$
+
+\verb$defint(x*abs(d(g,t)),t,0,1)$
+
+$$1\over2^{1/2}$$
+
+\medskip
+\noindent
+For $\int_Cx\,dx$ we have
+
+\medskip
+\verb$x=t$
+
+\verb$y=t$
+
+\verb$g=(x,y)$
+
+\verb$F=(x,0)$
+
+\verb$defint(dot(F,d(g,t)),t,0,1)$
+
+$$1\over2$$
+
+\newpage
+
+\noindent
+The following line integral problems are from
+{\it Advanced Calculus, Fifth Edition} by Wilfred Kaplan.
+
+\medskip
+\noindent
+Evaluate $\int y^2\,dx$ along the straight
+line from $(0,0)$ to $(2,2)$.
+
+\medskip
+\verb$x=2t$
+
+\verb$y=2t$
+
+\verb$g=(x,y)$
+
+\verb$F=(y^2,0)$
+
+\verb$defint(dot(F,d(g,t)),t,0,1)$
+
+$$8\over3$$
+
+\medskip
+\noindent
+Evaluate $\int z\,dx+x\,dy+y\,dz$
+along the path
+$x=2t+1$, $y=t^2$, $z=1+t^3$, $0\le t\le 1$.
+
+\medskip
+\verb$x=2t+1$
+
+\verb$y=t^2$
+
+\verb$z=1+t^3$
+
+\verb$g=(x,y,z)$
+
+\verb$F=(z,x,y)$
+
+\verb$defint(dot(F,d(g,t)),t,0,1)$
+
+$$163\over30$$
+

doc/manual/linear-algebra.tex

@@ -0,0 +1,101 @@
+\section{Linear algebra}
+\index{linear algebra}
+$dot$ is used to multiply vectors and matrices.
+The following example shows how to use $dot$ and $inv$ to solve for
+$\bf X$ in $\bf AX=B$.
+
+\medskip
+{\tt A=((3.8,7.2),(1.3,-0.9))}
+
+{\tt B=(16.5,-22.1)}
+
+{\tt X=dot(inv(A),B)}
+
+{\tt X}
+
+$$\left(\matrix{-11.2887\cr8.24961}\right)$$
+
+\medskip
+\noindent
+One might wonder why the $dot$ function is necessary.
+Why not simply use $X=inv(A)*B$ like scalar multiplication?
+The reason is that the software normally reorders factors internally to optimize processing.
+For example, $inv(A)*B$ in symbolic form is changed to $B*inv(A)$ internally.
+Since the dot product is not commutative, this reordering would give the wrong result.
+Using a function to do the multiply avoids the problem because
+function arguments are not reordered.
+
+\medskip
+\noindent
+It should be noted that $dot$ can have more than two arguments.
+For example, $dot(A,B,C)$ can be used for the dot product of three tensors.
+
+\bigskip
+\noindent
+The following example demonstrates the relation
+${\bf A}^{-1}=\mathop{\rm adj}{\bf A}/\mathop{\rm det}{\bf A}$.
+
+\medskip
+\verb$A=((a,b),(c,d))$
+
+\medskip
+\verb$inv(A)$
+$$\left(\matrix{
+\displaystyle{d\over ad-bc} & \displaystyle{-{b\over ad-bc}}\cr
+\cr
+\displaystyle{-{c\over ad-bc}} & \displaystyle{a\over ad-bc}\cr
+}\right)$$
+
+\medskip
+\verb$adj(A)$
+$$\left(\matrix{
+d & -b\cr
+-c & a\cr
+}\right)$$
+
+\medskip
+\verb$det(A)$
+$$ad-bc$$
+
+\medskip
+\verb$inv(A)-adj(A)/det(A)$
+$$\left(\matrix{
+0 & 0\cr
+0 & 0\cr
+}\right)$$
+
+\medskip
+\noindent
+Sometimes a calculation will be simpler if it can be reorganized to use $adj$ instead of $inv$.
+The main idea is to try to prevent the determinant from appearing as a divisor.
+For example, suppose for matrices $\bf A$ and $\bf B$ you want to check that
+$${\bf A}-{\bf B}^{-1}=0$$
+Depending on the complexity of $\mathop{\rm det}\bf B$, the software
+may not be able to find a simplification that yields zero.
+Should that occur, the following alternative can be tried.
+$$(\mathop{\rm det}{\bf B})\cdot{\bf A}-\mathop{\rm adj}{\bf B}=0$$
+
+\bigskip
+\noindent
+The adjunct of a matrix is related to the cofactors as follows.
+
+\medskip
+\verb$A=((a,b),(c,d))$
+
+\verb$C=((0,0),(0,0))$
+
+\verb$C[1,1]=cofactor(A,1,1)$
+
+\verb$C[1,2]=cofactor(A,1,2)$
+
+\verb$C[2,1]=cofactor(A,2,1)$
+
+\verb$C[2,2]=cofactor(A,2,2)$
+
+\verb$C$
+
+$$C=\left(\matrix{d&-c\cr -b&a}\right)$$
+
+\verb$adj(A)-transpose(C)$
+
+$$\left(\matrix{0&0\cr0&0\cr}\right)$$

doc/manual/list-of-functions.tex

@@ -0,0 +1,402 @@
+\section{Built-in functions}
+
+\section*{abs}
+abs($x$) returns the absolute value or vector length of $x$.
+The mag function should be used for complex $x$.
+
+\medskip
+{\tt P=(x,y)}
+
+{\tt abs(P)}
+
+$$(x^2+y^2)^{1/2}$$
+
+\section*{adj}
+adj($m$) returns the adjunct of matrix $m$.
+
+\section*{and}
+and($a,b,\ldots$) returns the logical ``and'' of predicate expressions.
+
+\section*{arccos}
+arccos($x$) returns the inverse cosine of $x$.
+
+\section*{arccosh}
+arccosh($x$) returns the inverse hyperbolic cosine of $x$.
+
+\section*{arcsin}
+arcsin($x$) returns the inverse sine of $x$.
+
+\section*{arcsinh}
+arcsinh($x$) returns the inverse hyperbolic sine of $x$.
+
+\section*{arctan}
+arcttan($x$) returns the inverse tangent of $x$.
+
+\section*{arctanh}
+arctanh($x$) returns the inverse hyperbolic tangent of $x$.
+
+\section*{arg}
+arg($z$) returns the angle of complex $z$.
+
+\section*{ceiling}
+ceiling($x$) returns the smallest integer not less than $x$.
+
+\section*{check}
+check($x$) In a script, if the predicate $x$ is true then continue, else stop.
+
+\section*{choose}
+choose($n,k$) returns $\displaystyle\left({n \atop k}\right)$
+
+\section*{circexp}
+circexp($x$) returns expression $x$ with circular functions converted
+to exponential forms.
+Sometimes this will simplify an expression.
+
+\section*{coeff}
+coeff($p,x,n$) returns the coefficient of $x^n$ in polynomial $p$.
+
+\section*{cofactor}
+cofactor($m,i,j$) returns of the cofactor of matrix $m$ with respect to row $i$ and column $j$.
+
+\section*{conj}
+conj($z$) returns the complex conjugate of $z$.
+
+\section*{contract}
+\index{trace}
+contract($a,i,j$) returns tensor $a$ summed over indices $i$ and $j$.
+If $i$ and $j$ are omitted then indices 1 and 2 are used.
+contract($m$) is equivalent to the trace of matrix $m$.
+
+\section*{cos}
+cos($x$) returns the cosine of $x$.
+%If $x$ is a floating point number then $\cos(x)$ is evaluated numerically.
+
+\section*{cosh}
+cosh($x$) returns the hyperbolic cosine of $x$.
+
+\section*{cross}
+cross($u,v$) returns the cross product of vectors $u$ and $v$.
+
+\section*{curl}
+curl($u$) returns the curl of vector $u$.
+
+\section*{d}
+d($f,x$) returns the derivative of $f$ with respect to $x$.
+
+\section*{defint}
+defint($f,x,a,b,\ldots$)
+returns the definite integral of $f$ with respect to $x$ evaluated from $a$ to $b$.
+The argument list can be extended for multiple integrals.
+For example, $d(f,x,a,b,y,c,d)$.
+
+\section*{deg}
+deg($p,x$) returns the degree of polynomial $p$ in $x$.
+
+\section*{denominator}
+denominator($x$) returns the denominator of expression $x$.
+
+\section*{det}
+det($m$) returns the determinant of matrix $m$.
+
+\section*{do}
+do($a,b,\ldots$) evaluates the argument list from left to right.
+Returns the result of the last argument.
+
+\section*{dot}
+dot($a,b,\ldots$) returns the dot product of tensors.
+
+\section*{draw}
+draw($f,x$) draws the function $f$ with respect to $x$.
+
+\section*{erf}
+erf($x$) returns the error function of $x$.
+
+\section*{erfc}
+erf($x$) returns the complementary error function of $x$.
+
+\section*{eval}
+eval($f,x,n$) returns $f$ evaluated at $x=n$.
+
+\section*{exp}
+exp($x$) returns $e^x$.
+
+\section*{expand}
+expand($r,x$) returns the partial fraction expansion of the ratio of
+polynomials $r$ in $x$.
+
+\medskip
+\verb$expand(1/(x^3+x^2),x)$
+
+$${1\over x^2}-{1\over x}+{1\over x+1}$$
+
+\section*{expcos}
+expcos($x$) returns the cosine of $x$ in exponential form.
+
+\medskip
+{\tt expcos(x)}
+
+$${1\over2}\exp(-ix)+{1\over2}\exp(ix)$$
+
+\section*{expsin}
+expsin($x$) returns the sine of $x$ in exponential form.
+
+\medskip
+{\tt expsin(x)}
+
+$${1\over2}i\exp(-ix)-{1\over2}i\exp(ix)$$
+
+\section*{factor}
+factor($n$) factors the integer $n$.
+
+\medskip
+{\tt factor(12345)}
+
+$$3\times 5\times 823$$
+
+\medskip
+\noindent
+factor($p,x$) factors polynomial $p$ in $x$.
+The last argument can be omitted for polynomials in $x$.
+The argument list can be extended for multivariate polynomials.
+For example, factor($p,x,y$) factors $p$ over $x$ and then over $y$.
+
+\medskip
+{\tt factor(125*x{\char94}3-1)}
+
+$$(5x-1)(25x^2+5x+1)$$
+
+\section*{factorial}
+Example:
+
+\medskip
+{\tt 10!}
+
+$$3628800$$
+
+\section*{filter}
+filter($f,a,b,\ldots$) returns $f$ with terms involving $a$, $b$, etc. removed.
+
+\medskip
+{\tt 1/a+1/b+1/c}
+
+$${1\over a}+{1\over b}+{1\over c}$$
+
+{\tt filter(last,a)}
+
+$${1\over b}+{1\over c}$$
+
+\section*{float}
+float($x$) converts $x$ to a floating point value.
+
+\medskip
+{\tt sum(n,0,20,(-1/2){\char94}n)}
+
+$$699051\over1048576$$
+
+{\tt float(last)}
+
+$$0.666667$$
+
+\section*{floor}
+floor($x$) returns the largest integer not greater than $x$.
+
+\section*{for}
+for($i,j,k,a,b,\ldots$) For $i$ equals $j$ through $k$ evaluate $a$, $b$, etc.
+
+\medskip
+{\tt x=0}
+
+{\tt y=2}
+
+{\tt for(k,1,9,x=sqrt(2+x),y=2*y/x)}
+
+{\tt float(y)}
+
+$$3.14159$$
+
+\section*{gcd}
+gcd($a,b,\ldots$) returns the greatest common divisor.
+
+\section*{hermite}
+hermite($x,n$) returns the $n$th Hermite polynomial in $x$.
+
+\section*{hilbert}
+hilbert($n$) returns a Hilbert matrix of order $n$.
+
+\section*{imag}
+imag($z$) returns the imaginary part of complex $z$.
+
+\section*{inner}
+inner($a,b,\ldots$) returns the inner product of tensors.
+Same as the dot product.
+
+\section*{integral}
+integral($f,x$) returns the integral of $f$ with respect to $x$.
+
+\section*{inv}
+inv($m$) returns the inverse of matrix $m$.
+
+\section*{isprime}
+isprime($n$) returns 1 if $n$ is prime, zero otherwise.
+
+\medskip
+{\tt isprime(2{\char94}53-111)}
+
+$$1$$
+
+\section*{laguerre}
+laguerre($x,n,a$) returns the $n$th Laguerre polynomial in $x$.
+If $a$ is omitted then $a=0$ is used.
+
+\section*{lcm}
+lcm($a,b,\ldots$) returns the least common multiple.
+
+\section*{leading}
+leading($p,x$) returns the leading coefficient of polynomial $p$ in $x$.
+
+\medskip
+\verb$leading(5x^2+x+1,x)$
+
+$$5$$
+
+\section*{legendre}
+legendre($x,n,m$) returns the $n$th Legendre polynomial in $x$.
+If $m$ is omitted then $m=0$ is used.
+
+\section*{log}
+log($x$) returns the natural logarithm of $x$.
+
+\section*{mag}
+mag($z$) returns the magnitude of complex $z$.
+
+\section*{mod}
+mod($a,b$) returns the remainder of $a$ divided by $b$.
+
+\section*{not}
+not($x$) negates the result of predicate expression $x$.
+
+\section*{nroots}
+nroots($p,x$) returns all of the roots, both real and complex, of
+polynomial $p$ in $x$.
+The roots are computed numerically.
+The coefficients of $p$ can be real or complex.
+
+\section*{numerator}
+numerator($x$) returns the numerator of expression $x$.
+%\begin{itemize}
+%\item[$\scriptstyle1$]{\tt numerator(a/b+b/a)}
+%\item[$\scriptstyle2$]\hspace{50pt} $a^2+b^2$
+%\end{itemize}
+
+\section*{or}
+or($a,b,\ldots$) returns the logical ``or'' of predicate expressions.
+
+\section*{outer}
+outer($a,b,\ldots$) returns the outer product of tensors.
+
+\section*{polar}
+polar($z$) converts complex $z$ to polar form.
+
+\section*{prime}
+prime($n$) returns the $n$th prime number, $1\le n\le10{,}000$.
+
+\section*{print}
+print($a,b,\ldots$) evaluates expressions and prints the results..
+Useful for printing from inside a ``for'' loop.
+
+\section*{product}
+product($i,j,k,f$) returns $\displaystyle\prod_{i=j}^k f$
+
+\section*{quote}
+quote($x$) returns expression $x$ unevaluated.
+
+\section*{quotient}
+quotient($p,q,x$) returns the quotient of polynomials in $x$.
+
+\section*{rank}
+rank($a$) returns the number of indices that tensor $a$ has.
+A scalar has no indices so its rank is zero.
+
+\section*{rationalize}
+rationalize($x$) puts everything over a common denominator.
+
+\medskip
+{\tt rationalize(a/b+b/a)}
+
+$$a^2+b^2\over ab$$
+
+\section*{real}
+real($z$) returns the real part of complex $z$.
+
+\section*{rect}
+rect($z$) returns complex $z$ in rectangular form.
+
+\section*{roots}
+roots($p,x$) returns the values of $x$ such that the polynomial $p(x)=0$.
+The polynomial should be factorable over integers.
+
+\section*{simplify}
+simplify($x$) returns $x$ in a simpler form.
+
+\section*{sin}
+sin($x$) returns the sine of $x$.
+
+\section*{sinh}
+sinh($x$) returns the hyperbolic sine of $x$.
+
+\section*{sqrt}
+sqrt($x$) returns the square root of $x$.
+
+\section*{stop}
+In a script, it does what it says.
+
+\section*{subst}
+subst($a,b,c$) substitutes $a$ for $b$ in $c$ and returns the result.
+
+\section*{sum}
+sum($i,j,k,f$) returns $\displaystyle\sum_{i=j}^k f$
+
+\section*{tan}
+tan($x$) returns the tangent of $x$.
+
+\section*{tanh}
+tanh($x$) returns the hyperbolic tangent of $x$.
+
+\section*{taylor}
+taylor($f,x,n,a$) returns the Taylor expansion of $f$ of $x$ at $a$.
+The argument $n$ is the degree of the expansion.
+If $a$ is omitted then $a=0$ is used.
+
+\medskip
+{\tt taylor(1/cos(x),x,4)}
+
+$${5\over24}x^4+{1\over2}x^2+1$$
+
+\section*{test}
+test($a,b,c,d,\ldots$)
+If $a$ is true then $b$ is returned else if $c$ is true then $d$ is returned, etc.
+If the number of arguments is odd then the last argument is returned when all else fails.
+
+\section*{transpose}
+transpose($a,i,j$) returns the transpose of tensor $a$ with respect to indices $i$ and $j$.
+If $i$ and $j$ are omitted then 1 and 2 are used.
+Hence a matrix can be transposed with a single argument.
+
+\medskip
+{\tt A=((a,b),(c,d))}
+
+{\tt transpose(A)}
+
+$$\left(\matrix{a & c\cr b & d\cr}\right)$$
+
+\section*{unit}
+unit($n$) returns an $n\times n$ identity matrix.
+
+\medskip
+{\tt unit(2)}
+
+$$\left(\matrix{1&0\cr0&1\cr}\right)$$
+
+\section*{zero}
+zero($i,j,\ldots$) returns a null tensor with dimensions $i$, $j$, etc.
+Useful for creating a tensor and then setting the component values.

doc/manual/nabokov.tex

@@ -0,0 +1,111 @@
+\section{Introduction}
+The following is an excerpt from Vladimir Nabokov's
+autobiography {\it Speak, Memory.}
+\begin{quote}
+A foolish tutor had explained logarithms to me much too early, and I had
+read (in a British publication, the {\it Boy's Own Paper}, I believe)
+about a certain Hindu calculator who in exactly two seconds could find the
+seventeenth root of, say,
+3529471145 760275132301897342055866171392
+(I am not sure I have got this right; anyway the root was 212).
+\end{quote}
+We can check Nabokov's arithmetic by typing the following into Eigenmath.
+
+\medskip
+\verb$212^17$
+
+\medskip
+\noindent
+After pressing the return key, Eigenmath displays the following result.
+$$3529471145760275132301897342055866171392$$
+So Nabokov did get it right after all.
+We can enter {\it float} or click on the float button to scale the number
+down to size.
+
+\medskip
+\verb$float$
+$$3.52947\times10^{39}$$
+
+\medskip
+\noindent
+Now let us see if Eigenmath can find the
+seventeenth root of this number, like the Hindu calculator could.
+
+\medskip
+\verb$N=212^17$
+
+\verb$N$
+$$N=3529471145760275132301897342055866171392$$
+
+\verb$N^(1/17)$
+$$212$$
+
+\medskip
+\noindent
+It is worth mentioning that when a symbol is assigned a value,
+no result is printed.
+To see the value of a symbol, just evaluate it by putting it on a line by
+itself.
+
+\medskip
+\verb$N$
+$$N=3529471145760275132301897342055866171392$$
+
+\newpage
+
+\subsection{Negative exponents}
+Eigenmath requires parentheses around negative exponents.
+For example,
+
+\medskip
+\verb$10^(-3)$
+
+\medskip
+\noindent
+instead of
+
+\medskip
+\verb$10^-3$
+
+\medskip
+\noindent
+The reason for this is that the binding of the negative sign is not always
+obvious.
+For example, consider
+
+\medskip
+\verb$x^-1/2$
+
+\medskip
+\noindent
+It is not clear whether the exponent should be $-1$ or $-1/2$.
+So Eigenmath requires
+
+\medskip
+\verb$x^(-1/2)$
+
+\medskip
+\noindent
+which is unambiguous.
+
+\medskip
+\noindent
+Now a new question arises.
+Never mind the minus sign, what is the binding of the caret symbol itself?
+The answer is, it binds to the first symbol that follows it and nothing else.
+For example, the following is parsed as $(x^1)/2$.
+
+\medskip
+\verb$x^1/2$
+
+$$\hbox{$1\over2$}x$$
+
+\medskip
+\noindent
+So in general, parentheses are needed when the exponent is an expression.
+
+\medskip
+\verb$x^(1/2)$
+
+$$x^{1/2}$$
+

doc/manual/qho.tex

@@ -0,0 +1,26 @@
+\subsection{Quantum harmonic oscillator}
+For total energy $E$, kinetic energy $K$ and potential energy $V$ we have
+$$E=K+V$$
+The corresponding formula for a quantum harmonic oscillator is
+$$(2n+1)\psi=-{d^2\psi\over dx^2}+x^2\psi$$
+where $n$ is an integer and represents the quantization of energy values.
+The solution to the above equation is
+$$\psi_n(x)=\exp(-x^2/2)H_n(x)$$
+where $H_n(x)$ is the $n$th Hermite polynomial in $x$.
+The following Eigenmath code checks $E=K+V$ for $n=7$.
+
+\medskip
+\verb$n=7$
+
+\verb$psi=exp(-x^2/2)*hermite(x,n)$
+
+\verb$E=(2*n+1)*psi$
+
+\verb$K=-d(psi,x,x)$
+
+\verb$V=x^2*psi$
+
+\verb$E-K-V$
+
+$$0$$
+

doc/manual/scripting.tex

@@ -0,0 +1,73 @@
+\section{Scripts}
+\index{scripts}
+Here is a simple example that draws the graph of $y=mx+b$.
+
+\medskip
+\verb$y=m*x+b$
+
+\verb$m=1/2$
+
+\verb$b=-3$
+
+\verb$draw(y)$
+
+\begin{center}
+\includegraphics[scale=0.4]{1.png}
+\end{center}
+
+\noindent
+Now suppose that we want to draw the graph
+with a different $m$.
+We could type in everything all over again, but it would be easier
+in the long run to write a script.
+Then we can go back and quickly change $m$ and $b$ as many times as we want.
+
+\medskip
+\noindent
+To prepare a script, click on the Edit Script button.
+Then enter the script commands, one per line, as shown above.
+Then click on the Run Script button to see the graph.
+
+\medskip
+\noindent
+Eigenmath runs a script by stepping through it line by line.
+Each line is evaluated just like a regular command.
+This continues until the end of the script is reached.
+After the script runs, you can click Edit Script and go back and change something.
+%By the way, Eigenmath automatically does a clear
+%running a script.
+
+\newpage
+
+\noindent
+Sometimes it is desirable to have a script print a few comments when it runs.
+This can be accomplished by placing the desired text in quotes
+on a single line.
+For example, the script
+
+\medskip
+\verb$"Here is the value of pi."$
+
+\verb$float(pi)$
+
+\medskip
+\noindent
+displays the following when run.
+
+\medskip
+\verb$Here is the value of pi.$
+
+$$3.14159$$
+
+\medskip
+\noindent
+Eigenmath includes a simple debug facility.
+Setting the variable $trace$ to 1 causes each line of the script to be
+printed as the script runs.
+Normally this setting would be the first line in the script.
+
+\medskip
+\verb$trace=1$
+
+\verb$--Now each line of the script is printed as it runs.$
+

doc/manual/space-shuttle-and-corvette.tex

@@ -0,0 +1,46 @@
+\subsection{Space shuttle and Corvette}
+The space shuttle accelerates from zero to 17{,}000 miles per hour
+in 8 minutes.
+A Corvette accelerates from zero to 60 miles per hour in 4.5 seconds.
+The following script compares the two.
+
+\medskip
+\verb$vs=17000*"mile"/"hr"$
+
+\verb$ts=8*"min"/(60*"min"/"hr")$
+
+\verb$as=vs/ts$
+
+\verb$as$
+
+\verb$vc=60*"mile"/"hr"$
+
+\verb$tc=4.5*"sec"/(3600*"sec"/"hr")$
+
+\verb$ac=vc/tc$
+
+\verb$ac$
+
+\verb$"Time for Corvette to reach orbital velocity:"$
+
+\verb$vs/ac$
+
+\verb$vs/ac*60*"min"/"hr"$
+
+\medskip
+\noindent
+Here is the result when the script runs.
+It turns out that the space shuttle accelerates more than twice as fast as a
+Corvette.
+
+\medskip
+$$a_s={\hbox{127500 mile}\over(\hbox{hr})^2}$$
+
+$$a_c={\hbox{48000 mile}\over(\hbox{hr})^2}$$
+
+\verb$Time for Corvette to reach orbital velocity:$
+
+$$\hbox{0.354167 hr}$$
+
+$$\hbox{21.25 min}$$
+

doc/manual/stokes-theorem-standalone.tex

@@ -0,0 +1,5 @@
+\documentclass{article}
+\pagestyle{empty}
+\begin{document}
+\include{stokes-theorem}
+\end{document}

doc/manual/stokes-theorem.tex

@@ -0,0 +1,81 @@
+\subsection{Stokes' theorem}
+\index{Stokes' theorem}
+Stokes' theorem proves the following equivalence of line and surface
+integrals.
+%\bigskip
+%\noindent
+%$\displaystyle{\oint_C P\,dx+Q\,dy+R\,dz}$
+%$$=
+%\int\!\!\!\int_S
+%\left({\partial Q\over\partial x}-{\partial P\over\partial y}\right)\,dx\,dy
+%+
+%\left({\partial R\over\partial y}-{\partial Q\over\partial z}\right)\,dy\,dz
+%+
+%\left({\partial P\over\partial z}-{\partial R\over\partial x}\right)\,dz\,dx
+%$$
+%
+%\noindent
+%Curve $C$ is the perimeter around $S$.
+%The theorem can be also be written as
+$$\oint P\,dx+Q\,dy+R\,dz
+=\int\!\!\!\int_S(\mathop{\rm curl}{\bf F})\cdot{\bf n}\,d\sigma
+$$
+where ${\bf F}=(P,Q,R)$.
+For $S$ parametrized by $x$ and $y$ we have
+$${\bf n}\,d\sigma=\left(
+{\partial S\over\partial x}\times{\partial S\over\partial y}
+\right)dx\,dy$$
+In many cases, converting an integral according to
+Stokes' theorem can turn a difficult problem into an easy one.
+
+\medskip
+\noindent
+Let ${\bf F}=(y,z,x)$ and let $S$ be the part of the paraboloid
+$z=4-x^2-y^2$
+that is above the $xy$ plane.
+The perimeter of the paraboloid is the circle $x^2+y^2=2$.
+Calculate both the line and surface integrals.
+It turns out that we need to use polar coordinates so that {\it defint} can
+succeed.
+
+\medskip
+\verb$--Surface integral$
+
+\verb$z=4-x^2-y^2$
+
+\verb$F=(y,z,x)$
+
+\verb$S=(x,y,z)$
+
+\verb$f=dot(curl(F),cross(d(S,x),d(S,y)))$
+
+\verb$x=r*cos(theta)$
+
+\verb$y=r*sin(theta)$
+
+\verb$defint(f*r,r,0,2,theta,0,2pi)$
+
+$$-4\pi$$
+
+\verb$--Line integral$
+
+\verb$x=2*cos(t)$
+
+\verb$y=2*sin(t)$
+
+\verb$z=4-x^2-y^2$
+
+\verb$P=y$
+
+\verb$Q=z$
+
+\verb$R=x$
+
+\verb$f=P*d(x,t)+Q*d(y,t)+R*d(z,t)$
+
+\verb$f=circexp(f)$
+
+\verb$defint(f,t,0,2pi)$
+
+$$-4\pi$$
+

doc/manual/surface-area.tex

@@ -0,0 +1,73 @@
+\subsection{Surface area}
+Let $S$ be a surface parameterized by $x$ and $y$.
+That is, let $S=(x,y,z)$ where $z=f(x,y)$.
+The tangent lines at a point on $S$ form a tiny parallelogram.
+The area $a$ of the parallelogram is given by the magnitude of the cross product.
+$$a=\left|{\partial S\over\partial x}\times{\partial S\over\partial y}\right|$$
+By summing over all the parallelograms we obtain the total surface area $A$.
+Hence
+$$A=\int\!\!\!\int dA=\int\!\!\!\int a\,dx\,dy$$
+The following example computes the surface area of a unit disk
+parallel to the $xy$ plane.
+
+\medskip
+\verb$z=2$
+
+\verb$S=(x,y,z)$
+
+\verb$a=abs(cross(d(S,x),d(S,y)))$
+
+\verb$defint(a,y,-sqrt(1-x^2),sqrt(1-x^2),x,-1,1)$
+
+$$\pi$$
+
+\medskip
+\noindent
+The result is $\pi$, the area of a unit circle, which is what we expect.
+The following example computes the surface area of $z=x^2+2y$ over
+a unit square.
+
+\medskip
+\verb$z=x^2+2y$
+
+\verb$S=(x,y,z)$
+
+\verb$a=abs(cross(d(S,x),d(S,y)))$
+
+\verb$defint(a,x,0,1,y,0,1)$
+
+$${3\over2}+{5\over8}\log(5)$$
+
+\medskip
+\noindent
+As a practical matter, $f(x,y)$ must be very simple in order
+for Eigenmath to solve the double integral.
+
+\newpage
+
+\noindent
+Find the area of the spiral ramp defined by\footnote{
+Williamson and Trotter, {\it Multivariable Mathematics,} p. 598.}
+$$S=\left(\matrix{u\cos v\cr u\sin v\cr v}\right),\qquad 0\le u\le1,\qquad 0\le v\le3\pi$$
+In this example, the coordinates $x$, $y$ and $z$ are all
+functions of an independent parameter space.
+
+\medskip
+\verb$x=u*cos(v)$
+
+\verb$y=u*sin(v)$
+
+\verb$z=v$
+
+\verb$S=(x,y,z)$
+
+\verb$a=abs(cross(d(S,u),d(S,v)))$
+
+\verb$defint(a,u,0,1,v,0,3pi)$
+
+$${3\over2}\pi\log(1+2^{1/2})+{3\pi\over2^{1/2}}$$
+
+\verb$float$
+
+$$10.8177$$
+

doc/manual/surface-integral.tex

@@ -0,0 +1,53 @@
+\subsection{Surface integrals}
+\index{surface integral}
+%\begin{center}
+%\includegraphics[scale=0.5]{sailboat.png}
+%\end{center}
+%\bigskip
+%\noindent
+A surface integral is like adding up all the wind on a sail.
+In other words, we want to compute
+$$\int\!\!\!\int{\bf F\cdot n}\,dA$$
+where ${\bf F\cdot n}$ is the amount of wind normal to a tiny parallelogram $dA$.
+The integral sums over the entire area of the sail.
+Let $S$ be the surface of the sail parameterized by $x$ and $y$.
+(In this model, the $z$ direction points downwind.)
+By the properties of the cross product we have the following for the unit normal $\bf n$
+and for $dA$.
+$${\bf n}={ {{\partial S\over\partial x}\times{\partial S\over\partial y}}\over
+ {\left|{\partial S\over\partial x}\times{\partial S\over\partial y}\right|}}\qquad
+dA=\left|{\partial S\over\partial x}\times{\partial S\over\partial y}\right|\,dx\,dy$$
+Hence
+$$\int\!\!\!\int{\bf F\cdot n}\,dA=\int\!\!\!\int{\bf F}\cdot
+\left({{\partial S\over\partial x}\times{\partial S\over\partial y}}\right)\,dx\,dy$$
+
+\noindent
+For example, evaluate the surface integral
+$$\int\!\!\!\int_S{\bf F\cdot n}\,d\sigma$$
+where ${\bf F}=xy^2z{\bf i}-2x^3{\bf j}+yz^2{\bf k}$, $S$ is the surface
+$z=1-x^2-y^2$, $x^2+y^2\le1$ and $\bf n$ is upper.\footnote{
+Kaplan, {\it Advanced Calculus,} p. 313.}
+
+\medskip
+\noindent
+Note that the surface intersects the $xy$ plane in a circle.
+By the right hand rule, crossing $x$ into $y$ yields $\bf n$ pointing upwards hence
+$${\bf n}\,d\sigma=\left({{\partial S\over\partial x}\times{\partial S\over\partial y}}\right)\,dx\,dy$$
+The following Eigenmath code computes the surface integral.
+The symbols $f$ and $h$ are used as temporary variables.
+
+\medskip
+\verb$z=1-x^2-y^2$
+
+\verb$F=(x*y^2*z,-2*x^3,y*z^2)$
+
+\verb$S=(x,y,z)$
+
+\verb$f=dot(F,cross(d(S,x),d(S,y)))$
+
+\verb$h=sqrt(1-x^2)$
+
+\verb$defint(f,y,-h,h,x,-1,1)$
+
+$${1\over48}\pi$$
+

doc/manual/symbols.tex

@@ -0,0 +1,151 @@
+
+\subsection{Defining symbols}
+As we saw earlier, Eigenmath uses the same syntax as dear old Fortran.
+
+\medskip
+\verb$N=212^17$
+
+\medskip
+\noindent
+No result is printed when a symbol is defined.
+To see a symbol's value, just evaluate it.
+
+\medskip
+\verb$N$
+
+$$3529471145760275132301897342055866171392$$
+
+\medskip
+\noindent
+Beyond its prosaic syntax, Eigenmath does have a few tricks up its sleeve.
+For example, a symbol can have a subscript.
+
+\medskip
+\verb$NA=6.02214*10^23$
+
+\verb$NA$
+
+$$N_A=6.02214\times10^{23}$$
+
+\medskip
+\noindent
+A symbol can be the name of a Greek letter.
+
+\medskip
+\verb$xi=1/2$
+
+\verb$xi$
+
+$$\xi=\hbox{$1\over2$}$$
+
+\medskip
+\noindent
+Since xi is $\xi$, how is $x_i$ entered?
+Well, that is an issue that may get resolved in the future.
+For now, xi is always $\xi$.
+
+\medskip
+\noindent
+Greek letters can appear in the subscript too.
+
+\medskip
+\verb$Amu=2.0$
+
+\verb$Amu$
+
+$$A_\mu=2.0$$
+
+\medskip
+\noindent
+The general rule is this.
+Eigenmath scans the entire symbol looking for Greek letters.
+
+\medskip
+\verb$alphamunu$
+
+$$\alpha_{\mu\nu}$$
+
+\newpage
+
+\medskip
+\noindent
+Let us turn now to what happens when a symbolic expression is evaluated.
+The most important point is that
+Eigenmath exhaustively evaluates symbolic subexpressions.
+
+\medskip
+\verb$A=B$
+
+\verb$B=C$
+
+\verb$C=D$
+
+\verb$sin(A)$
+
+$$\sin(D)$$
+
+\medskip
+\noindent
+In the above example, evaluating $\sin(A)$ yields $\sin(D)$ because Eigenmath
+resolves $A$ as far as it can, in this case down to $D$.
+However, internally the binding of $A$ is still $B$, as can be seen with
+the $binding$ function.
+
+\medskip
+\verb$binding(A)$
+
+$$B$$
+
+\medskip
+\noindent
+Let us return to symbolic definitions for a moment.
+It should be kept in mind that the right hand side of the definition is an
+expression that is evaluated before the binding is done.
+For example,
+
+\medskip
+\verb$B=1$
+
+\verb$A=B$
+
+\verb$binding(A)$
+
+$$1$$
+
+\medskip
+\noindent
+The binding of $A$ is 1 and not $B$ because $B$ was already defined before
+the $A=B$ occurred.
+The $quote$ function can be used to give a literal binding.
+
+\medskip
+\verb$A=quote(B)$
+
+\verb$binding(A)$
+
+$$B$$
+
+\newpage
+
+\noindent
+What this all means is that symbols have a dual nature.
+A symbol has a binding which may be different from its evaluation.
+Normally this difference is not important.
+The functions $quote$ and $binding$ are mentioned here mainly to provide insight
+into what is happening belowdecks.
+Normally you should not really need to use these functions.
+However, one notable exception is the use of $quote$ to clear a symbol.
+
+\medskip
+\verb$x=3$
+
+\verb$x$
+
+$$x=3$$
+
+\verb$x=quote(x)$
+
+\verb$x$
+
+$$x$$
+

doc/manual/syntax.tex

@@ -0,0 +1,41 @@
+\section{Syntax}
+
+%The symbol {\tt\char32} indicates a mandatory space.
+
+\begin{center}
+\begin{tabular}{clll}
+{\it Math} & & {\it Eigenmath} & {\it Alternate form and/or comment} \\
+\\
+$-a$ & & {\tt -a} \\
+\\
+$a+b$ & & {\tt a+b} \\
+\\
+$a-b$ & & {\tt a-b} \\
+\\
+$ab$ & & {\tt a*b} & \verb$a b$ \hspace{10pt}
+{\it with a space in between} \\
+\\
+$\displaystyle{a\over b}$ & & {\tt a/b} \\
+\\
+$\displaystyle{a\over bc}$ & & {\tt a/b/c} \\
+\\
+$a^2$ & & {\tt a{\char94}2} \\
+\\
+$\sqrt{a}$ & & {\tt a{\char94}(1/2)} & {\tt sqrt(a)} \\
+\\
+$\displaystyle{1\over\sqrt a}$ & & {\tt a{\char94}(-1/2)} & {\tt 1/sqrt(a)} \\
+\\
+$a(b+c)$ & & {\tt a*(b+c)} & \verb$a (b+c)$
+\hspace{10pt} {\it with a space in between} \\
+\\
+$f(a)$ & & {\tt f(a)} \\
+\\
+$\displaystyle{\left(\matrix{a\cr b\cr c\cr}\right)}$ & & {\tt (a,b,c)} \\
+\\
+$\displaystyle{\left(\matrix{a&b\cr c&d\cr}\right)}$ & & {\tt ((a,b),(c,d))} \\
+\\
+$T^{12}$ & & {\tt T[1,2]} & {\it tensor component access} \\
+\\
+$2\,\rm km$ & & {\tt 2*"km"} & {\it units of measure are quoted} \\
+\end{tabular}
+\end{center}

doc/manual/units-of-measure.tex

@@ -0,0 +1,15 @@
+\subsection{Units of measure}
+\index{units of measure}
+Quoted strings can be used to express units of measurement in a calculation.
+For example, the space shuttle accelerates from zero to
+17{,}000 miles per hour in 8 minutes.
+The average acceleration of the space shuttle is
+
+\medskip
+\verb$v=17000*"mile"/"hr"$
+
+\verb$t=8*"min"/(60*"min"/"hr")$
+
+\verb$v/t$
+$$127500\,\hbox{mile}\over(\hbox{hr})^2$$
+

doc/manual/zerozero.tex

@@ -0,0 +1,37 @@
+\subsection{Zero to the zero power}
+
+The following example draws a graph of the function $f(x)=|x^x|$.
+The graph shows why the convention $0^0=1$ makes sense.
+
+\medskip
+\verb$f(x)=abs(x^x)$
+
+\verb$xrange=(-2,2)$
+
+\verb$yrange=(-2,2)$
+
+\verb$draw(f)$
+
+\begin{center}
+\includegraphics[scale=0.4]{zerozero.png}
+\end{center}
+
+\medskip
+\noindent
+We can see how $0^0=1$ results in a continuous line through $x=0$.
+Now let us see how $x^x$ behaves in the complex plane.
+
+\medskip
+\verb$f(t)=(real(t^t),imag(t^t))$
+
+\verb$xrange=(-2,2)$
+
+\verb$yrange=(-2,2)$
+
+\verb$trange=(-4,2)$
+
+\verb$draw(f)$
+
+\begin{center}
+\includegraphics[scale=0.4]{zerozero2.png}
+\end{center}

doc/ross/Makefile

@@ -0,0 +1,2 @@
+ross.pdf : *.tex
+	pdftex ross

doc/ross/ross-1.1.tex

@@ -0,0 +1,19 @@
+\beginsection 1.1
+
+Prove $1^2+2^2+\cdots+n^2=n(n+1)(2n+1)/6$ for all natural numbers $n$.
+\medskip
+Use mathematical induction.
+First show that proposition $P_1$ is true:
+$$1^2=1(1+1)(2\cdot1+1)/6$$
+Next show that $P_{n+1}$ is true whenever $P_n$ is true:
+$$\eqalign{
+P_n+(n+1)^2&=P_{n+1}\cr
+n(n+1)(2n+1)+(n+1)^2&=(n+1)(n+2)(2n+3)/6\cr
+}$$
+Divide both sides by $(n+1)$ to get
+$$n(2n+1)/6+(n+1)=(n+2)(2n+3)/6$$
+Expand and simplify.
+$$\eqalign{
+(2n^2+n)/6+n+1&=(2n^2+3n+4n+6)/6\cr
+{1\over3}n^2+{7\over6}n+1&={1\over3}n^2+{7\over6}n+1\cr
+}$$

doc/ross/ross-1.2.tex

@@ -0,0 +1,18 @@
+\beginsection 1.2
+
+Prove $3+11+\cdots+(8n-5)=4n^2-n$ for all natural numbers $n$.
+
+\medskip
+
+Induction Step 1: Show that $P_1$ is true.
+$$P_1=(8\cdot1-5)=3$$
+Induction Step 2: Show that $P_n+(8(n+1)-5)=P_{n+1}$.
+$$\eqalign{
+P_n+(8(n+1)-5)&=4n^2-n+8n+3\cr
+&=4n^2+7n+3\cr
+\cr
+P_{n+1}&=4(n+1)^2-(n+1)\cr
+&=4(n^2+2n+1)-n-1\cr
+&=4n^2+8n+4-n-1\cr
+&=4n^2+7n+3\cr
+}$$

doc/ross/ross-1.3.tex

@@ -0,0 +1,20 @@
+\beginsection 1.3
+
+Prove $1^3+2^3+\cdots+n^3=(1+2+\cdots+n)^2$ for all natural numbers $n$.
+
+\medskip
+
+Induction Step 1: Show that $P_1$ is true.
+$$P_1=1^3=1^2$$
+Induction Step 2: Show that $P_n+(n+1)^3=P_{n+1}$.
+Note that $1+2+\cdots+n=n(n+1)/2$.
+$$\eqalign{
+P_n+(n+1)^3&=[n(n+1)/2]^2+(n+1)^3\cr
+&=(n+1)^2[(n/2)^2+(n+1)]\cr
+&=(n+1)^2(n^2/4+n+1)\cr
+\cr
+P_{n+1}&=[n(n+1)/2+(n+1)]^2\cr
+&=[(n+1)(n/2+1)]^2\cr
+&=(n+1)^2(n/2+1)^2\cr
+&=(n+1)^2(n^2/4+n+1)\cr
+}$$

doc/ross/ross-1.4.tex

@@ -0,0 +1,27 @@
+\beginsection{1.4}
+
+\medskip
+(a) Guess a formula for $1+3+\cdots+(2n-1)$ by evaluating the sum for
+$n=1$, 2, 3, and 4. [For $n=1$, the sum is simply 1.]
+
+\medskip
+For each $n$ we have
+$$\eqalign{
+1&=1\cr
+1+3&=4\cr
+1+3+5&=9\cr
+1+3+5+7&=16\cr
+}$$
+therefore a reasonable guess would be
+$$1+3+\cdots+(2n-1)=n^2$$
+
+\medskip
+(b) Prove your formula using mathematical induction.
+
+\medskip
+We already have $n^2=1$ for $n=1$.
+For $n+1$ we have
+$$1+3+\cdots+(2n-1)+(2(n+1)-1)=n^2+2n+1=(n+1)^2$$
+Therefore the formula is true for $n+1$ whenever it is true for $n$.
+Hence by induction the formula is true for all $n\ge1$.
+

doc/ross/ross-1.5.tex

@@ -0,0 +1,14 @@
+\beginsection{1.5}
+
+Prove $1+1/2+1/4+\cdots+1/2^n=2-1/2^n$ for all natural numbers $n$.
+
+\medskip
+For $n=1$ we have $1+1/2=3/2=2-1/2$ so the formula is true for $n=1$.
+Now we want to consider the expression
+$$1+1/2+1/4+\cdots+1/2^n+1/2^{n+1}$$
+All of the terms except the last can be replaced with $2-1/2^n$ therefore
+$$1+1/2+1/4+\cdots+1/2^n+1/2^{n+1}=2-1/2^n+1/2^{n+1}=2-(1/2^n)(1-1/2)
+=2-1/2^{n+1}$$
+Hence the formula is true for $n+1$ whenever it is true for $n$.
+Therefore by induction the formula is true for all $n\ge1$.
+

doc/ross/ross-14.1.tex

@@ -0,0 +1,64 @@
+\beginsection{14.1}
+
+Determine which of the following series converge.
+Justify your answers.
+
+\medskip
+(a) $\sum n^4/2^n$
+\medskip
+Try the ratio test.
+$$
+{a_{n+1}\over a_n}
+={(n+1)^4/2^{n+1}\over n^4/2^n}
+={(n+1)^4\over2^{n+1}}\times{2^n\over n^4}
+={(n+1)^4\over2n^4}\rightarrow{1\over2}<1
+$$
+So by the ratio test $\sum n^4/2^n$ converges.
+
+\medskip
+(b) $\sum2^n/n!$
+\medskip
+Try the ratio test.
+$$
+{a_{n+1}\over a_n}={2^{n+1}/(n+1)!\over 2^n/n!}
+={2^{n+1}\over(n+1)!}\times{n!\over2^n}={2\over n+1}<1
+$$
+So by the ratio test $\sum2^n/n!$ converges.
+
+\medskip
+(c) $\sum n^2/3^n$
+\medskip
+Try the ratio test.
+$$
+{a_{n+1}\over a_n}
+={(n+1)^2/3^{n+1}\over n^2/3^n}
+={(n+1)^2\over 3^{n+1}}\times{3^n\over n^2}
+={(n+1)^2\over3n^2}\rightarrow{1\over3}<1
+$$
+So by the ratio test $\sum n^2/3^n$ converges.
+
+\medskip
+(d) $\sum n!/(n^4+3)$
+\medskip
+Try the ratio test.
+$$
+{a_{n+1}\over a_n}
+={(n+1)!/((n+1)^4+3)\over n!/(n^4+3)}
+={(n+1)!\over((n+1)^4+3)}\times{n^4+3\over n!}
+\rightarrow n+1>1
+$$
+So by the ratio test $\sum n!/(n^4+3)$ diverges.
+
+\medskip
+(e) $\sum\cos^2 n/n^2$
+\medskip
+Use the comparison test.
+Since $\sum1/n^2$ converges and
+$|\cos^2 n/n^2|\le1/n^2$ for all $n$, $\sum\cos^2 n/n^2$ must also converge.
+
+\medskip
+(f) $\sum_{n=2}^\infty1/(\log n)$
+\medskip
+Use the comparison test.
+Since $\sum(1/n)$ diverges and $1/(\log n)>1/n$ for all $n>1$,
+$\sum_{n=2}^\infty1/(\log n)$ must also diverge.

doc/ross/ross-17.3.tex

@@ -0,0 +1,90 @@
+\beginsection 17.3
+
+Accept on faith that the following familiar functions are continuous on
+their domains: $\sin x$, $\cos x$, $e^x$, $2^x$, $\log_e x$ for
+$x>0$, $x^p$ for $x>0$ [$p$ any real number].
+Use these facts and theorems in this section to prove that the following
+functions are continuous.
+
+\medskip
+\noindent
+(a) $\log_e(1+\cos^4 x)$
+
+\medskip
+\itemitem{}
+The function $\cos^4 x$ is continuous by Theorem 17.4 (ii),
+product of continuous functions.
+The function $1+\cos^4 x$ is continuous by Theorem 17.4 (i),
+sum of continuous functions.
+The function $\log_e(1+\cos^4 x)$ is continuous by Theorem 17.5,
+composition of continuous functions.
+
+\medskip
+\noindent
+(b) $[\sin^2x+\cos^6x]^\pi$
+
+\medskip
+\itemitem{}
+This is equivalent to $\exp[\pi\log_e(\sin^2x+\cos^6x)]$ which is continuous by
+Theorems 17.4 and 17.5.
+Now show that $\sin^2x+\cos^6x>0$ to ensure the domain requirement of $\log_e$.
+We have $\sin^2x\ge0$, $\cos^6x\ge0$ by even exponents.
+By $\sin^2x+\cos^2x=1$, $\sin^2x$ and $\cos^2x$ cannot both be zero for the
+same $x$.
+If $\cos^2x\ne0$ then $\cos^6x\ne0$.
+Therefore $\sin^2x+\cos^6x>0$.
+
+\medskip
+\noindent
+(c) $2^{x^2}$
+
+\medskip
+\itemitem{}
+The function $x^2$ is continuous by $x^2=xx$ and
+Theorem 17.4 (ii), product of continuous functions.
+By Theorem 17.5, composition of continuous fuctions, $2^{x^2}$ is continuous.
+
+\medskip
+\noindent
+(d) $8^x$
+
+\medskip
+\itemitem{}
+This is equivalent to $\exp(x\log_e8)$ which is continuous by Theorems 17.4 and 17.5.
+
+\medskip
+\noindent
+(e) $\tan x$ for $x\ne$ odd multiple of $\pi/2$.
+
+\medskip
+\itemitem{}
+$\tan x=\sin x/\cos x$ is continuous by Theorem 17.4 (iii), ratio of
+continuous functions.
+The odd multiple restriction ensures $\cos x\ne0$.
+
+\medskip
+\noindent
+(f) $x\sin(1/x)$ for $x\ne0$
+
+\medskip
+\itemitem{}
+$1/x=\exp(-\log_ex)$ is continuous for $x>0$.
+For $x<0$, $1/x=-1/|x|$ so $1/x$ is also continous for $x<0$.
+Finally, $x\sin(1/x)$ is continuous by Theorems 17.4 and 17.5.
+
+\medskip
+\noindent
+(g) $x^2\sin(1/x)$ for $x\ne0$
+
+\medskip
+\itemitem{}
+$x^2\sin(1/x)=x[x\sin(1/x)]$, see (f).
+
+\medskip
+\noindent
+(h) $(1/x)\sin(1/x^2)$ for $x\ne0$
+
+\medskip
+\itemitem{}
+The only thing new here is $1/x^2$ which is continuous by $1/x^2=(1/x)(1/x)$
+and Theorem 17.4 (ii). See (f) for continuity of $1/x$.

doc/ross/ross-17.4.tex

@@ -0,0 +1,98 @@
+\beginsection 17.4
+
+Prove that the function $\sqrt x$ is continuous on its domain
+$[0,\infty)$. {\it Hint:} Apply Example 5 in \S8.
+
+\medskip
+
+Step 1. We want to convert $|f(x)-f(x_0)|$ to an expression involving $|x-x_0|$.
+The trick from \S8 is ``irrationalize the denominator.''
+
+\medskip
+
+$\displaystyle{
+|f(x)-f(x_0)|
+=\left|\sqrt x-\sqrt{x_0}\right|
+=\left|{(\sqrt x-\sqrt{x_0})(\sqrt x+\sqrt{x_0})\over\sqrt x+\sqrt{x_0}}\right|
+=\left|{x-x_0\over\sqrt x+\sqrt{x_0}}\right|
+={|x-x_0|\over\sqrt x+\sqrt{x_0}}
+}$
+
+\medskip
+
+Step 2. We want to get the $\sqrt x$ out of the denominator so that we can
+solve for $|x-x_0|$.
+One thing we can use to our advantage is that we don't necessarily have to use
+the exact representation of $|f(x)-f(x_0)|$.
+All we really need to say is that $|f(x)-f(x_0)|$ is less than something.
+We can do that if we can find an expression that is less than
+$\sqrt x+\sqrt{x_0}$ since making the denominator smaller makes the right
+side of the equation larger.
+We notice that $\sqrt{x_0}\le\sqrt x+\sqrt{x_0}$ so we can write
+
+\medskip
+
+$\displaystyle{
+|f(x)-f(x_0)|
+\le{|x-x_0|\over\sqrt{x_0}}
+}$
+
+\medskip
+
+Step 3. We want to arrange for $|f(x)-f(x_0)|$ to be less than some epsilon
+so we write
+
+\medskip
+
+$\displaystyle{
+|f(x)-f(x_0)|
+\le{|x-x_0|\over\sqrt{x_0}}<\epsilon
+}$
+
+
+\medskip
+
+Step 4. Solving for $|x-x_0|$ we have
+
+\medskip
+
+$\displaystyle{
+|x-x_0|<\epsilon\cdot\sqrt{x_0}
+}$
+
+\medskip
+
+So if we choose $\delta=\epsilon\cdot\sqrt x_0$ then $|x-x_0|<\delta$ implies
+that $|f(x)-f(x_0)|<\epsilon$.
+
+\medskip
+
+Step 5. Note that the above fails for $x_0=0$.
+We have to show by another means that $\sqrt x$ is continuous at zero.
+First we write
+
+\medskip
+
+$\displaystyle{
+|f(x)-f(0)|
+=\left|\sqrt x-\sqrt{0}\right|
+=\sqrt x<\epsilon
+}$
+
+\medskip
+
+Solving for $x$ we have
+
+\medskip
+
+$\displaystyle{
+x<\epsilon^2
+}$
+
+\medskip
+
+So if we choose $\delta=\epsilon^2$ then
+$|x-x_0|=x<\delta$ implies that $|f(x)-f(x_0)|=\sqrt x<\epsilon$.
+Note that we can say $|x-x_0|=x$ when $x_0=0$ because the domain we are using
+is $[0,\infty)$.
+

doc/ross/ross-17.5.tex

@@ -0,0 +1,18 @@
+\beginsection 17.5
+
+(a) Prove that if $m\in N$, then the function $f(x)=x^m$ is continuous
+on $R$.
+\medskip
+The function $f(x)=x$ is continuous. Theorem 17.4 (ii) tells us that the product
+of continous functions is continuous.
+Since $x^m$ is the product of $m$ continuous functions, $x^m$ is continuous.
+
+\medskip
+(b) Prove that every {\it polynomial function}
+$p(x)=a_0+a_1x+\cdots+a_nx^n$ is continuous on $R$.
+\medskip
+By Theorem 17.4 (ii) and Exercise 17.5 (a) above, each term in the polynomial
+is continuous.
+By Theorem 17.4 (i), the sum of the terms in the polynomial is continuous.
+Therefore, $p(x)$ is continuous.
+

doc/ross/ross-18.1.tex

@@ -0,0 +1,12 @@
+\beginsection 18.1
+
+Let $f$ be as in Theorem 18.1. Show that if $-f$ assumes its maximum at
+$x_0\in[a,b]$, then $f$ assumes its minimum at $x_0$.
+
+\medskip
+Since $x_0$ is a maximum for $-f$ we have
+$$-f(x_0)\ge-f(x).$$
+Multiply through by $-1$ to get
+$$f(x_0)\le f(x).$$
+We observe that $x_0$ is a minimum for $f$.
+

doc/ross/ross-18.2.tex

@@ -0,0 +1,25 @@
+\beginsection 18.2
+
+Reread the proof of Theorem 18.1 with $[a,b]$ replaced by $(a,b)$.
+Where does it break down? Discuss.
+
+\medskip
+
+Try going through the proof line by line replacing $[a,b]$ with $(a,b)$.
+
+1. Assume that $f$ is not bounded on $(a,b)$. {\it Ok.}
+
+2. Then to each $n\in N$ there corresponds an $x_n\in (a,b)$ such that
+$|f(x_n)|>n$. {\it Ok, because f is unbounded.}
+
+3. By the B-W Theorem 11.5, $(x_n)$ has a subsequence $(x_{n_k})$ that converges
+to some real number $x_0$.
+{\it Ok, because $x_n$ is bounded on $(a,b)$.}
+
+4. The number $x_0$ must also belong to the open interval $(a,b)$.
+{\it Uh-oh, here's a problem, $x_0$ may be equal to $a$ or $b$ in the limit.}
+
+5. Since $f$ is continuous at $x_0$, we have $\lim f(x_{n_k})=f(x_0)$,
+but we also have $\lim|f(x_{n_k})|=+\infty$ which is a contradiction.
+{\it If $x_0$ is equal to $a$ or $b$ then $f$ doesn't have to be continuous
+there so we cannot say that $\lim f(x_{n_k})=f(x_0)$.}

doc/ross/ross-18.3.tex

@@ -0,0 +1,48 @@
+\beginsection 18.3
+
+Use calculus to find the maximum and minimum of
+$f(x)=x^3-6x^2+9x+1$ on $[0,5)$.
+
+\medskip
+
+The maximum and minimum occur where $df/dx=0$.
+$$df/dx=3x^2-12x+9=0$$
+This equation has two solutions, $x=1$ and $x=3$.
+That was easy. Now we have to work harder.
+
+1. Use the ``first derivative test'' to see if $x=1$ is a minimum or maximum.
+Since $f'(0)>0$ and $f'(2)<0$ we conclude that $x=1$ is a local maximum.
+
+2. Repeat for $x=3$.
+Since $f'(2)<0$ and $f'(4)>0$ we conclude that $x=3$ is a local minimum.
+
+3. Now compute $f$ at $x=0,1,3,5$ and compare.
+
+$f(0)=1$
+
+$f(1)=5$
+
+$f(3)=1$
+
+$f(5)=21$.
+
+We observe that $f(5)=21$ is the maximum 
+but since the interval is open on the 5 side, we conclude that $f$
+has no maximum.
+The function $f$ has two minima, one at $f(0)$ and the other at $f(1)$.
+
+4. I just realized that I could have skipped steps 1 and 2.
+The computation of $f$ at $x=0,1,3,5$ is all that's needed.
+Also, the ``second derivative test'' is an easier way to find out if $x_0$ is a
+local minimum or maximum:
+
+$f''(x_0)>0$ means $x_0$ is a local minimum,
+
+$f''(x_0)<0$ means $x_0$ is a local maximum.
+
+
+\bigskip
+BTW, {\tt gnuplot} was a big help in solving the problem.
+
+\centerline{\tt plot [0:5] f(x)=x**3-6*x**2+9*x+1, f(x), g(x)=3*x**2-12*x+9, g(x)}
+

doc/ross/ross-18.5.tex

@@ -0,0 +1,23 @@
+\beginsection 18.5
+
+(a) Let $f$ and $g$ be continuous functions on $[a,b]$ such that
+$f(a)\ge g(a)$ and $f(b)\le g(b)$. Prove that $f(x_0)=g(x_0)$ for at least
+one $x_0$ in $[a,b]$.
+\medskip
+Define a function $h=f-g.$
+The function $h$ is continuous because $f$ and $g$ are continuous
+(theorem 17.4, p. 92).
+We have $h(a)\ge0$ and $h(b)\le0$.
+Now we can apply IVT and assert that some $x_0$ exists at which $h(x_0)=0$
+(simply replace $y$ with 0 in the theorem).
+Since $h(x_0)=f(x_0)-g(x_0)=0$, we have proved that $f(x_0)=g(x_0)$.
+
+\medskip
+(b) Show that Example 1 can be viewed as a special case of part (a).
+\medskip
+Example 1 includes a ``little trick'' if defining $g(x)=f(x)-x$.
+(There's the hint for the solution of 18.5 (a) above.)
+The $x$ can be considered a function so we have the difference of two
+functions as in part (a). The example is a special case because we have the
+specific function $x$.
+

doc/ross/ross-18.6.tex

@@ -0,0 +1,21 @@
+\beginsection 18.6
+
+Prove that $x=\cos x$ for some $x$ in $(0,\pi/2)$.
+
+\medskip
+
+Solution: We have $\cos(0)=1$ and $\cos(\pi/2)=0$.
+The problem is that the interval is open, not closed as required by IVT.
+However, since $\cos(0)\ne0$ and $\cos(\pi/2)\ne\pi/2$ we know that
+$\cos x=x$ does not occur at the endpoints.
+Therefore it is safe to use the closed interval $[0,\pi/2]$.
+%(Is this o.k. or b.s.?)
+Use the trick from Example 1: $f(x)=x-\cos x$.
+We have $f(0)=0-1=-1$ and $f(\pi/2)=\pi/2-0=\pi/2$ and
+therefore $f(0)<0<f(\pi/2)$.
+By IVT there exits an $x_0$ such that $f(x_0)=0$.
+Because $f(x_0)=x_0-\cos x_0=0$, we have a point $x_0$
+at which $x=\cos x$.
+We have already shown that $x\ne\cos x$ at either $0$ or $\pi/2$ so
+$x=\cos x$ must occur on the interval $(0,\pi/2)$.
+

doc/ross/ross-18.7.tex

@@ -0,0 +1,13 @@
+\beginsection 18.7
+
+Prove that $x2^x=1$ for some $x$ in $(0,1)$.
+
+\medskip
+
+We define the function $f(x)=x2^x$.
+For this function we have $f(0)=0$ and $f(1)=2$.
+Since $f(0)<1<f(1)$, by IVT $f(x)=1$ must occur somewhere in the
+closed interval $[0,1]$.
+We have $f(0)\ne1$ and $f(1)\ne1$ so $x2^x=1$ must occur somewhere in the
+open interval $(0,1)$.
+

doc/ross/ross-18.8.tex

@@ -0,0 +1,14 @@
+\beginsection 18.8
+
+Suppose that $f$ is a real-valued continuous function in $R$ and that
+$f(a)f(b)<0$ for some $a,b\in R$.
+Prove that there exists $x$ between $a$ and $b$ such that $f(x)=0$.
+
+\medskip
+
+From $f(a)f(b)<0$ we conclude that
+$f(a)$ and $f(b)$ have opposite signs and neither is zero.
+Therefore either $f(a)<0<f(b)$ or $f(b)<0<f(a)$.
+In both cases we have the existence of $f(x)=0$ and $x\ne a,b$ by
+the intermediate value theorem.
+

doc/ross/ross-19.2.tex

@@ -0,0 +1,53 @@
+\beginsection 19.2
+
+Prove that each of the following functions is uniformly continuous on the
+indicated set by directly verifying the $\epsilon$-$\delta$ property in
+Definition 19.1.
+\medskip
+To verify the $\epsilon$-$\delta$ property, start with the expression
+$|x-y|<\delta$ and then find a way to transform it into the expression
+$|f(x)-f(y)|<\epsilon$.
+
+\medskip
+(a) $f(x)=3x+11$ on $R$.
+\medskip
+We have
+$$|f(x)-f(y)|=|3x+11-3y-11|=3|x-y|$$
+Starting with
+$$|x-y|<\delta$$
+we can make the substitution
+$$|f(x)-f(y)|<3\delta$$
+Then for $\epsilon=3\delta$ we have
+$$|f(x)-f(y)|<\epsilon$$
+
+\medskip
+(b) $f(x)=x^2$ on $[0,3]$
+\medskip
+We have
+$$|f(x)-f(y)|=|x^2-y^2|=|x+y|\cdot|x-y|$$
+and
+$$|f(x)-f(y)|\le6|x-y|$$
+on the interval $[0,3]$.
+Starting with
+$$|x-y|<\delta$$
+we can make the substitution
+$$|f(x)-f(y)|<6\delta$$
+Then for $\epsilon=6\delta$ we have
+$$|f(x)-f(y)|<\epsilon$$
+
+\medskip
+(c) $f(x)=1/x$ on $[{1\over2},\infty)$
+\medskip
+We have
+$$|f(x)-f(y)|=\left|{1\over x}-{1\over y}\right|
+=\left|{y-x}\over xy\right|
+={|x-y|\over|xy|}$$
+and
+$$|f(x)-f(y)|\le4|x-y|$$
+on the interval $[{1\over2},\infty)$.
+Starting with
+$$|x-y|<\delta$$
+we can make the substitution
+$$|f(x)-f(y)|<4\delta$$
+Then for $\epsilon=4\delta$ we have
+$$|f(x)-f(y)|<\epsilon$$

doc/ross/ross-2.1.tex

@@ -0,0 +1,25 @@
+\beginsection 2.1
+
+Show that $\sqrt3$, $\sqrt5$, $\sqrt7$, $\sqrt{24}$, and
+$\sqrt{31}$ are not rational numbers.
+
+\medskip
+The solution to $x^2-3=0$ is $\sqrt3$.
+By the Rational Zeros Theorem, if $\sqrt3=p/q$ then
+$q$ divides $a_1=1$ and $p$ divides $a_0=3$.
+Since $q$ must be $\pm1$, the only possible values for
+$p/q$ are $\pm1$ and $\pm3$, none of which are
+solutions to $x^2-3=0$.
+Therefore $\sqrt3$ must be irrational.
+
+$\sqrt5$: $p/q=\pm1,\pm5$, none of which are
+solutions to $x^2-5=0$.
+
+$\sqrt7$: $p/q=\pm1,\pm7$, none of which are
+solutions to $x^2-7=0$.
+
+$\sqrt24$: $p/q=\pm1,\pm2,\pm3,\pm4,\pm6,\pm8,\pm12$,
+none of which are solutions to $x^2-24=0$.
+
+$\sqrt31$: $p/q=\pm1,\pm31$, none of which are
+solutions to $x^2-31=0$.

doc/ross/ross-2.2.tex

@@ -0,0 +1,23 @@
+\beginsection{2.2}
+
+Show that $2^{1/3}$, $5^{1/7}$, and $(13)^{1/4}$ do not represent
+rational numbers.
+
+\medskip
+$2^{1/3}$ is a solution to the polynomial $x^3-2=0$.
+By the Rational Zeroes Theorem, the only possible rational solutions are
+$\pm1$ and $\pm2$, neither of which is an actual solution.
+Therefore, $2^{1/3}$ is not a rational number.
+
+\medskip
+$5^{1/7}$ is a solution to the polynomial $x^7-5=0$.
+By the Rational Zeroes Theorem, the only possible rational solutions are
+$\pm1$ and $\pm5$, neither of which is an actual solution.
+Therefore, $5^{1/7}$ is not a rational number.
+
+\medskip
+$13^{1/4}$ is a solution to the polynomial $x^4-13=0$.
+By the Rational Zeroes Theorem, the only possible rational solutions are
+$\pm1$ and $\pm13$, neither of which is an actual a solution.
+Therefore, $13^{1/4}$ is not a rational number.
+

doc/ross/ross-2.3.tex

@@ -0,0 +1,17 @@
+\beginsection{2.3}
+
+Show that $(2+\sqrt2)^{1/2}$ does not represent a rational number.
+
+\medskip
+$(2+\sqrt2)^{1/2}$ is a solution to the equation
+$x^2-2=\sqrt2$.
+This can be rewritten as
+$$(x^2-2)^2=2$$
+which expands to
+$$x^4-4x^2+4=2$$
+Hence $(2+\sqrt2)^{1/2}$ is a solution to the polynomial
+$$x^4-4x^2+2=0$$
+By the Rational Zeroes Theorem the only possible rational solutions
+are $\pm1$ and $\pm2$, neither of which is an actual solution.
+Therefore $(2+\sqrt2)^{1/2}$ is not a rational number.
+

doc/ross/ross-28.15.tex

@@ -0,0 +1,62 @@
+\beginsection{28.15}
+
+Prove Leibniz' rule
+$$(fg)^{(n)}=\sum_{k=0}^n\left({n\atop k}\right)f^{(k)}(a)g^{(n-k)}(a)$$
+{\it Hint:} Use mathematical induction.
+For $n=1$, apply Theorem 28.3(iii).
+
+\medskip
+The following equalities will be used.
+$$\left({n\atop k}\right)+\left({n\atop k-1}\right)=\left({n+1\atop k}\right)
+\eqno\hbox{(A)}$$
+$$\left({n\atop0}\right)=1\eqno\hbox{(B)}$$
+$$\left({n\atop n}\right)=1\eqno\hbox{(C)}$$
+Theorem 28.3(iii) is not really needed.
+Just using the case $n=0$ for induction step 1.
+For induction step 2, show that Leibniz' rule is true for $n+1$ whenever it is
+true for $n$.
+$$\eqalign{
+(fg)^{(n+1)}&=((fg)^{(n)})^\prime\cr
+&=\left(\sum_{k=0}^n\left({n\atop k}\right)
+f^{(k)}g^{(n-k)}\right)^\prime\cr
+&=\sum_{k=0}^n\left({n\atop k}\right)
+\left[f^{(k)}\left(g^{(n-k)}\right)^\prime
++\left(f^{(k)}\right)^\prime g^{(n-k)}\right]\cr
+&=\sum_{k=0}^n\left({n\atop k}\right)
+\left[f^{(k)}g^{(n+1-k)}+f^{(k+1)}g^{(n-k)}\right]\cr
+&=\sum_{k=0}^n\left({n\atop k}\right)f^{(k)}g^{(n+1-k)}
++\sum_{k=0}^n\left({n\atop k}\right)f^{(k+1)}g^{(n-k)}
+}$$
+Use (B) to extract the $k=0$ term of the sum on the left and use (C)
+to extract the $k=n$ term of the sum on the right.
+$$\eqalign{
+(fg)^{(n+1)}&=
+f^{(0)}g^{(n+1)}
++\sum_{k=1}^n\left({n\atop k}\right)f^{(k)}g^{(n+1-k)}
++\sum_{k=0}^{n-1}\left({n\atop k}\right)f^{(k+1)}g^{(n-k)}
++f^{(n+1)}g^{(0)}
+}$$
+Here is the main trick.
+The index $k$ is just a dummy variable that takes on a range of values.
+We can shift the range of $k$ if we balance the math where $k$ is used.
+For the sum on the right, change $k$ so that it runs from
+1 to $n$.
+$$\eqalign{
+(fg)^{(n+1)}&=
+f^{(0)}g^{(n+1)}
++\sum_{k=1}^n\left({n\atop k}\right)f^{(k)}g^{(n+1-k)}
++\sum_{k=1}^n\left({n\atop k-1}\right)f^{(k)}g^{(n+1-k)}
++f^{(n+1)}g^{(0)}
+}$$
+Note that the $n$ in the binomial factor does not change because $n$ does
+not depend on $k$.
+Next use (A) to combine the two sums.
+$$\eqalign{
+(fg)^{(n+1)}&=
+f^{(0)}g^{(n+1)}
++\sum_{k=1}^n\left(n+1\atop k\right)f^{(k)}g^{(n+1-k)}
++f^{(n+1)}g^{(0)}
+}$$
+Finally, use (B) and (C) to combine all the terms.
+$$(fg)^{(n+1)}=\sum_{k=0}^{n+1}
+\left({n+1\atop k}\right)f^{(k)}g^{(n+1-k)}$$

doc/ross/ross-28.2.tex

@@ -0,0 +1,23 @@
+\beginsection 28.2
+
+Use the {\it definition} of derivative to calculate the derivatives of the
+following functions at the indicated points.
+
+\medskip
+(a) $f(x)=x^3$ at $x=2$
+\medskip
+$$\eqalign{
+f^\prime&=\lim_{x\rightarrow2}{f(x)-f(2)\over x-2}\cr
+&=\lim_{x\rightarrow2}{x^3-8\over x-2}\cr
+&=\lim_{x\rightarrow2}x^2+2x+4\cr
+&=12
+}$$
+
+\medskip
+(b) $g(x)=x+2$ at $x=a$
+\medskip
+$$\eqalign{
+g^\prime(a)&=\lim_{x\rightarrow a}{g(x)-g(a)\over x-a}\cr
+&=\lim_{x\rightarrow a}{x+2-a-2\over x-a}\cr
+&=1
+}$$

doc/ross/ross-29.11.tex

@@ -0,0 +1,13 @@
+\beginsection 29.11
+
+Show that $\sin x\le x$ for all $x\ge0$.
+{\it Hint:} Show that $f(x)=x-\sin x$ is increasing on $[0,\infty).$
+
+\medskip
+$$f^\prime(x)=1+\cos x$$
+Since the range of $\cos x$ is $[-1,1]$, the range of $f^\prime(x)$ is
+$[0,2]$. By corollary 29.7 we conclude that $f$ is an increasing function
+because $f^\prime(x)\ge0$ for all $x\in[0,\infty)$.
+Since $f(0)=0$ and $f$ is an increasing function we must have
+$f(x)\ge0$ for $x\in[0,\infty)$.
+Therefore $x\ge\sin x$.

doc/ross/ross-3.1.tex

@@ -0,0 +1,43 @@
+\beginsection{3.1}
+
+(a) Which of the properties A1--A4, M1--M4, DL, O1--O5 fail for the
+natural numbers $N$?
+
+\medskip
+A1. $a+(b+c)=(a+b+c)$
+
+A2. $a+b=b+a$
+
+A3. $a+0=a$
+
+A4. For each $a$ there is as element $-a$ such that $a+(-a)=0$.
+{\it Fails, no negative numbers.}
+
+M1. $a(bc)=(ab)c$
+
+M2. $ab=ba$
+
+M3. $a\cdot1=a$
+
+M4. For each $a\ne0$ there is an element $a^{-1}$ such that $aa^{-1}=1$.
+{\it Fails, no multiplicative inverse.}
+
+DL. $a(b+c)=ab+ac$
+
+O1. Either $a\le b$ or $b\le a$.
+
+O2. If $a\le b$ and $b\le a$ then $a=b$.
+
+O3. If $a\le b$ and $b\le c$ then $a\le c$.
+
+O4. If $a\le b$ then $a+c\le b+c$.
+
+O5. If $a\le b$ and $0\le c$ then $ac\le bc$.
+{\it Interesting since $0\not\in N$.}
+
+\medskip
+(b) Which of these properties fail for the set of integers $Z$?
+
+\medskip
+M4, no multiplicative inverse.
+

doc/ross/ross-3.5.tex

@@ -0,0 +1,21 @@
+\beginsection{3.5}
+
+(a) Show that $|b|\le a$ if and only if $-a\le b\le a$.
+
+\medskip
+First prove the implication.
+Let $|b|\le a$.
+Then by Theorem 3.2, $-a\le-|b|$.
+By $-|b|\le b\le|b|$ we have
+$$-a\le-|b|\le b\le|b|\le a$$
+Hence $|b|\le a$ implies $-a\le b\le a$.
+
+\medskip
+Now prove the converse.
+By hypothesis we have
+$$b\le a\eqno(1)$$
+Also by hypothesis we have $-a\le b$ which implies
+$$-b\le a\eqno(2)$$
+Therefore by (1) and (2) we have
+$$|b|\le a$$
+

doc/ross/ross-32.1.tex

@@ -0,0 +1,4 @@
+\beginsection 32.1
+
+Find the upper and lower Darboux integrals for $f(x)=x^3$ on the interval
+$[0,b]$. {\it Hint:} Exercise 1.3 and Example 1 in \S1 will be useful.

doc/ross/ross-7.1.tex

@@ -0,0 +1,39 @@
+\beginsection{7.1}
+
+Write out the first five terms of the following sequences.
+
+\medskip
+(a) $s_n=1/(3n+1)$
+\medskip
+$s_1={1\over4}$,
+$s_2={1\over7}$,
+$s_3={1\over10}$,
+$s_4={1\over13}$,
+$s_5={1\over16}$
+
+\medskip
+(b) $b_n=(3n+1)/(4n-1)$
+\medskip
+$b_1={4\over3}$,
+$b_2={7\over7}$,
+$b_3={10\over11}$,
+$b_4={13\over15}$,
+$b_5={16\over19}$
+
+\medskip
+(c) $c_n=n/3^n$
+\medskip
+$c_1={1\over3}$,
+$c_2={2\over9}$,
+$c_3={3\over27}$,
+$c_4={4\over81}$,
+$c_5={5\over243}$
+
+\medskip
+(d) $\sin(n\pi/4)$
+\medskip
+$\sqrt2\over2$,
+$1$,
+$\sqrt2\over2$,
+$0$,
+$-{\sqrt2\over2}$

doc/ross/ross-7.2.tex

@@ -0,0 +1,17 @@
+\beginsection{7.2}
+
+For each sequence in Exercise 7.1, determine whether it converges.
+If it converges, give its limit.
+No proofs are required.
+
+\medskip
+(a) $s_n=1/(3n+1)$ converges to 0.
+
+\medskip
+(b) $b_n=(3n+1)/(4n-1)$ converges to $3\over4$.
+
+\medskip
+(c) $c_n=n/3^n$ converges to 0.
+
+\medskip
+(d) $\sin(n\pi/4)$ does not converge.

doc/ross/ross.tex

@@ -0,0 +1,45 @@
+% export CVS_RSH=ssh
+% export CVSROOT=:ext:gweigt@cvs.sf.net:/cvsroot/eigenmath
+
+\parindent=0pt
+
+\input ross-1.1.tex
+\input ross-1.2.tex
+\input ross-1.3.tex
+\input ross-1.4.tex
+\input ross-1.5.tex
+
+\input ross-2.1.tex
+\input ross-2.2.tex
+\input ross-2.3.tex
+
+\input ross-3.1.tex
+\input ross-3.5.tex
+
+\input ross-7.1.tex
+\input ross-7.2.tex
+
+\input ross-14.1.tex
+
+\input ross-17.3.tex
+\input ross-17.4.tex
+\input ross-17.5.tex
+
+\input ross-18.1.tex
+\input ross-18.2.tex
+\input ross-18.3.tex
+\input ross-18.5.tex
+\input ross-18.6.tex
+\input ross-18.7.tex
+\input ross-18.8.tex
+
+\input ross-19.2.tex
+
+\input ross-28.2.tex
+\input ross-28.15.tex
+
+\input ross-29.11.tex
+
+\input ross-32.1.tex
+
+\end

lisp/example.lisp

@@ -0,0 +1,64 @@
+"Bondi metric..."
+
+; coordinate system
+
+(setq x0 u)
+(setq x1 r)
+(setq x2 theta)
+(setq x3 phi)
+
+; U, V, beta and gamma are functions of u, r and theta
+
+(setq g_uu (sum
+  (product
+    (V u r theta)
+    (power r -1)
+    (power e (product 2 (beta u r theta)))
+  )
+  (product
+    -1
+    (power (U u r theta) 2)
+    (power r 2)
+    (power e (product 2 (gamma u r theta)))
+  )
+))
+
+(setq g_ur (product 2 (power e (product 2 (beta u r theta)))))
+
+(setq g_utheta (product
+  2
+  (U u r theta)
+  (power r 2)
+  (power e (product 2 (gamma u r theta)))
+))
+
+(setq g_thetatheta (product
+  -1
+  (power r 2)
+  (power e (product 2 (gamma u r theta)))
+))
+
+(setq g_phiphi (product
+  -1
+  (power r 2)
+  (power e (product -2 (gamma u r theta)))
+  (power (sin theta) 2)
+))
+
+; metric tensor
+
+(setq gdd (sum
+  (product g_uu (tensor u u))
+  (product g_ur (tensor u r))
+  (product g_ur (tensor r u))
+  (product g_utheta (tensor u theta))
+  (product g_utheta (tensor theta u))
+  (product g_thetatheta (tensor theta theta))
+  (product g_phiphi (tensor phi phi))
+))
+
+(gr) ; compute g, guu, GAMUDD, RUDDD, RDD, R, GDD, GUD and GUU
+
+"Is the Einstein tensor GUU divergence-free?"
+
+(zerop (contract23 (covariant-derivative-of-up-up GUU)))

lisp/example2.lisp

@@ -0,0 +1,220 @@
+; Page references are for the book "Gravitation."
+; generic metric
+(setq gdd (sum
+  (product (g00) (tensor x0 x0))
+  (product (g11) (tensor x1 x1))
+  (product (g22) (tensor x2 x2))
+  (product (g33) (tensor x3 x3))
+))
+(gr) ; compute g, guu, GAMUDD, RUDDD, RDD, R, GDD, GUD and GUU
+; generic vectors
+(setq u (sum
+  (product (u0) (tensor x0))
+  (product (u1) (tensor x1))
+  (product (u2) (tensor x2))
+  (product (u3) (tensor x3))
+))
+(setq v (sum
+  (product (v0) (tensor x0))
+  (product (v1) (tensor x1))
+  (product (v2) (tensor x2))
+  (product (v3) (tensor x3))
+))
+(setq w (sum
+  (product (w0) (tensor x0))
+  (product (w1) (tensor x1))
+  (product (w2) (tensor x2))
+  (product (w3) (tensor x3))
+))
+; how to antisymmetrize three indices (p. 86)
+(setq V (sum
+  (product 1/6 V)
+  (product 1/6 (transpose12 (transpose23 V))) ; nu lambda mu -> mu nu lambda
+  (product 1/6 (transpose23 (transpose12 V))) ; lambda mu nu -> mu nu lambda
+  (product -1/6 (transpose12 V))              ; nu mu lambda -> mu nu lambda
+  (product -1/6 (transpose23 V))              ; mu lambda nu -> mu nu lambda
+  (product -1/6 (transpose13 V))              ; lambda nu mu -> mu nu lambda
+))
+; commutator (p. 206)
+(define commutator (sum
+  (contract13 (product +1 arg1 (gradient arg2)))
+  (contract13 (product -1 arg2 (gradient arg1)))
+))
+(print "connection coefficients (p. 210)")
+(setq temp (gradient gdd))
+(setq Gamma (contract23 (product 1/2 guu (sum
+  temp
+  (transpose23 temp)
+  (product -1 (transpose12 (transpose23 temp)))
+))))
+; check
+(print (equal Gamma GAMUDD))
+; covariant derivative of a vector (p. 211)
+(define covariant-derivative (sum
+  (gradient arg)
+  (contract13 (product arg GAMUDD))
+))
+(print "divergence of einstein is zero (p. 222)")
+; covariant-derivative-of-up-up is already defined in grlib
+(setq temp (covariant-derivative-of-up-up GUU))
+(setq temp (contract23 temp)) ; sum over 2nd and 3rd indices
+(print (equal temp 0))
+(print "computing riemann tensor (p. 219)")
+(setq temp1 (gradient GAMUDD))
+(setq temp2 (contract24 (product GAMUDD GAMUDD)))
+(setq riemann (sum
+  (transpose34 temp1)
+  (product -1 temp1)
+  (transpose23 temp2)
+  (product -1 (transpose34 (transpose23 temp2)))
+))
+; check
+(print (equal riemann RUDDD))
+; p. 259
+(define covariant-derivative-of-down-down (prog (temp)
+  (setq temp (product arg GAMUDD))
+  (return (sum
+    (gradient arg)
+    (product -1 (transpose12 (contract13 temp)))
+    (product -1 (contract23 temp))
+  ))
+))
+(define covariant-derivative-of-up-down (prog (temp)
+  (setq temp (product arg GAMUDD))
+  (return (sum
+    (gradient arg)
+    (transpose12 (contract14 temp))
+    (product -1 (contract23 temp))
+  ))
+))
+(define directed-covariant-derivative
+  (contract23 (product (covariant-derivative arg1) arg2))
+)
+(print "symmetry of covariant derivative (p. 252)")
+(setq temp1 (sum
+  (directed-covariant-derivative v u)
+  (product -1 (directed-covariant-derivative u v))
+))
+(setq temp2 (commutator u v))
+(equal temp1 temp2)
+(print "covariant derivative chain rule (p. 252)")
+(setq temp1 (directed-covariant-derivative (product (f) v) u))
+(setq temp2 (sum
+  (product (f) (directed-covariant-derivative v u))
+  (product v (contract12 (product (gradient (f)) u)))
+))
+(print (equal temp1 temp2))
+(print "additivity of covariant derivative (p. 257)")
+(setq temp1 (directed-covariant-derivative u (sum v w)))
+(setq temp2 (sum
+  (directed-covariant-derivative u v)
+  (directed-covariant-derivative u w)
+))
+(print (equal temp1 temp2))
+(print "riemann is antisymmetric on last two indices (p. 286)")
+(setq temp (sum
+  (product 1/2 RUDDD)
+  (product -1/2 (transpose34 RUDDD))
+))
+(print (equal RUDDD temp))
+(print "riemann vanishes when antisymmetrized on last three indices (p. 286)")
+(setq temp (sum
+  (product 1/6 RUDDD)
+  (product 1/6 (transpose34 (transpose24 RUDDD)))
+  (product 1/6 (transpose34 (transpose23 RUDDD)))
+  (product -1/6 (transpose23 RUDDD))
+  (product -1/6 (transpose34 RUDDD))
+  (product -1/6 (transpose24 RUDDD))
+))
+(print (equal temp 0))
+(print "double dual of riemann (p. 325)")
+(setq temp (contract23 (product gdd RUDDD))) ; lower 1st index
+(setq temp (transpose34 (contract35 (product temp guu)))) ; raise 3rd index
+(setq RDDUU (contract45 (product temp guu))) ; raise 4th index
+(setq temp (product epsilon RDDUU))
+(setq temp (contract35 temp)) ; sum over mu
+(setq temp (contract34 temp)) ; sum over nu
+(setq temp (product temp epsilon))
+(setq temp (contract35 temp)) ; sum over rho
+(setq temp (contract34 temp)) ; sum over sigma
+(setq GUUDD (product -1/4 temp)) ; negative due to levi-civita tensor
+; check
+(print (equal
+  (contract13 GUUDD)
+  GUD
+))
+(print "noncommutation of covariant derivatives (p. 389)")
+(setq B (sum
+  (product (B0) (tensor x0))
+  (product (B1) (tensor x1))
+  (product (B2) (tensor x2))
+  (product (B3) (tensor x3))
+))
+(setq temp (covariant-derivative-of-up-down (covariant-derivative B)))
+(setq temp1 (sum temp (product -1 (transpose23 temp))))
+(setq temp2 (transpose23 (contract25 (product RUDDD B))))
+(print (equal temp1 temp2))
+(print "bondi metric")
+; erase any definitions for U, V, beta and gamma
+(setq U (quote U))
+(setq V (quote V))
+(setq beta (quote beta))
+(setq gamma (quote gamma))
+; erase any definitions for u, r, theta and phi
+(setq u (quote u))
+(setq r (quote r))
+(setq theta (quote theta))
+(setq phi (quote phi))
+; new coordinate system
+(setq x0 u)
+(setq x1 r)
+(setq x2 theta)
+(setq x3 phi)
+; U, V, beta and gamma are functions of u, r and theta
+(setq g_uu (sum
+  (product
+    (V u r theta)
+    (power r -1)
+    (power e (product 2 (beta u r theta)))
+  )
+  (product
+    -1
+    (power (U u r theta) 2)
+    (power r 2)
+    (power e (product 2 (gamma u r theta)))
+  )
+))
+(setq g_ur (product 2 (power e (product 2 (beta u r theta)))))
+(setq g_utheta (product
+  2
+  (U u r theta)
+  (power r 2)
+  (power e (product 2 (gamma u r theta)))
+))
+(setq g_thetatheta (product
+  -1
+  (power r 2)
+  (power e (product 2 (gamma u r theta)))
+))
+(setq g_phiphi (product
+  -1
+  (power r 2)
+  (power e (product -2 (gamma u r theta)))
+  (power (sin theta) 2)
+))
+; metric tensor
+(setq gdd (sum
+  (product g_uu (tensor u u))
+  (product g_ur (tensor u r))
+  (product g_ur (tensor r u))
+  (product g_utheta (tensor u theta))
+  (product g_utheta (tensor theta u))
+  (product g_thetatheta (tensor theta theta))
+  (product g_phiphi (tensor phi phi))
+))
+(gr) ; compute g, guu, GAMUDD, RUDDD, RDD, R, GDD, GUD and GUU
+; is covariant derivative of metric zero?
+(print (equal (covariant-derivative-of-down-down gdd) 0))
+; is divergence of einstein zero?
+(setq temp (contract23 (covariant-derivative-of-up-up GUU)))
+(print (equal temp 0))

lisp/example2.tex

@@ -0,0 +1,487 @@
+\parindent=0pt
+{\tt ;\ Page\ references\ are\ for\ the\ book\ "Gravitation."}
+
+{\tt ;\ generic\ metric}
+
+{\tt (setq\ gdd\ (sum}
+
+{\tt \ \ (product\ (g00)\ (tensor\ x0\ x0))}
+
+{\tt \ \ (product\ (g11)\ (tensor\ x1\ x1))}
+
+{\tt \ \ (product\ (g22)\ (tensor\ x2\ x2))}
+
+{\tt \ \ (product\ (g33)\ (tensor\ x3\ x3))}
+
+{\tt ))}
+
+{\tt (gr)\ ;\ compute\ g,\ guu,\ GAMUDD,\ RUDDD,\ RDD,\ R,\ GDD,\ GUD\ and\ GUU}
+
+{\tt ;\ generic\ vectors}
+
+{\tt (setq\ u\ (sum}
+
+{\tt \ \ (product\ (u0)\ (tensor\ x0))}
+
+{\tt \ \ (product\ (u1)\ (tensor\ x1))}
+
+{\tt \ \ (product\ (u2)\ (tensor\ x2))}
+
+{\tt \ \ (product\ (u3)\ (tensor\ x3))}
+
+{\tt ))}
+
+{\tt (setq\ v\ (sum}
+
+{\tt \ \ (product\ (v0)\ (tensor\ x0))}
+
+{\tt \ \ (product\ (v1)\ (tensor\ x1))}
+
+{\tt \ \ (product\ (v2)\ (tensor\ x2))}
+
+{\tt \ \ (product\ (v3)\ (tensor\ x3))}
+
+{\tt ))}
+
+{\tt (setq\ w\ (sum}
+
+{\tt \ \ (product\ (w0)\ (tensor\ x0))}
+
+{\tt \ \ (product\ (w1)\ (tensor\ x1))}
+
+{\tt \ \ (product\ (w2)\ (tensor\ x2))}
+
+{\tt \ \ (product\ (w3)\ (tensor\ x3))}
+
+{\tt ))}
+
+$$V_{[\mu\nu\lambda]}={1\over3!}
+(V_{\mu\nu\lambda}
++V_{\nu\lambda\mu}
++V_{\lambda\mu\nu}
+-V_{\nu\mu\lambda}
+-V_{\mu\lambda\nu}
+-V_{\lambda\nu\mu}
+)$$
+{\tt ;\ how\ to\ antisymmetrize\ three\ indices\ (p.\ 86)}
+
+{\tt (setq\ V\ (sum}
+
+{\tt \ \ (product\ 1/6\ V)}
+
+{\tt \ \ (product\ 1/6\ (transpose12\ (transpose23\ V)))\ ;\ nu\ lambda\ mu\ ->\ mu\ nu\ lambda}
+
+{\tt \ \ (product\ 1/6\ (transpose23\ (transpose12\ V)))\ ;\ lambda\ mu\ nu\ ->\ mu\ nu\ lambda}
+
+{\tt \ \ (product\ -1/6\ (transpose12\ V))\ \ \ \ \ \ \ \ \ \ \ \ \ \ ;\ nu\ mu\ lambda\ ->\ mu\ nu\ lambda}
+
+{\tt \ \ (product\ -1/6\ (transpose23\ V))\ \ \ \ \ \ \ \ \ \ \ \ \ \ ;\ mu\ lambda\ nu\ ->\ mu\ nu\ lambda}
+
+{\tt \ \ (product\ -1/6\ (transpose13\ V))\ \ \ \ \ \ \ \ \ \ \ \ \ \ ;\ lambda\ nu\ mu\ ->\ mu\ nu\ lambda}
+
+{\tt ))}
+
+$$[{\bf u},{\bf v}]=
+(u^\beta{v^\alpha}_{,\beta}-v^\beta{u^\alpha}_{,\beta}){\bf e}_\alpha$$
+{\tt ;\ commutator\ (p.\ 206)}
+
+{\tt (define\ commutator\ (sum}
+
+{\tt \ \ (contract13\ (product\ +1\ arg1\ (gradient\ arg2)))}
+
+{\tt \ \ (contract13\ (product\ -1\ arg2\ (gradient\ arg1)))}
+
+{\tt ))}
+
+$${\Gamma^\alpha}_{\mu\nu}=\hbox{$1\over2$}g^{\alpha\beta}
+(g_{\beta\mu,\nu}+g_{\beta\nu,\mu}-g_{\mu\nu,\beta})$$
+{\tt (print\ "connection\ coefficients\ (p.\ 210)")}
+
+{\tt (setq\ temp\ (gradient\ gdd))}
+
+{\tt (setq\ Gamma\ (contract23\ (product\ 1/2\ guu\ (sum}
+
+{\tt \ \ temp}
+
+{\tt \ \ (transpose23\ temp)}
+
+{\tt \ \ (product\ -1\ (transpose12\ (transpose23\ temp)))}
+
+{\tt ))))}
+
+{\tt ;\ check}
+
+{\tt (print\ (equal\ Gamma\ GAMUDD))}
+
+$${T^\alpha}_{;\gamma}={T^\alpha}_{,\gamma}+
+T^\mu{\Gamma^\alpha}_{\mu\gamma}$$
+{\tt ;\ covariant\ derivative\ of\ a\ vector\ (p.\ 211)}
+
+{\tt (define\ covariant-derivative\ (sum}
+
+{\tt \ \ (gradient\ arg)}
+
+{\tt \ \ (contract13\ (product\ arg\ GAMUDD))}
+
+{\tt ))}
+
+$${G^{\mu\nu}}_{;\nu}=0$$
+{\tt (print\ "divergence\ of\ einstein\ is\ zero\ (p.\ 222)")}
+
+{\tt ;\ covariant-derivative-of-up-up\ is\ already\ defined\ in\ grlib}
+
+{\tt (setq\ temp\ (covariant-derivative-of-up-up\ GUU))}
+
+{\tt (setq\ temp\ (contract23\ temp))\ ;\ sum\ over\ 2nd\ and\ 3rd\ indices}
+
+{\tt (print\ (equal\ temp\ 0))}
+
+$${R^\alpha}_{\beta\gamma\delta}=
+{\partial{\Gamma^\alpha}_{\beta\delta}\over\partial x^\gamma}-
+{\partial{\Gamma^\alpha}_{\beta\gamma}\over\partial x^\delta}+
+{\Gamma^\alpha}_{\mu\gamma}{\Gamma^\mu}_{\beta\delta}-
+{\Gamma^\alpha}_{\mu\delta}{\Gamma^\mu}_{\beta\gamma}$$
+{\tt (print\ "computing\ riemann\ tensor\ (p.\ 219)")}
+
+{\tt (setq\ temp1\ (gradient\ GAMUDD))}
+
+{\tt (setq\ temp2\ (contract24\ (product\ GAMUDD\ GAMUDD)))}
+
+{\tt (setq\ riemann\ (sum}
+
+{\tt \ \ (transpose34\ temp1)}
+
+{\tt \ \ (product\ -1\ temp1)}
+
+{\tt \ \ (transpose23\ temp2)}
+
+{\tt \ \ (product\ -1\ (transpose34\ (transpose23\ temp2)))}
+
+{\tt ))}
+
+{\tt ;\ check}
+
+{\tt (print\ (equal\ riemann\ RUDDD))}
+
+$$T_{\alpha\beta;\gamma}=T_{\alpha\beta,\gamma}
+-T_{\mu\beta}{\Gamma^\mu}_{\alpha\gamma}
+-T_{\alpha\mu}{\Gamma^\mu}_{\beta\gamma}$$
+{\tt ;\ p.\ 259}
+
+{\tt (define\ covariant-derivative-of-down-down\ (prog\ (temp)}
+
+{\tt \ \ (setq\ temp\ (product\ arg\ GAMUDD))}
+
+{\tt \ \ (return\ (sum}
+
+{\tt \ \ \ \ (gradient\ arg)}
+
+{\tt \ \ \ \ (product\ -1\ (transpose12\ (contract13\ temp)))}
+
+{\tt \ \ \ \ (product\ -1\ (contract23\ temp))}
+
+{\tt \ \ ))}
+
+{\tt ))}
+
+$${T^\alpha}_{\beta;\gamma}={T^\alpha}_{\beta,\gamma}+
+{T^\mu}_\beta{\Gamma^\alpha}_{\mu\gamma}-
+{T^\alpha}_\mu{\Gamma^\mu}_{\beta\gamma}$$
+{\tt (define\ covariant-derivative-of-up-down\ (prog\ (temp)}
+
+{\tt \ \ (setq\ temp\ (product\ arg\ GAMUDD))}
+
+{\tt \ \ (return\ (sum}
+
+{\tt \ \ \ \ (gradient\ arg)}
+
+{\tt \ \ \ \ (transpose12\ (contract14\ temp))}
+
+{\tt \ \ \ \ (product\ -1\ (contract23\ temp))}
+
+{\tt \ \ ))}
+
+{\tt ))}
+
+$$\nabla_{\bf u}{\bf v}={v^\alpha}_{;\beta}u^{\beta}$$
+{\tt (define\ directed-covariant-derivative}
+
+{\tt \ \ (contract23\ (product\ (covariant-derivative\ arg1)\ arg2))}
+
+{\tt )}
+
+$$\nabla_{\bf u}{\bf v}-\nabla_{\bf v}{\bf u}=[{\bf u},{\bf v}]$$
+{\tt (print\ "symmetry\ of\ covariant\ derivative\ (p.\ 252)")}
+
+{\tt (setq\ temp1\ (sum}
+
+{\tt \ \ (directed-covariant-derivative\ v\ u)}
+
+{\tt \ \ (product\ -1\ (directed-covariant-derivative\ u\ v))}
+
+{\tt ))}
+
+{\tt (setq\ temp2\ (commutator\ u\ v))}
+
+{\tt (equal\ temp1\ temp2)}
+
+$$\nabla_{\bf u}(f{\bf v})=
+f\nabla_{\bf u}{\bf v}+{\bf v}\partial_{\bf u}f$$
+{\tt (print\ "covariant\ derivative\ chain\ rule\ (p.\ 252)")}
+
+{\tt (setq\ temp1\ (directed-covariant-derivative\ (product\ (f)\ v)\ u))}
+
+{\tt (setq\ temp2\ (sum}
+
+{\tt \ \ (product\ (f)\ (directed-covariant-derivative\ v\ u))}
+
+{\tt \ \ (product\ v\ (contract12\ (product\ (gradient\ (f))\ u)))}
+
+{\tt ))}
+
+{\tt (print\ (equal\ temp1\ temp2))}
+
+$$\nabla_{{\bf v}+{\bf w}}{\bf u}=
+\nabla_{\bf v}{\bf u}+\nabla_{\bf w}{\bf u}$$
+{\tt (print\ "additivity\ of\ covariant\ derivative\ (p.\ 257)")}
+
+{\tt (setq\ temp1\ (directed-covariant-derivative\ u\ (sum\ v\ w)))}
+
+{\tt (setq\ temp2\ (sum}
+
+{\tt \ \ (directed-covariant-derivative\ u\ v)}
+
+{\tt \ \ (directed-covariant-derivative\ u\ w)}
+
+{\tt ))}
+
+{\tt (print\ (equal\ temp1\ temp2))}
+
+$${R^\alpha}_{\beta\gamma\delta}={R^\alpha}_{\beta[\gamma\delta]}$$
+{\tt (print\ "riemann\ is\ antisymmetric\ on\ last\ two\ indices\ (p.\ 286)")}
+
+{\tt (setq\ temp\ (sum}
+
+{\tt \ \ (product\ 1/2\ RUDDD)}
+
+{\tt \ \ (product\ -1/2\ (transpose34\ RUDDD))}
+
+{\tt ))}
+
+{\tt (print\ (equal\ RUDDD\ temp))}
+
+$${R^\alpha}_{[\beta\gamma\delta]}=0$$
+{\tt (print\ "riemann\ vanishes\ when\ antisymmetrized\ on\ last\ three\ indices\ (p.\ 286)")}
+
+{\tt (setq\ temp\ (sum}
+
+{\tt \ \ (product\ 1/6\ RUDDD)}
+
+{\tt \ \ (product\ 1/6\ (transpose34\ (transpose24\ RUDDD)))}
+
+{\tt \ \ (product\ 1/6\ (transpose34\ (transpose23\ RUDDD)))}
+
+{\tt \ \ (product\ -1/6\ (transpose23\ RUDDD))}
+
+{\tt \ \ (product\ -1/6\ (transpose34\ RUDDD))}
+
+{\tt \ \ (product\ -1/6\ (transpose24\ RUDDD))}
+
+{\tt ))}
+
+{\tt (print\ (equal\ temp\ 0))}
+
+$${G^{\alpha\beta}}_{\gamma\delta}=
+\hbox{$1\over2$}\epsilon^{\alpha\beta\mu\nu}
+{R_{\mu\nu}}^{\rho\sigma}
+\hbox{$1\over2$}\epsilon_{\rho\sigma\gamma\delta}$$
+{\tt (print\ "double\ dual\ of\ riemann\ (p.\ 325)")}
+
+{\tt (setq\ temp\ (contract23\ (product\ gdd\ RUDDD)))\ ;\ lower\ 1st\ index}
+
+{\tt (setq\ temp\ (transpose34\ (contract35\ (product\ temp\ guu))))\ ;\ raise\ 3rd\ index}
+
+{\tt (setq\ RDDUU\ (contract45\ (product\ temp\ guu)))\ ;\ raise\ 4th\ index}
+
+{\tt (setq\ temp\ (product\ epsilon\ RDDUU))}
+
+{\tt (setq\ temp\ (contract35\ temp))\ ;\ sum\ over\ mu}
+
+{\tt (setq\ temp\ (contract34\ temp))\ ;\ sum\ over\ nu}
+
+{\tt (setq\ temp\ (product\ temp\ epsilon))}
+
+{\tt (setq\ temp\ (contract35\ temp))\ ;\ sum\ over\ rho}
+
+{\tt (setq\ temp\ (contract34\ temp))\ ;\ sum\ over\ sigma}
+
+{\tt (setq\ GUUDD\ (product\ -1/4\ temp))\ ;\ negative\ due\ to\ levi-civita\ tensor}
+
+{\tt ;\ check}
+
+{\tt (print\ (equal}
+
+{\tt \ \ (contract13\ GUUDD)}
+
+{\tt \ \ GUD}
+
+{\tt ))}
+
+$${B^\mu}_{;\alpha\beta}-{B^\mu}_{;\beta\alpha}=
+{R^\mu}_{\nu\beta\alpha}B^\nu$$
+{\tt (print\ "noncommutation\ of\ covariant\ derivatives\ (p.\ 389)")}
+
+{\tt (setq\ B\ (sum}
+
+{\tt \ \ (product\ (B0)\ (tensor\ x0))}
+
+{\tt \ \ (product\ (B1)\ (tensor\ x1))}
+
+{\tt \ \ (product\ (B2)\ (tensor\ x2))}
+
+{\tt \ \ (product\ (B3)\ (tensor\ x3))}
+
+{\tt ))}
+
+{\tt (setq\ temp\ (covariant-derivative-of-up-down\ (covariant-derivative\ B)))}
+
+{\tt (setq\ temp1\ (sum\ temp\ (product\ -1\ (transpose23\ temp))))}
+
+{\tt (setq\ temp2\ (transpose23\ (contract25\ (product\ RUDDD\ B))))}
+
+{\tt (print\ (equal\ temp1\ temp2))}
+
+{\tt (print\ "bondi\ metric")}
+
+{\tt ;\ erase\ any\ definitions\ for\ U,\ V,\ beta\ and\ gamma}
+
+{\tt (setq\ U\ (quote\ U))}
+
+{\tt (setq\ V\ (quote\ V))}
+
+{\tt (setq\ beta\ (quote\ beta))}
+
+{\tt (setq\ gamma\ (quote\ gamma))}
+
+{\tt ;\ erase\ any\ definitions\ for\ u,\ r,\ theta\ and\ phi}
+
+{\tt (setq\ u\ (quote\ u))}
+
+{\tt (setq\ r\ (quote\ r))}
+
+{\tt (setq\ theta\ (quote\ theta))}
+
+{\tt (setq\ phi\ (quote\ phi))}
+
+{\tt ;\ new\ coordinate\ system}
+
+{\tt (setq\ x0\ u)}
+
+{\tt (setq\ x1\ r)}
+
+{\tt (setq\ x2\ theta)}
+
+{\tt (setq\ x3\ phi)}
+
+{\tt ;\ U,\ V,\ beta\ and\ gamma\ are\ functions\ of\ u,\ r\ and\ theta}
+
+$$g_{uu}=Ve^{2\beta}/r-U^2r^2e^{2\gamma}$$
+{\tt (setq\ g\_uu\ (sum}
+
+{\tt \ \ (product}
+
+{\tt \ \ \ \ (V\ u\ r\ theta)}
+
+{\tt \ \ \ \ (power\ r\ -1)}
+
+{\tt \ \ \ \ (power\ e\ (product\ 2\ (beta\ u\ r\ theta)))}
+
+{\tt \ \ )}
+
+{\tt \ \ (product}
+
+{\tt \ \ \ \ -1}
+
+{\tt \ \ \ \ (power\ (U\ u\ r\ theta)\ 2)}
+
+{\tt \ \ \ \ (power\ r\ 2)}
+
+{\tt \ \ \ \ (power\ e\ (product\ 2\ (gamma\ u\ r\ theta)))}
+
+{\tt \ \ )}
+
+{\tt ))}
+
+$$g_{ur}=2e^{2\beta}$$
+{\tt (setq\ g\_ur\ (product\ 2\ (power\ e\ (product\ 2\ (beta\ u\ r\ theta)))))}
+
+$$g_{u\theta}=2Ur^2e^{2\gamma}$$
+{\tt (setq\ g\_utheta\ (product}
+
+{\tt \ \ 2}
+
+{\tt \ \ (U\ u\ r\ theta)}
+
+{\tt \ \ (power\ r\ 2)}
+
+{\tt \ \ (power\ e\ (product\ 2\ (gamma\ u\ r\ theta)))}
+
+{\tt ))}
+
+$$g_{\theta\theta}=-r^2e^{2\gamma}$$
+{\tt (setq\ g\_thetatheta\ (product}
+
+{\tt \ \ -1}
+
+{\tt \ \ (power\ r\ 2)}
+
+{\tt \ \ (power\ e\ (product\ 2\ (gamma\ u\ r\ theta)))}
+
+{\tt ))}
+
+$$g_{\phi\phi}=-r^2e^{-2\gamma}\sin^2\theta$$
+{\tt (setq\ g\_phiphi\ (product}
+
+{\tt \ \ -1}
+
+{\tt \ \ (power\ r\ 2)}
+
+{\tt \ \ (power\ e\ (product\ -2\ (gamma\ u\ r\ theta)))}
+
+{\tt \ \ (power\ (sin\ theta)\ 2)}
+
+{\tt ))}
+
+{\tt ;\ metric\ tensor}
+
+{\tt (setq\ gdd\ (sum}
+
+{\tt \ \ (product\ g\_uu\ (tensor\ u\ u))}
+
+{\tt \ \ (product\ g\_ur\ (tensor\ u\ r))}
+
+{\tt \ \ (product\ g\_ur\ (tensor\ r\ u))}
+
+{\tt \ \ (product\ g\_utheta\ (tensor\ u\ theta))}
+
+{\tt \ \ (product\ g\_utheta\ (tensor\ theta\ u))}
+
+{\tt \ \ (product\ g\_thetatheta\ (tensor\ theta\ theta))}
+
+{\tt \ \ (product\ g\_phiphi\ (tensor\ phi\ phi))}
+
+{\tt ))}
+
+{\tt (gr)\ ;\ compute\ g,\ guu,\ GAMUDD,\ RUDDD,\ RDD,\ R,\ GDD,\ GUD\ and\ GUU}
+
+{\tt ;\ is\ covariant\ derivative\ of\ metric\ zero?}
+
+{\tt (print\ (equal\ (covariant-derivative-of-down-down\ gdd)\ 0))}
+
+{\tt ;\ is\ divergence\ of\ einstein\ zero?}
+
+{\tt (setq\ temp\ (contract23\ (covariant-derivative-of-up-up\ GUU)))}
+
+{\tt (print\ (equal\ temp\ 0))}
+
+\end

lisp/index.html

@@ -0,0 +1,39 @@
+Simple Lisp
+
+<pre>
+	[roscoe ~/lisp]$ gcc -o lisp lisp.c
+	[roscoe ~/lisp]$ ./lisp example.lisp
+	Bondi metric...
+	Is the Einstein tensor GUU divergence-free?
+	t
+	&gt; (length GUU)
+	280
+	&gt; (length (covariant-derivative-of-up-up GUU))
+	6723
+	&gt; ^C
+	[roscoe ~/lisp]$ ./lisp example2.lisp
+	connection coefficients (p. 210)
+	t
+	divergence of einstein is zero (p. 222)
+	t
+	computing riemann tensor (p. 219)
+	t
+	symmetry of covariant derivative (p. 252)
+	t
+	covariant derivative chain rule (p. 252)
+	t
+	additivity of covariant derivative (p. 257)
+	t
+	riemann is antisymmetric on last two indices (p. 286)
+	t
+	riemann vanishes when antisymmetrized on last three indices (p. 286)
+	t
+	double dual of riemann (p. 325)
+	t
+	noncommutation of covariant derivatives (p. 389)
+	t
+	bondi metric
+	t
+	t
+	&gt;
+</pre>

lisp/qed.lisp

@@ -0,0 +1,222 @@
+; From "Quantum Electrodynamics" by Richard P. Feynman
+; pp. 40-43
+; generic spacetime vectors a, b and c
+(setq a (sum
+  (product at (tensor t))
+  (product ax (tensor x))
+  (product ay (tensor y))
+  (product az (tensor z))
+))
+(setq b (sum
+  (product bt (tensor t))
+  (product bx (tensor x))
+  (product by (tensor y))
+  (product bz (tensor z))
+))
+(setq c (sum
+  (product ct (tensor t))
+  (product cx (tensor x))
+  (product cy (tensor y))
+  (product cz (tensor z))
+))
+; define this function for multiplying spactime vectors
+; how it works: (dot arg1 (tensor t)) picks off the t'th element, etc.
+; the -1's are for the spacetime metric
+(define spacetime-dot (sum
+  (dot (dot arg1 (tensor t)) (dot arg2 (tensor t)))
+  (dot -1 (dot arg1 (tensor x)) (dot arg2 (tensor x)))
+  (dot -1 (dot arg1 (tensor y)) (dot arg2 (tensor y)))
+  (dot -1 (dot arg1 (tensor z)) (dot arg2 (tensor z)))
+))
+(setq temp1 (spacetime-dot a a))
+(setq temp2 (sum
+  (power at 2)
+  (product -1 (power ax 2))
+  (product -1 (power ay 2))
+  (product -1 (power az 2))
+))
+(equal temp1 temp2) ; print "t" if it's true
+(setq I (sum
+  (tensor t t)
+  (tensor x x)
+  (tensor y y)
+  (tensor z z)
+))
+(setq gammat (sum
+  (tensor t t)
+  (tensor x x)
+  (product -1 (tensor y y))
+  (product -1 (tensor z z))
+))
+(setq gammax (sum
+  (tensor t z)
+  (tensor x y)
+  (product -1 (tensor y x))
+  (product -1 (tensor z t))
+))
+(setq gammay (sum
+  (product -1 i (tensor t z))
+  (product i (tensor x y))
+  (product i (tensor y x))
+  (product -1 i (tensor z t))
+))
+(setq gammaz (sum
+  (tensor t y)
+  (product -1 (tensor x z))
+  (product -1 (tensor y t))
+  (tensor z x)
+))
+(equal (dot gammat gammat) I) ; print "t" if it's true
+(equal (dot gammax gammax) (dot -1 I)) ; print "t" if it's true
+(equal (dot gammay gammay) (dot -1 I)) ; print "t" if it's true
+(equal (dot gammaz gammaz) (dot -1 I)) ; print "t" if it's true
+(setq gamma5 (dot gammax gammay gammaz gammat))
+(equal (dot gamma5 gamma5) (dot -1 I)) ; print "t" if it's true
+; gamma is a "vector" of dirac matrices
+(setq gamma (sum
+  (product gammat (tensor t))
+  (product gammax (tensor x))
+  (product gammay (tensor y))
+  (product gammaz (tensor z))
+))
+(setq agamma (spacetime-dot a gamma))
+(setq bgamma (spacetime-dot b gamma))
+(setq cgamma (spacetime-dot c gamma))
+(setq temp1 agamma)
+(setq temp2 (sum
+  (product at gammat)
+  (product -1 ax gammax)
+  (product -1 ay gammay)
+  (product -1 az gammaz)
+))
+(equal temp1 temp2) ; print "t" if it's true
+; note: gammas are square matrices, use "dot" to multiply
+; use "spacetime-dot" to multiply spacetime vectors
+(setq temp1 (dot agamma bgamma))
+(setq temp2 (sum
+  (dot -1 bgamma agamma)
+  (dot 2 (spacetime-dot a b) I)
+))
+(equal temp1 temp2) ; print "t" if it's true
+(setq temp1 (dot agamma gamma5))
+(setq temp2 (dot -1 gamma5 agamma))
+(equal temp1 temp2) ; print "t" if it's true
+(setq temp1 (dot gammax agamma gammax))
+(setq temp2 (sum agamma (dot 2 ax gammax)))
+(equal temp1 temp2) ; print "t" if it's true
+(setq temp1 (spacetime-dot gamma gamma))
+(setq temp2 (dot 4 I))
+(equal temp1 temp2) ; print "t" if it's true
+(setq temp1 (spacetime-dot gamma (dot agamma gamma)))
+(setq temp2 (dot -2 agamma))
+(equal temp1 temp2) ; print "t" if it's true
+(setq temp1 (spacetime-dot gamma (dot agamma bgamma gamma)))
+(setq temp2 (dot 4 (spacetime-dot a b) I))
+(equal temp1 temp2) ; print "t" if it's true
+(setq temp1 (spacetime-dot gamma (dot agamma bgamma cgamma gamma)))
+(setq temp2 (dot -2 cgamma bgamma agamma))
+(equal temp1 temp2) ; print "t" if it's true
+; define series approximations for some transcendental functions
+; for 32-bit integers, overflow occurs for powers above 5
+(define order 5)
+(define yexp (prog temp count
+  (setq temp 0)
+  (setq count order)
+loop
+  (setq temp (product (power count -1) arg (sum 1 temp)))
+  (setq count (sum count -1))
+  (cond ((greaterp count 0) (goto loop)))
+  (return (sum 1 temp))
+))
+(define ysin (sum
+  (product -1/2 i (yexp (product i arg)))
+  (product 1/2 i (yexp (product -1 i arg)))
+))
+(define ycos (sum
+  (product 1/2 (yexp (product i arg)))
+  (product 1/2 (yexp (product -1 i arg)))
+))
+(define ysinh (sum
+  (product 1/2 (yexp arg))
+  (product -1/2 (yexp (product -1 arg)))
+))
+(define ycosh (sum
+  (product 1/2 (yexp arg))
+  (product 1/2 (yexp (product -1 arg)))
+))
+; same as above but for matrices
+(define YEXP (prog temp count
+  (setq temp 0)
+  (setq count order)
+loop
+  (setq temp (dot (power count -1) arg (sum I temp)))
+  (setq count (sum count -1))
+  (cond ((greaterp count 0) (goto loop)))
+  (return (sum I temp))
+))
+(define YSIN (sum
+  (product -1/2 i (YEXP (product i arg)))
+  (product 1/2 i (YEXP (product -1 i arg)))
+))
+(define YCOS (sum
+  (product 1/2 (YEXP (product i arg)))
+  (product 1/2 (YEXP (product -1 i arg)))
+))
+(define YSINH (sum
+  (product 1/2 (YEXP arg))
+  (product -1/2 (YEXP (product -1 arg)))
+))
+(define YCOSH (sum
+  (product 1/2 (YEXP arg))
+  (product 1/2 (YEXP (product -1 arg)))
+))
+; for truncating products of power series
+(define POWER (cond
+  ((greaterp arg2 order) 0)
+  (t (list 'power arg1 arg2))
+))
+(define truncate (eval (subst 'POWER 'power arg)))
+(setq temp1 (YEXP (dot 1/2 u gammat gammax)))
+(setq temp2 (sum
+  (product I (ycosh (product 1/2 u))) ; could use "dot" but not necessary
+  (dot gammat gammax (ysinh (product 1/2 u)))
+))
+(equal temp1 temp2) ; print t if it's true
+(setq temp1 (YEXP (dot 1/2 theta gammax gammay)))
+(setq temp2 (sum
+  (product I (ycos (product 1/2 theta))) ; could use "dot" but not necessary
+  (dot gammax gammay (ysin (product 1/2 theta)))
+))
+(equal temp1 temp2) ; print t if it's true
+(setq temp1 (truncate (dot
+  (YEXP (dot -1/2 u gammat gammaz))
+  gammat
+  (YEXP (dot 1/2 u gammat gammaz))
+)))
+(setq temp2 (sum
+  (product gammat (ycosh u)) ; could use "dot" but not necessary
+  (product gammaz (ysinh u)) ; could use "dot" but not necessary
+))
+(equal temp1 temp2) ; print t if it's true
+(setq temp1 (truncate (dot
+  (YEXP (dot -1/2 u gammat gammaz))
+  gammaz
+  (YEXP (dot 1/2 u gammat gammaz))
+)))
+(setq temp2 (sum
+  (product gammaz (ycosh u)) ; could use "dot" but not necessary
+  (product gammat (ysinh u)) ; could use "dot" but not necessary
+))
+(equal temp1 temp2) ; print t if it's true
+(setq temp1 (truncate (dot
+  (YEXP (dot -1/2 u gammat gammaz))
+  gammay
+  (YEXP (dot 1/2 u gammat gammaz))
+)))
+(equal temp1 gammay) ; print t if it's true
+(setq temp1 (truncate (dot
+  (YEXP (dot -1/2 u gammat gammaz))
+  gammax
+  (YEXP (dot 1/2 u gammat gammaz))
+)))
+(equal temp1 gammax) ; print t if it's true

lisp/qed.tex

@@ -0,0 +1,535 @@
+\parindent=0pt
+{\tt ;\ From\ "Quantum\ Electrodynamics"\ by\ Richard\ P.\ Feynman}
+
+{\tt ;\ pp.\ 40-43}
+
+$$
+a=\left(\matrix{
+a_t\cr
+a_x\cr
+a_y\cr
+a_z\cr
+}\right)
+\quad
+b=\left(\matrix{
+b_t\cr
+b_x\cr
+b_y\cr
+b_z\cr
+}\right)
+\quad
+c=\left(\matrix{
+c_t\cr
+c_x\cr
+c_y\cr
+c_z\cr
+}\right)
+$$
+{\tt ;\ generic\ spacetime\ vectors\ a,\ b\ and\ c}
+
+{\tt (setq\ a\ (sum}
+
+{\tt \ \ (product\ at\ (tensor\ t))}
+
+{\tt \ \ (product\ ax\ (tensor\ x))}
+
+{\tt \ \ (product\ ay\ (tensor\ y))}
+
+{\tt \ \ (product\ az\ (tensor\ z))}
+
+{\tt ))}
+
+{\tt (setq\ b\ (sum}
+
+{\tt \ \ (product\ bt\ (tensor\ t))}
+
+{\tt \ \ (product\ bx\ (tensor\ x))}
+
+{\tt \ \ (product\ by\ (tensor\ y))}
+
+{\tt \ \ (product\ bz\ (tensor\ z))}
+
+{\tt ))}
+
+{\tt (setq\ c\ (sum}
+
+{\tt \ \ (product\ ct\ (tensor\ t))}
+
+{\tt \ \ (product\ cx\ (tensor\ x))}
+
+{\tt \ \ (product\ cy\ (tensor\ y))}
+
+{\tt \ \ (product\ cz\ (tensor\ z))}
+
+{\tt ))}
+
+{\tt ;\ define\ this\ function\ for\ multiplying\ spactime\ vectors}
+
+{\tt ;\ how\ it\ works:\ (dot\ arg1\ (tensor\ t))\ picks\ off\ the\ t'th\ element,\ etc.}
+
+{\tt ;\ the\ -1's\ are\ for\ the\ spacetime\ metric}
+
+{\tt (define\ spacetime-dot\ (sum}
+
+{\tt \ \ (dot\ (dot\ arg1\ (tensor\ t))\ (dot\ arg2\ (tensor\ t)))}
+
+{\tt \ \ (dot\ -1\ (dot\ arg1\ (tensor\ x))\ (dot\ arg2\ (tensor\ x)))}
+
+{\tt \ \ (dot\ -1\ (dot\ arg1\ (tensor\ y))\ (dot\ arg2\ (tensor\ y)))}
+
+{\tt \ \ (dot\ -1\ (dot\ arg1\ (tensor\ z))\ (dot\ arg2\ (tensor\ z)))}
+
+{\tt ))}
+
+$$a^2\mathrel{\mathop=^?}a_t^2-a_x^2-a_y^2-a_z^2$$
+{\tt (setq\ temp1\ (spacetime-dot\ a\ a))}
+
+{\tt (setq\ temp2\ (sum}
+
+{\tt \ \ (power\ at\ 2)}
+
+{\tt \ \ (product\ -1\ (power\ ax\ 2))}
+
+{\tt \ \ (product\ -1\ (power\ ay\ 2))}
+
+{\tt \ \ (product\ -1\ (power\ az\ 2))}
+
+{\tt ))}
+
+{\tt (equal\ temp1\ temp2)\ ;\ print\ "t"\ if\ it's\ true}
+
+$$I=\left(\matrix{
+1&0&0&0\cr
+0&1&0&0\cr
+0&0&1&0\cr
+0&0&0&1\cr
+}\right)$$
+{\tt (setq\ I\ (sum}
+
+{\tt \ \ (tensor\ t\ t)}
+
+{\tt \ \ (tensor\ x\ x)}
+
+{\tt \ \ (tensor\ y\ y)}
+
+{\tt \ \ (tensor\ z\ z)}
+
+{\tt ))}
+
+$$\gamma_t=\left(\matrix{
+1&0&0&0\cr
+0&1&0&0\cr
+0&0&-1&0\cr
+0&0&0&-1\cr
+}\right)$$
+{\tt (setq\ gammat\ (sum}
+
+{\tt \ \ (tensor\ t\ t)}
+
+{\tt \ \ (tensor\ x\ x)}
+
+{\tt \ \ (product\ -1\ (tensor\ y\ y))}
+
+{\tt \ \ (product\ -1\ (tensor\ z\ z))}
+
+{\tt ))}
+
+$$\gamma_x=\left(\matrix{
+0&0&0&1\cr
+0&0&1&0\cr
+0&-1&0&0\cr
+-1&0&0&0\cr
+}\right)$$
+{\tt (setq\ gammax\ (sum}
+
+{\tt \ \ (tensor\ t\ z)}
+
+{\tt \ \ (tensor\ x\ y)}
+
+{\tt \ \ (product\ -1\ (tensor\ y\ x))}
+
+{\tt \ \ (product\ -1\ (tensor\ z\ t))}
+
+{\tt ))}
+
+$$\gamma_y=\left(\matrix{
+0&0&0&-i\cr
+0&0&i&0\cr
+0&i&0&0\cr
+-i&0&0&0\cr
+}\right)$$
+{\tt (setq\ gammay\ (sum}
+
+{\tt \ \ (product\ -1\ i\ (tensor\ t\ z))}
+
+{\tt \ \ (product\ i\ (tensor\ x\ y))}
+
+{\tt \ \ (product\ i\ (tensor\ y\ x))}
+
+{\tt \ \ (product\ -1\ i\ (tensor\ z\ t))}
+
+{\tt ))}
+
+$$\gamma_z=\left(\matrix{
+0&0&1&0\cr
+0&0&0&-1\cr
+-1&0&0&0\cr
+0&1&0&0\cr
+}\right)$$
+{\tt (setq\ gammaz\ (sum}
+
+{\tt \ \ (tensor\ t\ y)}
+
+{\tt \ \ (product\ -1\ (tensor\ x\ z))}
+
+{\tt \ \ (product\ -1\ (tensor\ y\ t))}
+
+{\tt \ \ (tensor\ z\ x)}
+
+{\tt ))}
+
+$$\gamma_t^2\mathrel{\mathop=^?}1$$
+{\tt (equal\ (dot\ gammat\ gammat)\ I)\ ;\ print\ "t"\ if\ it's\ true}
+
+$$\gamma_x^2=\gamma_y^2=\gamma_z^2\mathrel{\mathop=^?}-1$$
+{\tt (equal\ (dot\ gammax\ gammax)\ (dot\ -1\ I))\ ;\ print\ "t"\ if\ it's\ true}
+
+{\tt (equal\ (dot\ gammay\ gammay)\ (dot\ -1\ I))\ ;\ print\ "t"\ if\ it's\ true}
+
+{\tt (equal\ (dot\ gammaz\ gammaz)\ (dot\ -1\ I))\ ;\ print\ "t"\ if\ it's\ true}
+
+$$\gamma_5=\gamma_x\gamma_y\gamma_z\gamma_t$$
+{\tt (setq\ gamma5\ (dot\ gammax\ gammay\ gammaz\ gammat))}
+
+$$\gamma_5^2\mathrel{\mathop=^?}-1$$
+{\tt (equal\ (dot\ gamma5\ gamma5)\ (dot\ -1\ I))\ ;\ print\ "t"\ if\ it's\ true}
+
+$$\gamma=\left(\matrix{
+\gamma_t\cr
+\gamma_x\cr
+\gamma_y\cr
+\gamma_z\cr
+}\right)$$
+{\tt ;\ gamma\ is\ a\ "vector"\ of\ dirac\ matrices}
+
+{\tt (setq\ gamma\ (sum}
+
+{\tt \ \ (product\ gammat\ (tensor\ t))}
+
+{\tt \ \ (product\ gammax\ (tensor\ x))}
+
+{\tt \ \ (product\ gammay\ (tensor\ y))}
+
+{\tt \ \ (product\ gammaz\ (tensor\ z))}
+
+{\tt ))}
+
+$$a\!\!\!/=a\gamma\qquad b\!\!\!/=b\gamma\qquad c\!\!\!/=c\gamma$$
+{\tt (setq\ agamma\ (spacetime-dot\ a\ gamma))}
+
+{\tt (setq\ bgamma\ (spacetime-dot\ b\ gamma))}
+
+{\tt (setq\ cgamma\ (spacetime-dot\ c\ gamma))}
+
+$$a\!\!\!/\mathrel{\mathop=^?}a_t\gamma_t-a_x\gamma_x-a_y\gamma_y-a_x\gamma_x$$
+{\tt (setq\ temp1\ agamma)}
+
+{\tt (setq\ temp2\ (sum}
+
+{\tt \ \ (product\ at\ gammat)}
+
+{\tt \ \ (product\ -1\ ax\ gammax)}
+
+{\tt \ \ (product\ -1\ ay\ gammay)}
+
+{\tt \ \ (product\ -1\ az\ gammaz)}
+
+{\tt ))}
+
+{\tt (equal\ temp1\ temp2)\ ;\ print\ "t"\ if\ it's\ true}
+
+$$a\!\!\!/b\!\!\!/\mathrel{\mathop=^?}-b\!\!\!/a\!\!\!/ + 2ab$$
+{\tt ;\ note:\ gammas\ are\ square\ matrices,\ use\ "dot"\ to\ multiply}
+
+{\tt ;\ use\ "spacetime-dot"\ to\ multiply\ spacetime\ vectors}
+
+{\tt (setq\ temp1\ (dot\ agamma\ bgamma))}
+
+{\tt (setq\ temp2\ (sum}
+
+{\tt \ \ (dot\ -1\ bgamma\ agamma)}
+
+{\tt \ \ (dot\ 2\ (spacetime-dot\ a\ b)\ I)}
+
+{\tt ))}
+
+{\tt (equal\ temp1\ temp2)\ ;\ print\ "t"\ if\ it's\ true}
+
+$$a\!\!\!/\gamma_5\mathrel{\mathop=^?}-\gamma_5a\!\!\!/$$
+{\tt (setq\ temp1\ (dot\ agamma\ gamma5))}
+
+{\tt (setq\ temp2\ (dot\ -1\ gamma5\ agamma))}
+
+{\tt (equal\ temp1\ temp2)\ ;\ print\ "t"\ if\ it's\ true}
+
+$$\gamma_x a\!\!\!/\gamma_x\mathrel{\mathop=^?}a\!\!\!/+2a_x\gamma_x$$
+{\tt (setq\ temp1\ (dot\ gammax\ agamma\ gammax))}
+
+{\tt (setq\ temp2\ (sum\ agamma\ (dot\ 2\ ax\ gammax)))}
+
+{\tt (equal\ temp1\ temp2)\ ;\ print\ "t"\ if\ it's\ true}
+
+$$\gamma\gamma\mathrel{\mathop=^?}4$$
+{\tt (setq\ temp1\ (spacetime-dot\ gamma\ gamma))}
+
+{\tt (setq\ temp2\ (dot\ 4\ I))}
+
+{\tt (equal\ temp1\ temp2)\ ;\ print\ "t"\ if\ it's\ true}
+
+$$\gamma a\!\!\!/\gamma\mathrel{\mathop=^?}-2a\!\!\!/$$
+{\tt (setq\ temp1\ (spacetime-dot\ gamma\ (dot\ agamma\ gamma)))}
+
+{\tt (setq\ temp2\ (dot\ -2\ agamma))}
+
+{\tt (equal\ temp1\ temp2)\ ;\ print\ "t"\ if\ it's\ true}
+
+$$\gamma a\!\!\!/b\!\!\!/\gamma\mathrel{\mathop=^?}4ab$$
+{\tt (setq\ temp1\ (spacetime-dot\ gamma\ (dot\ agamma\ bgamma\ gamma)))}
+
+{\tt (setq\ temp2\ (dot\ 4\ (spacetime-dot\ a\ b)\ I))}
+
+{\tt (equal\ temp1\ temp2)\ ;\ print\ "t"\ if\ it's\ true}
+
+$$\gamma a\!\!\!/b\!\!\!/c\!\!\!/\gamma\mathrel{\mathop=^?}
+-2c\!\!\!/b\!\!\!/a\!\!\!/$$
+{\tt (setq\ temp1\ (spacetime-dot\ gamma\ (dot\ agamma\ bgamma\ cgamma\ gamma)))}
+
+{\tt (setq\ temp2\ (dot\ -2\ cgamma\ bgamma\ agamma))}
+
+{\tt (equal\ temp1\ temp2)\ ;\ print\ "t"\ if\ it's\ true}
+
+{\tt ;\ define\ series\ approximations\ for\ some\ transcendental\ functions}
+
+{\tt ;\ for\ 32-bit\ integers,\ overflow\ occurs\ for\ powers\ above\ 5}
+
+{\tt (define\ order\ 5)}
+
+{\tt (define\ yexp\ (prog\ temp\ count}
+
+{\tt \ \ (setq\ temp\ 0)}
+
+{\tt \ \ (setq\ count\ order)}
+
+{\tt loop}
+
+{\tt \ \ (setq\ temp\ (product\ (power\ count\ -1)\ arg\ (sum\ 1\ temp)))}
+
+{\tt \ \ (setq\ count\ (sum\ count\ -1))}
+
+{\tt \ \ (cond\ ((greaterp\ count\ 0)\ (goto\ loop)))}
+
+{\tt \ \ (return\ (sum\ 1\ temp))}
+
+{\tt ))}
+
+{\tt (define\ ysin\ (sum}
+
+{\tt \ \ (product\ -1/2\ i\ (yexp\ (product\ i\ arg)))}
+
+{\tt \ \ (product\ 1/2\ i\ (yexp\ (product\ -1\ i\ arg)))}
+
+{\tt ))}
+
+{\tt (define\ ycos\ (sum}
+
+{\tt \ \ (product\ 1/2\ (yexp\ (product\ i\ arg)))}
+
+{\tt \ \ (product\ 1/2\ (yexp\ (product\ -1\ i\ arg)))}
+
+{\tt ))}
+
+{\tt (define\ ysinh\ (sum}
+
+{\tt \ \ (product\ 1/2\ (yexp\ arg))}
+
+{\tt \ \ (product\ -1/2\ (yexp\ (product\ -1\ arg)))}
+
+{\tt ))}
+
+{\tt (define\ ycosh\ (sum}
+
+{\tt \ \ (product\ 1/2\ (yexp\ arg))}
+
+{\tt \ \ (product\ 1/2\ (yexp\ (product\ -1\ arg)))}
+
+{\tt ))}
+
+{\tt ;\ same\ as\ above\ but\ for\ matrices}
+
+{\tt (define\ YEXP\ (prog\ temp\ count}
+
+{\tt \ \ (setq\ temp\ 0)}
+
+{\tt \ \ (setq\ count\ order)}
+
+{\tt loop}
+
+{\tt \ \ (setq\ temp\ (dot\ (power\ count\ -1)\ arg\ (sum\ I\ temp)))}
+
+{\tt \ \ (setq\ count\ (sum\ count\ -1))}
+
+{\tt \ \ (cond\ ((greaterp\ count\ 0)\ (goto\ loop)))}
+
+{\tt \ \ (return\ (sum\ I\ temp))}
+
+{\tt ))}
+
+{\tt (define\ YSIN\ (sum}
+
+{\tt \ \ (product\ -1/2\ i\ (YEXP\ (product\ i\ arg)))}
+
+{\tt \ \ (product\ 1/2\ i\ (YEXP\ (product\ -1\ i\ arg)))}
+
+{\tt ))}
+
+{\tt (define\ YCOS\ (sum}
+
+{\tt \ \ (product\ 1/2\ (YEXP\ (product\ i\ arg)))}
+
+{\tt \ \ (product\ 1/2\ (YEXP\ (product\ -1\ i\ arg)))}
+
+{\tt ))}
+
+{\tt (define\ YSINH\ (sum}
+
+{\tt \ \ (product\ 1/2\ (YEXP\ arg))}
+
+{\tt \ \ (product\ -1/2\ (YEXP\ (product\ -1\ arg)))}
+
+{\tt ))}
+
+{\tt (define\ YCOSH\ (sum}
+
+{\tt \ \ (product\ 1/2\ (YEXP\ arg))}
+
+{\tt \ \ (product\ 1/2\ (YEXP\ (product\ -1\ arg)))}
+
+{\tt ))}
+
+{\tt ;\ for\ truncating\ products\ of\ power\ series}
+
+{\tt (define\ POWER\ (cond}
+
+{\tt \ \ ((greaterp\ arg2\ order)\ 0)}
+
+{\tt \ \ (t\ (list\ 'power\ arg1\ arg2))}
+
+{\tt ))}
+
+{\tt (define\ truncate\ (eval\ (subst\ 'POWER\ 'power\ arg)))}
+
+$$\exp[(u/2)\gamma_t\gamma_x]\mathrel{\mathop=^?}
+\cosh(u/2)+\gamma_t\gamma_x\sinh(u/2)$$
+{\tt (setq\ temp1\ (YEXP\ (dot\ 1/2\ u\ gammat\ gammax)))}
+
+{\tt (setq\ temp2\ (sum}
+
+{\tt \ \ (product\ I\ (ycosh\ (product\ 1/2\ u)))\ ;\ could\ use\ "dot"\ but\ not\ necessary}
+
+{\tt \ \ (dot\ gammat\ gammax\ (ysinh\ (product\ 1/2\ u)))}
+
+{\tt ))}
+
+{\tt (equal\ temp1\ temp2)\ ;\ print\ t\ if\ it's\ true}
+
+$$\exp[(\theta/2)\gamma_x\gamma_y]\mathrel{\mathop=^?}
+\cos(\theta/2)+\gamma_x\gamma_y\sin(\theta/2)$$
+{\tt (setq\ temp1\ (YEXP\ (dot\ 1/2\ theta\ gammax\ gammay)))}
+
+{\tt (setq\ temp2\ (sum}
+
+{\tt \ \ (product\ I\ (ycos\ (product\ 1/2\ theta)))\ ;\ could\ use\ "dot"\ but\ not\ necessary}
+
+{\tt \ \ (dot\ gammax\ gammay\ (ysin\ (product\ 1/2\ theta)))}
+
+{\tt ))}
+
+{\tt (equal\ temp1\ temp2)\ ;\ print\ t\ if\ it's\ true}
+
+$$\exp[-(u/2)\gamma_t\gamma_z]\gamma_t
+\exp[(u/2)\gamma_t\gamma_z]
+\mathrel{\mathop=^?}\gamma_t\cosh u+\gamma_z\sinh u$$
+{\tt (setq\ temp1\ (truncate\ (dot}
+
+{\tt \ \ (YEXP\ (dot\ -1/2\ u\ gammat\ gammaz))}
+
+{\tt \ \ gammat}
+
+{\tt \ \ (YEXP\ (dot\ 1/2\ u\ gammat\ gammaz))}
+
+{\tt )))}
+
+{\tt (setq\ temp2\ (sum}
+
+{\tt \ \ (product\ gammat\ (ycosh\ u))\ ;\ could\ use\ "dot"\ but\ not\ necessary}
+
+{\tt \ \ (product\ gammaz\ (ysinh\ u))\ ;\ could\ use\ "dot"\ but\ not\ necessary}
+
+{\tt ))}
+
+{\tt (equal\ temp1\ temp2)\ ;\ print\ t\ if\ it's\ true}
+
+$$\exp[-(u/2)\gamma_t\gamma_z]\gamma_z
+\exp[(u/2)\gamma_t\gamma_z]
+\mathrel{\mathop=^?}\gamma_z\cosh u+\gamma_t\sinh u$$
+{\tt (setq\ temp1\ (truncate\ (dot}
+
+{\tt \ \ (YEXP\ (dot\ -1/2\ u\ gammat\ gammaz))}
+
+{\tt \ \ gammaz}
+
+{\tt \ \ (YEXP\ (dot\ 1/2\ u\ gammat\ gammaz))}
+
+{\tt )))}
+
+{\tt (setq\ temp2\ (sum}
+
+{\tt \ \ (product\ gammaz\ (ycosh\ u))\ ;\ could\ use\ "dot"\ but\ not\ necessary}
+
+{\tt \ \ (product\ gammat\ (ysinh\ u))\ ;\ could\ use\ "dot"\ but\ not\ necessary}
+
+{\tt ))}
+
+{\tt (equal\ temp1\ temp2)\ ;\ print\ t\ if\ it's\ true}
+
+$$\exp[-(u/2)\gamma_t\gamma_z]\gamma_y
+\exp[(u/2)\gamma_t\gamma_z]
+\mathrel{\mathop=^?}\gamma_y$$
+{\tt (setq\ temp1\ (truncate\ (dot}
+
+{\tt \ \ (YEXP\ (dot\ -1/2\ u\ gammat\ gammaz))}
+
+{\tt \ \ gammay}
+
+{\tt \ \ (YEXP\ (dot\ 1/2\ u\ gammat\ gammaz))}
+
+{\tt )))}
+
+{\tt (equal\ temp1\ gammay)\ ;\ print\ t\ if\ it's\ true}
+
+$$\exp[-(u/2)\gamma_t\gamma_z]\gamma_x
+\exp[(u/2)\gamma_t\gamma_z]
+\mathrel{\mathop=^?}\gamma_x$$
+{\tt (setq\ temp1\ (truncate\ (dot}
+
+{\tt \ \ (YEXP\ (dot\ -1/2\ u\ gammat\ gammaz))}
+
+{\tt \ \ gammax}
+
+{\tt \ \ (YEXP\ (dot\ 1/2\ u\ gammat\ gammaz))}
+
+{\tt )))}
+
+{\tt (equal\ temp1\ gammax)\ ;\ print\ t\ if\ it's\ true}
+
+\end

lisp/reference.html

@@ -0,0 +1,627 @@
+<html>
+
+<p>
+<h1>adjunct</h1>
+
+<p>
+A adj A = det A I
+
+<pre>
+     > (setq A (sum
+     (product a (tensor x x))
+     (product b (tensor x y))
+     (product c (tensor y x))
+     (product d (tensor y y))
+     ))
+     > (setq A1 (adjunct A x y))
+     > (setq AA1 (contract23 (product A A1)))
+     > (setq I (sum (tensor x x) (tensor y y)))
+     > (equal AA1 (product (determinant A x y) I))
+     t
+</pre>
+
+<p>
+<h1>and</h1>
+
+<pre>
+     > (and t t t)
+     t
+     > (and t nil t)
+     nil
+</pre>
+
+<p>
+<h1>append</h1>
+
+<pre>
+     > (append (list a b) (list c d))
+     (a b c d)
+</pre>
+
+<p>
+<h1>atom</h1>
+
+<pre>
+     > (atom a)
+     t
+     > (atom (list a b c))
+     nil
+</pre>
+
+<p>
+<h1>car</h1>
+
+<pre>
+     > (car (list a b c))
+     a
+</pre>
+
+<p>
+<h1>cdr</h1>
+
+<pre>
+     > (cdr (list a b c))
+     (b c)
+</pre>
+
+<p>
+<h1>complex-conjugate</h1>
+
+<pre>
+     > (complex-conjugate (sum a (product b i)))
+     (sum a (product -1 b i))
+</pre>
+
+<p>
+<h1>component</h1>
+
+<pre>
+     > (setq A (sum
+     (product a (tensor x x))
+     (product b (tensor x y))
+     (product c (tensor y x))
+     (product d (tensor y y))
+     ))
+     > (component A x y)
+     b
+</pre>
+
+<p>
+<h1>cond</h1>
+
+<p>
+(cond (a1 b1) (a2 b2)... ) returns b of first non-nil a
+
+<pre>
+     > (define foo (cond ((zerop arg) "zero") (t "not zero")))
+     > (foo 0)
+     zero
+     > (foo 1)
+     not zero
+</pre>
+
+<p>
+<h1>cons</h1>
+
+<pre>
+     > (cons a b)
+     (a . b)
+</pre>
+
+<p>
+<h1>contract</h1>
+
+<p>
+Form: contract12, contract13, contract23, etc.
+For example (contract12 x) contracts, or sums over,
+the first and second indices.
+
+<pre>
+     > (setq A (sum
+     (product a (tensor x x))
+     (product b (tensor x y))
+     (product c (tensor y x))
+     (product d (tensor y y))
+     ))
+     > (contract12 A)
+     (sum a d)
+</pre>
+
+<p>
+<h1>cos</h1>
+
+<pre>
+     > (derivative (cos x) x)
+     (product -1 (sin x))
+</pre>
+
+<p>
+<h1>COS</h1>
+
+<p>
+This function returns the exponential form of the cosine function.
+
+<pre>
+     > (printcars (COS x))
+     sum
+     (product 1/2 (power e (product -1 i x)))
+     (product 1/2 (power e (product i x)))
+</pre>
+
+<p>
+<h1>define</h1>
+
+<pre>
+     > (define foo (cons arg2 arg1))
+     > (foo a b)
+     (b . a)
+     > foo
+     (cons arg2 arg1)
+</pre>
+
+<p>
+<h1>derivative</h1>
+
+<pre>
+     > (derivative (power x 2) x)
+     (product 2 x)
+</pre>
+
+<p>
+<h1>determinant</h1>
+
+<p>
+see adjunct
+
+<p>
+<h1>dot</h1>
+
+<p>
+The dot function returns the inner product of scalars and tensors.
+It is equivalent to an outer product followed by a contraction on
+the inner indices.
+
+<pre>
+     > (setq A (sum (product A11 (tensor 1 1)) (product A22 (tensor 2 2))))
+     > (setq X (sum (product X1 (tensor 1)) (product X2 (tensor 2))))
+     > (setq AX1 (dot A X))
+     > (setq AX2 (contract23 (product A X)))
+     > (equal AX1 AX2)
+     t
+     > (printcars AX1)
+     sum
+     (product A11 X1 (tensor 1))
+     (product A22 X2 (tensor 2))
+</pre>
+
+<p>
+The arguments are multiplied left to right,
+i.e. (dot a b c) = (dot (dot a b) c).
+
+<p>
+<h1>e</h1>
+
+<pre>
+     > (derivative (power e x) x)
+     (power e x)
+</pre>
+
+<p>
+<h1>equal</h1>
+
+<pre>
+     > (equal a a)
+     t
+</pre>
+
+<p>
+<h1>eval</h1>
+
+<pre>
+     > (setq foo (quote (sum a a)))
+     > foo
+     (sum a a)
+     > (eval foo)
+     (product 2 a)
+</pre>
+
+<p>
+<h1>exit</h1>
+
+<p>
+Exit simple lisp and return to the shell.
+Control-C also does this.
+
+<pre>
+     > (exit)
+     %
+</pre>
+
+<p>
+<h1>exp</h1>
+
+<pre>
+     > (exp x)
+     (power e x)
+</pre>
+
+<p>
+<h1>gc</h1>
+
+<p>
+Normally, garbage collection occurs when simple lisp runs
+out of memory.
+The gc function causes garbage collection to occur immediately.
+It returns the number of free cells, 1 cell = 12 bytes.
+
+<pre>
+     > (gc)
+     496053
+</pre>
+
+<p>
+<h1>goto</h1>
+
+<pre>
+     > (define foo (prog (k)
+     (setq k 1)
+     loop
+     (print k)
+     (setq k (sum k 1))
+     (cond ((lessp k 4) (goto loop)))
+     ))
+     > (foo)
+     1
+     2
+     3
+</pre>
+
+<p>
+<h1>greaterp</h1>
+
+<pre>
+     > (greaterp 1 0)
+     t
+     > (greaterp 0 1)
+     nil
+</pre>
+
+<p>
+<h1>i</h1>
+
+<p>
+The symbol "i" represents the square root of minus one.
+
+<pre>
+     > (product i i)
+     -1
+</pre>
+
+<p>
+<h1>integerp</h1>
+
+<p>
+Returns t if the argument is an integer, nil otherwise.
+
+<pre>
+     > (integerp 2)
+     t
+     > (integerp 2/3)
+     nil
+</pre>
+
+<p>
+<h1>integral</h1>
+
+<pre>
+     > (integral (power x 2) x)
+     (product 1/3 (power x 3))
+</pre>
+
+<p>
+The integral function uses a list of expression forms,
+much like a table of integrals.
+To extend the capability of the integral function,
+refer to the file "startup" and follow suit with "add-integral".
+
+<p>
+<h1>laguerre</h1>
+
+<p>
+See the "hydrogen atom wave functions" example.
+
+<p>
+<h1>laplacian</h1>
+
+<p>
+See the "hydrogen atom wave functions" example.
+
+<p>
+<h1>length</h1>
+
+<p>
+Returns the number of cars in a list.
+
+<pre>
+     > (length (list a b c))
+     3
+</pre>
+
+<p>
+<h1>lessp</h1>
+
+<pre>
+     > (lessp 0 1)
+     t
+     > (lessp 1 0)
+     nil
+</pre>
+
+<p>
+<h1>list</h1>
+
+<pre>
+     > (list a b c)
+     (a b c)
+</pre>
+
+<p>
+<h1>not</h1>
+
+<pre>
+     > (not nil)
+     t
+</pre>
+
+<p>
+<h1>null</h1>
+
+<pre>
+     > (null nil)
+     t
+</pre>
+
+<p>
+<h1>numberp</h1>
+
+<pre>
+     > (numberp 1)
+     t
+</pre>
+
+<p>
+<h1>or</h1>
+
+<pre>
+     > (or nil t nil)
+     t
+     > (or nil nil nil)
+     nil
+</pre>
+
+<p>
+<h1>power</h1>
+
+<pre>
+     > (power a (sum b c))
+     (product (power a b) (power a c))
+</pre>
+
+<p>
+<h1>print</h1>
+
+<pre>
+     > (print "hello")
+     hello
+</pre>
+
+<p>
+<h1>printcars</h1>
+
+<pre>
+     > (printcars (list a b c))
+     a
+     b
+     c
+</pre>
+
+<p>
+<h1>product</h1>
+
+<pre>
+     > (product a b a)
+     (product b (power a 2))
+</pre>
+
+<p>
+The "product" function expands products of sums.
+When applied to tensors, the result is an outer product.
+
+<pre>
+     > (setq A (sum
+     (product a11 (tensor 1 1))
+     (product a12 (tensor 1 2))
+     (product a21 (tensor 2 1))
+     (product a22 (tensor 2 2))
+     ))
+     > (setq X (sum
+     (product x1 (tensor 1))
+     (product x2 (tensor 2))
+     ))
+     > (printcars (product A X))
+     sum
+     (product a11 x1 (tensor 1 1 1))
+     (product a11 x2 (tensor 1 1 2))
+     (product a12 x1 (tensor 1 2 1))
+     (product a12 x2 (tensor 1 2 2))
+     (product a21 x1 (tensor 2 1 1))
+     (product a21 x2 (tensor 2 1 2))
+     (product a22 x1 (tensor 2 2 1))
+     (product a22 x2 (tensor 2 2 2))
+</pre>
+
+<p>
+The more familiar inner product is obtained by contraction.
+
+<pre>
+     > (printcars (contract23 (product A X)))
+     sum
+     (product a11 x1 (tensor 1))
+     (product a12 x2 (tensor 1))
+     (product a21 x1 (tensor 2))
+     (product a22 x2 (tensor 2))
+</pre>
+
+<p>
+<h1>prog</h1>
+
+<p>
+Use "prog" when more than one expression must be evaluated in a function.
+
+<pre>
+     > (define foo (prog ()
+     (print arg1)
+     (print arg2)
+     ))
+     > (foo a b)
+     a
+     b
+</pre>
+
+<p>
+<h1>quote</h1>
+
+<p>
+Quote means "don't evaluate."
+
+<pre>
+     > (quote (sum a a))
+     (sum a a)
+</pre>
+
+<p>
+An apostrophe does the same thing.
+
+<pre>
+     > '(sum a a)
+     (sum a a)
+</pre>
+
+<p>
+<h1>return</h1>
+
+<pre>
+     > (define foo (prog () (return "hello")))
+     > (foo)
+     hello
+</pre>
+
+<p>
+<h1>run</h1>
+
+<p>
+The "run" function evaluates the expressions in a file.
+
+<pre>
+     > (run "example1")
+     E=K+V for psi5?
+     t
+</pre>
+
+<p>
+<h1>setq</h1>
+
+<pre>
+     > (setq foo 1)
+     > foo
+     1
+</pre>
+
+<p>
+<h1>sin</h1>
+
+<pre>
+     > (derivative (sin x) x)
+     (cos x)
+</pre>
+
+<p>
+<h1>SIN</h1>
+
+<p>
+The SIN function returns the exponential form of the sine function.
+In some cases it is better to use SIN and COS because many
+trigonometric simplifications occur automatically.
+
+<pre>
+     > (printcars (SIN x))
+     sum
+     (product -1/2 i (power e (product i x)))
+     (product 1/2 i (power e (product -1 i x)))
+</pre>
+
+<pre>
+     > (sum (power (COS x) 2) (power (SIN x) 2))
+     1
+</pre>
+
+<p>
+<h1>sum</h1>
+
+<pre>
+     > (sum a b a)
+     (sum b (product 2 a))
+</pre>
+
+<p>
+<h1>subst</h1>
+
+<p>
+(subst a b c) means substitute a for b wherever it appears in c.
+
+<pre>
+     > (subst a b (list a b c))
+     (a a c)
+</pre>
+
+<p>
+<h1>tensor</h1>
+
+<pre>
+     > (product (tensor x) (tensor y))
+     (tensor x y)
+     > (product (tensor y) (tensor x))
+     (tensor y x)
+</pre>
+
+<p>
+<h1>transpose</h1>
+
+<p>
+Form: transpose12, transpose13, transpose23, etc.
+For example (transpose12 x) transposes the first and second indices.
+
+<pre>
+     > (setq A (sum
+     (product a (tensor x x))
+     (product b (tensor x y))
+     (product c (tensor y x))
+     (product d (tensor y y))
+     ))
+     > (printcars (transpose12 A))
+     sum
+     (product a (tensor x x))
+     (product b (tensor y x))
+     (product c (tensor x y))
+     (product d (tensor y y))
+</pre>
+
+<p>
+<h1>zerop</h1>
+
+<pre>
+     > (zerop 0)
+     t
+     > (zerop 1)
+     nil
+</pre>

lisp/robbins.lisp

@@ -0,0 +1,303 @@
+; December 7, 2001
+
+; See the paper "Solutions of the Robbins Problem" by William McCune.
+
+; Robbins equation (7)
+
+(setq E7 (eq (n (sum (n (sum (n x) y)) (n (sum x y)))) y))
+
+; From this derive the following:
+
+; n(n(3x)+x)+2x = 2x
+
+; STEP 10 ---------------------------------------------------------------------
+
+; with (7), let x be n(x)+y and y be n(x+y)
+
+(setq E10 (subst Y y E7))
+
+(setq E10 (subst (sum (n x) y) x E10))
+
+(setq E10 (eval (subst (n (sum x y)) Y E10)))
+
+; simplify with (7)
+
+(setq E10 (eval (subst (caddr E7) (cadr E7) E10)))
+
+; n(n(n(x+y)+n(x)+y)+y) = n(x+y) ?
+
+(setq tmp (eq (n (sum (n (sum (n (sum x y)) (n x) y)) y)) (n (sum x y))))
+
+(equal E10 tmp)
+
+; STEP 11 ---------------------------------------------------------------------
+
+; with (7), let y be n(n(x)+y) and x be x+y
+
+(setq E11 (subst X x E7))
+
+(setq E11 (subst (n (sum (n x) y)) y E11))
+
+(setq E11 (eval (subst (sum x y) X E11)))
+
+; simplify with (7)
+
+(setq E11 (eval (subst (caddr E7) (cadr E7) E11)))
+
+; n(n(n(n(x)+y)+x+y)+y) = n(n(x)+y) ?
+
+(setq tmp (eq (n (sum (n (sum (n (sum (n x) y)) x y)) y)) (n (sum (n x) y))))
+
+(equal E11 tmp)
+
+; STEP 29 ---------------------------------------------------------------------
+
+; with (7), let x be n(n(x)+y)+x+y and y be y
+
+(setq E29 (eval (subst (sum (n (sum (n x) y)) x y) x E7)))
+
+; simplify with (11)
+
+(setq E29 (eval (subst (caddr E11) (cadr E11) E29)))
+
+; n(n(n(n(x)+y)+x+2y)+n(n(x)+y)) = y ?
+
+(setq tmp (eq (n (sum (n (sum (n (sum (n x) y)) x y y)) (n (sum (n x) y)))) y))
+
+(equal E29 tmp)
+
+; STEP 54 ---------------------------------------------------------------------
+
+; with (7), let x be n(n(n(x)+y)+x+2y)+n(n(x)+y) and y be z
+
+(setq E54 (subst z y E7))
+
+(setq tmp (sum (n (sum x (n (sum y (n x))) (product 2 y))) (n (sum y (n x)))))
+
+(setq E54 (eval (subst tmp x E54)))
+
+; simplify with (29)
+
+(setq E54 (eval (subst (caddr E29) (cadr E29) E54)))
+
+; n(n(n(n(n(x)+y)+x+2y)+n(n(x)+y)+z)+n(y+z)) = z ?
+
+(setq tmp (eq (n (sum (n (sum y z)) (n (sum z (n (sum x (n (sum y (n x))) (product 2 y))) (n (sum y (n x))))))) z))
+
+(equal E54 tmp)
+
+; STEP 217 --------------------------------------------------------------------
+
+; with (7), let x be n(n(n(x)+y)+x+2y)+n(n(x)+y)+z and y be n(y+z)
+
+(setq E217 (subst Y y E7))
+
+(setq tmp (sum z (n (sum x (n (sum y (n x))) (product 2 y))) (n (sum y (n x)))))
+
+(setq E217 (subst tmp x E217))
+
+(setq E217 (eval (subst (n (sum y z)) Y E217)))
+
+; simplify with E54
+
+(setq E217 (eval (subst (caddr E54) (cadr E54) E217)))
+
+; n(n(n(n(n(x)+y)+x+2y)+n(n(x)+y)+n(y+z)+z)+z) = n(y+z) ?
+
+(setq tmp (eq (n (sum z (n (sum z (n (sum x (n (sum y (n x))) (product 2 y))) (n (sum y z)) (n (sum y (n x))))))) (n (sum y z))))
+
+(equal E217 tmp)
+
+; STEP 674 --------------------------------------------------------------------
+
+; with (7), let y be u and x be n(n(n(n(x)+y)+x+2y)+n(n(x)+y)+n(y+z)+z)+z
+
+(setq E674 (subst u y E7))
+
+(setq tmp (sum z (n (sum z (n (sum x (n (sum y (n x))) (product 2 y))) (n (sum y z)) (n (sum y (n x)))))))
+
+(setq E674 (eval (subst tmp x E674)))
+
+; simplify with E217
+
+(setq E674 (eval (subst (caddr E217) (cadr E217) E674)))
+
+; n(n(n(n(n(n(x)+y)+x+2y)+n(n(x)+y)+n(y+z)+z)+z+u)+n(n(y+z)+u)) = u ?
+
+(setq tmp (eq (n (sum (n (sum u z (n (sum z (n (sum x (n (sum y (n x))) (product 2 y))) (n (sum y z)) (n (sum y (n x))))))) (n (sum u (n (sum y z)))))) u))
+
+(equal E674 tmp)
+
+; STEP 6736 -------------------------------------------------------------------
+
+; with (10), let x be 3v and y be n(n(3v)+v)+2v
+
+(setq E10p (subst (product 3 v) x E10))
+
+(setq tmp (sum (n (sum v (n (product 3 v)))) (product 2 v)))
+
+(setq E10p (eval (subst tmp y E10p)))
+
+; with (674), let x be 3v, y be v, z be 2v, and u be n(n(3v)+v)
+
+(setq E674p (subst (product 3 v) x E674))
+
+(setq E674p (subst v y E674p))
+
+(setq E674p (subst (product 2 v) z E674p))
+
+(setq tmp (n (sum v (n (product 3 v)))))
+
+(setq E674p (eval (subst tmp u E674p)))
+
+; simplify with (10')
+
+(setq E6736 (eval (subst (caddr E10p) (cadr E10p) E674p)))
+
+; replace v with x
+
+(setq E6736 (eval (subst x v E6736)))
+
+; n(n(n(n(3x)+x)+n(3x))+n(n(n(3x)+x)+5x)) = n(n(3x)+x) ?
+
+(setq tmp (eq (n (sum (n (sum (n (product 3 x)) (n (sum x (n (product 3 x)))))) (n (sum (n (sum x (n (product 3 x)))) (product 5 x))))) (n (sum x (n (product 3 x))))))
+
+(equal E6736 tmp)
+
+; STEP 8855 -------------------------------------------------------------------
+
+; with (7), let x be n(n(3x)+x)+n(3x) and y be n(n(n(3x)+x)+5x)
+
+(setq E8855 (subst Y y E7))
+
+(setq tmp (sum (n (product 3 x)) (n (sum x (n (product 3 x))))))
+
+(setq E8855 (subst tmp x E8855))
+
+(setq tmp (n (sum (n (sum x (n (product 3 x)))) (product 5 x))))
+
+(setq E8855 (eval (subst tmp Y E8855)))
+
+; simplify with (6736)
+
+(setq E8855 (eval (subst (caddr E6736) (cadr E6736) E8855)))
+
+; n(n(n(3x)+x)+n(n(n(3x)+x)+n(3x)+n(n(n(3x)+x)+5x))) = n(n(n(3x)+x)+5x) ?
+
+(setq tmp (eq (n (sum (n (sum x (n (product 3 x)))) (n (sum (n (product 3 x)) (n (sum x (n (product 3 x)))) (n (sum (n (sum x (n (product 3 x)))) (product 5 x))))))) (n (sum (n (sum x (n (product 3 x)))) (product 5 x)))))
+
+(equal E8855 tmp)
+
+; with (54), let x be 3x, z be n(3x), and y be x
+
+(setq tmp (subst (product 3 x) x E54))
+
+(setq tmp (subst (n (product 3 x)) z tmp))
+
+(setq tmp (eval (subst x y tmp)))
+
+; simplify
+
+(setq E8855 (subst (caddr tmp) (cadr tmp) E8855))
+
+; flip
+
+(setq E8855 (eq (caddr E8855) (cadr E8855)))
+
+; n(n(n(3x)+x)+5x) = n(3x) ?
+
+(setq tmp (eq (n (sum (n (sum x (n (product 3 x)))) (product 5 x))) (n (product 3 x))))
+
+(equal E8855 tmp)
+
+; STEP 8865 -------------------------------------------------------------------
+
+; with (7), let y be n(n(3x)+x)+2x and x be 3x
+
+(setq E8865 (subst X x E7))
+
+(setq tmp (sum (n (sum x (n (product 3 x)))) (product 2 x)))
+
+(setq E8865 (subst tmp y E8865))
+
+(setq E8865 (eval (subst (product 3 x) X E8865)))
+
+; simplify with (8855)
+
+(setq E8865 (subst (caddr E8855) (cadr E8855) E8865))
+
+(setq E8865 (eval E8865))
+
+; n(n(n(n(3x)+x)+n(3x)+2x)+n(3x)) = n(n(3x)+x)+2x ?
+
+(setq tmp (eq (n (sum (n (product 3 x)) (n (sum (n (product 3 x)) (n (sum x (n (product 3 x)))) (product 2 x))))) (sum (n (sum x (n (product 3 x)))) (product 2 x))))
+
+(equal E8865 tmp)
+
+; STEP 8866 -------------------------------------------------------------------
+
+; with (7), let x be n(n(3x)+x)+4x and y be x
+
+(setq tmp (sum (n (sum x (n (product 3 x)))) (product 4 x)))
+
+(setq A2 (subst tmp x E7))
+
+(setq A2 (eval (subst x y A2)))
+
+; simplify with (8855)
+
+(setq A2 (eval (subst (caddr E8855) (cadr E8855) A2)))
+
+; with (11), let x be 3x and y be x
+
+(setq A3 (subst (product 3 x) x E11))
+
+(setq A3 (eval (subst x y A3)))
+
+; simplify A2 with A3
+
+(setq E8866 (eval (subst (caddr A3) (cadr A3) A2)))
+
+; n(n(n(3x)+x)+n(3x)) = x ?
+
+(setq tmp (eq (n (sum (n (product 3 x)) (n (sum x (n (product 3 x)))))) x))
+
+(equal E8866 tmp)
+
+; STEP 8870 -------------------------------------------------------------------
+
+; with (7), let x =  n(n(3x)+x)+n(3x)
+
+(setq tmp (sum (n (product 3 x)) (n (sum x (n (product 3 x))))))
+
+(setq E8870 (eval (subst tmp x E7)))
+
+; simplify with (8866)
+
+(setq E8870 (eval (subst (caddr E8866) (cadr E8866) E8870)))
+
+; n(n(n(n(3x)+x)+n(3x)+y)+n(x+y)) = y ?
+
+(setq tmp (eq (n (sum (n (sum x y)) (n (sum y (n (product 3 x)) (n (sum x (n (product 3 x)))))))) y))
+
+(equal E8870 tmp)
+
+; STEP 8871 -------------------------------------------------------------------
+
+; with (8870), let y be 2x
+
+(setq tmp (eval (subst (product 2 x) y E8870)))
+
+; use to simplify (8865)
+
+(setq E8871 (eval (subst (caddr tmp) (cadr tmp) E8865)))
+
+; flip
+
+(setq E8871 (eq (caddr E8871) (cadr E8871)))
+
+; n(n(3x)+x)+2x = 2x ?
+
+(setq tmp (eq (sum (n (sum x (n (product 3 x)))) (product 2 x)) (product 2 x)))
+
+(equal E8871 tmp)

lisp/startup.lisp

@@ -0,0 +1,612 @@
+; lisp reads this file on startup
+
+(define minus (product -1 arg))
+(define sqrt (power arg 1/2))
+(define exp (power e arg))
+(define complex (sum arg1 (product i arg2)))
+(define complex-conjugate (eval (subst (product -1 i) i arg)))
+(define magnitude2 (product arg (complex-conjugate arg)))
+(define magnitude (sqrt (magnitude2 arg)))
+(define zerop (equal arg 0))
+
+(define contract12 (contract 1 2 arg))
+(define contract13 (contract 1 3 arg))
+(define contract14 (contract 1 4 arg))
+(define contract15 (contract 1 5 arg))
+(define contract16 (contract 1 6 arg))
+(define contract23 (contract 2 3 arg))
+(define contract24 (contract 2 4 arg))
+(define contract25 (contract 2 5 arg))
+(define contract26 (contract 2 6 arg))
+(define contract34 (contract 3 4 arg))
+(define contract35 (contract 3 5 arg))
+(define contract36 (contract 3 6 arg))
+(define contract45 (contract 4 5 arg))
+(define contract46 (contract 4 6 arg))
+(define contract56 (contract 5 6 arg))
+
+(define transpose12 (transpose 1 2 arg))
+(define transpose13 (transpose 1 3 arg))
+(define transpose14 (transpose 1 4 arg))
+(define transpose15 (transpose 1 5 arg))
+(define transpose16 (transpose 1 6 arg))
+(define transpose23 (transpose 2 3 arg))
+(define transpose24 (transpose 2 4 arg))
+(define transpose25 (transpose 2 5 arg))
+(define transpose26 (transpose 2 6 arg))
+(define transpose34 (transpose 3 4 arg))
+(define transpose35 (transpose 3 5 arg))
+(define transpose36 (transpose 3 6 arg))
+(define transpose45 (transpose 4 5 arg))
+(define transpose46 (transpose 4 6 arg))
+(define transpose56 (transpose 5 6 arg))
+
+(define caar (car (car arg)))
+(define cadr (car (cdr arg)))
+(define cdar (cdr (car arg)))
+(define cddr (cdr (cdr arg)))
+
+(define caaar (car (car (car arg))))
+(define caadr (car (car (cdr arg))))
+(define cadar (car (cdr (car arg))))
+(define caddr (car (cdr (cdr arg))))
+(define cdaar (cdr (car (car arg))))
+(define cdadr (cdr (car (cdr arg))))
+(define cddar (cdr (cdr (car arg))))
+(define cdddr (cdr (cdr (cdr arg))))
+
+(define caaaar (car (car (car (car arg)))))
+(define caaadr (car (car (car (cdr arg)))))
+(define caadar (car (car (cdr (car arg)))))
+(define caaddr (car (car (cdr (cdr arg)))))
+(define cadaar (car (cdr (car (car arg)))))
+(define cadadr (car (cdr (car (cdr arg)))))
+(define caddar (car (cdr (cdr (car arg)))))
+(define cadddr (car (cdr (cdr (cdr arg)))))
+(define cdaaar (cdr (car (car (car arg)))))
+(define cdaadr (cdr (car (car (cdr arg)))))
+(define cdadar (cdr (car (cdr (car arg)))))
+(define cdaddr (cdr (car (cdr (cdr arg)))))
+(define cddaar (cdr (cdr (car (car arg)))))
+(define cddadr (cdr (cdr (car (cdr arg)))))
+(define cdddar (cdr (cdr (cdr (car arg)))))
+(define cddddr (cdr (cdr (cdr (cdr arg)))))
+
+(define COS (sum
+  (product 1/2 (exp (product i arg)))
+  (product 1/2 (exp (product -1 i arg)))
+))
+
+(define SIN (sum
+  (product -1/2 i (exp (product i arg)))
+  (product 1/2 i (exp (product -1 i arg)))
+))
+
+(define COSH (sum
+  (product 1/2 (exp arg))
+  (product 1/2 (exp (product -1 arg)))
+))
+
+(define SINH (sum
+  (product 1/2 (exp arg))
+  (product -1/2 (exp (product -1 arg)))
+))
+
+(define adjunct (cond
+  (arg5 (adjunct4 arg1 arg2 arg3 arg4 arg5))
+  (arg4 (adjunct3 arg1 arg2 arg3 arg4))
+  (arg3 (adjunct2 arg1 arg2 arg3))
+  (t nil)
+))
+
+(define adjunct2 (prog ()
+  (setq temp00 (component arg arg2 arg2))
+  (setq temp01 (component arg arg2 arg3))
+  (setq temp10 (component arg arg3 arg2))
+  (setq temp11 (component arg arg3 arg3))
+  (return (sum
+    (product +1 temp11 (tensor arg2 arg2))
+    (product -1 temp01 (tensor arg2 arg3))
+    (product -1 temp10 (tensor arg3 arg2))
+    (product +1 temp00 (tensor arg3 arg3))
+  ))
+))
+
+(define adjunct3 (prog ()
+  (setq temp00 (component arg arg2 arg2))
+  (setq temp01 (component arg arg2 arg3))
+  (setq temp02 (component arg arg2 arg4))
+  (setq temp10 (component arg arg3 arg2))
+  (setq temp11 (component arg arg3 arg3))
+  (setq temp12 (component arg arg3 arg4))
+  (setq temp20 (component arg arg4 arg2))
+  (setq temp21 (component arg arg4 arg3))
+  (setq temp22 (component arg arg4 arg4))
+  (return (sum
+    (product +1 temp11 temp22 (tensor arg2 arg2))
+    (product -1 temp21 temp12 (tensor arg2 arg2))
+    (product -1 temp10 temp22 (tensor arg3 arg2))
+    (product +1 temp20 temp12 (tensor arg3 arg2))
+    (product +1 temp10 temp21 (tensor arg4 arg2))
+    (product -1 temp20 temp11 (tensor arg4 arg2))
+    (product -1 temp01 temp22 (tensor arg2 arg3))
+    (product +1 temp21 temp02 (tensor arg2 arg3))
+    (product +1 temp00 temp22 (tensor arg3 arg3))
+    (product -1 temp20 temp02 (tensor arg3 arg3))
+    (product -1 temp00 temp21 (tensor arg4 arg3))
+    (product +1 temp20 temp01 (tensor arg4 arg3))
+    (product +1 temp01 temp12 (tensor arg2 arg4))
+    (product -1 temp11 temp02 (tensor arg2 arg4))
+    (product -1 temp00 temp12 (tensor arg3 arg4))
+    (product +1 temp10 temp02 (tensor arg3 arg4))
+    (product +1 temp00 temp11 (tensor arg4 arg4))
+    (product -1 temp10 temp01 (tensor arg4 arg4))
+  ))
+))
+
+(define adjunct4 (prog ()
+  (setq temp00 (component arg arg2 arg2))
+  (setq temp01 (component arg arg2 arg3))
+  (setq temp02 (component arg arg2 arg4))
+  (setq temp03 (component arg arg2 arg5))
+  (setq temp10 (component arg arg3 arg2))
+  (setq temp11 (component arg arg3 arg3))
+  (setq temp12 (component arg arg3 arg4))
+  (setq temp13 (component arg arg3 arg5))
+  (setq temp20 (component arg arg4 arg2))
+  (setq temp21 (component arg arg4 arg3))
+  (setq temp22 (component arg arg4 arg4))
+  (setq temp23 (component arg arg4 arg5))
+  (setq temp30 (component arg arg5 arg2))
+  (setq temp31 (component arg arg5 arg3))
+  (setq temp32 (component arg arg5 arg4))
+  (setq temp33 (component arg arg5 arg5))
+  (return (sum
+    (product +1 temp11 temp22 temp33 (tensor arg2 arg2))
+    (product -1 temp11 temp32 temp23 (tensor arg2 arg2))
+    (product -1 temp21 temp12 temp33 (tensor arg2 arg2))
+    (product +1 temp21 temp32 temp13 (tensor arg2 arg2))
+    (product +1 temp31 temp12 temp23 (tensor arg2 arg2))
+    (product -1 temp31 temp22 temp13 (tensor arg2 arg2))
+    (product -1 temp10 temp22 temp33 (tensor arg3 arg2))
+    (product +1 temp10 temp32 temp23 (tensor arg3 arg2))
+    (product +1 temp20 temp12 temp33 (tensor arg3 arg2))
+    (product -1 temp20 temp32 temp13 (tensor arg3 arg2))
+    (product -1 temp30 temp12 temp23 (tensor arg3 arg2))
+    (product +1 temp30 temp22 temp13 (tensor arg3 arg2))
+    (product +1 temp10 temp21 temp33 (tensor arg4 arg2))
+    (product -1 temp10 temp31 temp23 (tensor arg4 arg2))
+    (product -1 temp20 temp11 temp33 (tensor arg4 arg2))
+    (product +1 temp20 temp31 temp13 (tensor arg4 arg2))
+    (product +1 temp30 temp11 temp23 (tensor arg4 arg2))
+    (product -1 temp30 temp21 temp13 (tensor arg4 arg2))
+    (product -1 temp10 temp21 temp32 (tensor arg5 arg2))
+    (product +1 temp10 temp31 temp22 (tensor arg5 arg2))
+    (product +1 temp20 temp11 temp32 (tensor arg5 arg2))
+    (product -1 temp20 temp31 temp12 (tensor arg5 arg2))
+    (product -1 temp30 temp11 temp22 (tensor arg5 arg2))
+    (product +1 temp30 temp21 temp12 (tensor arg5 arg2))
+    (product -1 temp01 temp22 temp33 (tensor arg2 arg3))
+    (product +1 temp01 temp32 temp23 (tensor arg2 arg3))
+    (product +1 temp21 temp02 temp33 (tensor arg2 arg3))
+    (product -1 temp21 temp32 temp03 (tensor arg2 arg3))
+    (product -1 temp31 temp02 temp23 (tensor arg2 arg3))
+    (product +1 temp31 temp22 temp03 (tensor arg2 arg3))
+    (product +1 temp00 temp22 temp33 (tensor arg3 arg3))
+    (product -1 temp00 temp32 temp23 (tensor arg3 arg3))
+    (product -1 temp20 temp02 temp33 (tensor arg3 arg3))
+    (product +1 temp20 temp32 temp03 (tensor arg3 arg3))
+    (product +1 temp30 temp02 temp23 (tensor arg3 arg3))
+    (product -1 temp30 temp22 temp03 (tensor arg3 arg3))
+    (product -1 temp00 temp21 temp33 (tensor arg4 arg3))
+    (product +1 temp00 temp31 temp23 (tensor arg4 arg3))
+    (product +1 temp20 temp01 temp33 (tensor arg4 arg3))
+    (product -1 temp20 temp31 temp03 (tensor arg4 arg3))
+    (product -1 temp30 temp01 temp23 (tensor arg4 arg3))
+    (product +1 temp30 temp21 temp03 (tensor arg4 arg3))
+    (product +1 temp00 temp21 temp32 (tensor arg5 arg3))
+    (product -1 temp00 temp31 temp22 (tensor arg5 arg3))
+    (product -1 temp20 temp01 temp32 (tensor arg5 arg3))
+    (product +1 temp20 temp31 temp02 (tensor arg5 arg3))
+    (product +1 temp30 temp01 temp22 (tensor arg5 arg3))
+    (product -1 temp30 temp21 temp02 (tensor arg5 arg3))
+    (product +1 temp01 temp12 temp33 (tensor arg2 arg4))
+    (product -1 temp01 temp32 temp13 (tensor arg2 arg4))
+    (product -1 temp11 temp02 temp33 (tensor arg2 arg4))
+    (product +1 temp11 temp32 temp03 (tensor arg2 arg4))
+    (product +1 temp31 temp02 temp13 (tensor arg2 arg4))
+    (product -1 temp31 temp12 temp03 (tensor arg2 arg4))
+    (product -1 temp00 temp12 temp33 (tensor arg3 arg4))
+    (product +1 temp00 temp32 temp13 (tensor arg3 arg4))
+    (product +1 temp10 temp02 temp33 (tensor arg3 arg4))
+    (product -1 temp10 temp32 temp03 (tensor arg3 arg4))
+    (product -1 temp30 temp02 temp13 (tensor arg3 arg4))
+    (product +1 temp30 temp12 temp03 (tensor arg3 arg4))
+    (product +1 temp00 temp11 temp33 (tensor arg4 arg4))
+    (product -1 temp00 temp31 temp13 (tensor arg4 arg4))
+    (product -1 temp10 temp01 temp33 (tensor arg4 arg4))
+    (product +1 temp10 temp31 temp03 (tensor arg4 arg4))
+    (product +1 temp30 temp01 temp13 (tensor arg4 arg4))
+    (product -1 temp30 temp11 temp03 (tensor arg4 arg4))
+    (product -1 temp00 temp11 temp32 (tensor arg5 arg4))
+    (product +1 temp00 temp31 temp12 (tensor arg5 arg4))
+    (product +1 temp10 temp01 temp32 (tensor arg5 arg4))
+    (product -1 temp10 temp31 temp02 (tensor arg5 arg4))
+    (product -1 temp30 temp01 temp12 (tensor arg5 arg4))
+    (product +1 temp30 temp11 temp02 (tensor arg5 arg4))
+    (product -1 temp01 temp12 temp23 (tensor arg2 arg5))
+    (product +1 temp01 temp22 temp13 (tensor arg2 arg5))
+    (product +1 temp11 temp02 temp23 (tensor arg2 arg5))
+    (product -1 temp11 temp22 temp03 (tensor arg2 arg5))
+    (product -1 temp21 temp02 temp13 (tensor arg2 arg5))
+    (product +1 temp21 temp12 temp03 (tensor arg2 arg5))
+    (product +1 temp00 temp12 temp23 (tensor arg3 arg5))
+    (product -1 temp00 temp22 temp13 (tensor arg3 arg5))
+    (product -1 temp10 temp02 temp23 (tensor arg3 arg5))
+    (product +1 temp10 temp22 temp03 (tensor arg3 arg5))
+    (product +1 temp20 temp02 temp13 (tensor arg3 arg5))
+    (product -1 temp20 temp12 temp03 (tensor arg3 arg5))
+    (product -1 temp00 temp11 temp23 (tensor arg4 arg5))
+    (product +1 temp00 temp21 temp13 (tensor arg4 arg5))
+    (product +1 temp10 temp01 temp23 (tensor arg4 arg5))
+    (product -1 temp10 temp21 temp03 (tensor arg4 arg5))
+    (product -1 temp20 temp01 temp13 (tensor arg4 arg5))
+    (product +1 temp20 temp11 temp03 (tensor arg4 arg5))
+    (product +1 temp00 temp11 temp22 (tensor arg5 arg5))
+    (product -1 temp00 temp21 temp12 (tensor arg5 arg5))
+    (product -1 temp10 temp01 temp22 (tensor arg5 arg5))
+    (product +1 temp10 temp21 temp02 (tensor arg5 arg5))
+    (product +1 temp20 temp01 temp12 (tensor arg5 arg5))
+    (product -1 temp20 temp11 temp02 (tensor arg5 arg5))
+  ))
+))
+
+(define determinant (cond
+  (arg5 (determinant4 arg1 arg2 arg3 arg4 arg5))
+  (arg4 (determinant3 arg1 arg2 arg3 arg4))
+  (arg3 (determinant2 arg1 arg2 arg3))
+  (t nil)
+))
+
+(define determinant2 (prog ()
+  (setq temp00 (component arg arg2 arg2))
+  (setq temp01 (component arg arg2 arg3))
+  (setq temp10 (component arg arg3 arg2))
+  (setq temp11 (component arg arg3 arg3))
+  (return (sum
+    (product +1 temp00 temp11)
+    (product -1 temp01 temp10)
+  ))
+))
+
+(define determinant3 (prog ()
+  (setq temp00 (component arg arg2 arg2))
+  (setq temp01 (component arg arg2 arg3))
+  (setq temp02 (component arg arg2 arg4))
+  (setq temp10 (component arg arg3 arg2))
+  (setq temp11 (component arg arg3 arg3))
+  (setq temp12 (component arg arg3 arg4))
+  (setq temp20 (component arg arg4 arg2))
+  (setq temp21 (component arg arg4 arg3))
+  (setq temp22 (component arg arg4 arg4))
+  (return (sum
+    (product +1 temp00 temp11 temp22)
+    (product -1 temp00 temp21 temp12)
+    (product -1 temp10 temp01 temp22)
+    (product +1 temp10 temp21 temp02)
+    (product +1 temp20 temp01 temp12)
+    (product -1 temp20 temp11 temp02)
+  ))
+))
+
+(define determinant4 (prog ()
+  (setq temp00 (component arg arg2 arg2))
+  (setq temp01 (component arg arg2 arg3))
+  (setq temp02 (component arg arg2 arg4))
+  (setq temp03 (component arg arg2 arg5))
+  (setq temp10 (component arg arg3 arg2))
+  (setq temp11 (component arg arg3 arg3))
+  (setq temp12 (component arg arg3 arg4))
+  (setq temp13 (component arg arg3 arg5))
+  (setq temp20 (component arg arg4 arg2))
+  (setq temp21 (component arg arg4 arg3))
+  (setq temp22 (component arg arg4 arg4))
+  (setq temp23 (component arg arg4 arg5))
+  (setq temp30 (component arg arg5 arg2))
+  (setq temp31 (component arg arg5 arg3))
+  (setq temp32 (component arg arg5 arg4))
+  (setq temp33 (component arg arg5 arg5))
+  (return (sum
+    (product +1 temp00 temp11 temp22 temp33)
+    (product -1 temp00 temp11 temp32 temp23)
+    (product -1 temp00 temp21 temp12 temp33)
+    (product +1 temp00 temp21 temp32 temp13)
+    (product +1 temp00 temp31 temp12 temp23)
+    (product -1 temp00 temp31 temp22 temp13)
+    (product -1 temp10 temp01 temp22 temp33)
+    (product +1 temp10 temp01 temp32 temp23)
+    (product +1 temp10 temp21 temp02 temp33)
+    (product -1 temp10 temp21 temp32 temp03)
+    (product -1 temp10 temp31 temp02 temp23)
+    (product +1 temp10 temp31 temp22 temp03)
+    (product +1 temp20 temp01 temp12 temp33)
+    (product -1 temp20 temp01 temp32 temp13)
+    (product -1 temp20 temp11 temp02 temp33)
+    (product +1 temp20 temp11 temp32 temp03)
+    (product +1 temp20 temp31 temp02 temp13)
+    (product -1 temp20 temp31 temp12 temp03)
+    (product -1 temp30 temp01 temp12 temp23)
+    (product +1 temp30 temp01 temp22 temp13)
+    (product +1 temp30 temp11 temp02 temp23)
+    (product -1 temp30 temp11 temp22 temp03)
+    (product -1 temp30 temp21 temp02 temp13)
+    (product +1 temp30 temp21 temp12 temp03)
+  ))
+))
+
+; P (x) = 1
+;  0
+;
+; P (x) = x
+;  1
+;
+; n P (x) =  (2n - 1) x P   (x) - (n - 1) P   (x)
+;    n                   n-1               n-2
+;
+; (legendre x n)
+
+(define legendre (cond
+  ((equal arg2 0) 1)
+  ((equal arg2 1) arg1)
+  (t (product (power arg2 -1) (sum
+    (product (sum (product 2 arg2) -1) arg1 (legendre arg1 (sum arg2 -1)))
+    (product -1 (sum arg2 -1) (legendre arg1 (sum arg2 -2)))
+  )))
+))
+
+(define nth-derivative (cond
+  ((zerop arg3) arg1)
+  (t (derivative (nth-derivative arg1 arg2 (sum arg3 -1)) arg2))
+))
+
+; L (x, a) = 1
+;  0
+;
+; L (x, a) = 1 + a - x
+;  1
+;
+; n L (x, a) = (2n + a - 1 - x) L   (x, a) + (1 - n - a) L   (x, a)
+;    n                           n-1                      n-2
+;
+; (laguerre x n a)
+;
+; only works for n <= 6 due to 32-bit integers
+
+(define laguerre (cond
+  ((equal arg2 0) 1)
+  ((equal arg2 1) (sum 1 arg3 (product -1 arg1)))
+  (t (product (power arg2 -1) (sum
+    (product (sum (product 2 arg2) arg3 -1 (product -1 arg1)) (laguerre arg1 (sum arg2 -1) arg3))
+	(product (sum 1 (product -1 arg2) (product -1 arg3)) (laguerre arg1 (sum arg2 -2) arg3))
+  )))
+))
+
+;                      -a  x   n   n    n+a  -x
+; L(x, n, a) = (1/n!) x   e  (d /dx ) (x    e  )
+
+(define laguerre2 (product
+  (power (factorial arg2) -1)
+  (power arg1 (product -1 arg3))
+  (exp arg)
+  (nth-derivative
+    (product (power arg1 (sum arg2 arg3)) (exp (product -1 arg)))
+    arg1
+    arg2
+  )
+))
+
+(define absval (cond
+  ((not (numberp arg)) nil)
+  ((lessp arg 0) (product -1 arg))
+  (t arg)
+))
+
+(define factorial (cond
+  ((zerop arg) 1)
+  (t (product arg (factorial (sum arg -1))))
+))
+
+(define laplacian (sum
+  (product
+    (power r -2)
+    (derivative (product (power r 2) (derivative arg r)) r)
+  )
+  (product
+    (power r -2)
+    (power (sin theta) -1)
+    (derivative (product (sin theta) (derivative arg theta)) theta)
+  )
+  (product
+    (power r -2)
+    (power (sin theta) -2)
+    (derivative (derivative arg phi) phi)
+  )
+))
+
+; arg1=n arg2=x
+
+(define hermite (product
+  (power -1 arg1)
+  (power e (power arg2 2))
+  (nth-derivative (power e (product -1 (power arg2 2))) arg2 arg1)
+))
+
+; generic metric
+
+(setq gdd (sum
+  (product (g00) (tensor x0 x0))
+  (product (g11) (tensor x1 x1))
+  (product (g22) (tensor x2 x2))
+  (product (g33) (tensor x3 x3))
+))
+
+(setq eta (sum
+  (product -1 (tensor x0 x0))
+  (product +1 (tensor x1 x1))
+  (product +1 (tensor x2 x2))
+  (product +1 (tensor x3 x3))
+))
+
+(define epsilon (sum
+  (product +1 (tensor x0 x1 x2 x3)) ;  1
+  (product -1 (tensor x0 x1 x3 x2)) ;  2
+  (product -1 (tensor x0 x2 x1 x3)) ;  3
+  (product +1 (tensor x0 x2 x3 x1)) ;  4
+  (product +1 (tensor x0 x3 x1 x2)) ;  5
+  (product -1 (tensor x0 x3 x2 x1)) ;  6
+  (product -1 (tensor x1 x0 x2 x3)) ;  7
+  (product +1 (tensor x1 x0 x3 x2)) ;  8
+  (product +1 (tensor x1 x2 x0 x3)) ;  9
+  (product -1 (tensor x1 x2 x3 x0)) ; 10
+  (product -1 (tensor x1 x3 x0 x2)) ; 11
+  (product +1 (tensor x1 x3 x2 x0)) ; 12
+  (product +1 (tensor x2 x0 x1 x3)) ; 13
+  (product -1 (tensor x2 x0 x3 x1)) ; 14
+  (product -1 (tensor x2 x1 x0 x3)) ; 15
+  (product +1 (tensor x2 x1 x3 x0)) ; 16
+  (product +1 (tensor x2 x3 x0 x1)) ; 17
+  (product -1 (tensor x2 x3 x1 x0)) ; 18
+  (product -1 (tensor x3 x0 x1 x2)) ; 19
+  (product +1 (tensor x3 x0 x2 x1)) ; 20
+  (product +1 (tensor x3 x1 x0 x2)) ; 21
+  (product -1 (tensor x3 x1 x2 x0)) ; 22
+  (product -1 (tensor x3 x2 x0 x1)) ; 23
+  (product +1 (tensor x3 x2 x1 x0)) ; 24
+))
+
+(define gradient (sum
+  (product (derivative arg x0) (tensor x0))
+  (product (derivative arg x1) (tensor x1))
+  (product (derivative arg x2) (tensor x2))
+  (product (derivative arg x3) (tensor x3))
+))
+
+; compute g, guu, GAMUDD, RUDDD, RDD, R, GDD, GUD and GUU from gdd
+
+(define gr (prog (tmp tmp1 tmp2)
+  (setq g (determinant gdd x0 x1 x2 x3))
+  (setq guu (product (power g -1) (adjunct gdd x0 x1 x2 x3)))
+  ; connection coefficients
+  (setq tmp (gradient gdd))
+  (setq GAMUDD (contract23 (product 1/2 guu (sum
+    tmp
+    (transpose23 tmp)
+    (product -1 (transpose12 (transpose23 tmp)))
+  ))))
+  ; riemann tensor
+  (setq tmp1 (gradient GAMUDD))
+  (setq tmp2 (contract24 (product GAMUDD GAMUDD)))
+  (setq RUDDD (sum
+    (transpose34 tmp1)
+    (product -1 tmp1)
+    (transpose23 tmp2)
+    (product -1 (transpose34 (transpose23 tmp2)))
+  ))
+  (setq RDD (contract13 RUDDD)) ; ricci tensor
+  (setq R (contract12 (contract23 (product guu RDD)))) ; ricci scalar
+  (setq GDD (sum RDD (product -1/2 gdd R))) ; einstein tensor
+  (setq GUD (contract23 (product guu GDD))) ; raise 1st index
+  (setq GUU (contract23 (product GUD guu))) ; raise 2nd index
+))
+
+; covariant derivative of a vector
+
+(define covariant-derivative (sum
+  (gradient arg)
+  (contract13 (product arg GAMUDD))
+))
+
+(define covariant-derivative-of-up-up (prog (tmp)
+  (setq tmp (product arg GAMUDD))
+  (return (sum
+    (gradient arg)
+    (transpose12 (contract14 tmp))
+    (contract24 tmp)
+  ))
+))
+
+(define covariant-derivative-of-up-down (prog (tmp)
+  (setq tmp (product arg GAMUDD))
+  (return (sum
+    (gradient arg)
+    (transpose12 (contract14 tmp))
+    (product -1 (contract23 tmp))
+  ))
+))
+
+(define covariant-derivative-of-down-down (prog (tmp)
+  (setq tmp (product arg GAMUDD))
+  (return (sum
+    (gradient arg)
+    (product -1 (transpose12 (contract13 tmp)))
+    (product -1 (contract23 tmp))
+  ))
+))
+
+(define add-integral (prog (tmp)
+  (setq tmp arg)
+  (setq tmp (subst meta-a a tmp))
+  (setq tmp (subst meta-n n tmp))
+  (setq tmp (subst meta-x x tmp))
+  (setq integral-list (append integral-list (cons tmp nil)))
+))
+
+(add-integral (list
+  '(power x -1)
+  '(log x)
+))
+
+(add-integral (list
+  '(power x n)
+  '(product (power (sum n 1) -1) (power x (sum n 1)))
+  '(numberp n)
+  '(not (equal n -1))
+))
+
+(add-integral (list
+  '(power e (product a x))
+  '(product (power a -1) (power e (product a x)))
+))
+
+;  n  ax     -1  n  ax      -1            n-1  ax
+; x  e   -> a   x  e   - n a   (integral x    e  )
+;
+; integer n, n > 0
+
+(add-integral (list
+  '(product
+    (power x n)
+    (power e (product a x))
+  )
+  '(sum
+    (product
+      (power a -1)
+      (power x n)
+      (power e (product a x))
+    )
+    (product
+      -1
+      n
+      (power a -1)
+      (integral
+        (product
+          (power x (sum n -1))
+          (power e (product a x))
+        )
+        x
+      )
+    )
+  )
+  '(integerp n)
+  '(greaterp n 0)
+))