Thread Rating:
  • 0 Vote(s) - 0 Average
  • 1
  • 2
  • 3
  • 4
  • 5
A first hires look at tetration \(\lambda = 1\) and \(b = e\)
#21
(11/22/2021, 07:42 AM)Ember Edison Wrote: Your pictures of \( \mu = \lambda = 1+i \), The four slashes inside the image look like the function has crashed.

I also recently tested that in base=\(  e^10^24 \) or init(1,1E24,1000), the beta function would crash.

init(1,1E-24,1000) will also crash.


I also created images of a circle scanned in the complex plane with a radius of 1/1E16/1E-16. The images of the real axis still need some time.

I do not believe they crashed for \(\lambda = 1+i, \mu = 1+i\), I believe that's just what the branch cuts are going to look like. It is slanting at the right angle; that's the direction the branch cuts should go in; and the functional equation is still satisfied. You'll see this a lot for \(\lambda\) complex; it'll begin to branch in a chaotic manner. I like to categorize tetration functions as functions taking \(\mathcal{P} \to \mathbb{C}\) where \(\mathcal{B} = \mathbb{C}/\mathcal{P}\) are the set of discontinuities; and so long as:

$$
\int_{\mathcal{B}}\,dA = 0\\
$$

For \(dA\) the standard Lebesgue area measure of \(\mathbb{R}^2\); it's good enough to be a tetration in my boat. So the beta tetrations account for branch cuts and singularities; but are still holomorphic almost everywhere.

I am not surprised base \(e^{10^{24}}\) crashes; sadly I can't think of anyway to fix that. I'll take a look at \(e^{10^{-24}}\), though, there might be a way to fix that, I'll look.
Reply
#22
(11/23/2021, 03:21 AM)JmsNxn Wrote: I am not surprised base \(e^{10^{24}}\) crashes; sadly I can't think of anyway to fix that. I'll take a look at \(e^{10^{-24}}\), though, there might be a way to fix that, I'll look.


Do you have any opinion about \(  b \to \pm \infty \)? Like the folklore about \( tet_0(z) / tet_{\infty}(z) \) inside the forum, in the images they clearly show some kind of symmetry
Reply
#23
(11/23/2021, 03:43 AM)Ember Edison Wrote:
(11/23/2021, 03:21 AM)JmsNxn Wrote: I am not surprised base \(e^{10^{24}}\) crashes; sadly I can't think of anyway to fix that. I'll take a look at \(e^{10^{-24}}\), though, there might be a way to fix that, I'll look.


Do you have any opinion about \(  b \to \pm \infty \)? Like the folklore about \( tet_0(z) / tet_{\infty}(z) \) inside the forum, in the images they clearly show some kind of symmetry

I apologize, but my knowledge of this forum doesn't go as far back as I wished it did. There are many posts; even if I was a member; I was 16 at the time and was barely grasping calculus, lol. If you could explain what you mean exactly, I'd be happy to comment.

And good news. I have figured out how to make init(1,1E24) to work; I'm just thinking of the best way to program it at the moment, such that it works conducively with the rest of the program. Give me 5 days or so, and I think I can get it to work. It's a small fix, but I want to make sure it doesn't screw up anything major. But just so you know, the value \(\beta(1)\) will already be an overflow to the best of my calculations. The proposed solution will successfully grab the taylor series; but the recursive protocol may fail because it requires searching in a radius of \(|z| \le 10^{-24}\) and applying \(\exp\) a total of \(10^{24}\) times to get the next value; and this typically results in a recursion error on pari's part. I'm seeing if I can bypass this some how. Worst comes to worst. I can definitely produce the series:

$$
\beta_{1,10^{24}}(z) = \sum_{j=1}^\infty a_j e^{zj}\\
$$

But then trying to push forward I imagine may cause difficulties.  Again, give me five days, I'll see if I can integrate the code. I imagine you'll have to take a larger sample than a 1000 too, because this series will converge much slower than when \(b\) is a reasonable value.

EDIT: This is the final edit of this post.  I've got it to work at a beta level with previous code. Just let me debug at this point. We have to add another function which behaves like beta, but works for absurdly large values. I need to install more ram for the most part. Luckily I recently bought another 32 gb stick so I can let pari-gp get ridiculous. I'll update in 5 days or so. Again, I see how to make it work; it's just a matter of making it compatible with the previous programming. Also, there's a tad bit of arbitrary; so it fails for crazy values; but you can tweak it to make it work. I use the sqrt function but it isn't necessary; I just pulled it from a hat.
Reply
#24
We can mathematically recreate the function \(\beta_{\lambda,\mu}(z)\). Where as we'd use:

$$
\begin{align}
\beta_{\lambda,\mu}(s) &= \Omega_{j=1}^\infty \dfrac{e^{\mu z}}{1+e^{\lambda(j-s)}}\,\bullet z\\
\beta(s+1) &= \dfrac{e^{\mu\beta(s)}}{1+e^{-\lambda s}}\\
\end{align}
$$

To avoid a substantial amount of overflows we can use:

$$
\begin{align}
\widetilde{\beta}_{\lambda,\mu}(s) &= \Omega_{j=1}^\infty \dfrac{e^{\mu z}}{1+e^{\lambda(j-s)+\sqrt{|\mu|}}}\,\bullet z\\
\widetilde{\beta}(s+1) &= \dfrac{e^{\mu\widetilde{\beta}(s)}}{1+e^{-\lambda  s + \sqrt{|\mu|}}}\\
\end{align}
$$

Any function of the form \(g(|\mu|) = |\mu|^\delta\) for \(0 < \delta < 1\), should work fairly well. I tend to just choose \(\sqrt{\cdot}\). I've managed to get my code for the taylor series of beta to work for the theoretical memory limit of pari-gp. So if you've got the ram, I think I can make 1E24 work. Still scratching my head about how to code this perfectly though.

Regards, James.



Alright, so here's the code:

Code:
\p 1000
\ps 300

/*
Set digit precision to 1000, and series precision to 300
This is done a tad gratuitously but it's helpful for precision.
You can modify these constants and initialization will respect that.
*/



/*
This first function is an initialization function.

The value b = base, is the base of the superexponential;
with the convention that exp(b*z) is the exponential. Therein, our only bad value is b=0.

The value l = mult, is the multiplier of the beta function we choose.
Recall, that 2*Pi*I/mult is the period of the beta function.

count is telling us how deep we want to do the recursion.
By default I have set it to the seriesprecision you choose with \ps
And by default I've set \ps to 300, as this is already ridiculously accurate.
Also, by default it does at least 240 iterations.


Before any of the other functions can be run, we must run init.
Which requires a choice of mult and base given by the first and second argument, respectively.

This code is a hybridization between my code and Sheldon's code.
I had originally initialized with the multiplier as a free variable which had its benefits (graphing the multiplier, particularly).
Sheldon's method is significantly faster by requiring us to declare the multiplier; and a quick speed up with subst.
*/


offset = 0; /*this is a global variable that is unused in this protocol unless we call call large values of b (approx abs(b) = 1E6);
                at which; it prevents a good amount of overflows, and mostly just allows us to still grab a taylor series in the intialization process*/


init(l,b,{count=default(seriesprecision)}) = {
    base = b;
    mult = l;  /*declare these variables to work for all other functions; make them global*/
    
    if(count < 240  && default(seriesprecision) < 240,count=240);  /*This just makes sure we do at least 240 iterations regardless of series depth*/
    
    Cut = 0;  /*this is just a quick cut value, to limit how deep the recursion goes, otherwise we get deep recursion*/

    if(abs(real(mult)) < 1E-3, Cut = offset-1E3, Cut = offset-50 - 1/real(mult));

    beta_taylor = x/exp(mult) + O(x^2);  /*Sheldon is to credit for this polynomial translation*/
    for (n=1,count,
        beta_taylor=exp(base*subst(beta_taylor,x,x/exp(mult))) / (1 + exp(mult)/x) /*this is a degree faster than my old code; absolutely genius work by Sheldon*/
    );
    beta_taylor=Pol(beta_taylor); /*we treat these objects as polynomials*/
}

init_OFF(l,b,{count=default(seriesprecision)}) = {
    base = b;
    mult = l;  /*declare these variables to work for all other functions; make them global*/
    
    if(count < 240  && default(seriesprecision) < 240,count=240);  /*This just makes sure we do at least 240 iterations regardless of series depth*/
    
        
    Cut = 0;  /*this is just a quick cut value, to limit how deep the recursion goes, otherwise we get deep recursion*/

    if(abs(real(mult)) < 1E-3, Cut = offset-1E3, Cut = offset-50 - 1/real(mult));
    
    /*This variable offsets the nested exponentials so that the iterate will converge regardless how large b is; or how close it is to zero*/
    if(abs(b) > 1E6, offset = sqrt(abs(b)));
    if(abs(b) <1E-6, offset = sqrt(1/abs(b)));

    beta_taylor = x/exp(mult+offset) + O(x^2)/exp(offset);  /*Sheldon is to credit for this polynomial translation*/
    for (n=1,count,
        beta_taylor=exp(base*subst(beta_taylor,x,x/exp(mult+offset))) / (1 + exp(mult+offset)/x) /*this is a degree faster than my old code; absolutely genius work by Sheldon*/
    );
    beta_taylor=Pol(beta_taylor); /*we treat these objects as polynomials*/
}

/*
This truncates a series into the constant coefficient--needed for the recursion
We will use this to generate if conditions.
*/

Const(poly) = {
    my(vars = variables(poly));
    substvec(poly, vars, vector(#vars))
};

/*
A hybridization between mine and Sheldon's code. Sheldon optimized the exponential series defining beta.
But he didn't account for loss of accuracy in the functional equation.
But otherwise, mine and Sheldon's code are mathematically identical.

This function will satisfy the functional equation:

beta(z+1) = exp(base*beta(z))/(exp(-mult*z) + 1)

It will have singularities at the points mult*(z-j) = (2k+1)*Pi*I for j >= 1 and j,k integers.

It also has a period of 2*Pi*I/mult
*/



beta(z) = {
    
    /*
    The goal of this function: if we're in the radius of convergence of the exponential series (at an accurate level), just substitute the value.
    Otherwise, we pull backwards until we are in the radius of convergence,
    and iterate forwards with the functional equation.
    */
    
    /*
    The accurate radius of convergence (the radius where we are still sufficiently accurate) is about exp(-50 - 1/real(mult)) for most reasonable values.
    If mult is too small though, this can create recursion errors (Pari doesn't like 1E4 recursive calls)
    You may need to increase depth of iteration/series precision/digit precision for more anomalous values.
    */
    
    if(real(Const(z)) <= Cut,
        subst(beta_taylor,x,exp(mult*z)),
        exp(base*beta(z-1))/(1+exp(-mult*(z-1)))
    );
}

/*
This function is specifically for values of base greater than 1E6 or less than 1E-6; and satisfies a different functional equation than beta:
    
exp(base*beta_off(z))/(1+exp(mult*(offset-z))) = beta_off(z+1)

This effectively makes a shift to the argument so that we can still grab the Taylor series. It is related to beta by the relation:

beta_off(z+offset) = beta(z)
*/


beta_off(z) = {
    /*
    The goal of this function: if we're in the radius of convergence of the exponential series (at an accurate level), just substitute the value.
    Otherwise, we pull backwards until we are in the radius of convergence,
    and iterate forwards with the functional equation.
    */
    
    /*
    The accurate radius of convergence (the radius where we are still sufficiently accurate) is about exp(-50 - 1/real(mult)) for most reasonable values.
    If mult is too small though, this can create recursion errors (Pari doesn't like 1E4 recursive calls)
    You may need to increase depth of iteration/series precision/digit precision for more anomalous values.
    */

        if(real(Const(z)) <= Cut,
        subst(beta_taylor,x,exp(mult*(z))),
        exp(base*beta_off(z-1))/(1+exp(-mult*(z-1) + offset))
    );
}

if(offset != 0, beta(z) = beta_off(z+offset/mult));

****** Edited this code 27/11/2021 fixed a dumb arithmetic error*************



This can be added before we get to any of the super exponential stuff. I made a few cosmetic alterations. I didn't include any of the super exponential stuff; because I haven't found a good way to make the new functions into a super exponential. As you seem to be mostly dealing with the beta function, it should be fine for you, Ember.

To give a quick run through.

Now we have two different initialization methods: init, which is the same as before; and init_OFF.

If you use init_OFF(mult,base,N) if abs(base) > 1E6 or abs(base)<1E-6 it will initialize a slightly different polynomial than init(mult,base,N). It will create a new value offset = sqrt(abs(b)) or offset = sqrt(1/abs(b)) and effectively it will find the polynomial multiplied by a factor of exp(offset).

This means, instead of finding the exponential series:

$$
\beta(s) = \sum_{j=1}^N a_j e^{\lambda j s}\\
$$

It will now find

$$
\widetilde{\beta}(s) = \sum_{j=1}^N a_j e^{- j*\text{offset}}e^{\lambda j s} \\
$$

Where:

$$
\widetilde{\beta}(s+\text{offset}/\lambda) = \beta(s)\\
$$


From here, we have two different beta functions. We have beta, which runs precisely as before but will fail for the large values of base; and beta_off, which will work for large values of base; and basically just becomes \(\text{beta_off}(s) = \widetilde{\beta}(s)\). The trouble is if you plug in \(\text{beta_off}(\text{offset}/\text{mult}) = \beta(0)\); expect an overflow or a deep recursion error.

I have added in the last last line that if offset is non zero; make beta(z) = beta_off(z+offset/mult); and this seems to be working; but it could be glitchy. beta_off is much better; because again, beta(0) could be a huge overflow.




In such a sense, this is moreso intended to grab the taylor series near zero; if you add enough terms to the series (increase N); it is theoretically possible to avoid the deep recursion error; but it would require modifying how I defined the Cut variable; simply make this offset-100 instead of offset-1E3 for instance; but then you'll have to add plenty more terms to the series. And even then, if you avoid the deep recursion, for base = 1E24, expect beta(0) = beta_off(offset/mult) to overflow (no way around it, it is the pari-gp cap we're hitting here).

Also, as a disclaimer. I've only managed to make this work for 1E16, 1E-16, -1E16. On my computer the paristack overflows for 1E24--which is because I don't have enough ram... So beware, to get 1E24 and 1E-24 to work, you should be prepared to use a lot of ram. This may be fixable if we alter how offset is defined. Currently I just use the sqrt function; there may be a better way.

Notice, that there is a tad amount more of volatility. For example, for the 1E16 case, beta(-34) is about as large as you can make the argument. it starts to explode once you hit this point. But it will still satisfy the functional equation to the desired precision. It can be a little glitchy though. The quick plots I made are as expected; so everything seems to still work.

For example, here is:


Code:
\r beta.gp
init_OFF(1,1E16)
 *** _/_: Warning: increasing stack size to 16000000.
 *** _/_: Warning: increasing stack size to 32000000.
 *** _/_: Warning: increasing stack size to 64000000.
 *** _/_: Warning: increasing stack size to 128000000.
%50 = 0.E-16372977168*x^301 - 4.112449703997607243944723183055879593537065702677561890516738051002774944146442005932034696254832027178014659147185510860912437649631253473180366017745131974496148556063496698945081499280793553317117286374200962646778276149307902860132175848058942741198546835884827687707445019161970353324270976894561536876435325504346520062195662721201263069006487311834469140631510278354911651830730235400491998478473013106756680621618213362709708666862392995074448924602516592094860140297433878072233059290568135468112570516438784472729245674723826080752975799795024360547637834762835511377652173486940076203679070539549365621277780466162626283999483725232846435110256173135288788784475684127306956773019728068821482377119803137066761829627455956012377590538087545091292103536499004388482681377452123234857649463722133439645988430959120983464199326976403362598086958120718735149821693392785489608787216335580646463101544715527690645632902825400720869355169735258731902587401633871988526799444785050439347790602621 E-13028834572*x^299 - 7.300286911332839510693978307171106024501688280990622418844338820461971269384910811660295526946588535726541294209274370275989148544402753794643437568702842834806669355466390228312629263767867818322977780710587792169974643945450691176802318170955217579667716178753169703771213160892036599267409979999660242682178814557036094097403159657408687747771412156928163127766652586211829581727438889903975586385599555957384437558137107522755154650337781501696926373624421060219585134686438020115722441092052165343930231327713686477735412456671551916219670835831821257108000413748659748207918449760178808504542053345732883856079912256291412549255082986908751186773239784010466016989120768571232882685427464448416012570165346584330585507640125003807441481665893518137413656255711366961441062909214917483335987060773204271481607794704428076416142435230207117185995385561475708701185204468402019860187442147669684440399340419690890126612092981682395790855[+++]

\p9

ploth(X=-34,-33.5,real(beta(X)))
 *** exp: Warning: increasing stack size to 16000000.
 *** exp: Warning: increasing stack size to 32000000.
 *** exp: Warning: increasing stack size to 64000000.
%65 = [-34.0000000, -33.5000000, 1.57726759 E-12, 4.36942151 E8]

which spits out:

   

And after that it just blows up.

All in all, that's about it!

Regards.



Please note I've updated the code 27/11/2021; I made a really dumb arithmetic error.

I had coded with \(\exp(\text{offset})\) but I had written \(\exp(\lambda\text{offset})\) a couple of times (I mixed up a change of variables a couple of times). It had zero effect when \(\lambda = 1\); so I didn't notice it, but the old code caused errors for any \(\lambda \neq 1\). I apologize for the confusion. I was too focused on the \(2\pi i\) periodic case.
Reply
#25
Thanks, James! Your solution is valid in broad principle. 

But picking g(|u|)=sqrt(abs(u)) is too much of a burden on memory and cpu time: init_OFF(1,1E20) already requires >25GB RAM and several hours of cpu time.
The memory limit is less important for me, I have 128GB RAM and about 8TB HDD and 1TB SSD available for Virtual Memory, but it's just too slow. Not only is init_OFF slow, but beta_off is slow too. This means that for myself g(|u|) better not exceed 1E10.

sqrt(log(3E694127911065419641)) = 1.26E9, so I changed to g(|u|)=sqrt(log(abs(u))), which already solves init_OFF(1,1E25), but fail in init_OFF(1,1E50). A better form of the g function is still about a few weeks away from being determined.

setup g(|u|)=log(abs(u))^6, init_OFF(1E50) need >50GB, I don't think there is much time left before we get the final industrial upper limit of the beta function.

   
Reply
#26
Very cool to see!

I just chose sqrt because I figured it will definitely be large enough; but it definitely is nuking a mosquito.

I'm glad to see it's working!
Reply




Users browsing this thread: 1 Guest(s)