How can I create this function on MATLAB?

How can I write this function for the variable "n"? Assuming that I already specified the value of the constant theta.

[Star Strider]

The LHS seems to be missing ...

[Pablo Arboleya]

I'm new to MATLAB, I just want to know how to type this... Imagine LHS to be f(n).

[Star Strider]

What is θ?

[Pablo Arboleya]

θ = 1.306377883863080690468614492602605712916784585156713644368053759966434053766826598821501403701197395707296960938103086882238861447816353486887133922146194353457871100331881405093575355831932648017213832361522359062218601610856679057215197976095161992952797079925631721527841237130765849112456317518426331056521535131866841550790793723859233522084218420405320517689026025793443008695290636205698968726212274997876664385157661914387728449820775905648255609150041237885247936260880466881540643744253401310736114409413765036437930126767211713103026522838661546668804874760951441079075406984172603473107746775740640078109350834214374426542040853111654904209930908557470583487937577695233363648583054929273872814934167412502732669268404681540626763113223748823800118041206286013841914438857151609189388944789912125543384749359092744422082802260203323027106375022288131064778444817003723336406042118742608383328221769687812353049623008802672211104016065088809718347778314022490821844106377494000232824192701

[Guillaume]

And what precision do you want on your infinite sum? Also a 1000 digits of precision? Can you measure and angle with anywhere near that precision?

[Pablo Arboleya]

Sorry, what do you mean by a 1000 digits of precision?

[Guillaume]

The theta you wrote has 1000 digits. Can you measure an angle with that precison?

[Pablo Arboleya]

this theta is not an angle, is the mills constant... its used in number theory, the symbol has nothing to do with angles. It's used in mills theorem for prime numbers.

[Walter Roberson]

However, you are taking sin() of an expression involving theta, so it is acting as a rotation measurement.

In any case, you have theta to 1000 digits. How many digits do you need your formula to be accurate to?

[John D'Errico]

With only 1/k in the denominator, that sum will take forever to give you any degree of precisionin the sum. So that degree of precision theta is sort of wasted, because you would need literally trillions of terms to get even double precision accuracy.

[Pablo Arboleya]

I'm pretending to do a test to compare Wolframalpha's answer, and the website uses 1000 digits... so I want to use this 1000.

[Pablo Arboleya]

https://www.wolframalpha.com/input/?i=-1%2F2+%2B+(mills+constant%5E3%5E5)+%2B+(+sum_(k%3D1)%5E%E2%88%9E+sin(2+k+(mills+constant%5E3%5E5)+pi)%2Fk)%2F%CF%80 The result here is 16022236204009818131831320182, and it should be a prime number and end with "3" instead of a "2". I guess wolfram alpha is wrong, so I want to test this equation on MATLAB, all I wanna do right now is basically solve this "function" for n = 5, and check the results to compare.

[Walter Roberson]

What is known about the range of n ? Is it integer? Is it positive? Is it real-valued ?

[Pablo Arboleya]

n is any natural number

[Walter Roberson]

Interesting, I make it that the result is

theta^(3^n) - arccot(cot(theta^(3^n)*pi))/pi

[Pablo Arboleya]

wow, amazing... how can I define this function on MATLAB? I'm not familiar with the language, could you help me?

[Walter Roberson]

Using the intepretation that theta^3^n is theta^(3^n) then

theta =  sym('1306377883863080690468614492602605712916784585156713644368053759966434053766826598821501403701197395707296960938103086882238861447816353486887133922146194353457871100331881405093575355831932648017213832361522359062218601610856679057215197976095161992952797079925631721527841237130765849112456317518426331056521535131866841550790793723859233522084218420405320517689026025793443008695290636205698968726212274997876664385157661914387728449820775905648255609150041237885247936260880466881540643744253401310736114409413765036437930126767211713103026522838661546668804874760951441079075406984172603473107746775740640078109350834214374426542040853111654904209930908557470583487937577695233363648583054929273872814934167412502732669268404681540626763113223748823800118041206286013841914438857151609189388944789912125543384749359092744422082802260203323027106375022288131064778444817003723336406042118742608383328221769687812353049623008802672211104016065088809718347778314022490821844106377494000232824192701')/sym(10)^1000;
syms n positive
syms pi
result = theta^(3^n) - acot(cot(theta^(3^n)*pi))/pi

[Pablo Arboleya]

Thank you!

Answer 1

[John D'Errico]

With only 1/k in the denominator, that sum will take forever to give you any degree of precision in the sum. So that degree of precision theta is sort of wasted, because you would need literally trillions of terms to get even double precision accuracy.

Think of it like this: how many terms to do need to get 1000 digit accuracy in a sum with terms that decrease with 1/k?

When k is on the order of 1e16, you are finally approaching double precision, thus eps. But your computer is not fast enough to compute 1e16 terms in any reasonable amount of time. And don't forget that it takes a significant amount of computer time to compute anything in a thousand digits of precision.

But, just imagine that you can compute one term of that sum in a millionth of asecond. (That is being very gracious on my part, too. I just wish I had a computer that was that powerful, but I don't work for the NSA.)

So every second, one million terms. 31 million seconds per yar or so. So 3.1e13 terms per year. Or, roughly 300 years before you compute 1e16 terms. And that only gets you to roughly double precision. But you want to compute all this in 1000 digits. So you will need to compute on the order of 1e1000 terms, before you start to get close to that.

You should be able to do it, just around the time the universe suffers heat death. Of course, by then you will have a better computer.

However, you did ask how you might do this. You can use my HPF toolbox. Or you can use the symbolic toolbox. Done in HPF (you need to download it from the File Exchange first, and install it on your search path.)

DefaultNumberOfDigits 1000
theta = hpf('1.306377883863080690468614492602605712916784585156713644368053759966434053766826598821501403701197395707296960938103086882238861447816353486887133922146194353457871100331881405093575355831932648017213832361522359062218601610856679057215197976095161992952797079925631721527841237130765849112456317518426331056521535131866841550790793723859233522084218420405320517689026025793443008695290636205698968726212274997876664385157661914387728449820775905648255609150041237885247936260880466881540643744253401310736114409413765036437930126767211713103026522838661546668804874760951441079075406984172603473107746775740640078109350834214374426542040853111654904209930908557470583487937577695233363648583054929273872814934167412502732669268404681540626763113223748823800118041206286013841914438857151609189388944789912125543384749359092744422082802260203323027106375022288131064778444817003723336406042118742608383328221769687812353049623008802672211104016065088809718347778314022490821844106377494000232824192701');

Now, write a function to compute the sum. First I need to make sure that I've interpreted that iterated exponent properly.

Do you intend it as theta^(3^n), or or as (theta^3)^n? The difference is of course significant. By default, MATLAB would treat the computation of theta^3^n as the latter form.

Before I go too far, let me see how fast HPF is, working in 1000 digits.

timeit(@() theta^3^n)
ans =
    0.0015156
    
1/ans
ans =
       659.82       

So HPF can do 660 such computations per second, at least on my computer.

How fast is the symbolic toolboxhere?

digits 1000
theta = sym('1.306377883863080690468614492602605712916784585156713644368053759966434053766826598821501403701197395707296960938103086882238861447816353486887133922146194353457871100331881405093575355831932648017213832361522359062218601610856679057215197976095161992952797079925631721527841237130765849112456317518426331056521535131866841550790793723859233522084218420405320517689026025793443008695290636205698968726212274997876664385157661914387728449820775905648255609150041237885247936260880466881540643744253401310736114409413765036437930126767211713103026522838661546668804874760951441079075406984172603473107746775740640078109350834214374426542040853111654904209930908557470583487937577695233363648583054929273872814934167412502732669268404681540626763113223748823800118041206286013841914438857151609189388944789912125543384749359092744422082802260203323027106375022288131064778444817003723336406042118742608383328221769687812353049623008802672211104016065088809718347778314022490821844106377494000232824192701');
timeit(@() theta^3^n)
ans =
      0.00138
      
1/ans
ans =
       724.62

So the symbolic TB does 725 such computations per second, thus about the same as HPF. And, yes, my computer is getting a little long in the tooth, but it was quite fast when I bought it some years ago, so it is still decent. Not even close to a mllion per second. I think you are in deep trouble. Really deep. Massively deep. Humongously deep.

Sigh. anyway. I'll use the symbolic toolbox here.

function S = bigsum(n,kterms,theta)
  % precompute this number
  theta3n = theta^3^n;
  
  % initialize the sum. Note that 1/2 is representable exactly by
  % either HPF or by the symbolic toolbox, so I can just subtract 1/2.
  S = theta3n - 1/2;
  
  pie = sym('pi');
  
  for k = 1:kterms
      S = S + sin(2*k*pie*theta3n)/k;
  end
  S = vpa(S/pie);
end

It is not too bad, as long as I restrict the number of terms.

tic,S = bigsum(1,100,theta);toc
Elapsed time is 0.325404 seconds.

So 100 terms in 1/3 of a second. But do we even know the very first digit of that number yet?

for kt = 100:105
vpa(bigsum(1,kt,theta),20)
end
ans =
0.81880906421480201873
ans =
0.82165197521120346388
ans =
0.82334942686820709411
ans =
0.82099366721044812494
ans =
0.8187293340428497687
ans =
0.82046391141319398969

It seems to be converging. SLOWLY. The first digit might be an 8, that is, if this is a convergent series. After a fair amount of time,10000 terms gives:

vpa(bigsum(1,10000,theta),20)
ans =
0.82099823595744050401

Personally, I recommend buying a VERY fast computer. Better yet, you need to reconsider if you can do this computation in a finite amount of time, unless you are willing to wait for a few gazillion years.

[Walter Roberson]

>> syms k n theta pi positive

>> -1/2 + theta^(3^n) + symsum(sin(2*k*theta^(3^n)*pi)/k,k,1,inf)/pi

ans =

piecewise(~in(exp(exp(n*log(3))*log(theta)), 'integer'), theta^(3^n) - ((log(1 - exp(-pi*exp(exp(n*log(3))*log(theta))*2i))*1i)/2 - (log(1 - exp(pi*exp(exp(n*log(3))*log(theta))*2i))*1i)/2)/pi - 1/2)

With positive integer n and the known theta value, the condition drops out, leaving a finite expression that with sufficient prodding can be rewritten to something fairly simple.