# HG changeset patch
# User jwe
# Date 848438107 0
# Node ID 80b982e7f4b19a5f6d05314545a734066121ecd6
# Parent  1d63e820ee13135b1d70c6ad18efeb7e8ef13083
[project @ 1996-11-19 21:15:06 by jwe]

diff --git a/scripts/specfun/betai.m b/scripts/specfun/betai.m
--- a/scripts/specfun/betai.m
+++ b/scripts/specfun/betai.m
@@ -1,123 +1,93 @@
-## Copyright (C) 1996 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
+## Copyright (C) 1995, 1996  Kurt Hornik
+## 
+## This program 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
+## 
+## This program 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.
-##
+## 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, 59 Temple Place - Suite 330, Boston, MA
-## 02111-1307, USA.
+## along with this file.  If not, write to the Free Software Foundation,
+## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
 
 ## 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.
 
 ## Author: KH <Kurt.Hornik@ci.tuwien.ac.at>
-## Created: 2 August 1994.
+## Created: 2 August 1994
 ## Adapted-By: jwe
 
+## Computation is based on the series expansion
+##   betai(a, b, x) 
+##     = \frac{x^a}{B(a, b)}
+##         \sum_{l=0}^\infty \frac{(1-b)\cdots(l-b)}{a+l} \frac{x^l}{l!}
+## for x <= 1/2.  For x > 1/2, betai(a, b, x) = 1 - betai(b, a, 1-x).
+
 function y = betai (a, b, x)
-
-  ## 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)
-    usage (" betai (a, b, x)");
+  
+  if (nargin != 3)
+    usage ("betai (a, b, x)");
   endif
 
-  if (! (a > 0 && b > 0))
-    error ("betai: a and b must both be positive");
+  [retval, x, a, b] = common_size (x, a, b);
+  if (retval > 0)
+    error ("betai:  a, b and x must be of common size or scalar");
   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].");
+  
+  [r, c] = size (x);
+  s = r * c;
+  x = reshape (x, 1, s);
+  a = reshape (a, 1, s);
+  b = reshape (b, 1, s);
+  y = zeros (1, s);
+
+  k = find ((x < 0) | (x > 1) | !(a > 0) | !(b > 0) | isnan (x));
+  if any (k)
+    y(k) = NaN * ones (1, length (k));
   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);
-
+  k = find ((x == 1) & (a > 0) & (b > 0));
+  if any (k)
+    y(k) = ones (1, length (k));
   endif
 
+  ## Now do the series computations.  
+  ## The error when truncating at term L is always less than 2^(-L),
+  ## hence the following choice of L. 
+
+  L = ceil (-log (eps) / log (2));
+  h = ones (L, 1);
+  
+  k = find ((x >= 0) & (x <= 1/2) & (a > 0) & (b > 0));
+  if any (k)
+    l   = (1 : L)' * ones (1, length (k));
+    tmp = cumprod ((1 - (h * b(k)) ./ l) .* (h * x(k))) ...
+	./ ((h * a(k)) + l);
+    y(k) = exp (a(k) .* log (x(k))) .* (1 ./ a(k) + sum (tmp)) ...
+	./ beta (a(k), b(k));
+  endif
+  
+  k = find ((x > 1/2) & (x < 1) & (a > 0) & (b > 0));
+  if any (k)
+    l   = (1 : L)' * ones (1, length (k));
+    tmp = cumprod ((1 - (h * a(k)) ./ l) .* (h * (1 - x(k)))) ...
+	./ (h * b(k) + l);
+    y(k) = 1 - exp (b(k) .* log (1 - x(k))) ...
+	.* (1 ./ b(k) + sum (tmp)) ./ beta (a(k), b(k));
+  endif
+
+  y = reshape (y, r, c);
+  
 endfunction
+
+
+
+
+
diff --git a/scripts/specfun/erfinv.m b/scripts/specfun/erfinv.m
new file mode 100644
--- /dev/null
+++ b/scripts/specfun/erfinv.m
@@ -0,0 +1,56 @@
+## Copyright (C) 1995, 1996  Kurt Hornik
+## 
+## This program 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.
+## 
+## This program 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 this file.  If not, write to the Free Software Foundation,
+## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+
+## usage:  erfinv (x)
+##
+## Computes the inverse of the error function erf.
+
+## Author: KH <Kurt.Hornik@ci.tuwien.ac.at>
+## Created: 27 September 1994
+## Adapted-By: jwe
+
+function y = erfinv (x)
+  
+  if (nargin != 1)
+    usage ("erfinv (x)");
+  endif
+  
+  [m, n] = size (x);  
+  x = reshape (x, m * n, 1);
+  y = zeros (m * n, 1);
+  
+  i = find ((x < -1) | (x > 1));
+  if any (i)
+    y(i) = NaN * ones (length (i), 1);
+  endif
+
+  y (find (x == -1)) = (-Inf) * ones (sum (x == -1), 1);
+  y (find (x == 1)) = Inf * ones (sum (x == 1), 1);
+  
+  i = find ((x > -1) & (x < 1));
+  if any (i)
+    z_old = ones (length (i), 1);
+    z_new = zeros (length (i), 1);
+    while (any (any (abs (z_new - z_old) > eps)))
+      z_old = z_new;
+      z_new = z_old - (erf (z_old) - x(i)) .* exp (z_old.^2);
+    endwhile
+    y(i) = z_new;
+  endif
+  
+  y = reshape (y, m, n);
+    
+endfunction
diff --git a/scripts/specfun/gammai.m b/scripts/specfun/gammai.m
--- a/scripts/specfun/gammai.m
+++ b/scripts/specfun/gammai.m
@@ -1,28 +1,25 @@
-## Copyright (C) 1996 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
+## Copyright (C) 1995, 1996  Kurt Hornik
+## 
+## This program 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
+## 
+## This program 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.
-##
+## 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, 59 Temple Place - Suite 330, Boston, MA
-## 02111-1307, USA.
+## along with this file.  If not, write to the Free Software Foundation,
+## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
 
 ## 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.
@@ -35,95 +32,65 @@
 ## Adapted-By: jwe
 
 function y = gammai (a, x)
-
+  
   if (nargin != 2)
     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;
+  
+  [retval, a, x] = common_size (a, x);
+  if (retval > 0)
+    error ("gammai:  a and x must be of common size or scalar");
+  endif
+  
+  [r, c] = size (x);
+  s = r * c;
+  x = reshape (x, 1, s);
+  a = reshape (a, 1, s);
+  y = zeros (1, s);
 
-  ## 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");
+  k = find (!(a > 0) | isnan (x));
+  if any (k)
+    y(k) = NaN * ones (1, length (k));
   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");
+  
+  k = find ((x == Inf) & (a > 0));
+  if any (k)
+    y(k) = ones (1, length (k));
   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);
+  
+  ## For x < a + 1, use summation.  The below choice of L should ensure
+  ## that the overall error is less than eps ... 
+  k = find((x > 0) & (x < a + 1));
+  if any (k)
+    L = ceil (- max ([a(k), x(k)]) * log (eps));
+    A = cumprod ((ones (L, 1) * x(k)) ...
+	./ (ones (L, 1) * a(k) + (1 : L)' * ones (1, length (k))));
+    y(k) = exp (-x(k) + a(k) .* log (x(k))) ...
+	.* (1 + sum (A)) ./ gamma (a(k) + 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)
-    t1 = zeros (1, s);
-    t2 = ones (1, s);
-    u   = [t1; t2];
-    v   = [t2; x(S)];
+  k = find ((x >= a + 1) & (x < Inf) & (a > 0));
+  if any (k)
+    len = length (k);
+    u   = [zeros (1, len); ones (1, len)];
+    v   = [ones (1, len); x(k)];
     c_old = 0;
-    c_new = v(1,:) ./ v(2,:);
+    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,:);
+      u = v + u .* (ones (2, 1) * (n - a(k)));
+      v = u .* (ones (2, 1) * x(k)) + 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(k) = 1 - exp (-x(k) + a(k) .* log (x(k))) .* c_new ...
+	./ gamma (a(k));
   endif
+  
+  y = reshape (y, r, c);
 
-  y (find (x == Inf)) = ones (1, sum (x == Inf));
-
-  y = reshape (y, r_y, c_y);
-
-endfunction
+endfunction
\ No newline at end of file
diff --git a/scripts/specfun/log2.m b/scripts/specfun/log2.m
new file mode 100644
--- /dev/null
+++ b/scripts/specfun/log2.m
@@ -0,0 +1,49 @@
+## Copyright (C) 1995, 1996  Kurt Hornik
+## 
+## This program 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.
+## 
+## This program 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 this file.  If not, write to the Free Software Foundation,
+## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+
+## usage:  y = log2 (x) or [f, e] = log2 (x)
+##
+## y = log2 (x) returns the logarithm of base 2 of x.
+##
+## [f, e] = log2 (x) returns f and e with 1/2 <= abs(f) < 1 and
+## x = f * 2^e.
+
+## Author: AW <Andreas.Weingessel@ci.tuwien.ac.at>
+## Created: 17 October 1994
+## Adapted-By: jwe
+
+function [f, e] = log2 (x)
+
+  if (nargin != 1)
+    usage ("y = log2 (x) or [f, e] = log2 (x)");
+  endif
+
+  if (nargout < 2)
+    f = log (x) / log (2);
+  elseif (nargout == 2)
+    ## Only deal with the real parts ...
+    x = real (x);
+    ## Since log (0) gives problems, 0 entries are replaced by 1.  
+    ## This is corrected later by multiplication with the sign.
+    f = abs (x) + (x == 0);
+    e = (floor (log (f) / log (2)) + 1) .* (x != 0);
+    f = sign (x) .* f ./ (2 .^ e);
+  else
+    error ("log2 takes at most 2 output arguments");
+  endif
+
+endfunction
+
diff --git a/scripts/specfun/pow2.m b/scripts/specfun/pow2.m
new file mode 100644
--- /dev/null
+++ b/scripts/specfun/pow2.m
@@ -0,0 +1,36 @@
+## Copyright (C) 1995, 1996  Kurt Hornik
+## 
+## This program 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.
+## 
+## This program 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 this file.  If not, write to the Free Software Foundation,
+## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+
+## usage: y = pow2 (f [, e])
+##
+## y = pow2 (x) returns 2 .^ x.
+## y = pow2 (f, e) returns f .* (2 .^ e).
+
+## Author: AW <Andreas.Weingessel@ci.tuwien.ac.at>
+## Created: 17 October 1994
+## Adapted-By: jwe
+
+function y = pow2 (f, e)
+  
+  if (nargin == 1)
+    y = 2 .^ f;
+  elseif (nargin == 2)
+    y = f .* (2 .^ e);
+  else
+    usage ("y = pow2 (f [, e])");
+  endif
+
+endfunction