new file mode 100644
--- /dev/null
+++ b/scripts/elfun/gcd.m
@@ -0,0 +1,38 @@
+function [g, v] = gcd (a, ...)
+  
+  # [g [, v]] = gcd (a) returns the greatest common divisor g of the
+  # entries of the integer vector a, and an integer vector v such that
+  # g = v(1) * a(k) + ... + v(k) * a(k).
+  #
+  # [g [, v]] = gcd (a1, ..., ak) is the same with a = [a1, ..., ak].
+  
+  # Written by KH (Kurt.Hornik@ci.tuwien.ac.at) on Sep 16, 1994
+  # Copyright Dept of Statistics and Probability Theory TU Wien
+
+  if (nargin > 1)
+    va_start;
+    for k=2:nargin;
+      a = [a, va_arg()];
+    endfor
+  endif
+  
+  if (round(a) != a)
+    error("gcd:  all arguments must be integer");
+  endif
+  
+  g = abs(a(1));
+  v = sign(a(1));
+  for k=1:(length(a)-1)
+    x = [g, 1, 0];
+    y = [abs(a(k+1)), 0, 1];
+    while (y(1) > 0)
+      r = x - y * floor(x(1) / y(1));
+      x = y;
+      y = r;
+    endwhile
+    g = x(1);
+    v = [x(2) * v, x(3) * sign(a(k+1))];
+  endfor
+    
+endfunction
+    
\ No newline at end of file

new file mode 100644
--- /dev/null
+++ b/scripts/elfun/lcm.m
@@ -0,0 +1,32 @@
+function l = lcm (a, ...)
+  
+  # lcm (a) returns the least common multiple of the entries of the
+  # integer vector a.
+  # lcm (a1, ..., ak) is the same as lcm([a1, ..., ak]).
+  
+  # Written by KH (Kurt.Hornik@ci.tuwien.ac.at) on Sep 16, 1994
+  # Copyright Dept of Statistics and Probability Theory TU Wien
+
+  if (nargin > 1)
+    va_start;
+    for k=2:nargin;
+      a = [a, va_arg()];
+    endfor
+  endif
+  
+  if (round(a) != a)
+    error("lcm:  all arguments must be integer");
+  endif
+  
+  if (any(a) == 0)
+    l = 0;
+  else
+    a = abs(a);
+    l = a(1);
+    for k=1:(length(a)-1)
+      l = l * a(k+1) / gcd(l, a(k+1));
+    endfor
+  endif
+    
+endfunction
+    
\ No newline at end of file

new file mode 100644
--- /dev/null
+++ b/scripts/specfun/Makefile.in
@@ -0,0 +1,74 @@
+#
+# Makefile for octave's scripts/specfun directory
+#
+# John W. Eaton
+# jwe@che.utexas.edu
+# Department of Chemical Engineering
+# The University of Texas at Austin
+
+TOPDIR = ../..
+
+script_sub_dir = specfun
+
+srcdir = @srcdir@
+top_srcdir = @top_srcdir@
+VPATH = @srcdir@
+
+include $(TOPDIR)/Makeconf
+
+INSTALL = @INSTALL@
+INSTALL_PROGRAM = @INSTALL_PROGRAM@
+INSTALL_DATA = @INSTALL_DATA@
+
+SOURCES = *.m
+
+DISTFILES = Makefile.in $(SOURCES)
+
+FCN_FILES = $(wildcard $(srcdir)/*.m)
+FCN_FILES_NO_DIR = $(notdir $(FCN_FILES))
+
+all:
+.PHONY: all
+
+install: all
+	if test -d $(fcnfiledir)/$(script_sub_dir) ; then true ; \
+	else $(TOPDIR)/mkpath $(fcnfiledir)/$(script_sub_dir) ; fi
+	for f in $(FCN_FILES_NO_DIR) ; do \
+	  rm -f $(fcnfiledir)/$(script_sub_dir)/$$f ; \
+	  $(INSTALL_DATA) $(srcdir)/$$f $(fcnfiledir)/$(script_sub_dir)/$$f ; \
+	done
+.PHONY: install
+
+uninstall:
+	for f in $(FCN_FILES_NO_DIR) ; \
+	  do rm -f $(fcnfiledir)/$(script_sub_dir)/$$f ; \
+	done
+.PHONY: uninstall
+
+clean:
+.PHONY: clean
+
+tags: $(SOURCES)
+	ctags $(SOURCES)
+
+TAGS: $(SOURCES)
+	etags $(SOURCES)
+
+mostlyclean: clean
+.PHONY: mostlyclean
+
+distclean: clean
+	rm -f Makefile
+.PHONY: distclean
+
+realclean: distclean
+	rm -f tags TAGS
+.PHONY: realclean
+
+local-dist:
+	ln $(DISTFILES) ../../`cat ../../.fname`/scripts/specfun
+.PHONY: local-dist
+
+dist:
+	ln $(DISTFILES) ../../`cat ../../.fname`/scripts/specfun
+.PHONY: dist

new file mode 100644
--- /dev/null
+++ b/scripts/specfun/beta.m
@@ -0,0 +1,15 @@
+function retval = beta(a, b)
+  
+  # usage:  beta(a, b)
+  #
+  # Returns the beta function beta(a,b) = gamma(a) * gamma(b) / gamma(a+b)
+  # of a and b.
+
+  # Written by KH (Kurt.Hornik@ci.tuwien.ac.at) on Jun 13, 1993
+  # Updated by KH (Kurt.Hornik@ci.tuwien.ac.at) on Aug 13, 1994
+  # Copyright Dept of Probability Theory and Statistics TU Wien
+
+  retval = exp(lgamma(a) + lgamma(b) - lgamma(a+b));
+
+endfunction
+

new file mode 100644
--- /dev/null
+++ b/scripts/specfun/betai.m
@@ -0,0 +1,107 @@
+function y = betai(a, b, x)
+  
+  # usage:  betai(a, b, x)
+  #
+  # Returns the incomplete beta function
+  #   betai (a, b, x) = BETA(a,b)^(-1) INT_0^x t^(a-1) (1-t)^(b-1) dt.
+  # If x has more than one component, both a and b must be scalars.
+  # If x is a scalar, a and b must be of compatible dimensions.
+  
+  # Written by KH (Kurt.Hornik@ci.tuwien.ac.at) on Aug 2, 1994.
+  # Copyright Dept of Probability Theory and Statistics TU Wien.
+
+  # Computation is based on the series expansion
+  #   betai(a, b, x) 
+  #     = \frac{1}{B(a, b)} x^a 
+  #         \sum_{k=0}^\infty \frac{(1-b)\cdots(k-b)}{a+k} \frac{x^k}{k!}
+  # for x <= 1/2.  For x > 1/2, betai(a, b, x) = 1 - betai(b, a, 1-x).
+
+  if (nargin <> 3)
+    error("usage:  betai (a, b, x)");
+  endif
+
+  if !((a > 0) && (b > 0))
+    error("betai:  a and b must both be positive.");
+  endif
+  [nr, nc] = size(x);
+  if (min([nr nc]) == 0)
+    error("betai:  x must not be empty.");
+  endif
+  if (any(x < 0) || any(x > 1))
+    error("betai:  all entries of x must be in [0,1].");
+  endif
+
+  if ((nr > 1) || (nc > 1))
+    
+    if !(is_scalar(a) && is_scalar(b))
+      error("betai:  if x is not a scalar, a and b must be scalars.");
+    endif
+
+    n   = nr * nc;
+    x   = reshape(x, 1, n);
+    y   = zeros(1, n);
+    c   = exp(lgamma(a+b) - lgamma(a) - lgamma(b));
+
+    y(find(x == 1)) = ones(1, sum(x == 1));
+
+    # now do the series computation.  The error when truncating at term K
+    # is always less than 2^(-K), hence the following choice of K.
+    K   = ceil(-log(eps) / log(2));
+    k   = (1:K)';
+
+    ind = find((x > 0) & (x <= 1/2));
+    len = length(ind);
+    if (len > 0)
+      tmp    = cumprod((1 - b./k) * x(ind)) ./ ((a + k) * ones(1, len));
+      y(ind) = c * exp(a * log(x(ind))) .* (1/a + sum(tmp));
+    endif
+    
+    ind = find((x > 1/2) & (x < 1));
+    len = length(ind);
+    if (len > 0)
+      tmp    = cumprod((1 - a./k) * (1 - x(ind))) ./ ((b + k) * ones(1, len));
+      y(ind) = 1 - c * exp(b * log(1-x(ind))) .* (1/b + sum(tmp));
+    endif
+  
+    y = reshape(y, nr, nc);
+    
+  else
+    
+    [ra, ca] = size(a);
+    [rb, cb] = size(b);
+    if !((ra == rb) && (ca == cb))
+      error("betai:  a and b must have the same size.");
+    endif
+    
+    n   = ra * ca;
+    a   = reshape(a, 1, n);
+    b   = reshape(b, 1, n);
+    c   = exp(lgamma(a+b) - lgamma(a) - lgamma(b));
+    
+    if (x == 0)
+      y   = zeros(1, n);
+    elseif (x == 1)
+      y   = ones(1, n);
+    else
+      K   = ceil(-log(eps) / log(2));
+      k   = (1:K)' * ones(1, n);
+      h   = ones(K, 1);
+      if ((x > 0) && (x <= 1/2))
+	tmp = cumprod((1 - (h * b) ./ k) * x) ./ ((h * a) + k);
+	y   = c .* exp(a * log(x)) .* (1 ./ a + sum(tmp));
+      else
+	tmp = cumprod((1 - (h * a) ./ k) * (1-x)) ./ ((h * b) + k);
+	y   = 1 - c .* exp(b * log(1-x)) .* (1 ./ b + sum(tmp));
+      endif
+    endif
+  
+    y = reshape(y, ra, ca);
+    
+  endif
+
+endfunction
+
+
+
+
+

new file mode 100644
--- /dev/null
+++ b/scripts/specfun/gammai.m
@@ -0,0 +1,102 @@
+function y = gammai(a, x)
+  
+  # usage:  gammai(a, x)
+  #
+  # Computes the incomplete gamma function
+  #    gammai(a, x) 
+  #      = (integral from 0 to x of exp(-t) t^(a-1) dt) / gamma(a).
+  # If a is scalar, then gammai(a, x) is returned for each element of x
+  # and vice versa.
+  # If neither a nor x is scalar, the sizes of a and x must agree, and
+  # gammai is applied pointwise.
+  
+  # Written by KH (Kurt.Hornik@ci.tuwien.ac.at) on Aug 13, 1994
+  # Copyright Dept of Probability Theory and Statistics TU Wien
+
+  if (nargin != 2)
+    error("usage:  gammai(a, x)");
+  endif
+  
+  [r_a, c_a] = size(a);
+  [r_x, c_x] = size(x);
+  e_a = r_a * c_a;
+  e_x = r_x * c_x;
+  
+  # The following code is rather ugly.  We want the function to work
+  # whenever a and x have the same size or a or x is scalar.  
+  # We do this by reducing the latter cases to the former.
+  
+  if ((e_a == 0) || (e_x == 0))
+    error("gammai:  both a and x must be nonempty");
+  endif
+  if ((r_a == r_x) && (c_a == c_x))
+    n   = e_a;
+    a   = reshape(a, 1, n);
+    x   = reshape(x, 1, n);
+    r_y = r_a;
+    c_y = c_a;
+  elseif (e_a == 1)
+    n   = e_x;
+    a   = a * ones(1, n);
+    x   = reshape(x, 1, n);
+    r_y = r_x;
+    c_y = c_x;
+  elseif (e_x == 1)
+    n   = e_a;
+    a   = reshape(a, 1, n);
+    x   = x * ones(1, n);
+    r_y = r_a;
+    c_y = c_a;
+  else
+    error("gammai:  a and x must have the same size if neither is scalar"); 
+  endif
+
+  # Now we can do sanity checking ...
+  
+  if (any(a <= 0) || any(a == Inf))
+    error("gammai:  all entries of a must be positive anf finite");
+  endif
+  if any(x < 0)
+    error("gammai:  all entries of x must be nonnegative");
+  endif
+  
+  y = zeros(1, n);
+
+  # For x < a + 1, use summation.  The below choice of k should ensure
+  # that the overall error is less than eps ... 
+  S = find((x > 0) & (x < a + 1));
+  s = length(S);
+  if (s > 0)
+    k   = ceil(- max([a(S), x(S)]) * log(eps));
+    K   = (1:k)';
+    M   = ones(k, 1);
+    A   = cumprod((M * x(S)) ./ (M * a(S) + K * ones(1, s)));
+    y(S) = exp(-x(S) + a(S) .* log(x(S))) .* (1 + sum(A)) ./ gamma(a(S)+1);
+  endif
+
+  # For x >= a + 1, use the continued fraction.
+  # Note, however, that this converges MUCH slower than the series
+  # expansion for small a and x not too large!
+  S = find((x >= a + 1) & (x < Inf));
+  s = length(S);
+  if (s > 0)
+    u   = [zeros(1, s); ones(1, s)];
+    v   = [ones(1, s); x(S)];
+    c_old = 0;
+    c_new = v(1,:) ./ v(2,:);
+    n   = 1;
+    while (max(abs(c_old ./ c_new - 1)) > 10 * eps)
+      c_old = c_new;
+      u = v + u .* (ones(2, 1) * (n - a(S)));
+      v = u .* (ones(2, 1) * x(S)) + n * v;
+      c_new = v(1,:) ./ v(2,:);
+      n = n + 1;
+    endwhile
+    y(S) = 1 - exp(-x(S) + a(S) .* log(x(S))) .* c_new ./ gamma(a(S));
+  endif
+  
+  y(find(x == Inf)) = ones(1, sum(x == Inf));
+  
+  y = reshape(y, r_y, c_y);
+
+endfunction
\ No newline at end of file