MATLAB code for finding certain coefficients
2 ビュー (過去 30 日間)
古いコメントを表示
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.
9 件のコメント
Bruno Luong
2020 年 10 月 10 日
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?
採用された回答
Bruno Luong
2020 年 10 月 10 日
編集済み: Bruno Luong
2020 年 10 月 11 日
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
2 件のコメント
Bruno Luong
2020 年 10 月 11 日
編集済み: Bruno Luong
2020 年 10 月 11 日
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
その他の回答 (0 件)
参考
カテゴリ
Help Center および File Exchange で Logical についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!