## [nel, final] = xcolon (base, limit, inc = 1)
## Nth attempt to compute the same number of elements and final value
## as Matlab's colon operator. Does not handle Inf/NaN values or zero
## increments. That is expected to be done separately.
## I don't really know what to use for the tolerance, but chose
##
## eps (abs (final) + abs (inc)) ## INC is never zero; avoids eps (0)
##
## instead of just
##
## eps (final)
##
## because if FINAL is zero, then the tolerance would be extremely
## small (approximately 4.9e-324) which seems inappropriate for a
## range like -3:0.3:0, for example.
## NOTE: If the base and increment are integers and limit is not,
## this implementation allows the final value to be set to the
## specified limit. That appears to be compatible with Matlab but
## might not be desirable depending on the intended use of the range
## values. It would be bad for indexing arrays, but might be useful
## for plotting?
##
## We could:
##
## * Warn when any range contains floating point values. But that
## might be noisy and disabled by nearly everyone?
##
## * Warn when the base and increment are integers and the limit is
## not.
##
## * Warn only when base and increment are integers and a
## non-integer limit is selected because it falls within the
## tolerance? We would still be attempting compatible numerical
## behavior with this option.
##
## * Disregard the tolerance when the base and increment are
## integers so that the generated values are all integers when the
## base and increment are integers. This choice results in
## numerical results that are incompatible with (undocumented)
## Matlab behavior.
function [nel, final] = xcolon (base, limit, inc = 1)
nel = (limit - base) / inc + 1;
if (ceil (nel) - nel <= eps (nel))
nel = ceil (nel);
else
nel = floor (nel);
endif
final = base + (nel - 1) * inc;
if (inc > 0)
if (final > limit)
## WTF, no tolerance allowed when initial estimate is outside
## the specified range limit? If we do check against eps here,
## we get an incompatible result for 1:(1001/250)/(1/250)
[nel, final] = revise_nel_down (base, limit, inc, nel);
elseif (final < limit)
if (abs (final - limit) <= eps (abs (final) + abs (inc)))
final = limit;
else
[nel, final] = revise_nel_up (base, limit, inc, nel);
endif
endif
elseif (inc < 0)
if (final < limit)
## WTF, no tolerance allowed when initial estimate is outside
## the specified range limit?
[nel, final] = revise_nel_down (base, limit, inc, nel);
elseif (final > limit)
if (abs (final - limit) <= eps (abs (final) + abs (inc)))
final = limit;
else
[nel, final] = revise_nel_up (base, limit, inc, nel);
endif
endif
endif
endfunction
## REVISE_NEL_UP is only used if previous estimate of final value was
## inside the specified range limit by more than tolerance. One call
## should be sufficient to reach the limit but we allow subsequent
## revisions if needed.
function [nel, final] = revise_nel_up (base, limit, inc, nel)
nel = nel + 1;
final = base + (nel - 1) * inc;
if (inc > 0)
if (final > limit)
if (abs (final - limit) <= eps (abs (final) + abs (inc)))
final = limit;
else
## Previous value was not accepted and additional element
## pushed us outside the specified range limit by more than
## tolerance. Accept previous value.
nel = nel - 1;
final = base + (nel - 1) * inc;
endif
elseif (final < limit)
if (abs (final - limit) <= eps (abs (final) + abs (inc)))
final = limit;
else
[nel, final] = revise_nel_up (base, limit, inc, nel);
endif
endif
elseif (inc < 0)
if (final < limit)
if (abs (final - limit) <= eps (abs (final) + abs (inc)))
final = limit;
else
## Previous value was not accepted and additional element
## pushed us outside the specified range limit by more than
## tolerance. Accept previous value.
nel = nel - 1;
final = base + (nel - 1) * inc;
endif
elseif (final > limit)
if (abs (final - limit) <= eps (abs (final) + abs (inc)))
final = limit;
else
[nel, final] = revise_nel_up (base, limit, inc, nel);
endif
endif
endif
endfunction
## REVISE_NEL_DOWN is only used if previous estimate of final value was
## outside the specified range limit by more than tolerance. One call
## should bring us back inside but we allow subsequent revisions if
## necessary.
function [nel, final] = revise_nel_down (base, limit, inc, nel)
nel = nel - 1;
final = base + (nel - 1) * inc;
if (inc > 0)
if (final > limit)
if (abs (final - limit) <= eps (abs (final) + abs (inc)))
final = limit;
else
[nel, final] = revise_nel_down (base, limit, inc, nel);
endif
elseif (final < limit)
if (abs (final - limit) <= eps (abs (final) + abs (inc)))
final = limit;
else
## Previous estimate of final value was out of specified range
## limit and now we are inside by more than tolerance, so accept
## calculated value.
endif
endif
elseif (inc < 0)
if (final < limit)
if (abs (final - limit) <= eps (abs (final) + abs (inc)))
final = limit;
else
[nel, final] = revise_nel_down (base, limit, inc, nel);
endif
elseif (final > limit)
if (abs (final - limit) <= eps (abs (final) + abs (inc)))
final = limit;
else
## Previous estimate of final value was out of specified range
## limit and now we are inside by more than tolerance, so accept
## calculated value.
endif
endif
endif
endfunction