MATLAB code for finding certain coefficients
4 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Sehra Sahu
am 6 Okt. 2020
Bearbeitet: Bruno Luong
am 11 Okt. 2020
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 Kommentare
Bruno Luong
am 10 Okt. 2020
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?
Akzeptierte Antwort
Bruno Luong
am 10 Okt. 2020
Bearbeitet: Bruno Luong
am 11 Okt. 2020
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 Kommentare
Bruno Luong
am 11 Okt. 2020
Bearbeitet: Bruno Luong
am 11 Okt. 2020
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
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Logical finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!