(08/21/2014, 01:15 AM)jaydfox Wrote: Sorry to post two versions in one day, but the performance boost was too huge to ignore. In the table below, the top column shows the number of polynomials I calculated. Below is the time in seconds. As you can see, version 3 cuts the time by 30%-50%, compared to version 2. But version 4 blows them both out of the water, cutting the time by more than an order of magnitude.

Code:

Vers.   60      100     150     200
v2      0.7     9.1     91      --
v3      0.5     5.6     49
v4      0.15    0.5     3.2     15.3

How is this possible? Simple: I got rid of the rational coefficients. (...)

(08/07/2014, 01:19 AM)jaydfox Wrote:

(08/06/2014, 04:38 PM)Gottfried Wrote:

(08/05/2014, 05:57 PM)jaydfox Wrote: I'm working on a set of generator polynomials that would allow me to compute an arbitrary term in the sequence in log-polynomial time. For instance, the naive approach requires O(n) calculations to calculate the nth term.
(...)
I'll work on converting this to PARI/gp code, since I know that several of the active members of this forum seem to use that.

Wow, that's an ingenious approach. I was trying to arrive at a similar short-cutted computation-method using the idea of squaring the matrix M until sufficient few columns are relevant and only compute the required intermediate members of the sequence, but it seemed to become too complex and I left it for the time now. I'll be happy when I can try the Pari/GP-code....

I'll also be much interested if I can learn about the approach how to determine the function f(x) which prognoses/approximates the elements of the sequence. I could implement the function f(x) in Pari/Gp but just did not get the clue how one would arrive at it.
Gottfried

Alright, here's the first draft of PARI/gp code.

(Begin Edit: Adding link to attachment)
(End Edit)

Note the use of a constant V0 in my array indexes. Converting code from a typical language with 0-based arrays, to a language like PARI with 1-based arrays, can be very difficult to debug. So I took the easy way out and made all the code implicitly 0-based, using the V0 constant to add the 1 back in at the last moment. This allows me to port code back and forth between PARI/gp and SAGE/python or c++ or whatever.

Here are sample outputs. The first shows how to calculate elements from the sequence... sequentially.


The next shows how to initialize the generator polynomials and then print out the terms from the sequence with indexes as powers of 2 (i.e., A(1), A(2), A(4), A(, etc.).



The next couple pictures show me calculating the first 30 terms with power of 2 indexes. You can see that the calculation is very fast, a fraction of a second. You can validate the terms up to A(8192) against the list at the OEIS page:
http://oeis.org/A000123/b000123.txt

(08/08/2014, 12:55 AM)jaydfox Wrote: Deriving the continuous function is actually pretty straightforward.

(...)

In general:
\(
D\left[b_k x^k\right] = b_{k-1} \left(\frac{x}{2}\right)^{k-1} \\
\Rightarrow k b_k x^{k-1} = \frac{1}{2^{k-1}} b_{k-1} x^{k-1} \\
\Rightarrow b_k = \frac{1}{k} \left(\frac{1}{2^{k-1}}\right) b_{k-1}
\)

By induction:
\(
b_k = \frac{1}{2^{(k-1)k/2} k!}
\)

Then, as I posted previously:

\(
f(x) = \sum_{k=0}^{\infty}\frac{1}{2^{k(k-1)/2} k!} x^k
\)

(09/09/2014, 07:43 PM)jaydfox Wrote: Treating Gamma(k) at negative integers as infinity, and the reciprocal of such as zero, we can take the limit from negative to positive infinity. And we can replace k with (k+b), where b is zero in the original solution, but can now be treated as any real (well, any complex number, but the complex versions are less interesting).

\(
f_{\beta}(x) = \sum_{k=-\infty}^{\infty}\frac{2^{-(k+\beta)(k+\beta-1)/2}}{\Gamma(k+\beta+1)} x^{k+\beta}
\)

(07/31/2014, 01:51 AM)jaydfox Wrote: \( f(x) = \sum_{k=0}^{\infty}\frac{1}{2^{k(k-1)/2} k!} x^k \)

(08/20/2014, 02:11 PM)sheldonison Wrote: (...) I figured out how to initialize the SeqParts array for a power of 2, and then it is fast for numbers not to much bigger than that exact power of 2. What is the algorithm for initialization halfway to the next power of 2? For example (2^24+2^23)? I don't understand the polynomials well enough to write my own SeqPowers2 initialization starting with 3,6,12,24....

(08/20/2014, 07:57 PM)jaydfox Wrote: I'm still getting around to writing the code that calculates arbitrary terms of the sequence. It's actually a bit complicated, which is why I've been putting it off.

(08/04/2014, 11:49 PM)sheldonison Wrote: I made some numerical approximations, and for big enough values of x, using Jay's Taylor series, I get that \( \ln(f(x)) \approx \frac{(\ln(x))^2}{2\ln(2)} \), which match the equations in Gottfried's link. If ln(x)=1E100, then ln(f(x))~=7.213E199, ln(f(f(x)))~=3.753E399