MATLAB Answers

## How can I create this function on MATLAB?

Pablo Arboleya

### Pablo Arboleya (view profile)

さんによって質問されました 2019 年 3 月 30 日

### Walter Roberson (view profile)

さんによって コメントされました 2019 年 3 月 31 日
John D'Errico

### John D'Errico (view profile)

さんの 回答が採用されました How can I write this function for the variable "n"? Assuming that I already specified the value of the constant theta.

Pablo Arboleya

### Pablo Arboleya (view profile)

2019 年 3 月 30 日
wow, amazing... how can I define this function on MATLAB? I'm not familiar with the language, could you help me?
Walter Roberson

### Walter Roberson (view profile)

2019 年 3 月 30 日
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

### Pablo Arboleya (view profile)

2019 年 3 月 30 日
Thank you!

サインイン to comment.

R2018a

## 1 件の回答 2019 年 3 月 30 日

### John D'Errico (view profile)

2019 年 3 月 30 日
採用された回答

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

### Walter Roberson (view profile)

2019 年 3 月 31 日
>> 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.

サインイン to comment.