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 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
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 个评论
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
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.
0 个评论
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
1 个评论
Bruno Luong
2020-11-20
We note that numericaly
p(19) = 6167950454 == 1963319607*pi
but
>> sin(p(19))
ans =
1.497342914414640e-10
This is about 1.2e6 larger than
>> sin(pi)
ans =
1.224646799147353e-16
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Matrix Indexing 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!