diff --git a/example2 b/example2
index 2960fab..3089a97 100644
--- a/example2
+++ b/example2
@@ -1,52 +1,32 @@
 ; Page references are for the book "Gravitation."
-
 ; generic metric
-
 (setq gdd (sum
   (product (g00) (tensor x0 x0))
   (product (g11) (tensor x1 x1))
   (product (g22) (tensor x2 x2))
   (product (g33) (tensor x3 x3))
 ))
-
 (gr) ; compute g, guu, GAMUDD, RUDDD, RDD, R, GDD, GUD and GUU
-
 ; generic vectors
-
 (setq u (sum
   (product (u0) (tensor x0))
   (product (u1) (tensor x1))
   (product (u2) (tensor x2))
   (product (u3) (tensor x3))
 ))
-
 (setq v (sum
   (product (v0) (tensor x0))
   (product (v1) (tensor x1))
   (product (v2) (tensor x2))
   (product (v3) (tensor x3))
 ))
-
 (setq w (sum
   (product (w0) (tensor x0))
   (product (w1) (tensor x1))
   (product (w2) (tensor x2))
   (product (w3) (tensor x3))
 ))
-
-<hr>
-
-$$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}
-)$$
-
 ; how to antisymmetrize three indices (p. 86)
-
 (setq V (sum
   (product 1/6 V)
   (product 1/6 (transpose12 (transpose23 V))) ; nu lambda mu -> mu nu lambda
@@ -55,97 +35,42 @@ $$V_{[\mu\nu\lambda]}={1\over3!}
   (product -1/6 (transpose23 V))              ; mu lambda nu -> mu nu lambda
   (product -1/6 (transpose13 V))              ; lambda nu mu -> mu nu lambda
 ))
-
-<hr>
-
-$$[{\bf u},{\bf v}]=
-(u^\beta{v^\alpha}_{,\beta}-v^\beta{u^\alpha}_{,\beta}){\bf e}_\alpha$$
-
 ; commutator (p. 206)
-
 (define commutator (sum
   (contract13 (product +1 arg1 (gradient arg2)))
   (contract13 (product -1 arg2 (gradient arg1)))
 ))
-
-<hr>
-
-$${\Gamma^\alpha}_{\mu\nu}=\hbox{$1\over2$}g^{\alpha\beta}
-(g_{\beta\mu,\nu}+g_{\beta\nu,\mu}-g_{\mu\nu,\beta})$$
-
 (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))
-
-<hr>
-
-$${T^\alpha}_{;\gamma}={T^\alpha}_{,\gamma}+
-T^\mu{\Gamma^\alpha}_{\mu\gamma}$$
-
 ; covariant derivative of a vector (p. 211)
-
 (define covariant-derivative (sum
   (gradient arg)
   (contract13 (product arg GAMUDD))
 ))
-
-<hr>
-
-$${G^{\mu\nu}}_{;\nu}=0$$
-
 (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))
-
-<hr>
-
-$${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}$$
-
 (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))
-
-<hr>
-
-$$T_{\alpha\beta;\gamma}=T_{\alpha\beta,\gamma}
--T_{\mu\beta}{\Gamma^\mu}_{\alpha\gamma}
--T_{\alpha\mu}{\Gamma^\mu}_{\beta\gamma}$$
-
 ; p. 259
-
 (define covariant-derivative-of-down-down (prog (temp)
   (setq temp (product arg GAMUDD))
   (return (sum
@@ -154,13 +79,6 @@ $$T_{\alpha\beta;\gamma}=T_{\alpha\beta,\gamma}
     (product -1 (contract23 temp))
   ))
 ))
-
-<hr>
-
-$${T^\alpha}_{\beta;\gamma}={T^\alpha}_{\beta,\gamma}+
-{T^\mu}_\beta{\Gamma^\alpha}_{\mu\gamma}-
-{T^\alpha}_\mu{\Gamma^\mu}_{\beta\gamma}$$
-
 (define covariant-derivative-of-up-down (prog (temp)
   (setq temp (product arg GAMUDD))
   (return (sum
@@ -169,81 +87,37 @@ $${T^\alpha}_{\beta;\gamma}={T^\alpha}_{\beta,\gamma}+
     (product -1 (contract23 temp))
   ))
 ))
-
-<hr>
-
-$$\nabla_{\bf u}{\bf v}={v^\alpha}_{;\beta}u^{\beta}$$
-
 (define directed-covariant-derivative
   (contract23 (product (covariant-derivative arg1) arg2))
 )
-
-<hr>
-
-$$\nabla_{\bf u}{\bf v}-\nabla_{\bf v}{\bf u}=[{\bf u},{\bf v}]$$
-
 (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)
-
-<hr>
-
-$$\nabla_{\bf u}(f{\bf v})=
-f\nabla_{\bf u}{\bf v}+{\bf v}\partial_{\bf u}f$$
-
 (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))
-
-<hr>
-
-$$\nabla_{{\bf v}+{\bf w}}{\bf u}=
-\nabla_{\bf v}{\bf u}+\nabla_{\bf w}{\bf u}$$
-
 (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))
-
-<hr>
-
-$${R^\alpha}_{\beta\gamma\delta}={R^\alpha}_{\beta[\gamma\delta]}$$
-
 (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))
-
-<hr>
-
-$${R^\alpha}_{[\beta\gamma\delta]}=0$$
-
 (print "riemann vanishes when antisymmetrized on last three indices (p. 286)")
-
 (setq temp (sum
   (product 1/6 RUDDD)
   (product 1/6 (transpose34 (transpose24 RUDDD)))
@@ -252,92 +126,51 @@ $${R^\alpha}_{[\beta\gamma\delta]}=0$$
   (product -1/6 (transpose34 RUDDD))
   (product -1/6 (transpose24 RUDDD))
 ))
-
 (print (equal temp 0))
-
-<hr>
-
-$${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}$$
-
 (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
 ))
-
-<hr>
-
-$${B^\mu}_{;\alpha\beta}-{B^\mu}_{;\beta\alpha}=
-{R^\mu}_{\nu\beta\alpha}B^\nu$$
-
 (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))
-
-<hr>
-
 (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
-
-$$g_{uu}=Ve^{2\beta}/r-U^2r^2e^{2\gamma}$$
-
 (setq g_uu (sum
   (product
     (V u r theta)
@@ -351,39 +184,25 @@ $$g_{uu}=Ve^{2\beta}/r-U^2r^2e^{2\gamma}$$
     (power e (product 2 (gamma u r theta)))
   )
 ))
-
-$$g_{ur}=2e^{2\beta}$$
-
 (setq g_ur (product 2 (power e (product 2 (beta u r theta)))))
-
-$$g_{u\theta}=2Ur^2e^{2\gamma}$$
-
 (setq g_utheta (product
   2
   (U u r theta)
   (power r 2)
   (power e (product 2 (gamma u r theta)))
 ))
-
-$$g_{\theta\theta}=-r^2e^{2\gamma}$$
-
 (setq g_thetatheta (product
   -1
   (power r 2)
   (power e (product 2 (gamma u r theta)))
 ))
-
-$$g_{\phi\phi}=-r^2e^{-2\gamma}\sin^2\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))
@@ -393,15 +212,9 @@ $$g_{\phi\phi}=-r^2e^{-2\gamma}\sin^2\theta$$
   (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/qed.lisp b/qed.lisp
new file mode 100644
index 0000000..7fe3331
--- /dev/null
+++ b/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