diff --git a/doc/manual/1.png b/doc/manual/1.png new file mode 100644 index 0000000..0d279f3 Binary files /dev/null and b/doc/manual/1.png differ diff --git a/doc/manual/16.png b/doc/manual/16.png new file mode 100644 index 0000000..acdedd8 Binary files /dev/null and b/doc/manual/16.png differ diff --git a/doc/manual/16.tex b/doc/manual/16.tex new file mode 100644 index 0000000..6326cea --- /dev/null +++ b/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} +$$ + diff --git a/doc/manual/Eigenmath.tex b/doc/manual/Eigenmath.tex new file mode 100644 index 0000000..9d49dd9 --- /dev/null +++ b/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} diff --git a/doc/manual/JamesGregory.jpeg b/doc/manual/JamesGregory.jpeg new file mode 100644 index 0000000..b6f439c Binary files /dev/null and b/doc/manual/JamesGregory.jpeg differ diff --git a/doc/manual/Makefile b/doc/manual/Makefile new file mode 100644 index 0000000..60a6bab --- /dev/null +++ b/doc/manual/Makefile @@ -0,0 +1,6 @@ +all : + pdflatex Eigenmath + pdflatex Eigenmath + pdflatex Eigenmath + makeindex Eigenmath + pdflatex Eigenmath diff --git a/doc/manual/arc.png b/doc/manual/arc.png new file mode 100644 index 0000000..73e370f Binary files /dev/null and b/doc/manual/arc.png differ diff --git a/doc/manual/avogadro.tex b/doc/manual/avogadro.tex new file mode 100644 index 0000000..68b9ab0 --- /dev/null +++ b/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$$ + diff --git a/doc/manual/cardioid.png b/doc/manual/cardioid.png new file mode 100644 index 0000000..236c144 Binary files /dev/null and b/doc/manual/cardioid.png differ diff --git a/doc/manual/circle.png b/doc/manual/circle.png new file mode 100644 index 0000000..e01b81d Binary files /dev/null and b/doc/manual/circle.png differ diff --git a/doc/manual/circle2.png b/doc/manual/circle2.png new file mode 100644 index 0000000..1de17d9 Binary files /dev/null and b/doc/manual/circle2.png differ diff --git a/doc/manual/complex-numbers.tex b/doc/manual/complex-numbers.tex new file mode 100644 index 0000000..31893df --- /dev/null +++ b/doc/manual/complex-numbers.tex @@ -0,0 +1,5 @@ +\documentclass{article} +\pagestyle{empty} +\begin{document} +\include{complex} +\end{document} diff --git a/doc/manual/complex.tex b/doc/manual/complex.tex new file mode 100644 index 0000000..a89c3ed --- /dev/null +++ b/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}$$ + + diff --git a/doc/manual/curl-in-tensor-form.tex b/doc/manual/curl-in-tensor-form.tex new file mode 100644 index 0000000..09eedc5 --- /dev/null +++ b/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)$$ + diff --git a/doc/manual/derivative.tex b/doc/manual/derivative.tex new file mode 100644 index 0000000..3cacc4e --- /dev/null +++ b/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$$ diff --git a/doc/manual/draw.tex b/doc/manual/draw.tex new file mode 100644 index 0000000..32f078f --- /dev/null +++ b/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} + diff --git a/doc/manual/francois-viete.tex b/doc/manual/francois-viete.tex new file mode 100644 index 0000000..88b12f4 --- /dev/null +++ b/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$$ + diff --git a/doc/manual/fund-thm-of-calculus.tex b/doc/manual/fund-thm-of-calculus.tex new file mode 100644 index 0000000..b98cbc8 --- /dev/null +++ b/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$$ + diff --git a/doc/manual/geometric-series.tex b/doc/manual/geometric-series.tex new file mode 100644 index 0000000..594915e --- /dev/null +++ b/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$$ + diff --git a/doc/manual/greens-theorem.tex b/doc/manual/greens-theorem.tex new file mode 100644 index 0000000..c812fd1 --- /dev/null +++ b/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. diff --git a/doc/manual/how-it-works.tex b/doc/manual/how-it-works.tex new file mode 100644 index 0000000..2a10ac8 --- /dev/null +++ b/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. diff --git a/doc/manual/hydrogen-wavefunctions.tex b/doc/manual/hydrogen-wavefunctions.tex new file mode 100644 index 0000000..4007905 --- /dev/null +++ b/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$$ + diff --git a/doc/manual/integral-trick.tex b/doc/manual/integral-trick.tex new file mode 100644 index 0000000..2bcbae8 --- /dev/null +++ b/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$$ + diff --git a/doc/manual/integrals.tex b/doc/manual/integrals.tex new file mode 100644 index 0000000..a84518d --- /dev/null +++ b/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$$ diff --git a/doc/manual/lemniscate.png b/doc/manual/lemniscate.png new file mode 100644 index 0000000..69e0dff Binary files /dev/null and b/doc/manual/lemniscate.png differ diff --git a/doc/manual/line-integral.tex b/doc/manual/line-integral.tex new file mode 100644 index 0000000..82efccc --- /dev/null +++ b/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$$ + diff --git a/doc/manual/linear-algebra.tex b/doc/manual/linear-algebra.tex new file mode 100644 index 0000000..9dfd187 --- /dev/null +++ b/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)$$ diff --git a/doc/manual/list-of-functions.tex b/doc/manual/list-of-functions.tex new file mode 100644 index 0000000..5eb76f7 --- /dev/null +++ b/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. diff --git a/doc/manual/nabokov.tex b/doc/manual/nabokov.tex new file mode 100644 index 0000000..28278f3 --- /dev/null +++ b/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}$$ + diff --git a/doc/manual/parabola.png b/doc/manual/parabola.png new file mode 100644 index 0000000..adbd865 Binary files /dev/null and b/doc/manual/parabola.png differ diff --git a/doc/manual/parabola2.png b/doc/manual/parabola2.png new file mode 100644 index 0000000..1e00fb2 Binary files /dev/null and b/doc/manual/parabola2.png differ diff --git a/doc/manual/qho.tex b/doc/manual/qho.tex new file mode 100644 index 0000000..6701cd9 --- /dev/null +++ b/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$$ + diff --git a/doc/manual/sailboat.png b/doc/manual/sailboat.png new file mode 100644 index 0000000..94687be Binary files /dev/null and b/doc/manual/sailboat.png differ diff --git a/doc/manual/scripting.tex b/doc/manual/scripting.tex new file mode 100644 index 0000000..8cad659 --- /dev/null +++ b/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.$ + diff --git a/doc/manual/semicircle.png b/doc/manual/semicircle.png new file mode 100644 index 0000000..831019d Binary files /dev/null and b/doc/manual/semicircle.png differ diff --git a/doc/manual/space-shuttle-and-corvette.tex b/doc/manual/space-shuttle-and-corvette.tex new file mode 100644 index 0000000..ac5221e --- /dev/null +++ b/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}$$ + diff --git a/doc/manual/stokes-theorem-standalone.tex b/doc/manual/stokes-theorem-standalone.tex new file mode 100644 index 0000000..3624d8c --- /dev/null +++ b/doc/manual/stokes-theorem-standalone.tex @@ -0,0 +1,5 @@ +\documentclass{article} +\pagestyle{empty} +\begin{document} +\include{stokes-theorem} +\end{document} diff --git a/doc/manual/stokes-theorem.tex b/doc/manual/stokes-theorem.tex new file mode 100644 index 0000000..ecc501e --- /dev/null +++ b/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$$ + diff --git a/doc/manual/surface-area.tex b/doc/manual/surface-area.tex new file mode 100644 index 0000000..5b41ab2 --- /dev/null +++ b/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$$ + diff --git a/doc/manual/surface-integral.tex b/doc/manual/surface-integral.tex new file mode 100644 index 0000000..dac21d8 --- /dev/null +++ b/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$$ + diff --git a/doc/manual/symbols.tex b/doc/manual/symbols.tex new file mode 100644 index 0000000..bae9732 --- /dev/null +++ b/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$$ + diff --git a/doc/manual/syntax.tex b/doc/manual/syntax.tex new file mode 100644 index 0000000..e3f4c7c --- /dev/null +++ b/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} diff --git a/doc/manual/units-of-measure.tex b/doc/manual/units-of-measure.tex new file mode 100644 index 0000000..d3607ef --- /dev/null +++ b/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$$ + diff --git a/doc/manual/zerozero.png b/doc/manual/zerozero.png new file mode 100644 index 0000000..a138731 Binary files /dev/null and b/doc/manual/zerozero.png differ diff --git a/doc/manual/zerozero.tex b/doc/manual/zerozero.tex new file mode 100644 index 0000000..b4b717d --- /dev/null +++ b/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} diff --git a/doc/manual/zerozero2.png b/doc/manual/zerozero2.png new file mode 100644 index 0000000..9cb1058 Binary files /dev/null and b/doc/manual/zerozero2.png differ diff --git a/doc/ross/Makefile b/doc/ross/Makefile new file mode 100644 index 0000000..eb75f37 --- /dev/null +++ b/doc/ross/Makefile @@ -0,0 +1,2 @@ +ross.pdf : *.tex + pdftex ross diff --git a/doc/ross/ross-1.1.tex b/doc/ross/ross-1.1.tex new file mode 100644 index 0000000..43334d7 --- /dev/null +++ b/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 +}$$ diff --git a/doc/ross/ross-1.2.tex b/doc/ross/ross-1.2.tex new file mode 100644 index 0000000..31ad919 --- /dev/null +++ b/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 +}$$ diff --git a/doc/ross/ross-1.3.tex b/doc/ross/ross-1.3.tex new file mode 100644 index 0000000..089449d --- /dev/null +++ b/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 +}$$ diff --git a/doc/ross/ross-1.4.tex b/doc/ross/ross-1.4.tex new file mode 100644 index 0000000..b4891cc --- /dev/null +++ b/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$. + diff --git a/doc/ross/ross-1.5.tex b/doc/ross/ross-1.5.tex new file mode 100644 index 0000000..89e2c39 --- /dev/null +++ b/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$. + diff --git a/doc/ross/ross-14.1.tex b/doc/ross/ross-14.1.tex new file mode 100644 index 0000000..5d767a5 --- /dev/null +++ b/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. \ No newline at end of file diff --git a/doc/ross/ross-17.3.tex b/doc/ross/ross-17.3.tex new file mode 100644 index 0000000..8d9d1cc --- /dev/null +++ b/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$. diff --git a/doc/ross/ross-17.4.tex b/doc/ross/ross-17.4.tex new file mode 100644 index 0000000..0d3b652 --- /dev/null +++ b/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)$. + diff --git a/doc/ross/ross-17.5.tex b/doc/ross/ross-17.5.tex new file mode 100644 index 0000000..5907a33 --- /dev/null +++ b/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. + diff --git a/doc/ross/ross-18.1.tex b/doc/ross/ross-18.1.tex new file mode 100644 index 0000000..f8b9976 --- /dev/null +++ b/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$. + diff --git a/doc/ross/ross-18.2.tex b/doc/ross/ross-18.2.tex new file mode 100644 index 0000000..f7be884 --- /dev/null +++ b/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)$.} diff --git a/doc/ross/ross-18.3.tex b/doc/ross/ross-18.3.tex new file mode 100644 index 0000000..cf32983 --- /dev/null +++ b/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)} + diff --git a/doc/ross/ross-18.5.tex b/doc/ross/ross-18.5.tex new file mode 100644 index 0000000..50cd107 --- /dev/null +++ b/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$. + diff --git a/doc/ross/ross-18.6.tex b/doc/ross/ross-18.6.tex new file mode 100644 index 0000000..d87e34e --- /dev/null +++ b/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 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)) diff --git a/lisp/example2.tex b/lisp/example2.tex new file mode 100644 index 0000000..a8ccd99 --- /dev/null +++ b/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 diff --git a/lisp/index.html b/lisp/index.html new file mode 100644 index 0000000..eb7f86b --- /dev/null +++ b/lisp/index.html @@ -0,0 +1,39 @@ +Simple Lisp + +
+	[roscoe ~/lisp]$ gcc -o lisp lisp.c
+	[roscoe ~/lisp]$ ./lisp example.lisp
+	Bondi metric...
+	Is the Einstein tensor GUU divergence-free?
+	t
+	> (length GUU)
+	280
+	> (length (covariant-derivative-of-up-up GUU))
+	6723
+	> ^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
+	>
+
diff --git a/lisp/lisp.c b/lisp/lisp.c new file mode 100644 index 0000000..50c2c01 --- /dev/null +++ b/lisp/lisp.c @@ -0,0 +1,2846 @@ +/* 9-20-00 add "if" function */ + +#include +#include +#include +#include +#include + +#define MEM 3000000 /* (3,000,000 nodes) * (12 bytes/node) = 36 megabytes */ +#define TOS 100000 +#define BUF 1000 + +#define STRBUF 10000 + +#define CONS 0 +#define NUM 1 +#define STR 2 +#define SYM 3 +#define ZAND 4 +#define ZAPPEND 5 +#define ZARGLIST 6 +#define ZATOM 7 +#define ZCAR 8 +#define ZCDR 9 +#define ZCOMPONENT 10 +#define ZCOND 11 +#define ZCONS 12 +#define ZCONTRACT 13 +#define ZDEFINE 14 +#define ZDERIVATIVE 15 +#define ZDOT 16 +#define ZEQUAL 17 +#define ZEVAL 18 +#define ZEXIT 19 +#define ZGC 20 +#define ZGOTO 21 +#define ZGREATERP 22 +#define ZIF 23 +#define ZINTEGERP 24 +#define ZINTEGRAL 25 +#define ZLENGTH 26 +#define ZLESSP 27 +#define ZLIST 28 +#define ZLOAD 29 +#define ZNOT 30 +#define ZNULL 31 +#define ZNUMBERP 32 +#define ZOR 33 +#define ZPOWER 34 +#define ZPRINT 35 +#define ZPRINTCARS 36 +#define ZPRODUCT 37 +#define ZPROG 38 +#define ZQUOTE 39 +#define ZRETURN 40 +#define ZSETQ 41 +#define ZSUM 42 +#define ZSUBST 43 +#define ZTENSOR 44 +#define ZTRANSPOSE 45 + +#define iscons(p) ((p)->k == CONS) +#define isnum(p) ((p)->k == NUM) +#define isstr(p) ((p)->k == STR) +#define issym(p) ((p)->k >= SYM) + +#define car(p) (iscons(p) ? (p)->u.cons.car : nil) +#define cdr(p) (iscons(p) ? (p)->u.cons.cdr : nil) +#define caar(p) car(car(p)) +#define cadr(p) car(cdr(p)) +#define cdar(p) cdr(car(p)) +#define cddr(p) cdr(cdr(p)) +#define caddr(p) car(cdr(cdr(p))) +#define cadar(p) car(cdr(car(p))) +#define cddar(p) cdr(cdr(car(p))) +#define cdddr(p) cdr(cdr(cdr(p))) +#define caaddr(p) car(car(cdr(cdr(p)))) +#define cdaddr(p) cdr(car(cdr(cdr(p)))) +#define cadaddr(p) car(cdr(car(cdr(cdr(p))))) +#define cddaddr(p) cdr(cdr(car(cdr(cdr(p))))) +#define cdddaddr(p) cdr(cdr(cdr(car(cdr(cdr(p)))))) +#define caddaddr(p) car(cdr(cdr(car(cdr(cdr(p)))))) + +#define arg1(p) car(cdr(p)) +#define arg2(p) car(cdr(cdr(p))) +#define arg3(p) car(cdr(cdr(cdr(p)))) + +#define A u.num.a +#define B u.num.b + +typedef struct node { + char k, tag, gamma, pad; + union { + struct { struct node *car; struct node *cdr; } cons; + struct { struct node *binding; char *printname; } sym; + struct { int a, b; } num; + char *str; + } u; +} U; + +U _nil, mem[MEM]; + +typedef struct { + U *p; + int a, b; +} ST; + +ST stack[TOS]; + +U *_arg, *_arg1, *_arg2, *_arg3, *_arg4, *_arg5, *_arg6, *_arglist; +U *_cos; +U *derivative; +U *e; +U *eof; +U *freelist; +U *gamma[17]; +U *indexlist; +U *_integral; +U *nil; +U *nothing; +U *_power; +U *_product; +U *quote; +U *read1(); +U *read_expr(); +U *scann(); +U *_sin; +U *sum; +U *symtbl[64]; +U *t; +U *_tensor; +U *_integral_list; +U *_meta_a; +U *_meta_b; +U *_meta_f; +U *_meta_n; +U *_meta_x; + +FILE *infile; + +int tos; +int tensor_op; +int index_i; +int index_j; + +char buf[BUF]; + +int k; +char strbuf[STRBUF + 1]; + +/* gamma product matrix */ + +int gp[17][17] = { + 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, + 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, + 0,0,1,-6,-7,-8,-3,-4,-5,13,14,15,-16,9,10,11,-12, + 0,0,6,-1,-11,10,-2,-15,14,12,-5,4,-9,16,-8,7,-13, + 0,0,7,11,-1,-9,15,-2,-13,5,12,-3,-10,8,16,-6,-14, + 0,0,8,-10,9,-1,-14,13,-2,-4,3,12,-11,-7,6,16,-15, + 0,0,3,2,15,-14,1,11,-10,16,-8,7,13,12,-5,4,9, + 0,0,4,-15,2,13,-11,1,9,8,16,-6,14,5,12,-3,10, + 0,0,5,14,-13,2,10,-9,1,-7,6,16,15,-4,3,12,11, + 0,0,13,12,-5,4,16,-8,7,-1,-11,10,-3,-2,-15,14,-6, + 0,0,14,5,12,-3,8,16,-6,11,-1,-9,-4,15,-2,-13,-7, + 0,0,15,-4,3,12,-7,6,16,-10,9,-1,-5,-14,13,-2,-8, + 0,0,16,-9,-10,-11,-13,-14,-15,-3,-4,-5,1,-6,-7,-8,2, + 0,0,9,-16,8,-7,-12,5,-4,-2,-15,14,6,-1,-11,10,3, + 0,0,10,-8,-16,6,-5,-12,3,15,-2,-13,7,11,-1,-9,4, + 0,0,11,7,-6,-16,4,-3,-12,-14,13,-2,8,-10,9,-1,5, + 0,0,12,13,14,15,9,10,11,-6,-7,-8,-2,-3,-4,-5,-1, +}; + +U *ccons(); +U *num(); +U *str(); +U *sym(); +U *zand(); +U *zappend(); +U *zatom(); +U *zcar(); +U *zcdr(); +U *zcomponent(); +U *zcond(); +U *zcons(); +U *zcontract(); +U *zdefine(); +U *zderivative(); +U *zdot(); +U *zequal(); +U *zeval(); +U *zexit(); +U *zfixp(); +U *zgc(); +U *zgoto(); +U *zgreaterp(); +U *zif(); +U *zintegral(); +U *zlessp(); +U *zlist(); +U *zload(); +U *zlength(); +U *znot(); +U *znull(); +U *znumberp(); +U *zor(); +U *zpower(); +U *zprint(); +U *zprintcars(); +U *zproduct(); +U *zprog(); +U *zquote(); +U *zreturn(); +U *zsetq(); +U *zsum(); +U *zsubst(); +U *ztensor(); +U *ztranspose(); + +U *alloc(); +U *append(); +U *cons(); +U *contract(); +U *d(); +U *dcos(); +U *dd(); +U *dfunction(); +U *dpower(); +U *dproduct(); +U *dsin(); +U *dsum(); +U *eval(); +U *eval_user_function(); +U *eval_args(); +U *expand(int, int); +U *inner_expand(); +U *inner_product(U *, U *); +U *product(int, int); +U *ksum(); +U *number(); +U *pop(); +U *power(); +void push_factor(U *); +void push_term(U *); +U *subst(); +U *symbol(); +U *transpose(); + +void print(U *); +void print1(U *); +void untag(U *); +void add(int *, int *, int, int); +void mult(int *, int *, int, int); +void push(U *); +void pushvars(U *); +void popvars(U *); +void init(void); +void load(char *); +void stop(char *); +char *sdup(char *); + +main(argc, argv) +int argc; +char *argv[]; +{ + int i; + U *p; + infile = stdin; + init(); + load("startup.lisp"); + for (i = 1; i < argc; i++) + load(argv[i]); + for (;;) { + printf("> "); + p = read1(); + push(p); /* save from gc */ + p = eval(p); + pop(); + if (p != nothing) + print(p); + if (tos) { + tos = 0; + printf("stack error\n"); + } + } +} + +U * +eval(p) +U *p; +{ + if (p->k == SYM) + return p->u.sym.binding; /* symbol's binding */ + + if (p->k != CONS) + return p; /* numbers, strings, etc. */ + + switch (p->u.cons.car->k) { + case SYM: return eval_user_function(p); + case ZAND: return zand(p); + case ZAPPEND: return zappend(p); + case ZATOM: return zatom(p); + case ZCAR: return zcar(p); + case ZCDR: return zcdr(p); + case ZCOMPONENT: return zcomponent(p); + case ZCOND: return zcond(p); + case ZCONS: return zcons(p); + case ZCONTRACT: return zcontract(p); + case ZDEFINE: return zdefine(p); + case ZDERIVATIVE: return zderivative(p); + case ZDOT: return zdot(p); + case ZEQUAL: return zequal(p); + case ZEVAL: return zeval(p); + case ZEXIT: return zexit(p); + case ZGC: return zgc(p); + case ZGOTO: return zgoto(p); + case ZGREATERP: return zgreaterp(p); + case ZIF: return zif(p); + case ZINTEGERP: return zfixp(p); + case ZINTEGRAL: return zintegral(p); + case ZLENGTH: return zlength(p); + case ZLESSP: return zlessp(p); + case ZLIST: return zlist(p); + case ZLOAD: return zload(p); + case ZNOT: return znot(p); + case ZNULL: return znull(p); + case ZNUMBERP: return znumberp(p); + case ZOR: return zor(p); + case ZPOWER: return zpower(p); + case ZPRINT: return zprint(p); + case ZPRINTCARS: return zprintcars(p); + case ZPRODUCT: return zproduct(p); + case ZPROG: return zprog(p); + case ZQUOTE: return zquote(p); + case ZRETURN: return zreturn(p); + case ZSETQ: return zsetq(p); + case ZSUM: return zsum(p); + case ZSUBST: return zsubst(p); + case ZTENSOR: return ztensor(p); + case ZTRANSPOSE: return ztranspose(p); + default: return p; /* a list with ? head */ + } +} + +U * +eval_user_function(p) +U *p; +{ + U *q; + + q = eval_args(cdr(p)); + + push(_arg->u.sym.binding); + push(_arg1->u.sym.binding); + push(_arg2->u.sym.binding); + push(_arg3->u.sym.binding); + push(_arg4->u.sym.binding); + push(_arg5->u.sym.binding); + push(_arg6->u.sym.binding); + push(_arglist->u.sym.binding); + + _arglist->u.sym.binding = q; + + _arg->u.sym.binding = car(q); + _arg1->u.sym.binding = car(q); + + q = cdr(q); + _arg2->u.sym.binding = car(q); + + q = cdr(q); + _arg3->u.sym.binding = car(q); + + q = cdr(q); + _arg4->u.sym.binding = car(q); + + q = cdr(q); + _arg5->u.sym.binding = car(q); + + q = cdr(q); + _arg6->u.sym.binding = car(q); + + p = car(p); + if (p->u.sym.binding->k == CONS) + /* normal function */ + p = eval(p->u.sym.binding); + else + /* undefined function */ + p = cons(p->u.sym.binding, _arglist->u.sym.binding); + + _arglist->u.sym.binding = pop(); + _arg6->u.sym.binding = pop(); + _arg5->u.sym.binding = pop(); + _arg4->u.sym.binding = pop(); + _arg3->u.sym.binding = pop(); + _arg2->u.sym.binding = pop(); + _arg1->u.sym.binding = pop(); + _arg->u.sym.binding = pop(); + + return p; +} + +U * +eval_args(p) +U *p; +{ + int i, n = 0; + while (iscons(p)) { + push(eval(car(p))); + n++; + p = cdr(p); + } + p = nil; + for (i = 0; i < n; i++) + p = cons(stack[tos - i - 1].p, p); + tos -= n; + return p; +} + +U * +zand(p) +U *p; +{ + p = cdr(p); + while (iscons(p)) { + if (eval(car(p)) == nil) + return nil; + p = cdr(p); + } + return t; +} + +U * +zappend(p) +U *p; +{ + U *p1, *p2; + p1 = eval(arg1(p)); + push(p1); + p2 = eval(arg2(p)); + pop(); + p = append(p1, p2); + return p; +} + +U * +append(p1, p2) +U *p1, *p2; +{ + int i, n = 0; + while (iscons(p1)) { + push(car(p1)); + n++; + p1 = cdr(p1); + } + for (i = 0; i < n; i++) + p2 = cons(stack[tos - i - 1].p, p2); + tos -= n; + return p2; +} + +U * +zatom(p) +U *p; +{ + p = eval(arg1(p)); + if (iscons(p)) + return nil; + else + return t; +} + +U * +zcar(p) +U *p; +{ + return car(eval(arg1(p))); +} + +U * +zcdr(p) +U *p; +{ + return cdr(eval(arg1(p))); +} + +U * +zcomponent(p) +U *p; +{ + U *p1, *p2; + p1 = eval(arg1(p)); + push(p1); + p2 = eval_args(cddr(p)); + push(p2); + indexlist = p2; + tensor_op = 1; + p = eval(p1); + tensor_op = 0; + pop(); + pop(); + return p; +} + +U * +zcond(p) +U *p; +{ + p = cdr(p); + while (iscons(p)) { + if (eval(caar(p)) != nil) + return eval(cadar(p)); + p = cdr(p); + } + return nil; +} + +U * +zcons(p) +U *p; +{ + U *p1, *p2; + p1 = eval(arg1(p)); + push(p1); + p2 = eval(arg2(p)); + pop(); + return cons(p1, p2); +} + +U * +cons(p1, p2) +U *p1, *p2; +{ + U *p; + if (freelist == nil) { + push(p1); + push(p2); + p = alloc(); + pop(); + pop(); + } else { + p = freelist; + freelist = freelist->u.cons.cdr; + } + p->k = CONS; + p->u.cons.car = p1; + p->u.cons.cdr = p2; + return p; +} + +U * +zcontract(p) +U *p; +{ + int i, j; + i = tindex(eval(arg1(p))); + j = tindex(eval(arg2(p))); + p = eval(arg3(p)); + if (i == 0 || j == 0 || i == j) + return p; + index_i = i; + index_j = j; + tensor_op = 2; + push(p); + p = eval(p); + pop(); + tensor_op = 0; + return p; +} + +tindex(p) +U *p; +{ + if (isnum(p) && p->A > 0 && p->B == 1) + return p->A; + else + return 0; +} + +/* just like setq except arg2 is not evaluated */ + +U * +zdefine(p) +U *p; +{ + U *p1, *p2; + p1 = arg1(p); + p2 = arg2(p); + if (issym(p1)) + p1->u.sym.binding = p2; + return nothing; +} + +U * +zderivative(p) +U *p; +{ + U *p1, *p2; + p1 = eval(arg1(p)); + push(p1); + p2 = eval(arg2(p)); + push(p2); + p = d(p1, p2); + tos -= 2; + return p; +} + +U * +zdot(p) +U *p; +{ + int h = tos; + p = cdr(p); + while (iscons(p)) { + push_factor(eval(car(p))); + p = cdr(p); + } + return product(1, tos - h); +} + +U * +zequal(p) +U *p; +{ + U *p1, *p2; + p1 = eval(arg1(p)); + push(p1); + p2 = eval(arg2(p)); + pop(); + if (equal(p1, p2)) + return t; + else + return nil; +} + +equal(p1, p2) +U *p1, *p2; +{ + if (p1 == p2) + return 1; + + if (p1->k != p2->k) + return 0; + + switch (p1->k) { + case CONS: + while (iscons(p1) && iscons(p2)) { + if (!equal(car(p1), car(p2))) + return 0; + p1 = cdr(p1); + p2 = cdr(p2); + } + return equal(p1, p2); + case NUM: + if (p1->A == p2->A && p1->B == p2->B) + return 1; + else + return 0; + case STR: + if (strcmp(p1->u.str, p2->u.str) == 0) + return 1; + else + return 0; + default: + return 0; /* symbols would have matched p1 == p2 */ + } +} + +U * +zeval(p) +U *p; +{ + p = eval(arg1(p)); + push(p); + p = eval(p); + pop(); + return p; +} + +U * +zexit(p) +U *p; +{ + (void) p; + exit(1); + return nil; +} + +U * +zfixp(p) +U *p; +{ + p = eval(arg1(p)); + if (p->k == NUM && p->u.num.b == 1) + return t; + else + return nil; +} + +U * +zgc() +{ + return number(gc(), 1); +} + +gc() +{ + int i, n = 0; + + /* tag everything */ + + for (i = 0; i < MEM; i++) + mem[i].tag = 1; + + /* untag what's used */ + + for (i = 0; i < 64; i++) + untag(symtbl[i]); + + for (i = 0; i < tos; i++) + untag(stack[i].p); + + /* collect everything that's still tagged */ + + freelist = nil; + for (i = 0; i < MEM; i++) + if (mem[i].tag) { + mem[i].u.cons.cdr = freelist; + freelist = mem + i; + n++; + } + + return n; +} + +void +untag(p) +U *p; +{ + while (iscons(p) && p->tag) { + p->tag = 0; + untag(p->u.cons.car); + p = p->u.cons.cdr; + } + if (p->tag) { + p->tag = 0; + if (issym(p)) + untag(p->u.sym.binding); + } +} + +U * +zgoto(p) +U *p; +{ + return p; +} + +U * +zgreaterp(p) +U *p; +{ + + U *p1, *p2; + p1 = eval(arg1(p)); + push(p1); + p2 = eval(arg2(p)); + pop(); + if (lessp(p2, p1)) + return t; + else + return nil; +} + +U * +zif(p) +U *p; +{ + if (eval(arg1(p)) == nil) + return eval(arg3(p)); + else + return eval(arg2(p)); +} + +/* example: p = (integral (power x 2) x) */ + +U * +zintegral(p) +U *p; +{ + push(eval(arg1(p))); + push(eval(arg2(p))); + integral(); + return pop(); +} + +U * +zlessp(p) +U *p; +{ + U *p1, *p2; + p1 = eval(arg1(p)); + push(p1); + p2 = eval(arg2(p)); + pop(); + if (lessp(p1, p2)) + return t; + else + return nil; +} + +lessp(p1, p2) +U *p1, *p2; +{ + if (p1 == p2) + return 0; + + if (p1 == nil) + return 1; + + if (p2 == nil) + return 0; + + if (isnum(p1) && isnum(p2)) + if (p1->A * p2->B < p2->A * p1->B) + return 1; + else + return 0; + + if (isnum(p1)) + return 1; + + if (isnum(p2)) + return 0; + + if (isstr(p1) && isstr(p2)) + if (strcmp(p1->u.str, p2->u.str) < 0) + return 1; + else + return 0; + + if (isstr(p1)) + return 1; + + if (isstr(p2)) + return 0; + + if (issym(p1) && issym(p2)) + if (strcmp(p1->u.sym.printname, p2->u.sym.printname) < 0) + return 1; + else + return 0; + + if (issym(p1)) + return 1; + + if (issym(p2)) + return 0; + + while (iscons(p1) && iscons(p2) && equal(car(p1), car(p2))) { + p1 = cdr(p1); + p2 = cdr(p2); + } + + if (iscons(p1) && iscons(p2)) + return lessp(car(p1), car(p2)); + else + return lessp(p1, p2); +} + +U * +zlist(p) +U *p; +{ + int i, n = 0; + p = cdr(p); + while (iscons(p)) { + push(eval(car(p))); + n++; + p = cdr(p); + } + p = nil; + for (i = 0; i < n; i++) + p = cons(stack[tos - i - 1].p, p); + tos -= n; + return p; +} + +U * +zload(p) +U *p; +{ + p = eval(arg1(p)); + if (isstr(p)) + load(p->u.str); + else + printf("bad file name, string expected\n"); + return nothing; +} + +void +load(s) +char *s; +{ + U *p; + FILE *f = infile; + infile = fopen(s, "r"); + if (infile == NULL) { + printf("cannot open file %s\n", s); + infile = f; + return; + } + for (;;) { + p = read1(); + if (p == eof) + break; + push(p); + p = eval(p); + if (p != nothing) + print(p); + pop(); + } + fclose(infile); + infile = f; +} + +U * +zlength(p) +U *p; +{ + int n = 0; + p = eval(arg1(p)); + while (iscons(p)) { + n++; + p = cdr(p); + } + return number(n, 1); +} + +U * +znot(p) +U *p; +{ + if (eval(arg1(p)) == nil) + return t; + else + return nil; +} + +U * +znull(p) +U *p; +{ + if (eval(arg1(p)) == nil) + return t; + else + return nil; +} + +U * +znumberp(p) +U *p; +{ + p = eval(arg1(p)); + if (isnum(p)) + return t; + else + return nil; +} + +U * +zor(p) +U *p; +{ + p = cdr(p); + while (iscons(p)) { + if (eval(car(p)) != nil) + return t; + p = cdr(p); + } + return nil; +} + +U * +zprintcars(p) +U *p; +{ + p = eval(arg1(p)); + while (iscons(p)) { + print(car(p)); + p = cdr(p); + } + return nothing; +} + +void +add(pa, pb, a, b) +int *pa, *pb, a, b; +{ + int k; + a = *pa * b + *pb * a; + b = *pb * b; + k = gcd(a, b); + a /= k; + b /= k; + *pa = a; + *pb = b; +} + +U * +zpower(p) +U *p; +{ + push(eval(arg1(p))); + push(eval(arg2(p))); + return power(); +} + +U * +zprint(p) +U *p; +{ + p = cdr(p); + while (iscons(p)) { + print1(eval(car(p))); + p = cdr(p); + } + printf("\n"); + return nothing; +} + +void +print(p) +U *p; +{ + print1(p); + printf("\n"); +} + +void +print1(p) +U *p; +{ + static char buf[100]; + switch (p->k) { + case CONS: + printf("("); + print1(car(p)); + p = cdr(p); + while (iscons(p)) { + printf(" "); + print1(car(p)); + p = cdr(p); + } + if (p != nil) { + printf(" . "); + print1(p); + } + printf(")"); + break; + case STR: + printf("%s", p->u.str); + break; + case NUM: + if (p->B == 1) + printf("%d", p->A); + else + printf("%d/%d", p->A, p->B); + break; + default: + printf("%s", p->u.sym.printname); + break; + } +} + +U * +zproduct(p) +U *p; +{ + int h = tos; + p = cdr(p); + while (iscons(p)) { + push_factor(eval(car(p))); + p = cdr(p); + } + return product(0, tos - h); +} + +U * +zprog(p) +U *p; +{ + U *p1, *p2; + + pushvars(cadr(p)); + + p1 = cddr(p); + + while (iscons(p1)) { + p2 = eval(car(p1)); + if (car(p2)->k == ZRETURN) { + push(p2); + p2 = eval(cadr(p2)); + pop(); + popvars(cadr(p)); + return p2; + break; + } + if (car(p2)->k == ZGOTO) { + p1 = cddr(p); + p2 = cadr(p2); + while (iscons(p1) && car(p1) != p2) + p1 = cdr(p1); + } + p1 = cdr(p1); + } + + popvars(cadr(p)); + + return nothing; +} + +void +pushvars(p) +U *p; +{ + while (iscons(p)) { + if (issym(car(p))) + push(car(p)->u.sym.binding); + p = cdr(p); + } +} + +void +popvars(p) +U *p; +{ + if (iscons(p)) { + popvars(cdr(p)); + if (issym(car(p))) + car(p)->u.sym.binding = pop(); + } +} + +U * +zquote(p) +U *p; +{ + return arg1(p); +} + +U * +zreturn(p) +U *p; +{ + return p; +} + +U * +zsetq(p) +U *p; +{ + U *p1, *p2; + p1 = arg1(p); + p2 = eval(arg2(p)); + if (issym(p1)) + p1->u.sym.binding = p2; + return nothing; +} + +U * +zsubst(p) +U *p; +{ + U *p1, *p2, *p3; + + p1 = eval(arg1(p)); + push(p1); + + p2 = eval(arg2(p)); + push(p2); + + p3 = eval(arg3(p)); + push(p3); + + p = subst(p1, p2, p3); + + pop(); + pop(); + pop(); + + return p; +} + +/* substitute p1 for p2 in p3 */ + +U * +subst(p1, p2, p3) +U *p1, *p2, *p3; +{ + U *p4, *p5; + if (equal(p2, p3)) + return p1; + else if (iscons(p3)) { + p4 = subst(p1, p2, car(p3)); + push(p4); + p5 = subst(p1, p2, cdr(p3)); + pop(); + return cons(p4, p5); + } else + return p3; +} + +U * +zsum(p) +U *p; +{ + int h = tos; + p = cdr(p); + while (iscons(p)) { + push_term(eval(car(p))); + p = cdr(p); + } + return ksum(tos - h); +} + +U * +ztensor(p) +U *p; +{ + switch (tensor_op) { + default: + return cons(_tensor, eval_args(cdr(p))); + case 1: + if (equal(cdr(p), indexlist)) + return number(1, 1); + else + return number(0, 1); + case 2: + return contract(p, index_i, index_j); + case 3: + return transpose(p, index_i, index_j); + } +} + +U * +contract(p, i, j) +U *p; +int i, j; +{ + U *p1; + int k, n; + ST *s; + p1 = p; + n = 0; + s = stack + tos; + while (iscons(p1)) { + push(car(p1)); + n++; + p1 = cdr(p1); + } + if (i >= n || j >= n) { + tos -= n; + return p; + } + if (!equal(s[i].p, s[j].p)) { + tos -= n; + return number(0, 1); + } + if (n == 3) { + tos -= n; + return number(1, 1); + } + p = nil; + for (k = n - 1; k >= 0; k--) + if (k != i && k != j) + p = cons(s[k].p, p); + tos -= n; + return p; +} + +U * +transpose(p, i, j) +U *p; +int i, j; +{ + U *p1; + int k, n; + ST *s; + p1 = p; + n = 0; + s = stack + tos; + while (iscons(p1)) { + push(car(p1)); + n++; + p1 = cdr(p1); + } + if (i >= n || j >= n) { + tos -= n; + return p; + } + p = nil; + for (k = n - 1; k >= 0; k--) + if (k == i) + p = cons(s[j].p, p); + else if (k == j) + p = cons(s[i].p, p); + else + p = cons(s[k].p, p); + tos -= n; + return p; +} + +void +mult(pa, pb, a, b) +int *pa, *pb, a, b; +{ + int k; + a *= *pa; + b *= *pb; + k = gcd(a, b); + a /= k; + b /= k; + *pa = a; + *pb = b; +} + +U * +ztranspose(p) +U *p; +{ + int i, j; + i = tindex(eval(arg1(p))); + j = tindex(eval(arg2(p))); + p = eval(arg3(p)); + if (i == 0 || j == 0 || i == j) + return p; + index_i = i; + index_j = j; + tensor_op = 3; + push(p); + p = eval(p); + pop(); + tensor_op = 0; + return p; +} + +void +init() +{ + int i; + + nil = &_nil; + nil->k = SYM; + nil->u.sym.binding = nil; + nil->u.sym.printname = "nil"; + + for (i = 0; i < 64; i++) + symtbl[i] = nil; + + freelist = nil; + for (i = 0; i < MEM; i++) { + /* mem[i].gamma = 0; */ + mem[i].u.cons.cdr = freelist; + freelist = mem + i; + } + + symbol("and") ->k = ZAND; + symbol("append") ->k = ZAPPEND; + symbol("atom") ->k = ZATOM; + symbol("car") ->k = ZCAR; + symbol("cdr") ->k = ZCDR; + symbol("component") ->k = ZCOMPONENT; + symbol("cond") ->k = ZCOND; + symbol("cons") ->k = ZCONS; + symbol("contract") ->k = ZCONTRACT; + symbol("define") ->k = ZDEFINE; + symbol("derivative") ->k = ZDERIVATIVE; + symbol("equal") ->k = ZEQUAL; + symbol("eval") ->k = ZEVAL; + symbol("exit") ->k = ZEXIT; + symbol("gc") ->k = ZGC; + symbol("goto") ->k = ZGOTO; + symbol("greaterp") ->k = ZGREATERP; + symbol("if") ->k = ZIF; + symbol("integerp") ->k = ZINTEGERP; + symbol("integral") ->k = ZINTEGRAL; + symbol("lessp") ->k = ZLESSP; + symbol("list") ->k = ZLIST; + symbol("load") ->k = ZLOAD; + symbol("run") ->k = ZLOAD; + symbol("length") ->k = ZLENGTH; + symbol("not") ->k = ZNOT; + symbol("null") ->k = ZNULL; + symbol("numberp") ->k = ZNUMBERP; + symbol("or") ->k = ZOR; + symbol("power") ->k = ZPOWER; + symbol("print") ->k = ZPRINT; + symbol("printcars") ->k = ZPRINTCARS; + symbol("product") ->k = ZPRODUCT; + symbol("prog") ->k = ZPROG; + symbol("quote") ->k = ZQUOTE; + symbol("return") ->k = ZRETURN; + symbol("setq") ->k = ZSETQ; + symbol("sum") ->k = ZSUM; + symbol("subst") ->k = ZSUBST; + symbol("tensor") ->k = ZTENSOR; + symbol("transpose") ->k = ZTRANSPOSE; + symbol("dot") ->k = ZDOT; + + gamma[2] = symbol("gamma0"); + gamma[3] = symbol("gamma1"); + gamma[4] = symbol("gamma2"); + gamma[5] = symbol("gamma3"); + gamma[6] = symbol("sigma1"); /* (product gamma1 gamma0) */ + gamma[7] = symbol("sigma2"); /* (product gamma2 gamma0) */ + gamma[8] = symbol("sigma3"); /* (product gamma3 gamma0) */ + gamma[9] = symbol("isigma1"); /* (product gamma3 gamma2) */ + gamma[10] = symbol("isigma2"); /* (product gamma1 gamma3) */ + gamma[11] = symbol("isigma3"); /* (product gamma2 gamma1) */ + gamma[12] = symbol("igamma0"); /* (product gamma3 gamma2 gamma1) */ + gamma[13] = symbol("igamma1"); /* (product gamma0 gamma3 gamma2) */ + gamma[14] = symbol("igamma2"); /* (product gamma0 gamma1 gamma3) */ + gamma[15] = symbol("igamma3"); /* (product gamma0 gamma2 gamma1) */ + gamma[16] = symbol("i"); /* (product gamma0 gamma1 gamma2 gamma3) */ + + gamma[2]->gamma = 2; + gamma[3]->gamma = 3; + gamma[4]->gamma = 4; + gamma[5]->gamma = 5; + gamma[6]->gamma = 6; + gamma[7]->gamma = 7; + gamma[8]->gamma = 8; + gamma[9]->gamma = 9; + gamma[10]->gamma = 10; + gamma[11]->gamma = 11; + gamma[12]->gamma = 12; + gamma[13]->gamma = 13; + gamma[14]->gamma = 14; + gamma[15]->gamma = 15; + gamma[16]->gamma = 16; + + _arg = symbol("arg"); + _arg1 = symbol("arg1"); + _arg2 = symbol("arg2"); + _arg3 = symbol("arg3"); + _arg4 = symbol("arg4"); + _arg5 = symbol("arg5"); + _arg6 = symbol("arg6"); + _arglist = symbol("arglist"); + _cos = symbol("cos"); + derivative = symbol("derivative"); + e = symbol("e"); + eof = symbol("eof"); + _integral = symbol("integral"); + _integral_list = symbol("integral-list"); + _meta_a = symbol("meta-a"); + _meta_b = symbol("meta-b"); + _meta_f = symbol("meta-f"); + _meta_n = symbol("meta-n"); + _meta_x = symbol("meta-x"); + nothing = symbol("nothing"); + quote = symbol("quote"); + _power = symbol("power"); + _product = symbol("product"); + _sin = symbol("sin"); + sum = symbol("sum"); + t = symbol("t"); + _tensor = symbol("tensor"); +} + +void +stop(char *s) +{ + printf("\nERROR %s\n", s); + exit(1); +} + +/* the temp stack keeps intermediate results from garbage collection */ + +void +push(p) +U *p; +{ + if (tos == TOS) + stop("stack overflow"); + stack[tos++].p = p; +} + +U * +pop() +{ + if (tos == 0) + stop("stack underflow"); + return stack[--tos].p; +} + +void +push_term(p) +U *p; +{ + if (car(p) == sum) { + p = cdr(p); + while (iscons(p)) { + push(car(p)); + p = cdr(p); + } + } else + push(p); +} + +void +push_factor(p) +U *p; +{ + if (car(p) == _product) { + p = cdr(p); + while (iscons(p)) { + push(car(p)); + p = cdr(p); + } + } else + push(p); +} + +U * +alloc() +{ + U *p; + if (freelist == nil) { + gc(); + if (freelist == nil) + stop("out of memory"); + } + p = freelist; + freelist = freelist->u.cons.cdr; + return p; +} + +U * +number(a, b) +int a, b; +{ + int k; + U *p; + k = gcd(a, b); + a /= k; + b /= k; + p = alloc(); + p->k = NUM; + p->A = a; + p->B = b; + return p; +} + +gcd(a, b) +int a, b; +{ + int k, sign; + if (b == 0) + stop("divide by zero"); + if (a == 0) + return b; + if (a < 0) + a = -a; + if (b < 0) { + b = -b; + sign = -1; + } else + sign = 1; + while (b) { + k = a % b; + a = b; + b = k; + } + return sign * a; +} + +char * +sdup(char *s) +{ + int j = k; + while (*s) { + if (k >= STRBUF) + stop("string buffer full"); + strbuf[k++] = *s++; + } + strbuf[k++] = 0; + return strbuf + j; +} + +U * +symbol(s) +char *s; +{ + U *p; + int x; + if (strcmp(s, "nil") == 0) + return nil; + x = *s & 63; + p = symtbl[x]; + while (iscons(p)) { + if (strcmp(car(p)->u.sym.printname, s) == 0) + return car(p); + p = cdr(p); + } + p = alloc(); + p->k = SYM; + p->u.sym.printname = sdup(s); + p->u.sym.binding = p; + p = cons(p, symtbl[x]); + symtbl[x] = p; + return car(p); +} + +U * +read1() +{ + return read_expr(read_token()); +} + +U * +read_expr(c) +int c; +{ + int i, n; + U *p; + + while (c == ')' || c == '.') { + printf("discarding unexpected %c\n", c); + c = read_token(); + } + + switch (c) { + + case EOF: + return eof; + + case '(': + n = 0; + c = read_token(); + while (c != EOF && c != ')' && c != '.') { + p = read_expr(c); + push(p); + n++; + c = read_token(); + } + if (c == '.') { + c = read_token(); + p = read_expr(c); + c = read_token(); + } else + p = nil; + if (c != ')') + printf("missing )\n"); + for (i = 0; i < n; i++) + p = cons(stack[tos - i - 1].p, p); + tos -= n; + return p; + + case '\'': + c = read_token(); + p = read_expr(c); + p = cons(p, nil); + p = cons(quote, p); + return p; + + case 'n': + return scann(buf); + + case 'q': + p = alloc(); + p->k = STR; + p->u.str = sdup(buf); + return p; + + case 's': + return symbol(buf); + + default: + stop("bug in read_expr()"); + return nil; + } +} + +read_token() +{ + int c, i; + + /* skip spaces and comments */ + + for (;;) { + c = fgetc(infile); + if (c == ';') + while (c != EOF && c != '\n') + c = fgetc(infile); + if (c == EOF || !isspace(c)) + break; + } + + switch (c) { + + case EOF: + return EOF; + + case '(': + case ')': + case '.': + case '\'': + return c; + + case '\"': + for (i = 0; i < BUF; i++) { + c = fgetc(infile); + if (c == EOF) { + printf("unexpected end of file\n"); + break; + } + if (c == '\"') + break; + buf[i] = c; + } + if (i == BUF) + stop("input buffer overflow"); + buf[i] = 0; + return 'q'; /* quoted string */ + + default: + buf[0] = c; + for (i = 1; i < BUF; i++) { + c = fgetc(infile); + if (c == EOF) + break; + if (isspace(c) || c == '(' || c == ')' + || c == ';' || c == '.' || c == '\'') { + ungetc(c, infile); + break; + } + buf[i] = c; + } + if (i == BUF) + stop("input buffer overflow"); + buf[i] = 0; + if (*buf == '+' || *buf == '-' || isdigit(*buf)) + return 'n'; /* number */ + else + return 's'; /* symbol */ + } +} + +U * +scann(s) +char *s; +{ + int a, b, sgn; + + a = 0; + b = 1; + sgn = 1; + + if (*s == '+') + s++; + else if (*s == '-') { + s++; + sgn = -1; + } + + if (!isdigit(*s)) { + printf("syntax error in number\n"); + return nil; + } + + while (isdigit(*s)) { + a = 10 * a + *s - '0'; + s++; + } + + if (*s == '/') { + s++; + b = 0; + while (isdigit(*s)) { + b = 10 * b + *s - '0'; + s++; + } + } + + if (*s) { + printf("syntax error in number\n"); + return nil; + } + + return number(sgn * a, b); +} + +cmp(a, b) +ST *a, *b; +{ + if (equal(a->p, b->p)) + return 0; + else if (lessp(a->p, b->p)) + return -1; + else + return 1; +} + +U * +ksum(n) +int n; +{ + int a, b, i, j; + U *p, *q; + ST *s; + + s = stack + tos - n; + + /* add numbers */ + + a = 0; + b = 1; + + for (i = 0; i < n; i++) { + p = s[i].p; + if (isnum(p)) { + add(&a, &b, p->u.num.a, p->u.num.b); + s[i].p = nil; + } + } + + /* remove numeric coefficients */ + + for (i = 0; i < n; i++) { + p = s[i].p; + if (car(p) == _product && isnum(cadr(p))) { + s[i].a = cadr(p)->u.num.a; + s[i].b = cadr(p)->u.num.b; + if (iscons(cdddr(p))) + p = cons(_product, cddr(p)); + else + p = caddr(p); + s[i].p = p; + } else { + s[i].a = 1; + s[i].b = 1; + } + } + + /* sort terms */ + + qsort(s, n, sizeof (ST), cmp); + + /* combine terms */ + + for (i = 0; i < n - 1; i++) { + if (s[i].p == nil) + continue; + for (j = i + 1; j < n; j++) { + if (!equal(s[i].p, s[j].p)) + break; + add(&s[i].a, &s[i].b, s[j].a, s[j].b); + s[j].p = nil; + if (s[i].a == 0) { + s[i].p = nil; + break; + } + } + } + + /* put the coefficients back */ + + for (i = 0; i < n; i++) { + p = s[i].p; + if (p == nil) + continue; + if (s[i].a == 1 && s[i].b == 1) + continue; + q = number(s[i].a, s[i].b); + if (car(p) == _product) + p = cdr(p); + else { + push(q); /* save q from gc */ + p = cons(p, nil); + pop(); + } + p = cons(q, p); + p = cons(_product, p); + s[i].p = p; + } + + /* add number */ + + if (a != 0) { + for (i = 0; i < n; i++) + if (s[i].p == nil) + break; + if (i == n) + stop("bug in ksum()"); + s[i].p = number(a, b); + } + + /* remove nils */ + + j = 0; + for (i = 0; i < n; i++) + if (s[i].p != nil) + s[j++].p = s[i].p; + + /* sort terms */ + + qsort(s, j, sizeof (ST), cmp); + + /* result */ + + switch (j) { + case 0: + p = number(0, 1); + break; + case 1: + p = s[0].p; + break; + default: + p = nil; + for (i = j - 1; i >= 0; i--) + p = cons(s[i].p, p); + p = cons(sum, p); + break; + } + + tos -= n; + + return p; +} + +U * +product(inner, n) +int inner, n; +{ + int a, b, flag, g, i, j, x; + U *p, *q; + ST *s; + + s = stack + tos - n; + + /* check for product of sum or sums */ + + for (i = 0; i < n; i++) + if (car(s[i].p) == sum) + return expand(inner, n); + + /* multiply numbers, tensors and gammas */ + + a = 1; + b = 1; + g = -1; + + for (i = 0; i < n; i++) { + p = s[i].p; + if (isnum(p)) { + if (p->A == 0) { + tos -= n; + return number(0, 1); + } + mult(&a, &b, p->A, p->B); + s[i].p = nil; + } else if (p->gamma) { + if (g == -1) + g = i; + else { + x = gp[s[g].p->gamma][p->gamma]; + if (x < 0) { + x = -x; + a *= -1; + } + if (x == 1) { + s[g].p = nil; + s[i].p = nil; + g = -1; + } else { + s[g].p = gamma[x]; + s[i].p = nil; + } + } + } + } + + /* multiply tensors */ + + for (i = 0; i < n; i++) { + if (car(s[i].p) == _tensor) { + for (j = i + 1; j < n; j++) { + if (car(s[j].p) == _tensor) { + if (inner) { + s[i].p = inner_product(s[i].p, s[j].p); + s[j].p = nil; + if (isnum(s[i].p)) { + if (s[i].p->A == 0) { + tos -= n; + return number(0, 1); + } else { + s[i].p = nil; + i = j; + break; + } + } + } else { + s[i].p = append(s[i].p, cdr(s[j].p)); + s[j].p = nil; + } + } + } + } + } + + prep_exponents(n); + + /* sort factors */ + + qsort(s, n, sizeof (ST), cmp); + + /* combine factors */ + + for (i = 0; i < n - 1; i++) { + if (s[i].p == nil) + continue; + for (j = i + 1; j < n; j++) { + if (!equal(s[i].p, s[j].p)) + break; + add(&s[i].a, &s[i].b, s[j].a, s[j].b); + s[j].p = nil; + if (s[i].a == 0) { + s[i].p = nil; + break; + } +/* quick bug fix while working on general solution... */ +/* result is a number with exponent 1 */ + if (isnum(s[i].p) && s[i].a == 1 && s[i].b == 1) { + mult(&a, &b, s[i].p->A, s[i].p->B); + s[i].p = nil; + break; + } +/* result is a number with exponent -1 */ + if (isnum(s[i].p) && s[i].a == -1 && s[i].b == 1) { + mult(&a, &b, s[i].p->B, s[i].p->A); + s[i].p = nil; + break; + } + } + } + + /* restore exponents */ + + for (i = 0; i < n; i++) { + p = s[i].p; + if (p == nil) + continue; + if (s[i].a == 1 && s[i].b == 1) + continue; + q = number(s[i].a, s[i].b); + if (car(p) == _power) { + if (caaddr(p) == _product) + q = cons(q, cdaddr(p)); + else { + push(q); /* save from gc */ + q = cons(q, cons(caddr(p), nil)); + pop(); + } + q = cons(_product, q); + q = cons(q, nil); + q = cons(cadr(p), q); + } else + q = cons(p, cons(q, nil)); + q = cons(_power, q); + if (cadr(q) == sum && isnum(caddr(q))) + q = eval(q); + s[i].p = q; + } + + /* restore coefficient */ + + if (a != 1 || b != 1) { + for (i = 0; i < n; i++) + if (s[i].p == nil) + break; + if (i == n) + stop("bug in product()"); + s[i].p = number(a, b); + } + + /* remove nils */ + + j = 0; + flag = 0; + for (i = 0; i < n; i++) { + p = s[i].p; + if (p != nil) { + s[j++].p = p; + if (car(p) == sum) + flag = 1; + } + } + + if (flag) { + tos = tos - n + j; + return expand(inner, j); + } + + /* sort factors */ + + qsort(s, j, sizeof (ST), cmp); + + /* result */ + + switch (j) { + case 0: + p = number(1, 1); + break; + case 1: + p = s[0].p; + break; + default: + p = nil; + for (i = j - 1; i >= 0; i--) + p = cons(s[i].p, p); + p = cons(_product, p); + break; + } + + tos -= n; + + return p; +} + +/* Here are the cases that must be handled: + +input output + +s[i].p s[i].a s[i].b s[i].p +------ ------ ------ ------ + +a 1 1 a + +(power a b) 1 1 (power a b) + +(power a 1/2) 1 2 a + +(power a (product 1/2 b)) 1 2 (power a b) + +(power a (product 1/2 b c)) 1 2 (power a (product b c)) + + +Example for navigating cars and cdrs: + +x = (power a (product 3 b c)) + +(car x) -> power +(cdr x) -> (a (product 3 b c)) + +(cadr x) -> a +(cddr x) -> ((product 3 b c)) + +(caddr x) -> (product 3 b c) + +(caaddr x) -> product +(cdaddr x) -> (3 b c) + +(cadaddr x) -> 3 +(cddaddr x) -> (b c) + +(caddaddr x) -> b +(cdddaddr x) -> (c) + +*/ + +prep_exponents(n) +{ + int i; + U *p, *q; + for (i = tos - n; i < tos; i++) { + stack[i].a = 1; + stack[i].b = 1; + p = stack[i].p; + if (car(p) != _power) + continue; + if (isnum(caddr(p))) { + stack[i].p = cadr(p); + stack[i].a = caddr(p)->A; + stack[i].b = caddr(p)->B; + } else if (caaddr(p) == _product && isnum(cadaddr(p))) { + if (cdddaddr(p) == nil) + q = caddaddr(p); + else + q = cons(_product, cddaddr(p)); + q = cons(q, nil); + q = cons(cadr(p), q); + q = cons(_power, q); + stack[i].p = q; + stack[i].a = cadaddr(p)->A; + stack[i].b = cadaddr(p)->B; + } + } +} + +U * +expand(inner, n) +int inner, n; +{ + U *p; + int h, h1, h2, i, j, k, m; + + h = tos - n; + + /* stack[h + i].a = stack index of factor i */ + + /* stack[h + i].b = number of terms in factor i */ + + /* m = number of permutations */ + + m = 1; + for (i = 0; i < n; i++) { + p = stack[h + i].p; + if (car(p) == sum) { + stack[h + i].a = tos; + j = 0; + p = cdr(p); + while (iscons(p)) { + push(car(p)); + p = cdr(p); + j++; + } + } else { + stack[h + i].a = h + i; + j = 1; + } + stack[h + i].b = j; + m = j * m; + } + + h1 = tos; + + for (i = 0; i < m; i++) { + k = i; + h2 = tos; + for (j = 0; j < n; j++) { + p = stack[stack[h + j].a + k % stack[h + j].b].p; + k = k / stack[h + j].b; + push_factor(p); + } + push_term(product(inner, tos - h2)); + } + + p = ksum(tos - h1); + + tos = h; + + return p; +} + +U * +power() +{ + int a, b, h, i, n; + U *p1, *p2; + + p1 = stack[tos - 2].p; + p2 = stack[tos - 1].p; + + /* (power ? 0) -> 1 */ + + if (isnum(p2) && p2->A == 0) { + tos -= 2; + return number(1, 1); + } + + /* (power ? 1) -> ? */ + + if (isnum(p2) && p2->A == 1 && p2->B == 1) { + tos -= 2; + return p1; + } + + /* (power 1 ?) -> 1 */ + + if (isnum(p1) && p1->A == 1 && p1->B == 1) { + tos -= 2; + return p1; + } + + /* (power a (sum b c)) -> (product (power a b) (power a c)) */ + + if (car(p2) == sum) { + h = tos; + p2 = cdr(p2); + while (iscons(p2)) { + push(p1); + push(car(p2)); + push_factor(power()); + p2 = cdr(p2); + } + p1 = product(0, tos - h); + tos -= 2; + return p1; + } + + /* (power (product a b) c) -> (product (power a c) (power b c)) */ + + if (car(p1) == _product) { + h = tos; + p1 = cdr(p1); + while (iscons(p1)) { + push(car(p1)); + push(p2); + push_factor(power()); + p1 = cdr(p1); + } + p1 = product(0, tos - h); + tos -= 2; + return p1; + } + + /* (power (sum a b) 2) -> (sum (power a 2) (power b 2) (product 2 a b)) */ + + if (car(p1) == sum && isnum(p2) && p2->A > 1 && p2->B == 1) { + push(p1); + for (i = 1; i < p2->A; i++) { + push(p1); + push(expand(0, 2)); + } + p1 = pop(); + tos -= 2; + return p1; + } + + /* (power (power a b) c) -> (power a (product b c)) */ + + if (car(p1) == _power) { + push(cadr(p1)); + h = tos; + push_factor(caddr(p1)); + push_factor(p2); + push(product(0, tos - h)); + p1 = power(); + tos -= 2; + return p1; + } + + if (isnum(p1) && isnum(p2) && p2->B == 1) { + a = p1->A; + b = p1->B; + if (p2->A > 0) + n = p2->A; + else + n = -p2->A; + for (i = 1; i < n; i++) + mult(&a, &b, p1->A, p2->B); + if (p2->A > 0) + p1 = number(a, b); + else + p1 = number(b, a); + tos -= 2; + return p1; + } + + p1 = cons(_power, cons(p1, cons(p2, nil))); + tos -= 2; + return p1; +} + +U * +d(p, dx) +U *p, *dx; +{ + // for example, (derivative x x) is 1 + + if (equal(p, dx)) + return number(1, 1); + + // for example, (derivative a x) is 0 + + if (!iscons(p) || car(p) == _tensor) + return number(0, 1); + + if (car(p) == sum) + return dsum(p, dx); + + if (car(p) == _product) + return dproduct(p, dx); + + if (car(p) == _power) + return dpower(p, dx); + + if (car(p) == derivative) + return dd(p, dx); + + if (car(p) == _sin) + return dsin(p, dx); + + if (car(p) == _cos) + return dcos(p, dx); + + return dfunction(p, dx); +} + +U * +dsum(p, dx) +U *p, *dx; +{ + int h = tos; + p = cdr(p); + while (iscons(p)) { + push_term(d(car(p), dx)); /* result may be a sum */ + p = cdr(p); + } + return ksum(tos - h); +} + +U * +dproduct(p, dx) +U *p, *dx; +{ + int h1, h2, i, j, n = 0; + U *p1; + p1 = cdr(p); + while (iscons(p1)) { + n++; + p1 = cdr(p1); + } + h1 = tos; + for (i = 0; i < n; i++) { + h2 = tos; + p1 = cdr(p); + for (j = 0; j < n; j++) { + if (i == j) + push_factor(d(car(p1), dx)); + else + push(car(p1)); + p1 = cdr(p1); + } + push_term(product(0, tos - h2)); + } + return ksum(tos - h1); +} + +U * +dpower(p, dx) +U *p, *dx; +{ + int h = tos; + U *base, *exponent; + base = cadr(p); + exponent = caddr(p); + if (base == e) { + /* d exponent * e ^ exponent */ + push_factor(d(exponent, dx)); /* result may be a product */ + push(p); + return product(0, tos - h); + } else { + /* exponent * base ^ (exponent - 1) * d base */ + push_factor(exponent); /* exponent may be a product */ + push(base); + push(exponent); /* exponent is never a sum */ + push(number(-1, 1)); + push(ksum(2)); + push(power()); + push_factor(d(base, dx)); + return product(0, tos - h); + } +} + +/* derivative of derivative */ + +U * +dd(p, dx2) +U *p, *dx2; +{ + U *dx1; + dx1 = caddr(p); + p = d(cadr(p), dx2); + push(p); + if (car(p) == derivative) { + dx2 = caddr(p); + p = cadr(p); + if (lessp(dx1, dx2)) { + p = cons(derivative, cons(p, cons(dx1, nil))); + pop(); + push(p); + p = cons(derivative, cons(p, cons(dx2, nil))); + } else { + p = cons(derivative, cons(p, cons(dx2, nil))); + pop(); + push(p); + p = cons(derivative, cons(p, cons(dx1, nil))); + } + } else + p = d(p, dx1); + pop(); + return p; +} + +/* derivative of a generic function */ + +U * +dfunction(p, dx) +U *p, *dx; +{ + U *t; + + t = cdr(p); + + /* no arg list? */ + + if (t == nil) + return cons(derivative, cons(p, cons(dx, nil))); + + /* check arg list */ + + while (iscons(t)) { + if (equal(car(t), dx)) + return cons(derivative, cons(p, cons(dx, nil))); + t = cdr(t); + } + + return number(0, 1); +} + +U * +dsin(p, dx) +U *p, *dx; +{ + int h = tos; + push_factor(d(cadr(p), dx)); + push(cons(_cos, cons(cadr(p), nil))); + return product(0, tos - h); +} + +U * +dcos(p, dx) +U *p, *dx; +{ + int h = tos; + push(number(-1, 1)); + push_factor(d(cadr(p), dx)); + push(cons(_sin, cons(cadr(p), nil))); + return product(0, tos - h); +} + +/* returns the inner product of two tensor args + + p1 p2 return + + (tensor 1) (tensor 1) 1 + + (tensor 1) (tensor 2) 0 + + (tensor 1 2) (tensor 1) 0 + + (tensor 1 2) (tensor 2) (tensor 1) + + (tensor 1 2) (tensor 2 1) (tensor 1 1) + +*/ + +U * +inner_product(p1, p2) +U *p1, *p2; +{ + int i, n = 0; + while (iscons(p1)) { + push(car(p1)); + n++; + p1 = cdr(p1); + } + if (equal(stack[tos - 1].p, cadr(p2))) { + p2 = cddr(p2); + for (i = 1; i < n; i++) + p2 = cons(stack[tos - i - 1].p, p2); + if (cdr(p2) == nil) + p2 = number(1, 1); + } else + p2 = number(0, 1); + tos -= n; + return p2; +} + +#define f (_meta_f->u.sym.binding) +#define x (_meta_x->u.sym.binding) +#define a (_meta_a->u.sym.binding) +#define n (_meta_n->u.sym.binding) + +integral() +{ + U *tmp1, *tmp2; + + tmp2 = pop(); + tmp1 = pop(); + + push(f); + push(x); + push(a); + push(n); + + f = tmp1; + x = tmp2; + + if (does_not_depend_on_x(f)) + integral_of_constant(); + else if (car(f) == sum) + integral_of_sum(); + else if (car(f) == _product) + integral_of_product(); + else + integral_of_nub(); + + tmp1 = pop(); + + n = pop(); + a = pop(); + x = pop(); + f = pop(); + + push(tmp1); +} + +integral_of_constant() +{ + push(f); + push(x); + push(product(0, 2)); +} + +integral_of_sum() +{ + int h; + U *p; + + h = tos; + + p = cdr(f); + while (iscons(p)) { + push(car(p)); + push(x); + integral(); + push_term(pop()); + p = cdr(p); + } + + p = ksum(tos - h); + + push(p); +} + +integral_of_product() +{ + int h1, h2; + U *p; + + h1 = tos; + + p = cdr(f); + while (iscons(p)) { + if (does_not_depend_on_x(car(p))) + push(car(p)); + p = cdr(p); + } + + if (tos == h1) { + integral_of_nub(); + return; + } + + h2 = tos; + + p = cdr(f); + while (iscons(p)) { + if (depends_on_x(car(p))) + push(car(p)); + p = cdr(p); + } + + push(product(0, tos - h2)); + + push(x); + + integral(); + + push_factor(pop()); + + push(product(0, tos - h1)); +} + +integral_of_nub() +{ + int h; + U *p; + + h = tos; + + push(number(1, 1)); /* so meta variables can try the value 1 */ + + deconstruct(f); + + p = _integral_list->u.sym.binding; + + while (iscons(p)) { + if (match(caar(p), cddar(p), h)) { + p = eval(cadar(p)); + tos = h; + push(p); + return; + } + p = cdr(p); + } + + p = cons(_integral, cons(f, cons(x, nil))); + + tos = h; + + push(p); +} + +depends_on_x(p) +U *p; +{ + if (findx(p) || findf(p)) + return 1; + else + return 0; +} + +does_not_depend_on_x(p) +U *p; +{ + if (findx(p) || findf(p)) + return 0; + else + return 1; +} + +/* yes if dx somewhere in p */ + +findx(p) +U *p; +{ + if (p == x) + return 1; + while (iscons(p)) { + if (findx(car(p))) + return 1; + p = cdr(p); + } + return 0; +} + +/* yes if anonymous function somewhere in p, i.e. p = (sum (f) 1) */ + +findf(p) +U *p; +{ + if (iscons(p) && car(p)->k == SYM && cdar(p) == nil) + return 1; + while (iscons(p)) { + if (findf(car(p))) + return 1; + p = cdr(p); + } + return 0; +} + +/* push constant expressions */ + +deconstruct(p) +U *p; +{ + int h; + U *q; + + if (does_not_depend_on_x(p)) { + push(p); + return; + } + + if (car(p) != _product) { + p = cdr(p); + while (iscons(p)) { + deconstruct(car(p)); + p = cdr(p); + } + return; + } + + q = cdr(p); + while (iscons(q)) { + if (depends_on_x(car(q))) + deconstruct(car(q)); + q = cdr(q); + } + + h = tos; + + p = cdr(p); + while (iscons(p)) { + if (does_not_depend_on_x(car(p))) + push(car(p)); + p = cdr(p); + } + + if (tos - h) + push(product(0, tos - h)); +} + +/* p form + + l list of constraints + + h Where constant expressions start on the stack + + Juggle the constant expressions to try and match f and p. +*/ + +match(p, l, h) +U *p, *l; +int h; +{ + int j1, j2; + U *q; + + for (j1 = h; j1 < tos; j1++) + for (j2 = h; j2 < tos; j2++) + { + a = stack[j1].p; + n = stack[j2].p; + + /* check constraints */ + + q = l; + while (iscons(q)) { + if (eval(car(q)) == nil) + break; + q = cdr(q); + } + + if (iscons(q)) + continue; /* constraint failed */ + + if (equal(f, eval(p))) + return 1; + } + + return 0; +} diff --git a/lisp/qed.lisp b/lisp/qed.lisp new file mode 100644 index 0000000..7fe3331 --- /dev/null +++ b/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 diff --git a/lisp/qed.tex b/lisp/qed.tex new file mode 100644 index 0000000..aa1afaa --- /dev/null +++ b/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 diff --git a/lisp/reference.html b/lisp/reference.html new file mode 100644 index 0000000..7113fb4 --- /dev/null +++ b/lisp/reference.html @@ -0,0 +1,627 @@ + + +

+

adjunct

+ +

+A adj A = det A I + +

+     > (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
+
+ +

+

and

+ +
+     > (and t t t)
+     t
+     > (and t nil t)
+     nil
+
+ +

+

append

+ +
+     > (append (list a b) (list c d))
+     (a b c d)
+
+ +

+

atom

+ +
+     > (atom a)
+     t
+     > (atom (list a b c))
+     nil
+
+ +

+

car

+ +
+     > (car (list a b c))
+     a
+
+ +

+

cdr

+ +
+     > (cdr (list a b c))
+     (b c)
+
+ +

+

complex-conjugate

+ +
+     > (complex-conjugate (sum a (product b i)))
+     (sum a (product -1 b i))
+
+ +

+

component

+ +
+     > (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
+
+ +

+

cond

+ +

+(cond (a1 b1) (a2 b2)... ) returns b of first non-nil a + +

+     > (define foo (cond ((zerop arg) "zero") (t "not zero")))
+     > (foo 0)
+     zero
+     > (foo 1)
+     not zero
+
+ +

+

cons

+ +
+     > (cons a b)
+     (a . b)
+
+ +

+

contract

+ +

+Form: contract12, contract13, contract23, etc. +For example (contract12 x) contracts, or sums over, +the first and second indices. + +

+     > (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)
+
+ +

+

cos

+ +
+     > (derivative (cos x) x)
+     (product -1 (sin x))
+
+ +

+

COS

+ +

+This function returns the exponential form of the cosine function. + +

+     > (printcars (COS x))
+     sum
+     (product 1/2 (power e (product -1 i x)))
+     (product 1/2 (power e (product i x)))
+
+ +

+

define

+ +
+     > (define foo (cons arg2 arg1))
+     > (foo a b)
+     (b . a)
+     > foo
+     (cons arg2 arg1)
+
+ +

+

derivative

+ +
+     > (derivative (power x 2) x)
+     (product 2 x)
+
+ +

+

determinant

+ +

+see adjunct + +

+

dot

+ +

+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. + +

+     > (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))
+
+ +

+The arguments are multiplied left to right, +i.e. (dot a b c) = (dot (dot a b) c). + +

+

e

+ +
+     > (derivative (power e x) x)
+     (power e x)
+
+ +

+

equal

+ +
+     > (equal a a)
+     t
+
+ +

+

eval

+ +
+     > (setq foo (quote (sum a a)))
+     > foo
+     (sum a a)
+     > (eval foo)
+     (product 2 a)
+
+ +

+

exit

+ +

+Exit simple lisp and return to the shell. +Control-C also does this. + +

+     > (exit)
+     %
+
+ +

+

exp

+ +
+     > (exp x)
+     (power e x)
+
+ +

+

gc

+ +

+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. + +

+     > (gc)
+     496053
+
+ +

+

goto

+ +
+     > (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
+
+ +

+

greaterp

+ +
+     > (greaterp 1 0)
+     t
+     > (greaterp 0 1)
+     nil
+
+ +

+

i

+ +

+The symbol "i" represents the square root of minus one. + +

+     > (product i i)
+     -1
+
+ +

+

integerp

+ +

+Returns t if the argument is an integer, nil otherwise. + +

+     > (integerp 2)
+     t
+     > (integerp 2/3)
+     nil
+
+ +

+

integral

+ +
+     > (integral (power x 2) x)
+     (product 1/3 (power x 3))
+
+ +

+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". + +

+

laguerre

+ +

+See the "hydrogen atom wave functions" example. + +

+

laplacian

+ +

+See the "hydrogen atom wave functions" example. + +

+

length

+ +

+Returns the number of cars in a list. + +

+     > (length (list a b c))
+     3
+
+ +

+

lessp

+ +
+     > (lessp 0 1)
+     t
+     > (lessp 1 0)
+     nil
+
+ +

+

list

+ +
+     > (list a b c)
+     (a b c)
+
+ +

+

not

+ +
+     > (not nil)
+     t
+
+ +

+

null

+ +
+     > (null nil)
+     t
+
+ +

+

numberp

+ +
+     > (numberp 1)
+     t
+
+ +

+

or

+ +
+     > (or nil t nil)
+     t
+     > (or nil nil nil)
+     nil
+
+ +

+

power

+ +
+     > (power a (sum b c))
+     (product (power a b) (power a c))
+
+ +

+

print

+ +
+     > (print "hello")
+     hello
+
+ +

+

printcars

+ +
+     > (printcars (list a b c))
+     a
+     b
+     c
+
+ +

+

product

+ +
+     > (product a b a)
+     (product b (power a 2))
+
+ +

+The "product" function expands products of sums. +When applied to tensors, the result is an outer product. + +

+     > (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))
+
+ +

+The more familiar inner product is obtained by contraction. + +

+     > (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))
+
+ +

+

prog

+ +

+Use "prog" when more than one expression must be evaluated in a function. + +

+     > (define foo (prog ()
+     (print arg1)
+     (print arg2)
+     ))
+     > (foo a b)
+     a
+     b
+
+ +

+

quote

+ +

+Quote means "don't evaluate." + +

+     > (quote (sum a a))
+     (sum a a)
+
+ +

+An apostrophe does the same thing. + +

+     > '(sum a a)
+     (sum a a)
+
+ +

+

return

+ +
+     > (define foo (prog () (return "hello")))
+     > (foo)
+     hello
+
+ +

+

run

+ +

+The "run" function evaluates the expressions in a file. + +

+     > (run "example1")
+     E=K+V for psi5?
+     t
+
+ +

+

setq

+ +
+     > (setq foo 1)
+     > foo
+     1
+
+ +

+

sin

+ +
+     > (derivative (sin x) x)
+     (cos x)
+
+ +

+

SIN

+ +

+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. + +

+     > (printcars (SIN x))
+     sum
+     (product -1/2 i (power e (product i x)))
+     (product 1/2 i (power e (product -1 i x)))
+
+ +
+     > (sum (power (COS x) 2) (power (SIN x) 2))
+     1
+
+ +

+

sum

+ +
+     > (sum a b a)
+     (sum b (product 2 a))
+
+ +

+

subst

+ +

+(subst a b c) means substitute a for b wherever it appears in c. + +

+     > (subst a b (list a b c))
+     (a a c)
+
+ +

+

tensor

+ +
+     > (product (tensor x) (tensor y))
+     (tensor x y)
+     > (product (tensor y) (tensor x))
+     (tensor y x)
+
+ +

+

transpose

+ +

+Form: transpose12, transpose13, transpose23, etc. +For example (transpose12 x) transposes the first and second indices. + +

+     > (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))
+
+ +

+

zerop

+ +
+     > (zerop 0)
+     t
+     > (zerop 1)
+     nil
+
diff --git a/lisp/robbins.lisp b/lisp/robbins.lisp new file mode 100644 index 0000000..eaa07e5 --- /dev/null +++ b/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) diff --git a/lisp/startup.lisp b/lisp/startup.lisp new file mode 100644 index 0000000..e503f80 --- /dev/null +++ b/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) +)) diff --git a/MainOSX.cpp b/src/MainOSX.cpp similarity index 100% rename from MainOSX.cpp rename to src/MainOSX.cpp diff --git a/MainXP.cpp b/src/MainXP.cpp similarity index 100% rename from MainXP.cpp rename to src/MainXP.cpp diff --git a/Makefile b/src/Makefile similarity index 100% rename from Makefile rename to src/Makefile diff --git a/abs.cpp b/src/abs.cpp similarity index 100% rename from abs.cpp rename to src/abs.cpp diff --git a/add.cpp b/src/add.cpp similarity index 100% rename from add.cpp rename to src/add.cpp diff --git a/adj.cpp b/src/adj.cpp similarity index 100% rename from adj.cpp rename to src/adj.cpp diff --git a/alloc.cpp b/src/alloc.cpp similarity index 100% rename from alloc.cpp rename to src/alloc.cpp diff --git a/append.cpp b/src/append.cpp similarity index 100% rename from append.cpp rename to src/append.cpp diff --git a/arccos.cpp b/src/arccos.cpp similarity index 100% rename from arccos.cpp rename to src/arccos.cpp diff --git a/arccosh.cpp b/src/arccosh.cpp similarity index 100% rename from arccosh.cpp rename to src/arccosh.cpp diff --git a/arcsin.cpp b/src/arcsin.cpp similarity index 100% rename from arcsin.cpp rename to src/arcsin.cpp diff --git a/arcsinh.cpp b/src/arcsinh.cpp similarity index 100% rename from arcsinh.cpp rename to src/arcsinh.cpp diff --git a/arctan.cpp b/src/arctan.cpp similarity index 100% rename from arctan.cpp rename to src/arctan.cpp diff --git a/arctanh.cpp b/src/arctanh.cpp similarity index 100% rename from arctanh.cpp rename to src/arctanh.cpp diff --git a/arg.cpp b/src/arg.cpp similarity index 100% rename from arg.cpp rename to src/arg.cpp diff --git a/atomize.cpp b/src/atomize.cpp similarity index 100% rename from atomize.cpp rename to src/atomize.cpp diff --git a/bake.cpp b/src/bake.cpp similarity index 100% rename from bake.cpp rename to src/bake.cpp diff --git a/besselj.cpp b/src/besselj.cpp similarity index 100% rename from besselj.cpp rename to src/besselj.cpp diff --git a/bessely.cpp b/src/bessely.cpp similarity index 100% rename from bessely.cpp rename to src/bessely.cpp diff --git a/bignum.cpp b/src/bignum.cpp similarity index 100% rename from bignum.cpp rename to src/bignum.cpp diff --git a/binomial.cpp b/src/binomial.cpp similarity index 100% rename from binomial.cpp rename to src/binomial.cpp diff --git a/ceiling.cpp b/src/ceiling.cpp similarity index 100% rename from ceiling.cpp rename to src/ceiling.cpp diff --git a/choose.cpp b/src/choose.cpp similarity index 100% rename from choose.cpp rename to src/choose.cpp diff --git a/circexp.cpp b/src/circexp.cpp similarity index 100% rename from circexp.cpp rename to src/circexp.cpp diff --git a/clear.cpp b/src/clear.cpp similarity index 100% rename from clear.cpp rename to src/clear.cpp diff --git a/clock.cpp b/src/clock.cpp similarity index 100% rename from clock.cpp rename to src/clock.cpp diff --git a/cmdisplay.cpp b/src/cmdisplay.cpp similarity index 100% rename from cmdisplay.cpp rename to src/cmdisplay.cpp diff --git a/coeff.cpp b/src/coeff.cpp similarity index 100% rename from coeff.cpp rename to src/coeff.cpp diff --git a/cofactor.cpp b/src/cofactor.cpp similarity index 100% rename from cofactor.cpp rename to src/cofactor.cpp diff --git a/condense.cpp b/src/condense.cpp similarity index 100% rename from condense.cpp rename to src/condense.cpp diff --git a/conj.cpp b/src/conj.cpp similarity index 100% rename from conj.cpp rename to src/conj.cpp diff --git a/cons.cpp b/src/cons.cpp similarity index 100% rename from cons.cpp rename to src/cons.cpp diff --git a/contract.cpp b/src/contract.cpp similarity index 100% rename from contract.cpp rename to src/contract.cpp diff --git a/cos.cpp b/src/cos.cpp similarity index 100% rename from cos.cpp rename to src/cos.cpp diff --git a/cosh.cpp b/src/cosh.cpp similarity index 100% rename from cosh.cpp rename to src/cosh.cpp diff --git a/data.cpp b/src/data.cpp similarity index 100% rename from data.cpp rename to src/data.cpp diff --git a/decomp.cpp b/src/decomp.cpp similarity index 100% rename from decomp.cpp rename to src/decomp.cpp diff --git a/define.cpp b/src/define.cpp similarity index 100% rename from define.cpp rename to src/define.cpp diff --git a/defint.cpp b/src/defint.cpp similarity index 100% rename from defint.cpp rename to src/defint.cpp diff --git a/defs.h b/src/defs.h similarity index 100% rename from defs.h rename to src/defs.h diff --git a/degree.cpp b/src/degree.cpp similarity index 100% rename from degree.cpp rename to src/degree.cpp diff --git a/denominator.cpp b/src/denominator.cpp similarity index 100% rename from denominator.cpp rename to src/denominator.cpp diff --git a/derivative.cpp b/src/derivative.cpp similarity index 100% rename from derivative.cpp rename to src/derivative.cpp diff --git a/det.cpp b/src/det.cpp similarity index 100% rename from det.cpp rename to src/det.cpp diff --git a/dirac.cpp b/src/dirac.cpp similarity index 100% rename from dirac.cpp rename to src/dirac.cpp diff --git a/display.cpp b/src/display.cpp similarity index 100% rename from display.cpp rename to src/display.cpp diff --git a/distill.cpp b/src/distill.cpp similarity index 100% rename from distill.cpp rename to src/distill.cpp diff --git a/divisors.cpp b/src/divisors.cpp similarity index 100% rename from divisors.cpp rename to src/divisors.cpp diff --git a/dpow.cpp b/src/dpow.cpp similarity index 100% rename from dpow.cpp rename to src/dpow.cpp diff --git a/draw.cpp b/src/draw.cpp similarity index 100% rename from draw.cpp rename to src/draw.cpp diff --git a/dsolve.cpp b/src/dsolve.cpp similarity index 100% rename from dsolve.cpp rename to src/dsolve.cpp diff --git a/eigen.cpp b/src/eigen.cpp similarity index 100% rename from eigen.cpp rename to src/eigen.cpp diff --git a/erf.cpp b/src/erf.cpp similarity index 100% rename from erf.cpp rename to src/erf.cpp diff --git a/erfc.cpp b/src/erfc.cpp similarity index 100% rename from erfc.cpp rename to src/erfc.cpp diff --git a/eval.cpp b/src/eval.cpp similarity index 100% rename from eval.cpp rename to src/eval.cpp diff --git a/expand.cpp b/src/expand.cpp similarity index 100% rename from expand.cpp rename to src/expand.cpp diff --git a/expcos.cpp b/src/expcos.cpp similarity index 100% rename from expcos.cpp rename to src/expcos.cpp diff --git a/expsin.cpp b/src/expsin.cpp similarity index 100% rename from expsin.cpp rename to src/expsin.cpp diff --git a/factor.cpp b/src/factor.cpp similarity index 100% rename from factor.cpp rename to src/factor.cpp diff --git a/factorial.cpp b/src/factorial.cpp similarity index 100% rename from factorial.cpp rename to src/factorial.cpp diff --git a/factorpoly.cpp b/src/factorpoly.cpp similarity index 100% rename from factorpoly.cpp rename to src/factorpoly.cpp diff --git a/factors.cpp b/src/factors.cpp similarity index 100% rename from factors.cpp rename to src/factors.cpp diff --git a/filter.cpp b/src/filter.cpp similarity index 100% rename from filter.cpp rename to src/filter.cpp diff --git a/find.cpp b/src/find.cpp similarity index 100% rename from find.cpp rename to src/find.cpp diff --git a/float.cpp b/src/float.cpp similarity index 100% rename from float.cpp rename to src/float.cpp diff --git a/floor.cpp b/src/floor.cpp similarity index 100% rename from floor.cpp rename to src/floor.cpp diff --git a/for.cpp b/src/for.cpp similarity index 100% rename from for.cpp rename to src/for.cpp diff --git a/gamma.cpp b/src/gamma.cpp similarity index 100% rename from gamma.cpp rename to src/gamma.cpp diff --git a/gcd.cpp b/src/gcd.cpp similarity index 100% rename from gcd.cpp rename to src/gcd.cpp diff --git a/guess.cpp b/src/guess.cpp similarity index 100% rename from guess.cpp rename to src/guess.cpp diff --git a/help.h b/src/help.h similarity index 100% rename from help.h rename to src/help.h diff --git a/help.html b/src/help.html similarity index 100% rename from help.html rename to src/help.html diff --git a/hermite.cpp b/src/hermite.cpp similarity index 100% rename from hermite.cpp rename to src/hermite.cpp diff --git a/hilbert.cpp b/src/hilbert.cpp similarity index 100% rename from hilbert.cpp rename to src/hilbert.cpp diff --git a/history.cpp b/src/history.cpp similarity index 100% rename from history.cpp rename to src/history.cpp diff --git a/html-tool.c b/src/html-tool.c similarity index 100% rename from html-tool.c rename to src/html-tool.c diff --git a/icon.ico b/src/icon.ico old mode 100755 new mode 100644 similarity index 100% rename from icon.ico rename to src/icon.ico diff --git a/imag.cpp b/src/imag.cpp similarity index 100% rename from imag.cpp rename to src/imag.cpp diff --git a/index.cpp b/src/index.cpp similarity index 100% rename from index.cpp rename to src/index.cpp diff --git a/init.cpp b/src/init.cpp similarity index 100% rename from init.cpp rename to src/init.cpp diff --git a/inner.cpp b/src/inner.cpp similarity index 100% rename from inner.cpp rename to src/inner.cpp diff --git a/integral.cpp b/src/integral.cpp similarity index 100% rename from integral.cpp rename to src/integral.cpp diff --git a/inv.cpp b/src/inv.cpp similarity index 100% rename from inv.cpp rename to src/inv.cpp diff --git a/is.cpp b/src/is.cpp similarity index 100% rename from is.cpp rename to src/is.cpp diff --git a/isprime.cpp b/src/isprime.cpp similarity index 100% rename from isprime.cpp rename to src/isprime.cpp diff --git a/itab.cpp b/src/itab.cpp similarity index 100% rename from itab.cpp rename to src/itab.cpp diff --git a/itest.cpp b/src/itest.cpp similarity index 100% rename from itest.cpp rename to src/itest.cpp diff --git a/laguerre.cpp b/src/laguerre.cpp similarity index 100% rename from laguerre.cpp rename to src/laguerre.cpp diff --git a/laplace.cpp b/src/laplace.cpp similarity index 100% rename from laplace.cpp rename to src/laplace.cpp diff --git a/lcm.cpp b/src/lcm.cpp similarity index 100% rename from lcm.cpp rename to src/lcm.cpp diff --git a/leading.cpp b/src/leading.cpp similarity index 100% rename from leading.cpp rename to src/leading.cpp diff --git a/legendre.cpp b/src/legendre.cpp similarity index 100% rename from legendre.cpp rename to src/legendre.cpp diff --git a/list.cpp b/src/list.cpp similarity index 100% rename from list.cpp rename to src/list.cpp diff --git a/log.cpp b/src/log.cpp similarity index 100% rename from log.cpp rename to src/log.cpp diff --git a/logo.gif b/src/logo.gif similarity index 100% rename from logo.gif rename to src/logo.gif diff --git a/madd.cpp b/src/madd.cpp similarity index 100% rename from madd.cpp rename to src/madd.cpp diff --git a/mag.cpp b/src/mag.cpp similarity index 100% rename from mag.cpp rename to src/mag.cpp diff --git a/main.cpp b/src/main.cpp similarity index 100% rename from main.cpp rename to src/main.cpp diff --git a/mcmp.cpp b/src/mcmp.cpp similarity index 100% rename from mcmp.cpp rename to src/mcmp.cpp diff --git a/mfactor.cpp b/src/mfactor.cpp similarity index 100% rename from mfactor.cpp rename to src/mfactor.cpp diff --git a/mgcd.cpp b/src/mgcd.cpp similarity index 100% rename from mgcd.cpp rename to src/mgcd.cpp diff --git a/mini-test.cpp b/src/mini-test.cpp similarity index 100% rename from mini-test.cpp rename to src/mini-test.cpp diff --git a/misc.cpp b/src/misc.cpp similarity index 100% rename from misc.cpp rename to src/misc.cpp diff --git a/mmodpow.cpp b/src/mmodpow.cpp similarity index 100% rename from mmodpow.cpp rename to src/mmodpow.cpp diff --git a/mmul.cpp b/src/mmul.cpp similarity index 100% rename from mmul.cpp rename to src/mmul.cpp diff --git a/mod.cpp b/src/mod.cpp similarity index 100% rename from mod.cpp rename to src/mod.cpp diff --git a/mpow.cpp b/src/mpow.cpp similarity index 100% rename from mpow.cpp rename to src/mpow.cpp diff --git a/mprime.cpp b/src/mprime.cpp similarity index 100% rename from mprime.cpp rename to src/mprime.cpp diff --git a/mroot.cpp b/src/mroot.cpp similarity index 100% rename from mroot.cpp rename to src/mroot.cpp diff --git a/mscan.cpp b/src/mscan.cpp similarity index 100% rename from mscan.cpp rename to src/mscan.cpp diff --git a/msqrt.cpp b/src/msqrt.cpp similarity index 100% rename from msqrt.cpp rename to src/msqrt.cpp diff --git a/mstr.cpp b/src/mstr.cpp similarity index 100% rename from mstr.cpp rename to src/mstr.cpp diff --git a/multiply.cpp b/src/multiply.cpp similarity index 100% rename from multiply.cpp rename to src/multiply.cpp diff --git a/nroots.cpp b/src/nroots.cpp similarity index 100% rename from nroots.cpp rename to src/nroots.cpp diff --git a/numerator.cpp b/src/numerator.cpp similarity index 100% rename from numerator.cpp rename to src/numerator.cpp diff --git a/outer.cpp b/src/outer.cpp similarity index 100% rename from outer.cpp rename to src/outer.cpp diff --git a/partition.cpp b/src/partition.cpp similarity index 100% rename from partition.cpp rename to src/partition.cpp diff --git a/polar.cpp b/src/polar.cpp similarity index 100% rename from polar.cpp rename to src/polar.cpp diff --git a/pollard.cpp b/src/pollard.cpp similarity index 100% rename from pollard.cpp rename to src/pollard.cpp diff --git a/power.cpp b/src/power.cpp similarity index 100% rename from power.cpp rename to src/power.cpp diff --git a/prime.cpp b/src/prime.cpp similarity index 100% rename from prime.cpp rename to src/prime.cpp diff --git a/primetab.cpp b/src/primetab.cpp similarity index 100% rename from primetab.cpp rename to src/primetab.cpp diff --git a/print.cpp b/src/print.cpp similarity index 100% rename from print.cpp rename to src/print.cpp diff --git a/product.cpp b/src/product.cpp similarity index 100% rename from product.cpp rename to src/product.cpp diff --git a/prototype-tool.c b/src/prototype-tool.c similarity index 100% rename from prototype-tool.c rename to src/prototype-tool.c diff --git a/prototypes.h b/src/prototypes.h similarity index 100% rename from prototypes.h rename to src/prototypes.h diff --git a/qadd.cpp b/src/qadd.cpp similarity index 100% rename from qadd.cpp rename to src/qadd.cpp diff --git a/qdiv.cpp b/src/qdiv.cpp similarity index 100% rename from qdiv.cpp rename to src/qdiv.cpp diff --git a/qmul.cpp b/src/qmul.cpp similarity index 100% rename from qmul.cpp rename to src/qmul.cpp diff --git a/qpow.cpp b/src/qpow.cpp similarity index 100% rename from qpow.cpp rename to src/qpow.cpp diff --git a/qsub.cpp b/src/qsub.cpp similarity index 100% rename from qsub.cpp rename to src/qsub.cpp diff --git a/quickfactor.cpp b/src/quickfactor.cpp similarity index 100% rename from quickfactor.cpp rename to src/quickfactor.cpp diff --git a/quotient.cpp b/src/quotient.cpp similarity index 100% rename from quotient.cpp rename to src/quotient.cpp diff --git a/rationalize.cpp b/src/rationalize.cpp similarity index 100% rename from rationalize.cpp rename to src/rationalize.cpp diff --git a/real.cpp b/src/real.cpp similarity index 100% rename from real.cpp rename to src/real.cpp diff --git a/rect.cpp b/src/rect.cpp similarity index 100% rename from rect.cpp rename to src/rect.cpp diff --git a/rewrite.cpp b/src/rewrite.cpp similarity index 100% rename from rewrite.cpp rename to src/rewrite.cpp diff --git a/roots.cpp b/src/roots.cpp similarity index 100% rename from roots.cpp rename to src/roots.cpp diff --git a/run.cpp b/src/run.cpp similarity index 100% rename from run.cpp rename to src/run.cpp diff --git a/scan.cpp b/src/scan.cpp similarity index 100% rename from scan.cpp rename to src/scan.cpp diff --git a/selftest.cpp b/src/selftest.cpp similarity index 100% rename from selftest.cpp rename to src/selftest.cpp diff --git a/selftest.h b/src/selftest.h similarity index 100% rename from selftest.h rename to src/selftest.h diff --git a/sgn.cpp b/src/sgn.cpp similarity index 100% rename from sgn.cpp rename to src/sgn.cpp diff --git a/simfac.cpp b/src/simfac.cpp similarity index 100% rename from simfac.cpp rename to src/simfac.cpp diff --git a/simplify.cpp b/src/simplify.cpp similarity index 100% rename from simplify.cpp rename to src/simplify.cpp diff --git a/sin.cpp b/src/sin.cpp similarity index 100% rename from sin.cpp rename to src/sin.cpp diff --git a/sinh.cpp b/src/sinh.cpp similarity index 100% rename from sinh.cpp rename to src/sinh.cpp diff --git a/stack.cpp b/src/stack.cpp similarity index 100% rename from stack.cpp rename to src/stack.cpp diff --git a/stdafx.h b/src/stdafx.h similarity index 100% rename from stdafx.h rename to src/stdafx.h diff --git a/subst.cpp b/src/subst.cpp similarity index 100% rename from subst.cpp rename to src/subst.cpp diff --git a/sum.cpp b/src/sum.cpp similarity index 100% rename from sum.cpp rename to src/sum.cpp diff --git a/symbol.cpp b/src/symbol.cpp similarity index 100% rename from symbol.cpp rename to src/symbol.cpp diff --git a/tan.cpp b/src/tan.cpp similarity index 100% rename from tan.cpp rename to src/tan.cpp diff --git a/tanh.cpp b/src/tanh.cpp similarity index 100% rename from tanh.cpp rename to src/tanh.cpp diff --git a/taylor.cpp b/src/taylor.cpp similarity index 100% rename from taylor.cpp rename to src/taylor.cpp diff --git a/tensor.cpp b/src/tensor.cpp similarity index 100% rename from tensor.cpp rename to src/tensor.cpp diff --git a/test-script.txt b/src/test-script.txt similarity index 100% rename from test-script.txt rename to src/test-script.txt diff --git a/test.cpp b/src/test.cpp similarity index 100% rename from test.cpp rename to src/test.cpp diff --git a/transform.cpp b/src/transform.cpp similarity index 100% rename from transform.cpp rename to src/transform.cpp diff --git a/transpose.cpp b/src/transpose.cpp similarity index 100% rename from transpose.cpp rename to src/transpose.cpp diff --git a/userfunc.cpp b/src/userfunc.cpp similarity index 100% rename from userfunc.cpp rename to src/userfunc.cpp diff --git a/variables.cpp b/src/variables.cpp similarity index 100% rename from variables.cpp rename to src/variables.cpp diff --git a/vectorize.cpp b/src/vectorize.cpp similarity index 100% rename from vectorize.cpp rename to src/vectorize.cpp diff --git a/window.cpp b/src/window.cpp similarity index 100% rename from window.cpp rename to src/window.cpp diff --git a/zero.cpp b/src/zero.cpp similarity index 100% rename from zero.cpp rename to src/zero.cpp