MATLAB code for finding certain coefficients

I am trying to input

In the summation, s are distinct odd primes such that

How can it be written in MATLAB? Here values can be 0 also. All the other symbols involved in the expression are natural numbers and r can be any suitable value. I couldn't get any idea.

[Walter Roberson]

I do not understand what you mean by that "such that" ?

[Walter Roberson]

The summation is over m, but m is not used in the terms ??

[John D'Errico]

I found it immensely confusing. I even considered trying to answer the question until I tried to read it.

I assume that maybe the sum is over all partitions of the number m, so all possible sums of integers that sum to m? So perhaps for a fixed m, we compute the sum?

Is the denominator in that sum composed of binomial coefficients? Or are they Legendre symbols? (That would seem to make little sense, but the syntax could work.) Is that a factorial symbol in there?

A mathematical formula written completely out of context is sometimes just a set of meaningless characters.

[Walter Roberson]

Integer partions is an interesting idea -- I would not have thought of that!

[Bruno Luong]

I guess (pi|li) meaning pi is one of the (odd) prime factor of li. So li/pi is just standard division.

The "distinct" part I don't know, what if there is no such distincts r-tuplets odd primes? Remove the term from the sum?

This formula begs for some explanation.

[Sehra Sahu]

Yes, @Bruno, you are right. By 'distinct' I mean that if then we do not consider the same prime to divide any other . For example, if we are able to split , the summand concerning this partition will be . Even if you consider other primes in this figure, they will be like in the denominator, which is equal to 1. But we cannot split . If you consider like that, only 3 can divide the first two partitioned terms 3 and 6, since 2 is not involved. But here 3 should not simultaneously divide 3 and 6, or as a matter of fact, none of the primes should divide as such. So, when I split , then naturally we should be taking and , and hence the summand would be . I hope this is understandable.

[Sehra Sahu]

In this case m is fixed.

[Bruno Luong]

One more question, does the order of (l1, l2, ... lr) matter? Meaning do you take all the permutations of {li} in the sum or they are selected said once, in the increasing order?

[Sehra Sahu]

Order is not the matter. Just to remove the possible repetitions, we can consider them in increasing or decreasing order.

Answer 1

[Bruno Luong]

Some detail of the formula is not totally clear to me (do we consider li>=0 or li>=1), but it goes like this

(EDIT, assume now the constraints li>=1)

m = 20;
r = 2;
verbose = true;

l = Partition(m, r);
% list all the odd primes
ptab = arrayfun(@(x) oddprime(x), 1:m, 'unif', 0);

% do the sum
s = 0;
for i=1:size(l,1)
    [p, li] = alldistinctp(l(i,:), ptab);
    if ~isempty(p) && size(p,2)==r
         q = li./p;
         if verbose
            T = array2table([li; p; q]);
            T.Properties.RowNames = {'l' 'p' 'l/p'};
            disp(T)
         end
         s = s + 1/prod(factorial(q));
    end
end
fprintf('s(m=%d,r=%d) = %f\n', m, r, s)

function [p, l] = alldistinctp(l, ptab)
l(l==0) = [];
n = size(l,2);
p = cell(1,n);
[p{:}] = ndgrid(ptab{l});
p = cat(n+1,p{:});
p = reshape(p, [], n);
% discard n-tuplets with repeated p
b = any(diff(sort(p,2),1,2)==0,2);
p(b,:) = [];
end

function p = oddprime(x)
p = unique(factor(x));
p(p==2) = [];
end

function v = Partition(n, lgt)
% v = Partition(n)
% INPUT
%   n: non negative integer
%   lgt: optional non negative integer
% OUTPUT:
%   v: (m x lgt) non-negative integer array such as sum(v,2)==n
%       each row of v is descending sorted
%       v contains all possible combinations
%       m = p(n) in case lgt == n, where p is the partition function
%       v is (dictionnary) sorted
% Algorithm:
%    Recursive
% Example:
% >> Partition(5)
% 
% ans =
% 
%      5     0     0     0     0
%      4     1     0     0     0
%      3     2     0     0     0
%      3     1     1     0     0
%      2     2     1     0     0
%      2     1     1     1     0
%      1     1     1     1     1

if nargin < 2
    lgt = n;
end
v = PartR(lgt+1, n, Inf);

end % Partition

%% Recursive engine of integer partition
function v = PartR(n, L1, head)
rcall = isfinite(head);
if rcall
    L1 = L1-head;
end
if n <= 2
    if ~rcall
        v = L1;
    elseif head >= L1
        v = [head, L1];
    else
        v = zeros(0, n, class(L1));
    end
else % recursive call
    j = min(L1,head):-1:0;
    v = arrayfun(@(j) PartR(n-1, L1, j), j(:), 'UniformOutput', false);
    v = cat(1,v{:});
    if rcall
        v = [head+zeros([size(v,1),1], class(head)), v];
    end
end
end % PartR

[Sehra Sahu]

You may consider for time being.

Apparently the code shows some error in line 19, says function definition not permitted.

[Bruno Luong]

Easy to fix: If you run on older version MATLAB, save the three functions in separate mfiles from the script.

It works for me it give out this for m=20, r=2

           Var1    Var2
           ____    ____

    l       19      1  
    p       19      1  
    l/p      1      1  

           Var1    Var2
           ____    ____

    l       17      3  
    p       17      3  
    l/p      1      1  

           Var1    Var2
           ____    ____

    l       15      5  
    p        3      5  
    l/p      5      1  

           Var1    Var2
           ____    ____

    l       14      6  
    p        7      3  
    l/p      2      2  

           Var1    Var2
           ____    ____

    l       13      7  
    p       13      7  
    l/p      1      1  

           Var1    Var2
           ____    ____

    l       11      9  
    p       11      3  
    l/p      1      3  

s(m=20,r=2) = 3.425000