MATLAB 的冒号运算符是如何工作的？

Question

如 this answer by Sam Roberts and this other answer by gnovice 中所述，MATLAB 的冒号运算符 (start:step:stop) 创建值向量的方式与 linspace 不同。 Sam Roberts 特别指出：

The colon operator adds increments to the starting point, and subtracts decrements from the end point to reach a middle point. In this way, it ensures that the output vector is as symmetric as possible.

但是，The MathWorks 关于此的官方文档已从其网站上删除。

如果 Sam 的描述是正确的，那么步长的误差不是对称的吗？

>> step = 1/3;
>> C = 0:step:5;
>> diff(C) - step
ans =
   1.0e-15 *
  Columns 1 through 10
         0         0    0.0555   -0.0555   -0.0555    0.1665   -0.2776    0.6106   -0.2776    0.1665
  Columns 11 through 15
    0.1665   -0.2776   -0.2776    0.6106   -0.2776

关于冒号运算符的有趣注意事项：

其值取决于其长度：

>> step = 1/3;
>> C = 0:step:5;
>> X = 0:step:3;
>> C(1:10) - X
ans =
   1.0e-15 *
         0         0         0         0         0   -0.2220         0   -0.4441    0.4441         0

四舍五入会产生重复值:

>> E = 1-eps : eps/4 : 1+eps;
>> E-1
ans =
   1.0e-15 *
   -0.2220   -0.2220   -0.1110         0         0         0         0    0.2220    0.2220

最后一个值是有公差的，如果步长创建了一个刚好超过结束值的值，仍然使用这个结束值：

>> A = 0 : step : 5-2*eps(5)
A =
  Columns 1 through 10
         0    0.3333    0.6667    1.0000    1.3333    1.6667    2.0000    2.3333    2.6667    3.0000
  Columns 11 through 16
    3.3333    3.6667    4.0000    4.3333    4.6667    5.0000
>> A(end) == 5 - 2*eps(5)
ans =
  logical
   1
>> step*15 - 5
ans =
     0

Answer 1

Sam's answer is still archived by the Way Back Machine 引用的已删除页面。幸运的是，即使是附加的 M-file colonop 也在那里。而且这个函数似乎仍然与 MATLAB 的功能相匹配（我使用的是 R2017a）：

>> all(0:step:5 == colonop(0,step,5))
ans =
  logical
   1
>> all(-pi:pi/21:pi == colonop(-pi,pi/21,pi))
ans =
  logical
   1

我将在这里复制函数对一般情况的作用（有一些生成整数向量和处理特殊情况的快捷方式）。我正在用更有意义的变量名替换函数的变量名。输入为 start、step 和 stop。

首先计算 start 和 stop 之间有多少步。如果最后一步超过 stop 超过容差，则不执行：

n = round((stop-start)/step);
tol = 2.0*eps*max(abs(start),abs(stop));
sig = sign(step);
if sig*(start+n*step - stop) > tol
  n = n - 1;
end

这解释了问题中提到的最后一个观察结果。

接下来，它计算最后一个元素的值，并确保它不超过 stop 值，即使它在之前的计算中允许超过它。

last = start + n*step;
if sig*(last-stop) > -tol
   last = stop;
end

这就是为什么问题中向量 A 中的 lasat 值实际上具有 stop 值作为最后一个值的原因。

接下来，它分两部分计算输出数组，如广告所示：数组的左半部分和右半部分独立填充：

out = zeros(1,n+1);
k = 0:floor(n/2);
out(1+k) = start + k*step;
out(n+1-k) = last - k*step;

请注意，它们不是通过递增填充的，而是通过计算一个整数数组并将其乘以步长，就像linspace一样。这解释了问题中关于数组 E 的观察。不同之处在于，数组的右半部分是通过从 last 值中减去这些值来填充的。

作为最后一步，对于 odd-sized 数组，单独计算中间值以确保它恰好位于 half-way 两个端点：

if mod(n,2) == 0
   out(n/2+1) = (start+last)/2;
end

完整功能colonop复制在底部

请注意，分别填充数组的左侧和右侧并不意味着步长的误差应该完全对称。这些误差由舍入误差给出。但是，在 stop 点没有完全达到步长的情况下，它确实有所不同，就像问题中数组 A 的情况一样。在这种情况下，略短的步长是在数组的中间，而不是在末尾：

>> step=1/3;
>> A = 0 : step : 5-2*eps(5);
>> A/step-(0:15)
ans =
   1.0e-14 *
  Columns 1 through 10
         0         0         0         0         0         0         0   -0.0888   -0.4441   -0.5329
  Columns 11 through 16
   -0.3553   -0.3553   -0.5329   -0.5329   -0.3553   -0.5329

但即使在准确到达 stop 点的情况下，中间也会累积一些额外的错误。以问题中的数组 C 为例。 linspace:

不会发生此错误累积

C = 0:1/3:5;
lims = eps(C);
subplot(2,1,1)
plot(diff(C)-1/3,'o-')
hold on
plot(lims,'k:')
plot(-lims,'k:')
plot([1,15],[0,0],'k:')
ylabel('error')
title('0:1/3:5')
L = linspace(0,5,16);
subplot(2,1,2)
plot(diff(L)-1/3,'x-')
hold on
plot(lims,'k:')
plot(-lims,'k:')
plot([1,15],[0,0],'k:')
title('linspace(0,5,16)')
ylabel('error')

colonop:

function out = colonop(start,step,stop)
% COLONOP  Demonstrate how the built-in a:d:b is constructed.
%
%   v = colonop(a,b) constructs v = a:1:b.
%   v = colonop(a,d,b) constructs v = a:d:b.
%
%   v = a:d:b is not constructed using repeated addition.  If the
%   textual representation of d in the source code cannot be
%   exactly represented in binary floating point, then repeated
%   addition will appear to have accumlated roundoff error.  In
%   some cases, d may be so small that the floating point number
%   nearest a+d is actually a.  Here are two imporant examples.
%
%   v = 1-eps : eps/4 : 1+eps is the nine floating point numbers
%   closest to v = 1 + (-4:1:4)*eps/4.  Since the spacing of the
%   floating point numbers between 1-eps and 1 is eps/2 and the
%   spacing between 1 and 1+eps is eps,
%   v = [1-eps 1-eps 1-eps/2 1 1 1 1 1+eps 1+eps].
%
%   Even though 0.01 is not exactly represented in binary,
%   v = -1 : 0.01 : 1 consists of 201 floating points numbers
%   centered symmetrically about zero.
%
%   Ideally, in exact arithmetic, for b > a and d > 0,
%   v = a:d:b should be the vector of length n+1 generated by
%   v = a + (0:n)*d where n = floor((b-a)/d).
%   In floating point arithmetic, the delicate computatations
%   are the value of n, the value of the right hand end point,
%   c = a+n*d, and symmetry about the mid-point.

if nargin < 3
    stop = step;
    step = 1;
end

tol = 2.0*eps*max(abs(start),abs(stop));
sig = sign(step);

% Exceptional cases.

if ~isfinite(start) || ~isfinite(step) || ~isfinite(stop)
   out = NaN;
   return
elseif step == 0 || start < stop && step < 0 || stop < start && step > 0
   % Result is empty.
   out = zeros(1,0);
   return
end

% n = number of intervals = length(v) - 1.

if start == floor(start) && step == 1
   % Consecutive integers.
   n = floor(stop) - start;
elseif start == floor(start) && step == floor(step)
   % Integers with spacing > 1.
   q = floor(start/step);
   r = start - q*step;
   n = floor((stop-r)/step) - q;
else
   % General case.
   n = round((stop-start)/step);
   if sig*(start+n*step - stop) > tol
      n = n - 1;
   end
end

% last = right hand end point.

last = start + n*step;
if sig*(last-stop) > -tol
   last = stop;
end

% out should be symmetric about the mid-point.

out = zeros(1,n+1);
k = 0:floor(n/2);
out(1+k) = start + k*step;
out(n+1-k) = last - k*step;
if mod(n,2) == 0
   out(n/2+1) = (start+last)/2;
end

MATLAB 的冒号运算符是如何工作的？

How does MATLAB's colon operator work?

matlab

vector

colon