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

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.

[Bruno Luong]

p = zeros(1:100)

That is a 100-dimension array. I don't believe your PC can handle that statement.

[Stephen23]

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.

[David Goodmanson]

It's fairly clear that zeros(1,100) was intended.

[Thomas Pfeifer]

As David said, I intended to mean zeros(1,100).

I did a correction. Thx for the comments :)

Answer 1

[David Goodmanson]

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

[Thomas Pfeifer]

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]

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.

Answer 2

[Sai Veeramachaneni]

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.

Answer 3

[Bruno Luong]

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

[Bruno Luong]

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