MATLAB Answers

## Generating a particular sequnce of numbers

Asked by Adi gahlawat

### Adi gahlawat (view profile)

on 27 Jul 2013
Accepted Answer by Azzi Abdelmalek

### Azzi Abdelmalek (view profile)

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

#### 0 Comments

Sign in to comment.

## 7 Answers

### Azzi Abdelmalek (view profile)

Answer by Azzi Abdelmalek

### Azzi Abdelmalek (view profile)

on 27 Jul 2013
Edited by Azzi Abdelmalek

### Azzi Abdelmalek (view profile)

on 27 Jul 2013
Accepted Answer

d=4
cell2mat(arrayfun(@(x) x:-1:1,1:d,'un',0))

Azzi Abdelmalek

### Azzi Abdelmalek (view profile)

on 27 Jul 2013
Adi gahlawat, you can not compare by puting df=3, look at this:
df=1000;
tic
for k=1:500
A1=cell2mat(arrayfun(@(x) x:-1:1,1:df+1,'un',0))';
end
toc
tic
for k=1:500
A2=zeros(df+1,df+1);
for i=1:df+1
A2(i,i:df+1)=1:df+2-i;
end
A2=A2(:); % Converting matrix to a vector
A2=A2(A2~=0); % Removing zeros
end
toc
tic
for k=1:500
N = df*(df+1)/2;
A = zeros(1,N);
n = 1:df;
A((n.^2-n+2)/2) = n;
A = cumsum(A)-(1:N)+1;
end
toc
Elapsed time is 5.617201 seconds. % Azzi's answer
Elapsed time is 10.643052 seconds. % Adi's answer
Elapsed time is 4.755004 seconds. % Stafford's answer
Youssef Khmou

### Youssef Khmou (view profile)

on 27 Jul 2013
ok thank you for the explanation .
Jan

### Jan (view profile)

on 28 Jul 2013
Azzi's for loop approach is faster.

Sign in to comment.

### Roger Stafford (view profile)

Answer by Roger Stafford

### Roger Stafford (view profile)

on 27 Jul 2013

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;

Adi gahlawat

### Adi gahlawat (view profile)

on 27 Jul 2013
Hi Roger,
your method is excellent. It's about 2 times faster than my for loop based code. Much obliged.
Adi

Sign in to comment.

### Azzi Abdelmalek (view profile)

Answer by Azzi Abdelmalek

### Azzi Abdelmalek (view profile)

on 28 Jul 2013
Edited by Azzi Abdelmalek

### Azzi Abdelmalek (view profile)

on 28 Jul 2013

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

Jan

### Jan (view profile)

on 28 Jul 2013
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.

Sign in to comment.

### Richard Brown (view profile)

Answer by Richard Brown

### Richard Brown (view profile)

on 29 Jul 2013

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

Azzi Abdelmalek

### Azzi Abdelmalek (view profile)

on 29 Jul 2013
Almost, the same result
Elapsed time is 22.940850 seconds.
Elapsed time is 16.967270 seconds.
Jan

### Jan (view profile)

on 29 Jul 2013
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 (view profile)

on 29 Jul 2013
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 :)

Sign in to comment.

Answer by Jan

### Jan (view profile)

on 29 Jul 2013

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.

Richard Brown

### Richard Brown (view profile)

on 29 Jul 2013
Nice. I imagine d would be limited to less than 65535, that's a pretty huge vector otherwise

Sign in to comment.

### Richard Brown (view profile)

Answer by Richard Brown

### Richard Brown (view profile)

on 29 Jul 2013
Edited by Richard Brown

### Richard Brown (view profile)

on 29 Jul 2013

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;

Richard Brown

### Richard Brown (view profile)

on 29 Jul 2013
Potentially suffers from floating point errors, but I checked it up to d = 10000 :)
Jan

### Jan (view profile)

on 29 Jul 2013
@Richard: How did you find this formula?
Richard Brown

### Richard Brown (view profile)

on 29 Jul 2013
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.

Sign in to comment.

### Andrei Bobrov (view profile)

Answer by Andrei Bobrov

### Andrei Bobrov (view profile)

on 27 Jul 2013
Edited by Andrei Bobrov

### Andrei Bobrov (view profile)

on 30 Jul 2013

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));

#### 0 Comments

Sign in to comment.