(08/20/2014, 02:11 PM)sheldonison Wrote: I downloaded your memory efficient code, jdf_seq_v2.gp, and it works great 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....

Okay, here's the latest version:

jdf_seq_v3.gp (Size: 8.45 KB / Downloads: 450)

Note that for SeqPowers2 and SeqPowers2m, it converts large numbers to floating point for display purposes. To turn this off, set jPrec=0

You may have noticed that the graphs I've been posting in this discussion use matplotlib. I prefer graphing in SAGE/python with matplotlib. But my sequence code is in PARI/gp. I would actually prefer to write all my code in SAGE/python, but I'm trying to keep my code accessible to the group here, which means PARI/gp.

My workaround is to export the data from gp to a text file. I do this in my code by setting jWrite=1 when I'm ready to export data. I then import the data into SAGE, then save it to a SAGE native format for reuse later. From there, I can make all my graphs. (Sheldon, if you know of a better way to write stuff to disk from gp, please let me know.)

Anyway, updates in this code:

1. I created a slightly modified version of the NextSeq() function, aptly named NextSeq2(). It accomplishes the same effect as the previous version, but the bit checking is faster.

2. InitPolys will reuse previously calculated polynomials. It will only add what it needs to. If you run InitPolys(20), then InitPolys(25), it will only calculate 21..25 in the second call. Because of this, you can safely call InitPolys from other functions, to ensure that all required polynomials are available. See SeqPowers2, for example.

3. The new function SeqPowers2m, which takes two main arguments: n and m. The n parameter is the max exponent; the m parameter is also an exponent, but controls a sequential counter. The function calculates k*2^j, where [2^(m-1)+1 <= k <= 2^m], and [1 <= j <= n]. There is an optional parameter start, which limits j to the range [start <= j <= n].

Effectively, this calculates a grid of values, from about 2^(start+m-1) to 2^(n+m).

4. I made a bug fix and a couple performance enhancements to InitPolys:

4.1: The previous version scaled the polynomials incorrectly (they were reduced by one too many powers of 2). For example, the following code should always return 1, for any value of k:

Code:

`subst(SeqPol[k+1], x, 2^(k-1))`

The defect wasn't noticeable, because the code had a compensating "bug". In this case, two wrongs made a right. Anyway, both bugs are fixed now.

4.2: I realized that when calculating the points for polynomial interpolation, I had a running sum of terms, only to take the differences as the first step of interpolation. This seemed inefficient. Therefore, I don't calculate the values directly, but the first differences, and the interpolation code skips the first differences and starts on the second differences.

4.3: This is the big performance boost. I realized that I could reuse almost half the points from the previous polynomial grid. This saves some of the very expensive computations (high-degree polynomials with large rational coefficients). Speedup is roughly 40%. (Begin Edit) I.e., it takes 40% less time. Depending on your point of view, that means the speedup is roughly 70% ~= 1/(1-.4) (End Edit)

(Begin Edit)

5. I forgot to mentioned another big change in SeqNum: while I was writing SeqPowers2m, I found a pattern for calculating A(2^n) for integer n, which requires only one polynomial evaluation (as opposed to n evaluations). I haven't done a comparison between the old and new version, but I assume the speedup is significant, especially when n is very large.

(End Edit)

There are probably still a few bugs in the code, so if you run into any, please let me know.

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.

~ Jay Daniel Fox