/* This program is a test of the infinite compositions base b = 1/2 This program includes the file INIT_HALF.dat which stores a 100 x 100 matrix of Taylor coefficients which are enough for 100 digit precision. If the INIT.DAT file is corrupted or lost; run the command, beta_init(100); This will take approximately an hour or two to compile. If it takes too long; you can redownload the .zip file this code was packaged in. This code is almost entirely written by James David Nixon (excluding the graphing program); and is an attachment to his construction of the beta tetration. This code is coupled with the paper The Limits Of A Family; Of Asymptotic Solutions to the Tetration Equation. It can be found on the author's arxiv page. */ /* set series precision/numerical precision to 100/100 */ \p 100 \ps 100 /* A brief note on the function protocol for grabbing taylor series or not. For example, if you were to write, Abel_N(1+z,log(2)), this will print out the Taylor series in z about the point 1 with multiplier log(2). Or beta(3+z,log(2)) will print out the Taylor series in z about the point 3. It's fairly self explanatory. Or mu(1,1+y) will print out the Taylor Series in the multiplier y about the point 1. So we can do this for all function arguments. This code is not really based off of Taylor series; so it runs much slower than just numerical evaluation. But it's surprisingly effective. */ /*******************************************************************************************/ /* This is a construction of a faster running beta function, but it requires more data; and avoids much of the errors. This program should have come with the file INIT_HALF.dat which includes all the Taylor series Coefficients needed for 100 digit accuracy. To recreate this file, simply run beta_init(100); or beta_init(x) for any larger x; but 100 is enough for 100 digits of accuracy. I've coded everything for 100 digit accuracy or less; and if you want to meddle around there could be errors We first use the function Phi_Inv to discover a taylorseries for beta(-log(w)/log(2)) about w=0, and then use this Taylor series and the functional equation to extend beta. */ base = log(0.00001); Phi_Inv(w,l,{n=100}) = { my(out = 0); for(i=0,n, out = w*exp(base*out)/(exp(l*(n+1-i))+w) ); out; } /* this is the initialization protocol; you do not need to run it if you already have the INIT_HALF.dat file. */ beta_init(n) = { beta_taylor = Phi_Inv(w,l,n); writebin("INIT_TINY.DAT",beta_taylor); print("Done initializing beta."); } /*This just reads initialization data; just make sure it's in the same folder as the source. And pari-gp is reading from the source's folder.*/ beta_taylor = read("INIT_TINY.DAT"); /* This truncates a series into the constant coefficient--needed for the recursion We will use this to generate our if conditions */ Const(poly) = { my(vars = variables(poly)); substvec(poly, vars, vector(#vars)) }; /* This is the beta function we need. It loads the Taylor values for small values, and then we push forward using the functional equation. It satisfies the equation beta(z+1,y) = exp(base*beta(z,y))/(exp(-y*z)+1) to at least 100 digits */ beta(z,y) = { if(real(Const(z)) < -50, sum(j=1,99,subst(Pol(polcoef(beta_taylor,j,w),l),l,y)*exp(y*z*j)), exp(base*beta(z-1,y))/(exp(-y*(z-1)) + 1) ); } /* This is the error function between beta and the tetration function base associated with the multiplier y. 100 iterations managed to be enough to get 20 digit precision. We don't overflow like we usually do with b = e, so 100 iterations is possible. */ mu(z,y,{count=100})={ if(count>0, count--; log(1 + mu(z+1,y,count)/beta(z+1,y))/base - log(1+exp(-y*z))/base, 0 ); } /* This is the first Abel function we construct. It is associated to the multiplier y; which we will vary in the final result. This Abel function will not be normalized as Abel_N(0,y) = 1; */ Abel_N(z,y,{count=100}) = { beta(z,y) + mu(z,y,count); } /* HERE ENDS ALL THE CODE I'VE WRITTEN--the rest was procured from mike3 of the tetration forum*/ /* ============================================================================= ** Color graphing system - mike3 - 20.11.10 04:07 ** Hi. ** I thought I'd post the code I use to generate the color graphs from Pari/GP. ** Here it is. ** ** Note: the output is in .PPM format. ** You'll need something else to convert that to .PNG. (I use GIMP.) ** ** Also, I might warn you: it takes a LONG time to graph a complex function ** with a significantly complicated calculation procedure, as it must be ** evaluated at every pixel of the graph. ** ** (updated 12/16/2010 -- program was not all good: ** * spurious "func" parameter in MakeGraph and "safetyarg" was missing.) ** ------------------------------------------------------------------------------------------ ** ** ============================================================================= */ /* =============================== Code: ==================================================== */ /* Complex function magnitude/phase plotter. */ /* To use: * 1. Define function to graph as func(z). * 2. Load this program. * 3. Execute MakeGraph(width, height, x0, y0, x1, y1, filename) with the parameters given as follows: * width, height = width/height of image in pixels * x0, y0, x1, y1 = rectangle of complex plane to graph: x0 + y0i in upper-left corner to x1 + y1i in lower-right corner * filename = name of file to save as. * Output is in .PPM format. */ /* Color conversion (HSB to RGB). */ HSB2RGB(hsb) = { local(H=hsb[1]); local(S=hsb[2]); local(B=hsb[3]); local(HH); local(F); local(P); local(Q); local(T); HH = floor(6*H)%6; F = (6*H) - floor(6*H); P = B*(1 - S); Q = B*(1 - (F*S)); T = B*(1 - (1-F)*S); if(B > 1.0, B = 1.0); if(B < 0.0, B = 0.0); if(P > 1.0, P = 1.0); if(P < 0.0, P = 0.0); if(Q > 1.0, Q = 1.0); if(Q < 0.0, Q = 0.0); if(T > 1.0, T = 1.0); if(T < 0.0, T = 0.0); if(HH == 0, return([B, T, P])); if(HH == 1, return([Q, B, P])); if(HH == 2, return([P, B, T])); if(HH == 3, return([P, Q, B])); if(HH == 4, return([T, P, B])); if(HH == 5, return([B, P, Q])); } /* Safe argument. */ safetyarg(z) = if(z == 0, 0, arg(z)); /* Make graph. */ MakeGraph(width, height, x0, y0, x1, y1, filename) = { xstep = (x1 - x0)/width; ystep = (y1 - y0)/height; write(filename, "P3"); write(filename, "# ", filename); write(filename, width, " ", height); write(filename, "255"); for(y=0, height-1, for(x=0, width-1, xx = x0+(xstep*x); yy = y0+(ystep*y); z = xx+yy*I; funcvalue = func(z); mag = abs(funcvalue); phase = safetyarg(funcvalue); H = phase/(2*Pi); S = 1/(1 + 0.3*log(mag + 1)); B = 1 - 1/(1.1 + 5*log(mag + 1)); RGB = HSB2RGB([H, S, B]); Red = floor(RGB[1]*255.0); Green = floor(RGB[2]*255.0); Blue = floor(RGB[3]*255.0); write1(filename, Red, " ", Green, " ", Blue, " "); ); write(filename, ""); ); print("Done."); }