# HG changeset patch
# User jwe
# Date 789859021 0
# Node ID 9fc405c8c06c378cf718aff2b3c41924218f0cde
# Parent  f558749713f1fb9fceccbbe38578d3ed61e8b52e
[project @ 1995-01-11 21:17:01 by jwe]

diff --git a/scripts/elfun/gcd.m b/scripts/elfun/gcd.m
--- a/scripts/elfun/gcd.m
+++ b/scripts/elfun/gcd.m
@@ -1,13 +1,32 @@
+# Copyright (C) 1995 John W. Eaton
+# 
+# This file is part of Octave.
+# 
+# Octave is free software; you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published by the
+# Free Software Foundation; either version 2, or (at your option) any
+# later version.
+# 
+# Octave is distributed in the hope that it will be useful, but WITHOUT
+# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+# for more details.
+# 
+# You should have received a copy of the GNU General Public License
+# along with Octave; see the file COPYING.  If not, write to the Free
+# Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+
 function [g, v] = gcd (a, ...)
-  
+
+# usage: 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
+# Written by KH (Kurt.Hornik@ci.tuwien.ac.at) on Sep 16, 1994.
 
   if (nargin > 1)
     va_start;
@@ -17,7 +36,7 @@
   endif
   
   if (round (a) != a)
-    error ("gcd:  all arguments must be integer");
+    error ("gcd: all arguments must be integer");
   endif
   
   g = abs (a(1));
diff --git a/scripts/elfun/lcm.m b/scripts/elfun/lcm.m
--- a/scripts/elfun/lcm.m
+++ b/scripts/elfun/lcm.m
@@ -1,11 +1,30 @@
+# Copyright (C) 1995 John W. Eaton
+# 
+# This file is part of Octave.
+# 
+# Octave is free software; you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published by the
+# Free Software Foundation; either version 2, or (at your option) any
+# later version.
+# 
+# Octave is distributed in the hope that it will be useful, but WITHOUT
+# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+# for more details.
+# 
+# You should have received a copy of the GNU General Public License
+# along with Octave; see the file COPYING.  If not, write to the Free
+# Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+
 function l = lcm (a, ...)
 
+# usage: 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
+# Written by KH (Kurt.Hornik@ci.tuwien.ac.at) on Sep 16, 1994.
 
   if (nargin > 1)
     va_start;
@@ -24,7 +43,7 @@
     a = abs (a);
     l = a (1);
     for k = 1:(length (a) - 1)
-      l = l * a(k+1) / gcd(l, a(k+1));
+      l = l * a(k+1) / gcd (l, a(k+1));
     endfor
   endif
     
diff --git a/scripts/linear-algebra/null.m b/scripts/linear-algebra/null.m
--- a/scripts/linear-algebra/null.m
+++ b/scripts/linear-algebra/null.m
@@ -1,3 +1,21 @@
+# Copyright (C) 1995 John W. Eaton
+# 
+# This file is part of Octave.
+# 
+# Octave is free software; you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published by the
+# Free Software Foundation; either version 2, or (at your option) any
+# later version.
+# 
+# Octave is distributed in the hope that it will be useful, but WITHOUT
+# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+# for more details.
+# 
+# You should have received a copy of the GNU General Public License
+# along with Octave; see the file COPYING.  If not, write to the Free
+# Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+
 function retval = null (A, tol)
 
 # usage: null (A, tol)
@@ -10,8 +28,7 @@
 # max (size (A)) * sigma_max (A) * eps, where sigma_max (A) is the
 # maximal singular value of A. 
 
-# written by KH (Kurt.Hornik@ci.tuwien.ac.at) on Dec 24, 1993
-# copyright Dept of Probability Theory and Statistics TU Wien
+# Written by KH (Kurt.Hornik@ci.tuwien.ac.at) on Dec 24, 1993.
 
   [U, S, V] = svd (A);
 
@@ -22,7 +39,7 @@
   if (nargin == 1)
     tol = max (size (A)) * s (1) * eps;
   else (nargin != 2)
-    usage ("null(A [, tol])"); 
+    usage ("null (A [, tol])"); 
   endif
 
   rank = sum (s > tol);
diff --git a/scripts/linear-algebra/orth.m b/scripts/linear-algebra/orth.m
--- a/scripts/linear-algebra/orth.m
+++ b/scripts/linear-algebra/orth.m
@@ -1,3 +1,21 @@
+# Copyright (C) 1995 John W. Eaton
+# 
+# This file is part of Octave.
+# 
+# Octave is free software; you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published by the
+# Free Software Foundation; either version 2, or (at your option) any
+# later version.
+# 
+# Octave is distributed in the hope that it will be useful, but WITHOUT
+# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+# for more details.
+# 
+# You should have received a copy of the GNU General Public License
+# along with Octave; see the file COPYING.  If not, write to the Free
+# Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+
 function retval = orth (A, tol)
 
 # usage: orth (A, tol)
@@ -10,8 +28,7 @@
 # max (size (A)) * sigma_max (A) * eps, where sigma_max (A) is the
 # maximal singular value of A.
 
-# written by KH (Kurt.Hornik@ci.tuwien.ac.at) on Dec 24, 1993
-# copyright Dept of Probability Theory and Statistics TU Wien
+# Written by KH (Kurt.Hornik@ci.tuwien.ac.at) on Dec 24, 1993.
 
   [U, S, V] = svd (A);
 
diff --git a/scripts/signal/fftconv.m b/scripts/signal/fftconv.m
--- a/scripts/signal/fftconv.m
+++ b/scripts/signal/fftconv.m
@@ -1,6 +1,24 @@
+# Copyright (C) 1995 John W. Eaton
+# 
+# This file is part of Octave.
+# 
+# Octave is free software; you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published by the
+# Free Software Foundation; either version 2, or (at your option) any
+# later version.
+# 
+# Octave is distributed in the hope that it will be useful, but WITHOUT
+# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+# for more details.
+# 
+# You should have received a copy of the GNU General Public License
+# along with Octave; see the file COPYING.  If not, write to the Free
+# Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+
 function c = fftconv (a, b, N)
 
-# usage:  fftconv (a, b [, N])
+# usage: fftconv (a, b [, N])
 #
 # c = fftconv (a, b) returns the convolution of the vectors a and b,
 # a vector with length equal to length (a) + length (b) - 1.  
@@ -10,11 +28,10 @@
 # The computation uses the FFT by calling fftfilt.  If the optional
 # argument N is specified, an N-point FFT is used.
 
-# Written by KH (Kurt.Hornik@ci.tuwien.ac.at) on Sep 3, 1994
-# Copyright Dept of Statistics and Probability Theory TU Wien
+# Written by KH (Kurt.Hornik@ci.tuwien.ac.at) on Sep 3, 1994.
   
   if (nargin < 2 || nargin > 3)
-    usage (" fftconv (b, x [, N])");
+    usage ("fftconv (b, x [, N])");
   endif
   
   if (is_matrix (a) || is_matrix (b))
@@ -32,7 +49,7 @@
       c = fftfilt (a, b);
     else
       if !(is_scalar (N))
-	error ("fftconv:  N has to be a scalar");
+	error ("fftconv: N has to be a scalar");
       endif
       c = fftfilt (a, b, N);
     endif
diff --git a/scripts/signal/fftfilt.m b/scripts/signal/fftfilt.m
--- a/scripts/signal/fftfilt.m
+++ b/scripts/signal/fftfilt.m
@@ -1,3 +1,21 @@
+# Copyright (C) 1995 John W. Eaton
+# 
+# This file is part of Octave.
+# 
+# Octave is free software; you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published by the
+# Free Software Foundation; either version 2, or (at your option) any
+# later version.
+# 
+# Octave is distributed in the hope that it will be useful, but WITHOUT
+# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+# for more details.
+# 
+# You should have received a copy of the GNU General Public License
+# along with Octave; see the file COPYING.  If not, write to the Free
+# Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+
 function y = fftfilt (b, x, N)
 
 # usage:  fftfilt (b, x [, N])
@@ -7,7 +25,6 @@
 # b using an N-point FFT.
 
 # Written by KH (Kurt.Hornik@ci.tuwien.ac.at) on Sep 3, 1994
-# Copyright Dept of Statistics and Probability Theory TU Wien
 
 # Reference:  Oppenheim & Schafer (1989).  Discrete-time Signal
 # Processing (Chapter 8).  Prentice-Hall.
@@ -24,8 +41,8 @@
   
   [r_x, c_x] = size (x);
   [r_b, c_b] = size (b);
-  if !( (min ([r_x c_x]) == 1) || (min ([r_b c_b]) == 1) )
-    error ("fftfilt:  both x and b should be vectors.");
+  if (! (min ([r_x c_x]) == 1 || min ([r_b c_b]) == 1))
+    error ("fftfilt: both x and b should be vectors");
   endif
   l_x  = r_x * c_x;
   l_b  = r_b * c_b;
@@ -35,25 +52,25 @@
     return;
   endif
   
-  x    = reshape (x, 1, l_x);
-  b    = reshape (b, 1, l_b);
+  x = reshape (x, 1, l_x);
+  b = reshape (b, 1, l_b);
   
   if (nargin == 2)
-    # use FFT with the smallest power of 2 which is >= length (x) +
-    # length (b) - 1 as number of points ...
+# Use FFT with the smallest power of 2 which is >= length (x) +
+# length (b) - 1 as number of points ...
     N    = 2^(ceil (log (l_x + l_b - 1) / log(2)));
     y    = ifft (fft (x, N) .* fft(b, N));
   else
-    # use overlap-add method ...
+# Use overlap-add method ...
     if !(is_scalar (N))
-      error ("fftfilt:  N has to be a scalar");
+      error ("fftfilt: N has to be a scalar");
     endif
-    N    = 2^(ceil (log (max ([N l_b])) / log(2)));
-    L    = N - l_b + 1;
-    B    = fft (b, N);
-    R    = ceil (l_x / L);
-    y    = zeros (1, l_x);
-    for r=1:R;
+    N = 2^(ceil (log (max ([N l_b])) / log(2)));
+    L = N - l_b + 1;
+    B = fft (b, N);
+    R = ceil (l_x / L);
+    y = zeros (1, l_x);
+    for r = 1:R;
       lo  = (r - 1) * L + 1;
       hi  = min (r * L, l_x);
       tmp = ifft (fft (x(lo:hi), N) .* B);
@@ -62,14 +79,15 @@
     endfor  
   endif
     
-  y    = reshape (y(1:l_x), r_x, c_x);
+  y = reshape (y(1:l_x), r_x, c_x);
 
-  # final cleanups:  if both x and b are real respectively integer, y
-  # should also be so
-  if !(any (imag (x)) || any (imag (b)))
+# Final cleanups:  if both x and b are real respectively integer, y
+# should also be
+
+  if (! (any (imag (x)) || any (imag (b))))
     y = real (y);
   endif
-  if !(any (x - round (x)) || any (b - round (b)))
+  if (! (any (x - round (x)) || any (b - round (b))))
     y = round (y);
   endif
 
diff --git a/scripts/specfun/beta.m b/scripts/specfun/beta.m
--- a/scripts/specfun/beta.m
+++ b/scripts/specfun/beta.m
@@ -1,15 +1,35 @@
-function retval = beta(a, b)
+# Copyright (C) 1995 John W. Eaton
+# 
+# This file is part of Octave.
+# 
+# Octave is free software; you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published by the
+# Free Software Foundation; either version 2, or (at your option) any
+# later version.
+# 
+# Octave is distributed in the hope that it will be useful, but WITHOUT
+# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+# for more details.
+# 
+# You should have received a copy of the GNU General Public License
+# along with Octave; see the file COPYING.  If not, write to the Free
+# Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+
+function retval = beta (a, b)
   
-# usage:  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));
+  if (nargin != 2)
+    usage ("beta (a, b)");
+  endif
+
+  retval = exp (lgamma (a) + lgamma (b) - lgamma (a+b));
 
 endfunction
 
diff --git a/scripts/specfun/betai.m b/scripts/specfun/betai.m
--- a/scripts/specfun/betai.m
+++ b/scripts/specfun/betai.m
@@ -1,6 +1,24 @@
+# Copyright (C) 1995 John W. Eaton
+# 
+# This file is part of Octave.
+# 
+# Octave is free software; you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published by the
+# Free Software Foundation; either version 2, or (at your option) any
+# later version.
+# 
+# Octave is distributed in the hope that it will be useful, but WITHOUT
+# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+# for more details.
+# 
+# You should have received a copy of the GNU General Public License
+# along with Octave; see the file COPYING.  If not, write to the Free
+# Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+
 function y = betai(a, b, x)
   
-# usage:  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.
@@ -8,7 +26,6 @@
 # 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) 
@@ -24,78 +41,79 @@
     error("betai:  a and b must both be positive.");
   endif
   [nr, nc] = size(x);
-  if (min([nr nc]) == 0)
+  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].");
+  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.");
+    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));
+    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));
+    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)';
+# Now do the series computation.  The error when truncating at term K
+# is always less than 2^(-K), hence the following choice of K.
 
-    ind = find((x > 0) & (x <= 1/2));
-    len = length(ind);
+    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));
+    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));
+      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);
+    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.");
+    [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));
+
+    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);
+      y   = zeros (1, n);
     elseif (x == 1)
-      y   = ones(1, n);
+      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));
+      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));
+	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);
+    y = reshape (y, ra, ca);
     
   endif
 
diff --git a/scripts/specfun/gammai.m b/scripts/specfun/gammai.m
--- a/scripts/specfun/gammai.m
+++ b/scripts/specfun/gammai.m
@@ -1,24 +1,44 @@
-function y = gammai(a, x)
+# Copyright (C) 1995 John W. Eaton
+# 
+# This file is part of Octave.
+# 
+# Octave is free software; you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published by the
+# Free Software Foundation; either version 2, or (at your option) any
+# later version.
+# 
+# Octave is distributed in the hope that it will be useful, but WITHOUT
+# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+# for more details.
+# 
+# You should have received a copy of the GNU General Public License
+# along with Octave; see the file COPYING.  If not, write to the Free
+# Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+
+function y = gammai (a, x)
   
-# usage:  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).
+#
+#   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)
-    usage (" gammai(a, x)");
+    usage ("gammai (a, x)");
   endif
   
-  [r_a, c_a] = size(a);
-  [r_x, c_x] = size(x);
+  [r_a, c_a] = size (a);
+  [r_x, c_x] = size (x);
   e_a = r_a * c_a;
   e_x = r_x * c_x;
   
@@ -26,38 +46,38 @@
   # 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");
+  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))
+  if (r_a == r_x && c_a == c_x)
     n   = e_a;
-    a   = reshape(a, 1, n);
-    x   = reshape(x, 1, n);
+    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);
+    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);
+    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"); 
+    error ("gammai: a and x must have the same size if neither is scalar"); 
   endif
 
-  # Now we can do sanity checking ...
+# Now we can do sanity checking ...
   
   if (any (a <= 0) || any (a == Inf))
-    error ("gammai:  all entries of a must be positive anf finite");
+    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");
+    error ("gammai: all entries of x must be nonnegative");
   endif
   
   y = zeros(1, n);
@@ -65,14 +85,14 @@
 # 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);
+  S = find ((x > 0) & (x < a + 1));
+  s = length (S);
   if (s > 0)
-    k   = ceil(- max([a(S), x(S)]) * log(eps));
+    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);
+    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.
@@ -87,18 +107,18 @@
     c_old = 0;
     c_new = v(1,:) ./ v(2,:);
     n   = 1;
-    while (max(abs(c_old ./ c_new - 1)) > 10 * eps)
+    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;
+      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));
+    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 (find (x == Inf)) = ones (1, sum (x == Inf));
   
-  y = reshape(y, r_y, c_y);
+  y = reshape (y, r_y, c_y);
 
 endfunction
\ No newline at end of file
diff --git a/scripts/statistics/kurtosis.m b/scripts/statistics/kurtosis.m
--- a/scripts/statistics/kurtosis.m
+++ b/scripts/statistics/kurtosis.m
@@ -22,7 +22,7 @@
 #
 # If x is a vector of length N, return the kurtosis
 #
-# 	kurtosis(x) = N^(-1) std(x)^(-4) SUM_i (x(i)-mean(x))^4 - 3
+#   kurtosis(x) = N^(-1) std(x)^(-4) SUM_i (x(i)-mean(x))^4 - 3
 #
 # of x.
 #
@@ -30,7 +30,6 @@
 # column.
 
 # Written by KH (Kurt.Hornik@ci.tuwien.ac.at) on Jul 29, 1994.
-# Copyright Dept of Probability Theory and Statistics TU Wien, Austria.
 
   if (nargin != 1)
     usage ("kurtosis (x)");
diff --git a/scripts/statistics/skewness.m b/scripts/statistics/skewness.m
--- a/scripts/statistics/skewness.m
+++ b/scripts/statistics/skewness.m
@@ -30,7 +30,6 @@
 # column.
 
 # Written by KH (Kurt.Hornik@ci.tuwien.ac.at) on Jul 29, 1994.
-# Copyright Dept of Probability Theory and Statistics TU Wien, Austria.
 
   if (nargin != 1)
     usage ("skewness (x)");