# HG changeset patch
# User David Bateman <dbateman@free.fr>
# Date 1210625831 -7200
# Node ID df9519e9990c2b04db5cdc32d0953e0d9f4970ce
# Parent  2b458dfe31ae7e645cbdd0167b7d5d7bb44bd30b
Handle single precision eps values

diff -r 2b458dfe31ae -r df9519e9990c scripts/ChangeLog
--- a/scripts/ChangeLog	Mon May 12 01:35:30 2008 +0200
+++ b/scripts/ChangeLog	Mon May 12 22:57:11 2008 +0200
@@ -24,6 +24,26 @@
 
 2008-05-12  David Bateman  <dbateman@free.fr>
 
+	* control/base/__stepimp__.m, control/base/bode_bounds.m,
+	control/base/damp.m, control/base/dlqr.m, control/base/lsim.m,
+	control/base/tzero.m, control/hinf/hinfsyn.m,
+	control/hinf/is_dgkf.m, control/system/d2c.m,
+	control/system/is_controllable.m, control/system/is_detectable.m,
+	control/system/is_stabilizable.m, control/system/is_stable.m,
+	control/system/sysconnect.m, general/bicubic.m,
+	general/cplxpair.m, general/isdefinite.m, general/issymmetric.m,
+	general/quadgk.m, general/quadl.m, general/quadv.m,
+	geometry/delaunayn.m, linear-algebra/krylov.m,
+	linear-algebra/null.m, linear-algebra/onenormest.m,
+	linear-algebra/orth.m, linear-algebra/rank.m,
+	linear-algebra/rref.m, optimization/qp.m, optimization/sqp.m,
+	polynomial/polygcd.m, polynomial/residue.m, sparse/normest.m,
+	specfun/erfinv.m, statistics/distributions/betainv.m,
+	statistics/distributions/gaminv.m,
+	statistics/distributions/kolmogorov_smirnov_cdf.m,
+	statistics/tests/manova.m: Modify calls to eps to allow for single
+	precision types.
+
 	* general/isa.m: Also treat "float: and "numeric" as the class
 	argument.
 
diff -r 2b458dfe31ae -r df9519e9990c scripts/control/base/__stepimp__.m
--- a/scripts/control/base/__stepimp__.m	Mon May 12 01:35:30 2008 +0200
+++ b/scripts/control/base/__stepimp__.m	Mon May 12 22:57:11 2008 +0200
@@ -66,7 +66,8 @@
   DIGITAL = is_digital (sys);
   if (DIGITAL)
     NSTATES = ndstates;
-    if (TSAMPLE < eps)
+    if (isa (TSAMPLE, "single") && TSAMPLE < eps ("single") ||
+	!isa (TSAMPLE, "single") && TSAMPLE < eps)
       error ("__stepimp__: sampling time of discrete system too small")
     endif
   else
diff -r 2b458dfe31ae -r df9519e9990c scripts/control/base/bode_bounds.m
--- a/scripts/control/base/bode_bounds.m	Mon May 12 01:35:30 2008 +0200
+++ b/scripts/control/base/bode_bounds.m	Mon May 12 22:57:11 2008 +0200
@@ -48,11 +48,17 @@
     zer = reshape (zer, 1, length (zer));
   endif
 
+  if (isa (zer, "single") || isa (pol, "single"))
+    myeps = eps ("single");
+  else
+    myeps = eps;
+  endif
+
   ## check for natural frequencies away from omega = 0
   if (DIGITAL)
     ## The 2nd conditions prevents log(0) in the next log command
-    iiz = find (abs(zer-1) > norm(zer)*eps && abs(zer) > norm(zer)*eps);
-    iip = find (abs(pol-1) > norm(pol)*eps && abs(pol) > norm(pol)*eps);
+    iiz = find (abs(zer-1) > norm(zer)*myeps && abs(zer) > norm(zer)*myeps);
+    iip = find (abs(pol-1) > norm(pol)*myeps && abs(pol) > norm(pol)*myeps);
 
     ## avoid dividing empty matrices, it would work but looks nasty
     if (! isempty (iiz))
@@ -68,8 +74,8 @@
     endif
   else
     ## continuous
-    iip = find (abs(pol) > norm(pol)*eps);
-    iiz = find (abs(zer) > norm(zer)*eps);
+    iip = find (abs(pol) > norm(pol)*myeps);
+    iiz = find (abs(zer) > norm(zer)*myeps);
 
     if (! isempty (zer))
       czer = zer(iiz);
diff -r 2b458dfe31ae -r df9519e9990c scripts/control/base/damp.m
--- a/scripts/control/base/damp.m	Mon May 12 01:35:30 2008 +0200
+++ b/scripts/control/base/damp.m	Mon May 12 22:57:11 2008 +0200
@@ -76,7 +76,8 @@
     endif
     d0 = -cos (atan2 (imag (cpole), real (cpole)));
     f0 = 0.5 / pi * abs (cpole);
-    if (abs (imag (cpole)) < eps)
+    if (isa (cpole, "single") && abs(imag (cpole)) < eps ("single") || 
+	! isa (cpole, "single") && abs (imag (cpole)) < eps)
       printf ("%12f         ---  %12f  %10f  %12f\n",
              real (pole), abs (pole), d0, f0);
     else
diff -r 2b458dfe31ae -r df9519e9990c scripts/control/base/dlqr.m
--- a/scripts/control/base/dlqr.m	Mon May 12 01:35:30 2008 +0200
+++ b/scripts/control/base/dlqr.m	Mon May 12 22:57:11 2008 +0200
@@ -130,7 +130,11 @@
 
   ## Checking stabilizability and detectability (dimensions are checked
   ## inside these calls).
-  tol = 200*eps;
+  if (isa (a, "single") || isa (b, "single") || isa (q, "single") || isa (r, "single"))
+    tol = 200 * eps ("single");
+  else
+    tol = 200 * eps;
+  endif
   if (is_stabilizable (ao, b, tol, 1) == 0)
     error ("dlqr: (a,b) not stabilizable");
   endif
diff -r 2b458dfe31ae -r df9519e9990c scripts/control/base/lsim.m
--- a/scripts/control/base/lsim.m	Mon May 12 01:35:30 2008 +0200
+++ b/scripts/control/base/lsim.m	Mon May 12 22:57:11 2008 +0200
@@ -74,13 +74,17 @@
   t(2)-t(1);
   u=u';
   n = max (size (t));
-
+  
+  if (isa (t, "single"))
+    tol = 10 * eps ("single");
+  else
+    tol = 10 * eps;
+  endif
   for ii = 1:(n-1)
-
     ## check if step size changed
     ## FIXME -- this is probably not the best test, but it is
     ## better than a test for exact equality.
-    if (abs (t(ii+1) - t(ii) - Ts) > 10 * eps)
+    if (abs (t(ii+1) - t(ii) - Ts) > tol)
       Ts = t(ii+1) - t(ii);
       ## [F,G] = c2d(a,b,Ts);
       dsys = c2d (sys, Ts);
diff -r 2b458dfe31ae -r df9519e9990c scripts/control/base/tzero.m
--- a/scripts/control/base/tzero.m	Mon May 12 01:35:30 2008 +0200
+++ b/scripts/control/base/tzero.m	Mon May 12 22:57:11 2008 +0200
@@ -98,7 +98,11 @@
   ## balance coefficients
   Asys = __zgpbal__ (Asys);
   [A, B, C, D] = sys2ss (Asys);
-  meps = 2*eps*norm ([A, B; C, D], "fro");
+  if (isa ([A, B; C, D], "single"))
+    meps = 2*eps("single")*norm ([A, B; C, D], "fro");
+  else
+    meps = 2*eps*norm ([A, B; C, D], "fro");
+  endif
   ## ENVD algorithm
   Asys = zgreduce (Asys, meps);
   [A, B, C, D] = sys2ss (Asys);
diff -r 2b458dfe31ae -r df9519e9990c scripts/control/hinf/hinfsyn.m
--- a/scripts/control/hinf/hinfsyn.m	Mon May 12 01:35:30 2008 +0200
+++ b/scripts/control/hinf/hinfsyn.m	Mon May 12 22:57:11 2008 +0200
@@ -136,7 +136,12 @@
   endif
   ## set default arguments
   if (nargin < 8)
-    tol = 200*eps;
+    if (isa (Asys.a, "single") || isa (Asys.b, "single") || isa (Asys.c, "single") ||
+	isa (Asys.d, "single"))
+      tol = 200*eps("single");
+    else
+      tol = 200*eps;
+    endif
   elseif (! is_sample (tol))
     error ("tol must be a positive scalar.")
   endif
diff -r 2b458dfe31ae -r df9519e9990c scripts/control/hinf/is_dgkf.m
--- a/scripts/control/hinf/is_dgkf.m	Mon May 12 01:35:30 2008 +0200
+++ b/scripts/control/hinf/is_dgkf.m	Mon May 12 22:57:11 2008 +0200
@@ -139,8 +139,13 @@
     error ("Argument 1 must be a system data structure");
   endif
   if (nargin < 4)
-    tol = 200*eps;
-  elseif (! is_sample (tol))
+    if (isa (Asys.a, "single") || isa (Asys.b, "single") || isa (Asys.c, "single") ||
+	isa (Asys.d, "single"))
+      tol = 200*eps("single");
+    else
+      tol = 200*eps;
+    endif
+      elseif (! is_sample (tol))
     error ("is_dgkf: tol must be a positive scalar")
   endif
 
diff -r 2b458dfe31ae -r df9519e9990c scripts/control/system/d2c.m
--- a/scripts/control/system/d2c.m	Mon May 12 01:35:30 2008 +0200
+++ b/scripts/control/system/d2c.m	Mon May 12 22:57:11 2008 +0200
@@ -89,6 +89,13 @@
   endif
   T = sysgettsam (sys);
 
+  if (isa (sys.a, "single") || isa (sys.b, "single") || isa (sys.c, "single") ||
+      isa (sys.d, "single"))
+    myeps = eps ("single");
+  else
+    myeps = eps;
+  endif
+
   if (strcmp (opt, "bi"))
     ## bilinear transform
     ## convert with bilinear transform
@@ -98,7 +105,7 @@
     [a, b, c, d, tsam, n, nz, stname, inname, outname, yd] = sys2ss (sys);
 
     poles = eig (a);
-    if (find (abs (poles-1) < 200*(n+nz)*eps))
+    if (find (abs (poles-1) < 200*(n+nz)*myeps))
       warning ("d2c: some poles very close to one.  May get bad results.");
     endif
 
@@ -136,10 +143,10 @@
     endif
 
     poles = eig (a);
-    if (find (abs (poles) < 200*(n+nz)*eps))
+    if (find (abs (poles) < 200*(n+nz)*myeps))
       warning ("d2c: some poles very close to zero.  logm not performed");
       Mtop = zeros (ma, na+nb);
-    elseif (find (abs (poles-1) < 200*(n+nz)*eps))
+    elseif (find (abs (poles-1) < 200*(n+nz)*myeps))
       warning ("d2c: some poles very close to one.  May get bad results.");
       logmat = real (logm (Amat) / T);
       Mtop = logmat(1:na,:);
diff -r 2b458dfe31ae -r df9519e9990c scripts/control/system/is_controllable.m
--- a/scripts/control/system/is_controllable.m	Mon May 12 01:35:30 2008 +0200
+++ b/scripts/control/system/is_controllable.m	Mon May 12 22:57:11 2008 +0200
@@ -90,7 +90,11 @@
 
   ## check for default tolerance
   if (deftol)
-    tol = 1000*eps;
+    if (isa (a, "single") || isa (b, "single"))
+      tol = 1000 * eps("single");
+    else
+      tol = 1000*eps;
+    endif
   endif
 
   ## check tol dimensions
diff -r 2b458dfe31ae -r df9519e9990c scripts/control/system/is_detectable.m
--- a/scripts/control/system/is_detectable.m	Mon May 12 01:35:30 2008 +0200
+++ b/scripts/control/system/is_detectable.m	Mon May 12 22:57:11 2008 +0200
@@ -59,7 +59,11 @@
   endif
 
   if (! exist ("tol"))
-    tol = 200*eps;
+    if (isa (a, "single") || isa (c, "single"))
+      tol = 200 * eps("single");
+    else
+      tol = 200 * eps;
+    endif
   endif
 
   retval = is_stabilizable (a', c', tol, dflg);
diff -r 2b458dfe31ae -r df9519e9990c scripts/control/system/is_stabilizable.m
--- a/scripts/control/system/is_stabilizable.m	Mon May 12 01:35:30 2008 +0200
+++ b/scripts/control/system/is_stabilizable.m	Mon May 12 22:57:11 2008 +0200
@@ -70,10 +70,13 @@
   endif
 
   if (! exist ("tol"))
-    tol = 200*eps;
+    if (isa (a, "single") || isa (b, "single"))
+      tol = 200 * eps("single");
+    else
+      tol = 200 * eps;
+    endif
   endif
 
-
   ## Checking dimensions
   n = is_square (a);
   if (n == 0)
diff -r 2b458dfe31ae -r df9519e9990c scripts/control/system/is_stable.m
--- a/scripts/control/system/is_stable.m	Mon May 12 01:35:30 2008 +0200
+++ b/scripts/control/system/is_stable.m	Mon May 12 22:57:11 2008 +0200
@@ -69,7 +69,11 @@
   endif
 
   if (nargin < 2)
-    tol = 200*eps;
+    if (isa (a, "single"))
+      tol = 200 * eps("single");
+    else
+      tol = 200 * eps;
+    endif
   elseif (! isscalar (tol))
     error ("is_stable: tol(%dx%d) must be a scalar", rows (tol),
 	   columns (tol));
diff -r 2b458dfe31ae -r df9519e9990c scripts/control/system/sysconnect.m
--- a/scripts/control/system/sysconnect.m	Mon May 12 01:35:30 2008 +0200
+++ b/scripts/control/system/sysconnect.m	Mon May 12 22:57:11 2008 +0200
@@ -86,12 +86,19 @@
     error ("sysconnect: order must be either 0 or 1")
   endif
 
+  if (isa (sys.a, "single") || isa (sys.b, "single") || isa (sys.c, "single") ||
+      isa (sys.d, "single"))
+    myeps = eps ("single");
+  else
+    myeps = eps;
+  endif
+
   if (nargin < 5)
-    tol = 200*eps;
+    tol = 200*myeps;
   elseif (! is_sample (tol))
     error ("sysconnect: tol must be a positive scalar");
-  elseif (tol > 1e2*sqrt(eps))
-    warning ("sysconnect: tol set to large value=%g, eps=%g", tol, eps);
+  elseif (tol > 1e2*sqrt(myeps))
+    warning ("sysconnect: tol set to large value=%g, eps=%g", tol, myeps);
   endif
 
   ## convert signal names to indices
diff -r 2b458dfe31ae -r df9519e9990c scripts/general/bicubic.m
--- a/scripts/general/bicubic.m	Mon May 12 01:35:30 2008 +0200
+++ b/scripts/general/bicubic.m	Mon May 12 22:57:11 2008 +0200
@@ -48,6 +48,13 @@
     extrapval = NaN;
   endif
 
+  if (isa (X, "single") || isa (Y, "single") || isa (Z, "single") || 
+      isa (XI, "single") || isa (YI, "single"))
+    myeps = eps("single");
+  else
+    myeps = eps;
+  endif
+
   if (nargin <= 2)
     ## bicubic (Z) or bicubic (Z, 2)
     if (nargin == 1) 
@@ -94,9 +101,9 @@
 
 
     X = reshape (X, 1, cz);
-    X(cz) *= 1 + sign (X(cz))*eps;
+    X(cz) *= 1 + sign (X(cz))*myeps;
     if (X(cz) == 0) 
-      X(cz) = eps;
+      X(cz) = myeps;
     endif; 
     XI = reshape (XI, 1, length (XI));
     [m, i] = sort ([X, XI]);
@@ -104,9 +111,9 @@
     xidx = o(find (i > cz));
     
     Y = reshape (Y, rz, 1);
-    Y(rz) *= 1 + sign (Y(rz))*eps;
+    Y(rz) *= 1 + sign (Y(rz))*myeps;
     if (Y(rz) == 0) 
-      Y(rz) = eps;
+      Y(rz) = myeps;
     endif; 
     YI = reshape (YI, length (YI), 1);
     [m, i] = sort ([Y; YI]);
diff -r 2b458dfe31ae -r df9519e9990c scripts/general/cplxpair.m
--- a/scripts/general/cplxpair.m	Mon May 12 01:35:30 2008 +0200
+++ b/scripts/general/cplxpair.m	Mon May 12 22:57:11 2008 +0200
@@ -60,7 +60,11 @@
   endif
 
   if (nargin < 2 || isempty (tol))
-    tol = 100*eps; 
+    if (isa (z, "single"))
+      tol = 100 * eps("single");
+    else
+      tol = 100*eps; 
+    endif
   endif
 
   nd = ndims (z);
@@ -95,7 +99,11 @@
   z = z(idx + n * ones (n, 1) * [0:m-1]);
 
   ## Put the purely real values at the end of the returned list
-  [idxi, idxj] = find (abs (imag (z)) ./ (abs (z) + realmin) < tol );
+  cls = "double";
+  if (isa (z, "single"))
+    cls = "single";
+  endif
+  [idxi, idxj] = find (abs (imag (z)) ./ (abs (z) + realmin(cls)) < tol );
   q = sparse (idxi, idxj, 1, n, m);
   nr = sum (q, 1);
   [q, idx] = sort (q, 1);
diff -r 2b458dfe31ae -r df9519e9990c scripts/general/isdefinite.m
--- a/scripts/general/isdefinite.m	Mon May 12 01:35:30 2008 +0200
+++ b/scripts/general/isdefinite.m	Mon May 12 22:57:11 2008 +0200
@@ -33,7 +33,11 @@
 
   if (nargin == 1 || nargin == 2)
     if (nargin == 1)
-      tol = 100*eps;
+      if (isa (x, "single"))
+	tol = 100 * eps("single");
+      else
+	tol = 100*eps; 
+      endif
     endif
     sym = issymmetric (x, tol);
     if (sym > 0)
diff -r 2b458dfe31ae -r df9519e9990c scripts/general/issymmetric.m
--- a/scripts/general/issymmetric.m	Mon May 12 01:35:30 2008 +0200
+++ b/scripts/general/issymmetric.m	Mon May 12 22:57:11 2008 +0200
@@ -38,7 +38,11 @@
     retval = issquare (x);
     if (retval != 0)
       if (nargin == 1)
-        tol = eps;
+	if (isa (x, "single"))
+	  tol = eps("single");
+	else
+	  tol = eps; 
+	endif
       endif
       norm_x = norm (x, inf);
       if (norm_x != 0 && norm (x - x', inf) / norm_x > tol)
diff -r 2b458dfe31ae -r df9519e9990c scripts/general/quadgk.m
--- a/scripts/general/quadgk.m	Mon May 12 01:35:30 2008 +0200
+++ b/scripts/general/quadgk.m	Mon May 12 22:57:11 2008 +0200
@@ -291,12 +291,18 @@
       q0 = sum (q_subs);
       err0 = sum (q_errs);
     
+      if (isa (a, "single") || isa (b, "single") || isa (waypoints, "single"))
+	myeps = eps ("single");
+      else
+	myeps = eps;
+      endif
+
       first = true;
       while (true)
 	## Check for sub-intervals that are too small. Test must be
 	## performed in untransformed sub-intervals. What is a good
 	## value for this test. Shampine suggests 100*eps
-	if (any (diff (trans (subs), [], 2) / (b - a) < 100 * eps))
+	if (any (diff (trans (subs), [], 2) / (b - a) < 100 * myeps))
 	  q = q0;
 	  err = err0;
 	  break;
diff -r 2b458dfe31ae -r df9519e9990c scripts/general/quadl.m
--- a/scripts/general/quadl.m	Mon May 12 01:35:30 2008 +0200
+++ b/scripts/general/quadl.m	Mon May 12 22:57:11 2008 +0200
@@ -63,14 +63,19 @@
   if (nargin < 5)
     trace = []; 
   endif
+  if (isa (a, "single") || isa (b, "single"))
+    myeps = eps ("single");
+  else
+    myeps = eps;
+  endif
   if (isempty (tol))
-    tol = eps; 
+    tol = myeps; 
   endif
   if (isempty (trace))
     trace = 0; 
   endif
-  if (tol < eps)
-    tol = eps;
+  if (tol < myeps)
+    tol = myeps;
   endif
 
   m = (a+b)/2; 
@@ -119,7 +124,7 @@
   if (R > 0 && R < 1)
     tol = tol/R; 
   endif
-  is = s*abs(is)*tol/eps;
+  is = s*abs(is)*tol/myeps;
   if (is == 0)
     is = b-a;
   endif
diff -r 2b458dfe31ae -r df9519e9990c scripts/general/quadv.m
--- a/scripts/general/quadv.m	Mon May 12 01:35:30 2008 +0200
+++ b/scripts/general/quadv.m	Mon May 12 22:57:11 2008 +0200
@@ -52,6 +52,11 @@
   if (nargin < 5)
     trace = []; 
   endif
+  if (isa (a, "single") || isa (b, "single"))
+    myeps = eps ("single");
+  else
+    myeps = eps;
+  endif
   if (isempty (tol))
     tol = 1e-6; 
   endif
@@ -69,10 +74,10 @@
   ## If have edge singularities, move edge point by eps*(b-a) as
   ## discussed in Shampine paper used to implement quadgk
   if (isinf (fa))
-    fa = feval (f, a + eps * (b-a), varargin{:});
+    fa = feval (f, a + myeps * (b-a), varargin{:});
   endif
   if (isinf (fb))
-    fb = feval (f, b - eps * (b-a), varargin{:});
+    fb = feval (f, b - myeps * (b-a), varargin{:});
   endif
 
   h = (b - a) / 2;
@@ -85,7 +90,7 @@
     warning("Maximum iteration count reached");
   elseif (isnan(Q) || isinf (Q))
     warning ("Infinite or NaN function evaluations were returned");
-  elseif (hmin < (b - a) * eps)
+  elseif (hmin < (b - a) * myeps)
     warning ("Minimum step size reached. Possibly singular integral");
   endif
 endfunction
diff -r 2b458dfe31ae -r df9519e9990c scripts/geometry/delaunayn.m
--- a/scripts/geometry/delaunayn.m	Mon May 12 01:35:30 2008 +0200
+++ b/scripts/geometry/delaunayn.m	Mon May 12 22:57:11 2008 +0200
@@ -54,6 +54,12 @@
 
   t = __delaunayn__ (x, varargin{:});
 
+  if (isa (x, "single"))
+    myeps = eps ("single");
+  else
+    myeps = eps;
+  endif
+
   ## Try to remove the zero volume simplices. The volume of the i-th simplex is
   ## given by abs(det(x(t(i,1:end-1),:)-x(t(i,2:end),:)))/prod(1:n) 
   ## (reference http://en.wikipedia.org/wiki/Simplex). Any simplex with a 
@@ -67,7 +73,7 @@
   [nt, n] = size (t);
   for i = 1:nt
     X = x(t(i,1:end-1),:) - x(t(i,2:end),:);
-    if (abs (det (X)) /  sqrt (sum (X .^ 2, 2)) < 1e3 * eps)
+    if (abs (det (X)) /  sqrt (sum (X .^ 2, 2)) < 1e3 * myeps)
      idx = [idx, i];
     endif
   endfor
diff -r 2b458dfe31ae -r df9519e9990c scripts/linear-algebra/krylov.m
--- a/scripts/linear-algebra/krylov.m	Mon May 12 01:35:30 2008 +0200
+++ b/scripts/linear-algebra/krylov.m	Mon May 12 22:57:11 2008 +0200
@@ -59,7 +59,11 @@
 
 function [Uret, H, nu] = krylov (A, V, k, eps1, pflg);
 
-  defeps = 1e-12;
+  if (isa (A, "single") || isa (V, "single"))
+    defeps = 1e-6
+  else
+    defeps = 1e-12;
+  endif
 
   if (nargin < 3 || nargin > 5)
     print_usage ();
diff -r 2b458dfe31ae -r df9519e9990c scripts/linear-algebra/null.m
--- a/scripts/linear-algebra/null.m	Mon May 12 01:35:30 2008 +0200
+++ b/scripts/linear-algebra/null.m	Mon May 12 22:57:11 2008 +0200
@@ -52,7 +52,11 @@
     endif
 
     if (nargin == 1)
-      tol = max (size (A)) * s (1) * eps;
+      if (isa (A, "single"))
+	tol = max (size (A)) * s (1) * eps ("single");
+      else
+	tol = max (size (A)) * s (1) * eps;
+      endif
     elseif (nargin != 2)
       print_usage ();
     endif
@@ -61,7 +65,11 @@
 
     if (rank < cols)
       retval = V (:, rank+1:cols);
-      retval(abs (retval) < eps) = 0;
+      if (isa (A, "single"))
+	retval(abs (retval) < eps ("single")) = 0;
+      else
+	retval(abs (retval) < eps) = 0;
+      endif
     else
       retval = zeros (cols, 0);
     endif
diff -r 2b458dfe31ae -r df9519e9990c scripts/linear-algebra/onenormest.m
--- a/scripts/linear-algebra/onenormest.m	Mon May 12 01:35:30 2008 +0200
+++ b/scripts/linear-algebra/onenormest.m	Mon May 12 22:57:11 2008 +0200
@@ -111,6 +111,7 @@
     else
       t = min (n, default_t);
     endif
+    issing = isa (varargin {1}, "single");
   else
     if (size (varargin, 2) < 3)
       print_usage();
@@ -123,6 +124,7 @@
     else
       t = default_t;
     endif
+    issing = isa (varargin {3}, "single");
   endif
 
   ## Initial test vectors X.
@@ -133,6 +135,13 @@
   est_old = 0; # To check if the estimate has increased.
   S = zeros (n, t); # Normalized vector of signs.  The normalization is 
 
+  if (issing)
+    myeps = eps ("single");
+    X = single (X);
+  else
+    myeps = eps;
+  endif
+
   for iter = 1 : itmax + 1
     Y = feval (apply, X);
 
diff -r 2b458dfe31ae -r df9519e9990c scripts/linear-algebra/orth.m
--- a/scripts/linear-algebra/orth.m	Mon May 12 01:35:30 2008 +0200
+++ b/scripts/linear-algebra/orth.m	Mon May 12 22:57:11 2008 +0200
@@ -51,7 +51,11 @@
     endif
 
     if (nargin == 1)
-      tol = max (size (A)) * s (1) * eps;
+      if (isa (A, "single"))
+	tol = max (size (A)) * s (1) * eps ("single");
+      else
+	tol = max (size (A)) * s (1) * eps;
+      endif
     endif
 
     rank = sum (s > tol);
diff -r 2b458dfe31ae -r df9519e9990c scripts/linear-algebra/rank.m
--- a/scripts/linear-algebra/rank.m	Mon May 12 01:35:30 2008 +0200
+++ b/scripts/linear-algebra/rank.m	Mon May 12 22:57:11 2008 +0200
@@ -42,7 +42,11 @@
     if (isempty (sigma))
       tolerance = 0;
     else
-      tolerance = max (size (A)) * sigma (1) * eps;
+      if (isa (A, "single"))
+	tolerance = max (size (A)) * sigma (1) * eps ("single");
+      else
+	tolerance = max (size (A)) * sigma (1) * eps;
+      endif
     endif
   elseif (nargin == 2)
     sigma = svd (A);
diff -r 2b458dfe31ae -r df9519e9990c scripts/linear-algebra/rref.m
--- a/scripts/linear-algebra/rref.m	Mon May 12 01:35:30 2008 +0200
+++ b/scripts/linear-algebra/rref.m	Mon May 12 22:57:11 2008 +0200
@@ -44,7 +44,11 @@
   [rows, cols] = size (A);
 
   if (nargin < 2)
-    tolerance = eps * max (rows, cols) * norm (A, inf);
+    if (isa (A, "single"))
+      tolerance = eps ("single") * max (rows, cols) * norm (A, inf ("single"));
+    else
+      tolerance = eps * max (rows, cols) * norm (A, inf);
+    endif
   endif
 
   used = zeros (1, cols);
diff -r 2b458dfe31ae -r df9519e9990c scripts/optimization/qp.m
--- a/scripts/optimization/qp.m	Mon May 12 01:35:30 2008 +0200
+++ b/scripts/optimization/qp.m	Mon May 12 22:57:11 2008 +0200
@@ -188,7 +188,12 @@
     n_in = length (bin);
 
     ## Check if the initial guess is feasible.
-    rtol = sqrt (eps);
+    if (isa (x0, "single") || isa (H, "single") || isa (q, "single") || isa (A, "single")
+	|| isa (b, "single"))
+      rtol = sqrt (eps ("single"));
+    else
+      rtol = sqrt (eps);
+    endif
 
     eq_infeasible = (n_eq > 0 && norm (A*x0-b) > rtol*(1+norm (b)));
     in_infeasible = (n_in > 0 && any (Ain*x0-bin < -rtol*(1+norm (bin))));
diff -r 2b458dfe31ae -r df9519e9990c scripts/optimization/sqp.m
--- a/scripts/optimization/sqp.m	Mon May 12 01:35:30 2008 +0200
+++ b/scripts/optimization/sqp.m	Mon May 12 22:57:11 2008 +0200
@@ -281,7 +281,11 @@
 	if (isvector (lb))
 	  __sqp_lb__ = lb;
 	elseif (isempty (lb))
-	  __sqp_lb__ = -realmax;
+	  if (isa (x, "single"))
+	    __sqp_lb__ = -realmax ("single");
+	  else
+	    __sqp_lb__ = -realmax;
+	  endif
 	else
 	  error ("sqp: invalid lower bound");
 	endif
@@ -290,7 +294,11 @@
 	if (isvector (ub))
 	  __sqp_ub__ = ub;
 	elseif (isempty (lb))
-	  __sqp_ub__ = realmax;
+	  if (isa (x, "single"))
+	    __sqp_ub__ = realmax ("single");
+	  else
+	    __sqp_ub__ = realmax;
+	  endif
 	else
 	  error ("sqp: invalid upper bound");
 	endif
diff -r 2b458dfe31ae -r df9519e9990c scripts/polynomial/polygcd.m
--- a/scripts/polynomial/polygcd.m	Mon May 12 01:35:30 2008 +0200
+++ b/scripts/polynomial/polygcd.m	Mon May 12 22:57:11 2008 +0200
@@ -45,7 +45,11 @@
 
   if (nargin == 2 || nargin == 3)
     if (nargin == 2)
-      tol = sqrt (eps);
+      if (isa (a, "single") || isa (b, "single"))
+	tol = sqrt (eps ("single"));
+      else
+	tol = sqrt (eps);
+      endif
     endif
     if (length (a) == 1 || length (b) == 1)
       if (a == 0)
diff -r 2b458dfe31ae -r df9519e9990c scripts/polynomial/residue.m
--- a/scripts/polynomial/residue.m	Mon May 12 01:35:30 2008 +0200
+++ b/scripts/polynomial/residue.m	Mon May 12 22:57:11 2008 +0200
@@ -219,7 +219,11 @@
   ## Determine if the poles are (effectively) zero.
 
   small = max (abs (p));
-  small = max ([small, 1]) * eps*1e4 * (1 + numel (p))^2;
+  if (isa (a, "single") || isa (b, "single"))
+    small = max ([small, 1]) * eps ("single") * 1e4 * (1 + numel (p))^2;
+  else
+    small = max ([small, 1]) * eps * 1e4 * (1 + numel (p))^2;
+  endif
   p(abs (p) < small) = 0;
 
   ## Determine if the poles are (effectively) real, or imaginary.
@@ -334,8 +338,11 @@
   endif
 
   ## Check for leading zeros and trim the polynomial coefficients.
-
-  small = max ([max(abs(pden)), max(abs(pnum)), 1]) * eps;
+  if (isa (r, "single") || isa (p, "single") || isa (k, "single"))
+    small = max ([max(abs(pden)), max(abs(pnum)), 1]) * eps ("single");
+  else
+    small = max ([max(abs(pden)), max(abs(pnum)), 1]) * eps;
+  endif
 
   pnum(abs (pnum) < small) = 0;
   pden(abs (pden) < small) = 0;
diff -r 2b458dfe31ae -r df9519e9990c scripts/sparse/normest.m
--- a/scripts/sparse/normest.m	Mon May 12 01:35:30 2008 +0200
+++ b/scripts/sparse/normest.m	Mon May 12 22:57:11 2008 +0200
@@ -32,8 +32,14 @@
   if (nargin < 2)
     tol = 1e-6;
   endif
-  if (tol < eps)
-    tol = eps
+  if (isa (A, "single"))
+    if (tol < eps ("single"))
+      tol = eps ("single");
+    endif
+  else
+    if (tol < eps)
+      tol = eps
+    endif
   endif
   if (ndims(A) != 2)
     error ("normest: A must be a matrix");
diff -r 2b458dfe31ae -r df9519e9990c scripts/specfun/erfinv.m
--- a/scripts/specfun/erfinv.m	Mon May 12 01:35:30 2008 +0200
+++ b/scripts/specfun/erfinv.m	Mon May 12 22:57:11 2008 +0200
@@ -35,7 +35,11 @@
   endif
 
   maxit = 100;
-  tol = eps;
+  if (isa (x, "single"))
+    tol = eps ("single");
+  else
+    tol = eps;
+  endif
 
   iterations = 0;
 
diff -r 2b458dfe31ae -r df9519e9990c scripts/statistics/distributions/betainv.m
--- a/scripts/statistics/distributions/betainv.m	Mon May 12 01:35:30 2008 +0200
+++ b/scripts/statistics/distributions/betainv.m	Mon May 12 22:57:11 2008 +0200
@@ -62,29 +62,36 @@
       y = a / (a + b) * ones (size (k));
     endif
     x = x (k);
-    l = find (y < eps);
-    if (any (l))
-      y(l) = sqrt (eps) * ones (length (l), 1);
+
+    if (isa (y, "single"))
+      myeps = eps ("single");
+    else
+      myeps = eps;
     endif
-    l = find (y > 1 - eps);
+
+    l = find (y < myeps);
     if (any (l))
-      y(l) = 1 - sqrt (eps) * ones (length (l), 1);
+      y(l) = sqrt (myeps) * ones (length (l), 1);
+    endif
+    l = find (y > 1 - myeps);
+    if (any (l))
+      y(l) = 1 - sqrt (myeps) * ones (length (l), 1);
     endif
 
     y_old = y;
     for i = 1 : 10000
       h     = (betacdf (y_old, a, b) - x) ./ betapdf (y_old, a, b);
       y_new = y_old - h;
-      ind   = find (y_new <= eps);
+      ind   = find (y_new <= myeps);
       if (any (ind))
         y_new (ind) = y_old (ind) / 10;
       endif
-      ind = find (y_new >= 1 - eps);
+      ind = find (y_new >= 1 - myeps);
       if (any (ind))
         y_new (ind) = 1 - (1 - y_old (ind)) / 10;
       endif
       h = y_old - y_new;
-      if (max (abs (h)) < sqrt (eps))
+      if (max (abs (h)) < sqrt (myeps))
         break;
       endif
       y_old = y_new;
diff -r 2b458dfe31ae -r df9519e9990c scripts/statistics/distributions/gaminv.m
--- a/scripts/statistics/distributions/gaminv.m	Mon May 12 01:35:30 2008 +0200
+++ b/scripts/statistics/distributions/gaminv.m	Mon May 12 22:57:11 2008 +0200
@@ -63,21 +63,28 @@
       y = a * b * ones (size (k));
     endif
     x = x (k);
-    l = find (x < eps);
+
+    if (isa (x, "single"))
+      myeps = eps ("single");
+    else
+      myeps = eps;
+    endif
+
+    l = find (x < myeps);
     if (any (l))
-      y(l) = sqrt (eps) * ones (length (l), 1);
+      y(l) = sqrt (myeps) * ones (length (l), 1);
     endif
 
     y_old = y;
     for i = 1 : 100
       h     = (gamcdf (y_old, a, b) - x) ./ gampdf (y_old, a, b);
       y_new = y_old - h;
-      ind   = find (y_new <= eps);
+      ind   = find (y_new <= myeps);
       if (any (ind))
         y_new (ind) = y_old (ind) / 10;
         h = y_old - y_new;
       endif
-      if (max (abs (h)) < sqrt (eps))
+      if (max (abs (h)) < sqrt (myeps))
         break;
       endif
       y_old = y_new;
diff -r 2b458dfe31ae -r df9519e9990c scripts/statistics/distributions/kolmogorov_smirnov_cdf.m
--- a/scripts/statistics/distributions/kolmogorov_smirnov_cdf.m	Mon May 12 01:35:30 2008 +0200
+++ b/scripts/statistics/distributions/kolmogorov_smirnov_cdf.m	Mon May 12 22:57:11 2008 +0200
@@ -50,7 +50,11 @@
   endif
 
   if (nargin == 1)
-    tol = eps;
+    if (isa (x, "single"))
+      tol = eps ("single");
+    else
+      tol = eps;
+    endif
   else
     if (! isscalar (tol) || ! (tol > 0))
       error ("kolmogorov_smirnov_cdf: tol has to be a positive scalar");
diff -r 2b458dfe31ae -r df9519e9990c scripts/statistics/tests/manova.m
--- a/scripts/statistics/tests/manova.m	Mon May 12 01:35:30 2008 +0200
+++ b/scripts/statistics/tests/manova.m	Mon May 12 22:57:11 2008 +0200
@@ -83,7 +83,12 @@
   n_w = n - k;
 
   l = real (eig (SSB / SSW));
-  l (l < eps) = 0;
+
+  if (isa (l, "single"))
+    l (l < eps ("single")) = 0;
+  else
+    l (l < eps) = 0;
+  endif
 
   ## Wilks' Lambda
   ## =============
diff -r 2b458dfe31ae -r df9519e9990c src/ChangeLog
--- a/src/ChangeLog	Mon May 12 01:35:30 2008 +0200
+++ b/src/ChangeLog	Mon May 12 22:57:11 2008 +0200
@@ -1,3 +1,15 @@
+2008-05-21  David Bateman  <dbateman@free.fr>
+
+	* DLD-FUNCTIONS/sqrt.m: Replace DBL_* with FLT_* for single
+	precision types.
+	* data.cc (static octave_value fill_matrix (const
+	octave_value_list&, double, float, const char *)): Add function
+	with additional argument to allow for different valid for double
+	and single precision.
+	(Finf, FNaN, FNA, Frealmax, Frealmin): Use it here.
+	(Feps): Modify behavior for a single numerical argument to give
+	difference to next largest value in the class of the type passed.
+
 2008-05-21  John W. Eaton  <jwe@octave.org>
 
 	* pt-idx.h (tree_index_expression::tree_index_expression (int, int)): 
diff -r 2b458dfe31ae -r df9519e9990c src/DLD-FUNCTIONS/sqrtm.cc
--- a/src/DLD-FUNCTIONS/sqrtm.cc	Mon May 12 01:35:30 2008 +0200
+++ b/src/DLD-FUNCTIONS/sqrtm.cc	Mon May 12 22:57:11 2008 +0200
@@ -270,7 +270,7 @@
 		    normX = getmax (normX, abs (X(i,j)));
 		  }
 
-	      if (imagX < normX * 100 * DBL_EPSILON)
+	      if (imagX < normX * 100 * FLT_EPSILON)
 		retval(0) = real (X);
 	      else
 		retval(0) = X;
@@ -321,11 +321,11 @@
 
 	  if (nargout < 2)
 	    {
-	      if (err > 100*(minT+DBL_EPSILON)*n)
+	      if (err > 100*(minT+FLT_EPSILON)*n)
 		{
 		  if (minT == 0.0)
 		    error ("sqrtm: A is singular, sqrt may not exist");
-		  else if (minT <= sqrt (DBL_MIN))
+		  else if (minT <= sqrt (FLT_MIN))
 		    error ("sqrtm: A is nearly singular, failed to find sqrt");
 		  else
 		    error ("sqrtm: failed to find sqrt");
diff -r 2b458dfe31ae -r df9519e9990c src/data.cc
--- a/src/data.cc	Mon May 12 01:35:30 2008 +0200
+++ b/src/data.cc	Mon May 12 22:57:11 2008 +0200
@@ -2662,6 +2662,87 @@
 }
 
 static octave_value
+fill_matrix (const octave_value_list& args, double val, float fval, 
+	     const char *fcn)
+{
+  octave_value retval;
+
+  int nargin = args.length ();
+
+  oct_data_conv::data_type dt = oct_data_conv::dt_double;
+
+  dim_vector dims (1, 1);
+  
+  if (nargin > 0 && args(nargin-1).is_string ())
+    {
+      std::string nm = args(nargin-1).string_value ();
+      nargin--;
+
+      dt = oct_data_conv::string_to_data_type (nm);
+
+      if (error_state)
+	return retval;
+    }
+
+  switch (nargin)
+    {
+    case 0:
+      break;
+
+    case 1:
+      get_dimensions (args(0), fcn, dims);
+      break;
+
+    default:
+      {
+	dims.resize (nargin);
+
+	for (int i = 0; i < nargin; i++)
+	  {
+	    dims(i) = args(i).is_empty () ? 0 : args(i).idx_type_value ();
+
+	    if (error_state)
+	      {
+		error ("%s: expecting scalar integer arguments", fcn);
+		break;
+	      }
+	  }
+      }
+      break;
+    }
+
+  if (! error_state)
+    {
+      dims.chop_trailing_singletons ();
+
+      check_dimensions (dims, fcn);
+
+      // Note that automatic narrowing will handle conversion from
+      // NDArray to scalar.
+
+      if (! error_state)
+	{
+	  switch (dt)
+	    {
+	    case oct_data_conv::dt_single:
+	      retval = FloatNDArray (dims, fval);
+	      break;
+
+	    case oct_data_conv::dt_double:
+	      retval = NDArray (dims, val);
+	      break;
+
+	    default:
+	      error ("%s: invalid class name", fcn);
+	      break;
+	    }
+	}
+    }
+
+  return retval;
+}
+
+static octave_value
 fill_matrix (const octave_value_list& args, double val, const char *fcn)
 {
   octave_value retval;
@@ -2724,7 +2805,7 @@
 	  switch (dt)
 	    {
 	    case oct_data_conv::dt_single:
-	      retval = FloatNDArray (dims, val);
+	      retval = FloatNDArray (dims, static_cast <float> (val));
 	      break;
 
 	    case oct_data_conv::dt_double:
@@ -2805,7 +2886,7 @@
 	  switch (dt)
 	    {
 	    case oct_data_conv::dt_single:
-	      retval = FloatComplexNDArray (dims, static_cast <FloatComplex> (val));
+	      retval = FloatComplexNDArray (dims, static_cast<FloatComplex> (val));
 	      break;
 
 	    case oct_data_conv::dt_double:
@@ -2933,7 +3014,8 @@
 @samp{\"double\"}.  The default is @samp{\"double\"}.\n\
 @end deftypefn")
 {
-  return fill_matrix (args, lo_ieee_inf_value (), "Inf");
+  return fill_matrix (args, lo_ieee_inf_value (), 
+		      lo_ieee_float_inf_value (), "Inf");
 }
 
 DEFALIAS (inf, Inf);
@@ -2965,7 +3047,8 @@
 @samp{\"double\"}.  The default is @samp{\"double\"}.\n\
 @end deftypefn")
 {
-  return fill_matrix (args, lo_ieee_nan_value (), "NaN");
+  return fill_matrix (args, lo_ieee_nan_value (), 
+		      lo_ieee_float_nan_value (), "NaN");
 }
 
 DEFALIAS (nan, NaN);
@@ -3026,11 +3109,105 @@
  $2.2204\\times10^{-16}$.\n\
 @end tex\n\
 @end iftex\n\
+for double precision and\n\
+@ifinfo\n\
+ 1.1921e-07.\n\
+@end ifinfo\n\
+@iftex\n\
+@tex\n\
+ $1.1921\\times10^{-7}$.\n\
+@end tex\n\
+@end iftex\n\
+for single precision. Given a single argument @var{x}, return the\n\
+distance between @var{x} and the next largest value.\n\
 @end deftypefn")
 {
-  return fill_matrix (args, DBL_EPSILON, "eps");
+  int nargin = args.length ();
+  octave_value retval;
+
+  if (nargin == 1 && ! args(0).is_string ())
+    {
+      if (args(0).is_single_type ())
+	{
+	  float val = args(0).float_value ();
+
+	  if (! error_state)
+	    {
+	      val  = ::fabsf(val);
+	      if (xisnan (val) || xisinf (val))
+		retval = fill_matrix (octave_value ("single"), 
+				      lo_ieee_nan_value (), 
+				      lo_ieee_float_nan_value (), "eps");
+	      else if (val < FLT_MIN)
+		retval = fill_matrix (octave_value ("single"), 0e0, 
+				      powf (2.0, -149e0), "eps");
+	      else
+		{
+		  int expon;
+		  frexpf (val, &expon);
+		  val = std::pow (static_cast <float> (2.0), 
+				  static_cast <float> (expon - 24));
+		  retval = fill_matrix (octave_value ("single"), DBL_EPSILON, 
+					val, "eps");
+		}
+	    }
+	}
+      else
+	{
+	  double val = args(0).double_value ();
+
+	  if (! error_state)
+	    {
+	      val  = ::fabs(val);
+	      if (xisnan (val) || xisinf (val))
+		retval = fill_matrix (octave_value_list (), 
+				      lo_ieee_nan_value (), 
+				      lo_ieee_float_nan_value (), "eps");
+	      else if (val < DBL_MIN)
+		retval = fill_matrix (octave_value_list (),
+				      pow (2.0, -1074e0), 0e0, "eps");
+	      else
+		{
+		  int expon;
+		  frexp (val, &expon);
+		  val = std::pow (static_cast <double> (2.0), 
+				  static_cast <double> (expon - 53));
+		  retval = fill_matrix (octave_value_list (), val, 
+					FLT_EPSILON, "eps");
+		}
+	    }
+	}
+    }
+  else
+    retval = fill_matrix (args, DBL_EPSILON, FLT_EPSILON, "eps");
+
+  return retval;
 }
 
+/*
+
+%!assert(eps(1/2),2^(-53))
+%!assert(eps(1),2^(-52))
+%!assert(eps(2),2^(-51))
+%!assert(eps(realmax),2^971)
+%!assert(eps(0),2^(-1074))
+%!assert(eps(realmin/2),2^(-1074))
+%!assert(eps(realmin/16),2^(-1074))
+%!assert(eps(Inf),NaN)
+%!assert(eps(NaN),NaN)
+%!assert(eps(single(1/2)),single(2^(-24)))
+%!assert(eps(single(1)),single(2^(-23)))
+%!assert(eps(single(2)),single(2^(-22)))
+%!assert(eps(realmax('single')),single(2^104))
+%!assert(eps(single(0)),single(2^(-149)))
+%!assert(eps(realmin('single')/2),single(2^(-149)))
+%!assert(eps(realmin('single')/16),single(2^(-149)))
+%!assert(eps(single(Inf)),single(NaN))
+%!assert(eps(single(NaN)),single(NaN))
+
+*/
+
+
 DEFUN (pi, args, ,
   "-*- texinfo -*-\n\
 @deftypefn {Built-in Function} {} pi (@var{x})\n\
@@ -3072,7 +3249,7 @@
 @seealso{realmin}\n\
 @end deftypefn")
 {
-  return fill_matrix (args, DBL_MAX, "realmax");
+  return fill_matrix (args, DBL_MAX, FLT_MAX, "realmax");
 }
 
 DEFUN (realmin, args, ,
@@ -3096,7 +3273,7 @@
 @seealso{realmax}\n\
 @end deftypefn")
 {
-  return fill_matrix (args, DBL_MIN, "realmin");
+  return fill_matrix (args, DBL_MIN, FLT_MIN, "realmin");
 }
 
 DEFUN (I, args, ,
@@ -3136,7 +3313,8 @@
 to the special constant used to designate missing values.\n\
 @end deftypefn")
 {
-  return fill_matrix (args, lo_ieee_na_value (), "NA");
+  return fill_matrix (args, lo_ieee_na_value (), 
+		      lo_ieee_float_na_value (), "NA");
 }
 
 DEFUN (false, args, ,