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.