Find first item with a certain condition and use it, instead of finding all items with that condition

1 次查看(过去 30 天)
First of all, this is my first time on this site. So excuse me if i did something wrong.
My following code does what it should do. But for larger n, like n=10^10 it takes a very long time. If you have some advice to my code so feel free to tell me.
i = 1;
k = 1;
p = zeros(1,100);
p(1) = 3;
n = 10^5;
while i <= n
if abs(sin(i))<abs(sin(p(k)))
p(k+1) = i;
k = k+1;
end
i = i+1;
end
p = p(1:k);
For that reason i tried to use find.
For example for p(2):
p(2) = find(abs(sin(1:n)) < abs(sin(3)) , 1 , 'first');
An error is shown for n^10:
Requested 10000000000x1 (74.5GB) array exceeds maximum array size preference. Creation of arrays greater than this limit may
take a long time and cause MATLAB to become unresponsive.
My Question:
Is there an option to find the first number and use it, instead of creating a very long array, like 'stop after you found first item' .
I guess
while i <= n
if abs(sin(i))<abs(sin(p(k)))
p(k+1) = i;
k = k+1;
end
i = i+1;
end
isn't optimally.
  4 个评论

请先登录,再进行评论。

回答(3 个)

David Goodmanson
David Goodmanson 2020-11-20
编辑:David Goodmanson 2020-11-20
Hi Thomas,
[After reading comments from Bruno and Rik I realized I forgot to mention that in order to get your code to run I changed zeros(1:100) to the intended zeros(1,100)].
As you can see, when the number of integers goes up, checking them all will become untenable. I don’t know if your goal is or has to be extending that same basic approach, but other methods are available.
For large i, you are looking for sin(i) to be as near to 0 as possible. That means that i has to be as near as possible to some multiple of pi. So
i = m*pi +r
where m is an integer and r should be small. This means you are looking for better and better rational approximations to pi,
i/m ~~ pi
for large values of i. This suggests the method of continued fractions which can give a series of increasingly accurate rational approximations i/m for pi. The numerators of those fractions are the values p that you have calculated.
Running your code for N = 1e9 gave [saving small n for use later]
p = [3 22 333 355 103993 104348 208341 312689 833719 1146408 ...
4272943 5419351 80143857 165707065 245850922 411557987]
and took 5 minutes, which seemed long enough so that I was not keen on N = 1e10.
The code below does require the symbolic toolbox and variable precision arithmetic. It uses vpa with default precision (32 digits), and first calculates the coefficients in the continued fraction representation of pi, of which the first 33 (according to OEIS A001203) are accurate. That’s about what you would expect. Using those 33 then gives a series of rational n/d approximations to pi, which are contained in the 2xn matrix nd, where n is the top row and d is the bottom row. The early values of n agree with your calculations for p. The fourth column gives the famous approximation to pi, 355/113.
I could only easily check the first 25 instances of n/d [OEIS A002485, 2486] so I stopped there. In that case the largest value for n (or p) is approximately 8.9e12 which is going to take a bit of time by the check-them-all method.
pi32 = vpa(pi)
a = nd2cfrac([pi32, 1])
nd = [];
for k = 1:33;
ndk = cfrac2nd(a(1:k));
nd = [nd ndk'];
end
nd = nd(:,1:25);
disp(nd')
function a = nd2cfrac(nd)
% continued fraction representation of rational number n/d in lowest terms
% (n/d is automatically reduced to lowest terms on input).
% where a = [a0 a1 ... an] and nd is the 1x2 vector [n d]
%
% a = nd2cfrac(nd)
n = nd(1); d = nd(2);
g = gcd(n,d);
n = n/g; d = d/g;
a = [];
while d~=0
ndi = floor(n/d);
a = [a ndi];
nrem = n - ndi*d;
n = d;
d = nrem;
end
end
function nd = cfrac2nd(a)
% continued fraction evauation to find the rational fraction n/d
% in lowest terms
% where a = [a0 a1 ... an] and nd is the 1x2 vector [n d]
%
% nd = cfrac2nd(a)
n = 1;
d = 0;
for k = length(a):-1:1
napd = n*a(k) + d;
g = gcd(napd,n);
d = n/g;
n = napd/g;
end
nd = [n d];
end
  2 个评论
Thomas Pfeifer
Thomas Pfeifer 2020-11-20
编辑:Thomas Pfeifer 2020-11-20
Hi David,
I appreciate your answer. I like your Code.
But isn't there a possibility to use the sinus-method to get the numerators p in a faster they i did?
I couldn't find a possibility to break the search after i found the first item.
Like 22 is the first number >3 that fulfills
abs(sin(22)) < abs(sin(3)) (1)
--> don't look for more items, just take the 22 and break the search at that point
I guess the find-method searches for all numbers that fulfills (1) or?
Bruno Luong
Bruno Luong 2020-11-20
编辑:Bruno Luong 2020-11-20
Then the result outcome depends entirely on the implementation of SIN and storage of PI (the one used internally by SIN, that might or might not be the same as MATLAB PI) as double.
Just David rightly points out, using the rational fraction of pi is the correct way to go, unless if you want to know how good/bad sin(x) is implemented for large n.

请先登录,再进行评论。


Sai Veeramachaneni
Sai Veeramachaneni 2020-11-20
Your second example using find is not feasible due to memory limitations.
To my understanding you are doing linear search and for that your first example is optimal one.

Bruno Luong
Bruno Luong 2020-11-20
编辑:Bruno Luong 2020-11-20
Process by block, eventually you get the result after a couple of minutes
p = zeros(1,100);
p(1) = 3;
n = 10^10;
m = 1e7;
sk = abs(sin(p(1)));
k = 1;
for j=1:ceil(n/m)
fprintf('j=%d/%d\n', j, ceil(n/m));
i0 = (j-1)*m;
s = abs(sin(i0+(1:m)));
for i=1:m
if s(i)<sk
k = k+1;
p(k) = i0+i;
sk = s(i);
end
end
end
p = p(1:k);
abs(sin(p))
fprintf([repmat('%1.0f, ',1,6) '\n'], p(1:k))
mod(p/pi+0.5,1)-0.5
Result
>> abs(sin(p))
ans =
Columns 1 through 6
0.141120008059867 0.008851309290404 0.008821166113886 0.000030144353359 0.000019129335778 0.000011015017584
Columns 7 through 12
0.000008114318195 0.000002900699389 0.000002312919416 0.000000587779973 0.000000549579498 0.000000038200475
Columns 13 through 18
0.000000014772847 0.000000008654781 0.000000006118065 0.000000002536716 0.000000001044633 0.000000000447449
Column 19
0.000000000149734
>> fprintf([repmat('%1.0f, ',1,6) '\n'], p(1:k))
3, 22, 333, 355, 103993, 104348,
208341, 312689, 833719, 1146408, 4272943, 5419351,
80143857, 165707065, 245850922, 411557987, 1068966896, 2549491779,
6167950454,
>> mod(p/pi+0.5,1)-0.5
ans =
Columns 1 through 6
-0.045070341448628 0.002817496043395 -0.002807900797706 0.000009595245686 -0.000006089052476 0.000003506189387
Columns 7 through 12
-0.000002582863090 0.000000923319021 -0.000000736210495 0.000000187137630 -0.000000174855813 0.000000012340024
Columns 13 through 18
-0.000000003725290 0.000000007450581 0 0 0 0
Column 19
0

类别

Help CenterFile Exchange 中查找有关 Matrix Indexing 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by