Make for loop more efficient

I have a function with the following for-loop that seems to be pretty efficient. However, I'm wondering if this is the most efficient implementation, or is there a faster way to do this computation. There is a helper function too "Pnm_index" that finds the proper index to store the computations in the final output variable Pnms. I've included that function as well below. Since this helper function is so simple I'm also wondering if eliminating it will provide any speed up. Any insight is appreciated!

tic
s = sin(3*pi/7);
c = cos(3*pi/7);
p = 5; % input variable
Pnms = zeros(1,((p+1)*(p+2))*(2^-1)); % preallocate the output
Pn0 = 1.0; % initialize
fact = 1.0; % initialize
for m = 0:p
    Pnm = Pn0;
    i = Pnm_index(m, m); % find the proper index, see function below
    Pnms(i) = Pnm; % store the output
    tmp1 = Pnm;
    Pnm = c * (2 * m + 1) * tmp1;
    for n = m+1:p
        i = Pnm_index(n, m); % find the proper index, see function below
        Pnms(i) = Pnm; % store the output
        tmp2 = tmp1;
        tmp1 = Pnm;
        Pnm = (c * (2 * n + 1) * tmp1 - (n + m) * tmp2) / (n - m + 1);
    end
    Pn0 = -Pn0 * fact * s;
    fact = fact + 2;
end
toc

Elapsed time is 0.042007 seconds.

function out = Pnm_index(n,m)

out = (n*(n+1))*(2^-1) + m + 1;

end

[Epsilon]

Hi Christopher,

Eliminating the helper function Pnm_index(n,m) does indeed increase the speed.

%Original:
function out = Pnm_index(n,m)

out = (n*(n+1))*(2^-1) + m + 1;

end

tic
s = sin(3*pi/7);
c = cos(3*pi/7);
p = 5; % input variable
Pnms = zeros(1,((p+1)*(p+2))*(2^-1)); % preallocate the output
Pn0 = 1.0; % initialize
fact = 1.0; % initialize
for m = 0:p
    Pnm = Pn0;
    i = Pnm_index(m, m); % find the proper index, see function below
    Pnms(i) = Pnm; % store the output
    tmp1 = Pnm;
    Pnm = c * (2 * m + 1) * tmp1;
    for n = m+1:p
        i = Pnm_index(n, m); % find the proper index, see function below
        Pnms(i) = Pnm; % store the output
        tmp2 = tmp1;
        tmp1 = Pnm;
        Pnm = (c * (2 * n + 1) * tmp1 - (n + m) * tmp2) / (n - m + 1);
    end
    Pn0 = -Pn0 * fact * s;
    fact = fact + 2;
end
toc

Elapsed time is 0.017910 seconds.

tic
s = sin(3*pi/7);
c = cos(3*pi/7);
p = 5; % input variable
Pnms = zeros(1, ((p+1)*(p+2))*(2^-1)); % preallocate the output
Pn0 = 1.0; % initialize
fact = 1.0; % initialize

for m = 0:p
    Pnm = Pn0;
    % Calculate index directly
    i = (m*(m+1))/2 + m + 1;
    Pnms(i) = Pnm; % store the output
    tmp1 = Pnm;
    Pnm = c * (2 * m + 1) * tmp1;
    for n = m+1:p
        % Calculate index directly
        i = (n*(n+1))/2 + m + 1;
        Pnms(i) = Pnm; % store the output
        tmp2 = tmp1;
        tmp1 = Pnm;
        Pnm = (c * (2 * n + 1) * tmp1 - (n + m) * tmp2) / (n - m + 1);
    end
    Pn0 = -Pn0 * fact * s;
    fact = fact + 2;
end
toc

Elapsed time is 0.005703 seconds.

[Walter Roberson]

We can see from the below that the values calculated by the form loop are nonlinear, and so the inner for loop cannot be vectorized.

Pi = sym(pi);
s = sin(3*Pi/7);
c = cos(3*Pi/7);
syms Pnm
syms Pnms
m = 1;
tmp1 = Pnm;
for n = 2:5
    i = (n*(n+1))/2 + m + 1;
    Pnms(i) = Pnm; % store the output
    tmp2 = tmp1;
    tmp1 = Pnm;
    Pnm = (c * (2 * n + 1) * tmp1 - (n + m) * tmp2) / (n - m + 1);
end
Pnms(:)

ans = 

●