所引用的已删除页面萨姆的回答 https://stackoverflow.com/a/26297279/7328782还是由 Wayback Machine 存档 https://web.archive.org/web/20121120071523/http://www.mathworks.com/support/solutions/en/data/1-4FLI96/index.html?solution=1-4FLI96。幸运的是，即使是附带的 M 文件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 and stop.
首先它计算中间有多少步start and 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
这就是向量中的 lasat 值的原因A问题中实际上有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 value.
最后一步，对于奇数大小的数组，单独计算中间值以确保它恰好位于两个端点的中间：
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
