diff --git a/gsll.asd b/gsll.asd
index 1fc4d0320a38cfc196820a6dd2e68a01aed7705b..6fe0ffacb514ea1fac1297e52f58febe0f84388f 100644
--- a/gsll.asd
+++ b/gsll.asd
@@ -1,6 +1,6 @@
 ;; Definition of GSLL system 
 ;; Liam Healy
-;; Time-stamp: <2008-08-10 17:01:08EDT gsll.asd>
+;; Time-stamp: <2008-08-10 18:27:45EDT gsll.asd>
 ;; $Id$
 
 (asdf:defsystem "gsll"
@@ -82,6 +82,7 @@
 	    ((:file "blas1")
 	     (:file "blas2")
 	     (:file "blas3" :depends-on (blas2))
+	     (:file "exponential")
 	     ;(:file "lu")
 	     ;(:file "qr")
 	     ;(:file "qrpt")
diff --git a/linear-algebra/exponential.lisp b/linear-algebra/exponential.lisp
new file mode 100644
index 0000000000000000000000000000000000000000..ecd54695ea27b38d208114e9d9849162f2282796
--- /dev/null
+++ b/linear-algebra/exponential.lisp
@@ -0,0 +1,33 @@
+;; Exponential of a matrix
+;; Liam Healy 2008-08-10 17:25:35EDT exponential.lisp
+;; Time-stamp: <2008-08-10 18:27:04EDT exponential.lisp>
+;; $Id: $
+
+(in-package :gsl)
+
+;;; No documentation
+(defmfun matrix-exponential (matrix exponential &optional (mode :double))
+  "gsl_linalg_exponential_ss"
+  (((mpointer matrix) :pointer) ((mpointer exponential) :pointer)
+   (mode sf-mode))
+  :inputs (matrix)
+  :outputs (exponential)
+  :return (exponential)
+  :documentation
+  "Calculate the matrix exponential by the scaling and
+   squaring method described in Moler + Van Loan,
+   SIAM Rev 20, 801 (1978).  The matrix to be exponentiated
+   is matrix, the returned exponential is exponential.
+   The mode argument allows
+   choosing an optimal strategy, from the table
+   given in the paper, for a given precision.")
+
+#|
+;;; A matrix of the form [[0, 1], [-1, 0]]
+;;; when exponentiated gives [[cos x,  sin x], [-sin x,  cos x]]
+(letm ((mat (matrix-double-float (a (0.0d0 1.0d0) (-1.0d0 0.0d0))))
+       (exp (matrix-double-float '(2 2))))
+  (cl-array (matrix-exponential mat exp)))
+#2A((0.5403023058681384d0 0.841470984807895d0)
+    (-0.841470984807895d0 0.5403023058681384d0))
+|#