Are there alternative ways to write codes to calculate a_p ?

Here

given

m= 0:n+1

n is the degree of Bspline which can be any value

p = 0:n
k=k0:k1;
k0 =ceil(d-((n+1)/2));
k1 =floor(d+((n+1)/2));
num_points = 1; % Change this if you need more points
% Generate random variable d uniformly in the range (-0.5, 0.5)
d = -0.5 + (0.5 - (-0.5)) * rand(num_points, 1);

mu is the unit step function or heaviside function

I tried to develop the function a_P

function a_p = calculate_ap(n, d)
% n: Degree of the polynomial
% d: Interpolation point as a random variable with |d| < 0.5

% Define the range for p (0 to n)
p_values = 0:n;

% Calculate the lower and upper bounds of k
k0 = ceil(d - (n + 1) / 2);
k1 = floor(d + (n + 1) / 2);

% Define the range for k
k_values = k0:k1;

% Initialize a matrix to store a_p(k) values
a_p = zeros(length(p_values), length(k_values));

% Loop through each value of p to compute a_p for each k
for p_idx = 1:length(p_values)
    p = p_values(p_idx);
    
    % Calculate the binomial coefficient (n choose p)
    binomial_np = nchoosek(n, p);
    
    % Calculate the factorial term 1/n!
    n_factorial_term = 1 / factorial(n);
    
    % Loop through each value of k
    for k_idx = 1:length(k_values)
        k = k_values(k_idx);
        
        % Initialize the sum for this combination of p and k
        sum_result = 0;
        
        % Sum over m from 0 to n+1
        for m = 0:n+1
            % Calculate the binomial coefficient (n+1 choose m)
            binomial_n1m = nchoosek(n+1, m);
            
            % Calculate the term (-1)^m
            sign_term = (-1)^m;
            
            % Calculate the term (n+1)/2 - m
            difference_term = (n+1) / 2 - m;
            
            % Calculate the term (difference_term)^(n-p)
            power_term = difference_term^(n - p);
            
            % Calculate the Heaviside function mu(p - m + (n+1)/2)
            mu_term = heaviside(p - m + (n+1)/2);
            
            % Add the current term to the sum
            sum_result = sum_result + binomial_n1m * sign_term * power_term * mu_term;
        end
        
        % Calculate a_p(k) for the current p and k
        a_p(p_idx, k_idx) = n_factorial_term * binomial_np * sum_result;
    end
end

For solving this function

n= input("Enter the degree:"); %Degree
p=0:n;
a_p = calculate_ap(n,p);
disp(a_p)

solution is :

for n= 3

0.6667 0.6667 0.6667 0.6667 0.6667

-1.0000 -1.0000 -1.0000 -1.0000 -1.0000

0.5000 0.5000 0.5000 0.5000 0.5000

0 0 0 0 0

>> Is this function written correctly does it need any change ? Please could you suggest?

[Torsten]

I wonder why a_p depends on k, but the variable k does not appear in its definition.

Answer 1

[Subhajyoti]

Hi @Rohitashya,

You can try precomputing values and vectorized operations to efficiently calculate the given equation. Additionally, dynamic programming paradigm helps in avoiding repeated calculations and optimizes code efficiency.

Here, in the following implementation, I have utilized precomputed factorial and binomials to compute the equation.

MAX_N = 10000;

% preallocate memory for factorial values upto MAX_N
factorial_values = cumprod(1:MAX_N);
function y = fact(n)
    global MAX_N, factorial_values;
    if n==0
        y = 1;
    else 
        y = factorial_values(n);
    end
end

% preallocate 2-dimentions array to precompute binomial coefficients
binomial_values = zeros(MAX_N, MAX_N);

% precompute binomial coefficients using dynamic programming
for n = 1:MAX_N
    for k = 0:n
        if k == 0 || k == n
            binomial_values(n+1, k+1) = 1;
        else
            binomial_values(n+1, k+1) = binomial_values(n, k) + binomial_values(n, k+1);
        end
    end
end

function nCr(n, r)
    global binomial_values;
    binomial_values(n+1, r+1)
end

function a_p = calculate_ap(n, p)
    ans = 0;
    for m = 0:n+2
         sign_term = (-1)^m;
    
         base = ((n+1)/2 - m);
         power_term = (base)^(n-p);
    
         % Calculate the Heaviside function mu(p - m + (n+1)/2)
         mu = heaviside(p - m + (n+1)/2);
    
         ans = ans + nCr(n+1, m) * sign_term * power_term * mu;
    end

   ans = ans / (fact(n-p) * fact(p));
end

For repetitive calculations, you can try to further optimize the code by memorizing the “heaviside” function and utilising binary exponentiation algorithm.

Additionally, you can refer to the following resources to know more about using common techniques used to improve code run-time:

• https://www.mathworks.com/help/matlab/ref/matlab.lang.memoizedfunction.html
• https://www.mathworks.com/help/matlab/matlab_prog/vectorization.html

[Rohitashya]

Thank you so much

[Rohitashya]

MAX_N = 10000;

% preallocate memory for factorial values upto MAX_N

factorial_values = cumprod(1:MAX_N);

function y = fact(n)

global MAX_N, factorial_values;

if n==0

y = 1;

else

y = factorial_values(n);

end

end could you explain this part a bit ?

[Rohitashya]

Please explain in details about the frst two parts