zerodiv

 function f = zerodiv(x,y,val,wantcaution)
 
 <x>,<y> are matrices of the same size or either or both can be scalars.
 <val> (optional) is the value to use when <y> is 0.  default: 0.
 <wantcaution> (optional) is whether to perform special handling of weird
   cases (see below).  default: 1.

 calculate x./y but use <val> when y is 0.
 if <wantcaution>, then if the absolute value of one or more elements of y is 
   less than 1e-5 (but not exactly 0), we issue a warning and then treat these 
   elements as if they are exactly 0.
 if not <wantcaution>, then we do nothing special.

 note some weird cases:
   if either x or y is [], we return [].
   NaNs in x and y are handled in the usual way.

 history:
 2011/02/02 - in the case that y is not a scalar and wantcaution is set to 0, 
              we were allowing division by 0 to result in Inf and *not* replaced 
              with val as desired.  big mistake.  we have now fixed this.

 example:
 isequalwithequalnans(zerodiv([1 2 3],[1 0 NaN]),[1 0 NaN])

0001 function f = zerodiv(x,y,val,wantcaution)
0002 
0003 % function f = zerodiv(x,y,val,wantcaution)
0004 %
0005 % <x>,<y> are matrices of the same size or either or both can be scalars.
0006 % <val> (optional) is the value to use when <y> is 0.  default: 0.
0007 % <wantcaution> (optional) is whether to perform special handling of weird
0008 %   cases (see below).  default: 1.
0009 %
0010 % calculate x./y but use <val> when y is 0.
0011 % if <wantcaution>, then if the absolute value of one or more elements of y is
0012 %   less than 1e-5 (but not exactly 0), we issue a warning and then treat these
0013 %   elements as if they are exactly 0.
0014 % if not <wantcaution>, then we do nothing special.
0015 %
0016 % note some weird cases:
0017 %   if either x or y is [], we return [].
0018 %   NaNs in x and y are handled in the usual way.
0019 %
0020 % history:
0021 % 2011/02/02 - in the case that y is not a scalar and wantcaution is set to 0,
0022 %              we were allowing division by 0 to result in Inf and *not* replaced
0023 %              with val as desired.  big mistake.  we have now fixed this.
0024 %
0025 % example:
0026 % isequalwithequalnans(zerodiv([1 2 3],[1 0 NaN]),[1 0 NaN])
0027 
0028 % input
0029 if nargin < 4  % need optimal speed so try to bypass in the fully specified case if we can
0030   if ~exist('val','var') || isempty(val)
0031     val = 0;
0032   end
0033   if ~exist('wantcaution','var') || isempty(wantcaution)
0034     wantcaution = 1;
0035   end
0036 else
0037   if isempty(val)
0038     val = 0;
0039   end
0040   if isempty(wantcaution)
0041     wantcaution = 1;
0042   end
0043 end
0044 
0045 % handle special case of y being scalar
0046 if isscalar(y)
0047   if y==0
0048     f = repmat(val,size(x));
0049   else
0050     if wantcaution && abs(y) < 1e-5   % see allzero.m
0051       warning('abs value of divisor is less than 1e-5. we are treating the divisor as 0.');
0052       f = repmat(val,size(x));
0053     else
0054 %REMOVED:
0055 %      if abs(y) < 1e-5 && ~wantcaution
0056 %        warning('abs value of divisor is less than 1e-5. we are treating the divisor as-is.');
0057 %      end
0058       f = x./y;
0059     end
0060   end
0061 else
0062   % do it
0063   bad = y==0;
0064   bad2 = abs(y) < 1e-5;  % see allzero.m
0065   if wantcaution && any(bad2(:) & ~bad(:))
0066     warning('abs value of one or more divisors is less than 1e-5. we are treating these divisors as 0.');
0067   end
0068 %REMOVED:
0069 %  if any(bad2 & ~bad) && ~wantcaution
0070 %    warning('abs value of one or more divisors is less than 1e-5. we are treating the divisors as-is.');
0071 %  end
0072   if wantcaution
0073     y(bad2) = 1;
0074     f = x./y;
0075     f(bad2) = val;
0076   else
0077     y(bad) = 1;
0078     f = x./y;
0079     f(bad) = val;
0080   end
0081 end