5 views (last 30 days)

Show older comments

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.

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'];

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

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.

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.

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);

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

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

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

Start Hunting!