--- a/scripts/specfun/ellipke.m	Wed Sep 04 10:13:53 2013 -0700
+++ b/scripts/specfun/ellipke.m	Wed Sep 04 11:44:18 2013 -0700
@@ -25,7 +25,7 @@
 ## Compute complete elliptic integrals of the first K(@var{m}) and second
 ## E(@var{m}) kind.
 ##
-## @var{m} must be a scalar or real array with 0 @leq{} @var{m} @leq{} 1.
+## @var{m} must be a scalar or real array with -Inf @leq{} @var{m} @leq{} 1.
 ##
 ## The optional input @var{tol} is currently ignored (@sc{matlab} uses this
 ## to allow a faster, less accurate approximation).
@@ -48,55 +48,54 @@
     print_usage ();
   endif
 
-  k = e = zeros (size (m));
   m = m(:);
-  if (any (!isreal (m)))
-    error ("ellipke must have real m");
-  endif
-  if (any (m > 1))
-    error ("ellipke must have m <= 1");
+  if (! isreal (m))
+    error ("ellipke: M must be real");
+  elseif (any (m > 1))
+    error ("ellipke: M must be <= 1");
   endif
 
-  Nmax = 16;
-  idx = find (m == 1);
-  if (!isempty (idx))
-    k(idx) = Inf;
-    e(idx) = 1;
-  endif
+  k = e = zeros (size (m));
 
-  idx = find (m == -Inf);
-  if (!isempty (idx))
-    k(idx) = 0;
-    e(idx) = Inf;
-  endif
+  ## Handle extreme values
+  idx_1 = (m == 1);
+  k(idx_1) = Inf;
+  e(idx_1) = 1;
+
+  idx_neginf = (m == -Inf);
+  k(idx_neginf) = 0;
+  e(idx_neginf) = Inf;
 
   ## Arithmetic-Geometric Mean (AGM) algorithm
   ## ( Abramowitz and Stegun, Section 17.6 )
-  idx = find (m != 1 & m != -Inf);
-  if (!isempty (idx))
-    idx_neg = find (m < 0 & m != -Inf);
+  Nmax = 16;
+  idx = !idx_1 & !idx_neginf;
+  if (any (idx))
+    idx_neg = find (m < 0 & !idx_neginf);
     mult_k = 1./sqrt (1 - m(idx_neg));
     mult_e = sqrt (1 - m(idx_neg));
-    m(idx_neg) = -m(idx_neg)./(1 - m(idx_neg));
-    a = ones (length (idx), 1);
+    m(idx_neg) = -m(idx_neg) ./ (1 - m(idx_neg));
+    a = ones (sum (idx), 1);
     b = sqrt (1 - m(idx));
     c = sqrt (m(idx));
     f = 0.5;
-    sum = f*c.*c;
-    for n = 2:Nmax
+    sum = f*c.^2;
+    n = 2;
+    do
       t = (a + b)/2;
       c = (a - b)/2;
       b = sqrt (a.*b);
       a = t;
-      f = f * 2;
-      sum = sum + f*c.*c;
-      if (all (c./a < eps)) break; endif
-    endfor
-    if (n >= Nmax) error ("ellipke: not enough workspace"); endif
-    k(idx) = 0.5*pi./a;
-    e(idx) = 0.5*pi.*(1 - sum)./a;
-    k(idx_neg) = mult_k.*k(idx_neg);
-    e(idx_neg) = mult_e.*e(idx_neg);
+      f *= 2;
+      sum += f*c.^2;
+    until (all (c./a < eps) || (++n > Nmax))
+    if (n >= Nmax)
+      error ("ellipke: algorithm did not converge in %d iterations", Nmax);
+    endif
+    k(idx) = 0.5*pi ./ a;
+    e(idx) = 0.5*pi*(1 - sum) ./ a;
+    k(idx_neg) = mult_k .* k(idx_neg);
+    e(idx_neg) = mult_e .* e(idx_neg);
   endif
 
 endfunction
@@ -105,17 +104,16 @@
 ## Test complete elliptic functions of first and second kind
 ## against "exact" solution from Mathematica 3.0
 %!test
-%! m = [0.0; 0.01; 0.1; 0.5; 0.9; 0.99; 1.0 ];
+%! m = [0.0; 0.01; 0.1; 0.5; 0.9; 0.99; 1.0];
 %! [k,e] = ellipke (m);
 %!
-%! ## K(1.0) is really infinity - see below
 %! k_exp = [1.5707963267948966192;
 %!          1.5747455615173559527;
 %!          1.6124413487202193982;
 %!          1.8540746773013719184;
 %!          2.5780921133481731882;
 %!          3.6956373629898746778;
-%!          0.0 ];
+%!          Inf ];
 %! e_exp = [1.5707963267948966192;
 %!          1.5668619420216682912;
 %!          1.5307576368977632025;
@@ -123,7 +121,6 @@
 %!          1.1047747327040733261;
 %!          1.0159935450252239356;
 %!          1.0 ];
-%! if (k(7)==Inf), k(7)=0; endif;
 %! assert (k, k_exp, 8*eps);
 %! assert (e, e_exp, 8*eps);
 
@@ -156,6 +153,25 @@
 %! assert (k, k_exp, 1e-15);
 %! assert (e, e_exp, 1e-8);
 
+## Test negative values against "exact" solution from Mathematica.
+%! m = [-0.01; -1; -5; -100; -1000; -Inf];
+%! [k,e] = ellipke (m);
+%!
+%! k_exp = [1.5668912730681963584;
+%!          1.3110287771460599052;
+%!          0.9555039270640439337;
+%!          0.3682192486091410329;
+%!          0.1530293349884987857;
+%!          0];
+%! e_exp = [1.5747159850169884130;
+%!          1.9100988945138560089;
+%!          2.8301982463458773125;
+%!          10.209260919814572009;
+%!          31.707204053711259719;
+%!          Inf ];
+%! assert (k, k_exp, 8*eps);
+%! assert (e, e_exp, 8*eps (e_exp));
+
 ## Test input validation
 %!error ellipke ()
 %!error ellipke (1,2,3)