MATLAB Answers

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

2 views (last 30 days)
Thomas Pfeifer
Thomas Pfeifer on 17 Nov 2020
Commented: Bruno Luong on 20 Nov 2020
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;
i = i+1;
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;
i = i+1;
isn't optimally.


Show 1 older comment
Stephen Cobeldick
Stephen Cobeldick on 20 Nov 2020
p = zeros(1:100);
will produce an array of size 1x2x3x4x5x ... x99xx100, containing 100! elements (in full: 93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000 elements). If every atom of the known universe could store one bit of information you would require 64 * 100! / 1e80 universes to hold all of that information (approximately 6e79 universes). You could reduce that requirement by using single or integer class, but otherwise if you are not in posession of 6e79 universes to store your data, I suggest that you create a smaller array.

Sign in to comment.

Answers (3)

David Goodmanson
David Goodmanson on 20 Nov 2020
Edited: David Goodmanson on 20 Nov 2020
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'];
nd = nd(:,1:25);
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;
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;
nd = [n d];


Thomas Pfeifer
Thomas Pfeifer on 20 Nov 2020
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 on 20 Nov 2020
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.

Sign in to comment.

Sai Veeramachaneni
Sai Veeramachaneni on 20 Nov 2020
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.


Sign in to comment.

Bruno Luong
Bruno Luong on 20 Nov 2020
Edited: Bruno Luong on 20 Nov 2020
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);
p = p(1:k);
fprintf([repmat('%1.0f, ',1,6) '\n'], p(1:k))
>> 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
>> 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,
>> 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

  1 Comment

Bruno Luong
Bruno Luong on 20 Nov 2020
We note that numericaly
p(19) = 6167950454 == 1963319607*pi
>> sin(p(19))
ans =
This is about 1.2e6 larger than
>> sin(pi)
ans =

Sign in to comment.


Community Treasure Hunt

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

Start Hunting!