Generating a particular sequnce of numbers

2 次查看（过去 30 天）
Adi gahlawat
Adi gahlawat 2013-7-27
Hi,
given a variable natural number d, I'm trying to generate a sequence of the form:
[1 2 1 3 2 1 4 3 2 1.......d d-1 d-2......3 2 1].
I don't want to use for loop for this process, does anyone know a better (faster) method. I tried the colon operator without any success.
Thank you.
Adi

请先登录，再回答此问题。

采纳的回答

Azzi Abdelmalek
Azzi Abdelmalek 2013-7-27
编辑：Azzi Abdelmalek 2013-7-27
d=4
cell2mat(arrayfun(@(x) x:-1:1,1:d,'un',0))
  7 个评论
显示 5更早的评论

请先登录，再进行评论。

更多回答（6 个）

Roger Stafford
Roger Stafford 2013-7-27
Here's another method to try:
N = d*(d+1)/2;
A = zeros(1,N);
n = 1:d;
A((n.^2-n+2)/2) = n;
A = cumsum(A)-(1:N)+1;
  1 个评论
Adi gahlawat
Adi gahlawat 2013-7-27
编辑：Adi gahlawat 2013-7-27
Hi Roger,
your method is excellent. It's about 2 times faster than my for loop based code. Much obliged.
Adi

请先登录，再进行评论。

Azzi Abdelmalek
Azzi Abdelmalek 2013-7-28
编辑：Azzi Abdelmalek 2013-7-28
Edit
This is twice faster then Stafford's answer
A4=zeros(1,d*(d+1)/2); % Pre-allocate
c=0;
for k=1:d
A4(c+1:c+k)=k:-1:1;
c=c+k;
end
  1 个评论
Jan
Jan 2013-7-28
编辑：Jan 2013-7-28
Yes, this is exactly the kind of simplicity, which runs fast. While the one-liners with anonymous functions processed by cellfun or arrayfun look sophisticated, such basic loops hit the point. +1
I'd replace sum(1:d) by: d*(d+1)/2 . Anbd you can omit idx.

请先登录，再进行评论。

Richard Brown
Richard Brown 2013-7-29
Even faster:
k = 1;
n = d*(d+1)/2;
out = zeros(n, 1);
for i = 1:d
for j = i:-1:1
out(k) = j;
k = k + 1;
end
end
  7 个评论
显示 5更早的评论
Jan
Jan 2013-7-29
编辑：Jan 2013-7-29
Under R2011b I get for d=1000 and 500 repetitions:
Elapsed time is 3.466296 seconds. Azzi's loop
Elapsed time is 3.765340 seconds. Richard's double loop
Elapsed time is 1.897343 seconds. C-Mex (see my answer)
Richard Brown
Richard Brown 2013-7-29
编辑：Richard Brown 2013-7-29
I checked again, and I agree with Azzi. My method was running faster because of another case I had in between his and mine. The JIT was doing some kind of unanticipated optimisation between cases.
I get similar orders of magnitude results to Azzi for R2012a if I remove that case, and if I run in R2013a (Linux), his method is twice as fast.
Shame, I like it when JIT brings performance of completely naive loops up to vectorised speed :)

请先登录，再进行评论。

Jan
Jan 2013-7-29
An finally the C-Mex:
#include "mex.h"
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray*prhs[]) {
mwSize d, i, j;
double *r;
d = (mwSize) mxGetScalar(prhs[0]);
plhs[0] = mxCreateDoubleMatrix(1, d * (d + 1) / 2, mxREAL);
r = mxGetPr(plhs[0]);
for (i = 1; i <= d; i++) {
for (j = i; j != 0; *r++ = j--) ;
}
}
And if your number d can be limited to 65535, the times shrink from 1.9 to 0.34 seconds:
#include "mex.h"
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray*prhs[]) {
uint16_T d, i, j, *r;
d = (uint16_T) mxGetScalar(prhs[0]);
plhs[0] = mxCreateNumericMatrix(1, d * (d + 1) / 2, mxUINT16_CLASS, mxREAL);
r = (uint16_T *) mxGetData(plhs[0]);
for (i = 1; i <= d; i++) {
for (j = i; j != 0; *r++ = j--) ;
}
}
For UINT32 0.89 seconds are required.
  1 个评论
Richard Brown
Richard Brown 2013-7-29
Nice. I imagine d would be limited to less than 65535, that's a pretty huge vector otherwise

请先登录，再进行评论。

Richard Brown
Richard Brown 2013-7-29
编辑：Richard Brown 2013-7-29
Also comparable, but not (quite) faster
n = 1:(d*(d+1)/2);
a = ceil(0.5*(-1 + sqrt(1 + 8*n)));
out = a.*(a + 1)/2 - n + 1;
  3 个评论
显示 1更早的评论
Jan
Jan 2013-7-29
@Richard: How did you find this formula?
Richard Brown
Richard Brown 2013-7-29
If you look at the sequence, and add 0, 1, 2, 3, 4 ... you get
n: 1 2 3 4 5 6 7 8 9 10
1 3 3 6 6 6 10 10 10 10
Note that these are the triangular numbers, and that the triangular numbers 1, 3, 6, 10 appear in their corresponding positions, The a-th triangular number is given by
n = a (a + 1) / 2
So if you solve this quadratic for a where n is a triangular number, you get the index of the triangular number. If you do this for a value of n in between two triangular numbers, you can round this up, and invert the formula to get the nearest triangular number above (which is what the sequence is). Finally, you just subtract the sequence 0, 1, 2, ... to recover the original one.

请先登录，再进行评论。

Andrei Bobrov
Andrei Bobrov 2013-7-27
编辑：Andrei Bobrov 2013-7-30
out = nonzeros(triu(toeplitz(1:d)));
or
out = bsxfun(@minus,1:d,(0:d-1)');
out = out(out>0);
or
z = 1:d;
z2 = cumsum(z);
z1 = z2 - z + 1;
for jj = d:-1:1
out(z1(jj):z2(jj)) = jj:-1:1;
end
or
out = ones(d*(d+1)/2,1);
ii = cumsum(d:-1:1) - (d:-1:1) + 1;
out(ii(2:end)) = 1-d : -1;
out = flipud(cumsum(out));

类别

在 Help Center 和 File Exchange 中查找有关 Performance and Memory 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by