Mercurial > octave-nkf
changeset 7795:df9519e9990c
Handle single precision eps values
line wrap: on
line diff
--- 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.
--- 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
--- 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);
--- 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
--- 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
--- 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);
--- 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);
--- 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
--- 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
--- 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,:);
--- 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
--- 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);
--- 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)
--- 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));
--- 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
--- 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]);
--- 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);
--- 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)
--- 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)
--- 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;
--- 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
--- 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
--- 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
--- 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 ();
--- 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
--- 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);
--- 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);
--- 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);
--- 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);
--- 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))));
--- 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
--- 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)
--- 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;
--- 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");
--- 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;
--- 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;
--- 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;
--- 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");
--- 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 ## =============
--- 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)):
--- 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");
--- 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, ,