09/11/2014, 01:33 AM

As promised previously, here's a version that initializes the SeqParts array to an arbitrary non-negative index in the sequence:

jdf_seq_v6.gp (Size: 13.35 KB / Downloads: 451)

Okay, I've updated quite a few functions, mostly to make them compatible with the new InitSeq() functions (which replace InitBinPartSeq() ).

There's a new PrevSeq() function, which runs the sequence backwards. As such, the Seq() function can now take negative arguments, and will handle them accordingly.

The InitSeq() function should demonstrate how the polynomials are used to calculate terms in the sequence. It is not optimized, but the sub-calls are, so it should still be reasonably fast.

The InitSeq2() function is optimized. There are two main optimizations. First, the SeqPowers2() function is used to initialize the SeqParts array for the highest power of 2 in the provided index. Second, "negative" bits are used to reduce the number of calculations required. Between these two optimizations, you can quickly calculate terms which have low information.

For example, InitSeq2(2^100) takes about 10 ms on my computer. InitSeq2(2^100 - 2^95) takes about 90 ms. However, InitSeq2(2^100\3) takes about 640 ms. Note that the binary representation of 2^100\3 has alternating 1's and 0's, and the use of negative bits cannot simplify this.

For the InitSeq2() funciton, set the "show" parameter to 2 to view the bits. (Set the show parameter to "3" to get a printout of the SeqParts array at completion of the function call.) The output will look something like this:

In the example above, the first size/bits printout is using only positive bits. The second size/bits printout shows the same number, using negative bits. The left-most bit is the least significant.

Once the sequence is initialized with InitSeq() or InitSeq2(), you can then use Seq() with positive or negative counts, or use NextSeq() and PrevSeq() directly. For example, to calculate the terms from 2^13-5 to 2^13+5 (which you can validate against this list at OEIS):

Now let's check out 2^60-5 through 2^60+5:

Final disclaimer: this was a complicated set of logic to pull together, and it doesn't help that I'm still getting used to indexing arrays from 1 instead of from 0 (PARI/gp is an outlier in this regard). So, I expect there are still a few bugs. Indeed, I had to re-edit the file several times while writing this post, as I found bugs while making screenshots or validating what I was writing. If I find any more bugs in the coming days, I'll post a small update, perhaps version "6b".

jdf_seq_v6.gp (Size: 13.35 KB / Downloads: 451)

(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.

Okay, I've updated quite a few functions, mostly to make them compatible with the new InitSeq() functions (which replace InitBinPartSeq() ).

There's a new PrevSeq() function, which runs the sequence backwards. As such, the Seq() function can now take negative arguments, and will handle them accordingly.

The InitSeq() function should demonstrate how the polynomials are used to calculate terms in the sequence. It is not optimized, but the sub-calls are, so it should still be reasonably fast.

The InitSeq2() function is optimized. There are two main optimizations. First, the SeqPowers2() function is used to initialize the SeqParts array for the highest power of 2 in the provided index. Second, "negative" bits are used to reduce the number of calculations required. Between these two optimizations, you can quickly calculate terms which have low information.

For example, InitSeq2(2^100) takes about 10 ms on my computer. InitSeq2(2^100 - 2^95) takes about 90 ms. However, InitSeq2(2^100\3) takes about 640 ms. Note that the binary representation of 2^100\3 has alternating 1's and 0's, and the use of negative bits cannot simplify this.

For the InitSeq2() funciton, set the "show" parameter to 2 to view the bits. (Set the show parameter to "3" to get a printout of the SeqParts array at completion of the function call.) The output will look something like this:

Code:

`(15:02) gp > InitSeq2(2^27-2^23-2^17+2^11-2^7, 2)`

size=27, bits=List([0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 0])

size=28, bits=List([0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 1, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, -1, 0, 0, 0, 1])

A(125699968)=13989020196114224701143563476000978180153807452722418367111885874553950782754824115506668

(15:03) gp >

In the example above, the first size/bits printout is using only positive bits. The second size/bits printout shows the same number, using negative bits. The left-most bit is the least significant.

Once the sequence is initialized with InitSeq() or InitSeq2(), you can then use Seq() with positive or negative counts, or use NextSeq() and PrevSeq() directly. For example, to calculate the terms from 2^13-5 to 2^13+5 (which you can validate against this list at OEIS):

Now let's check out 2^60-5 through 2^60+5:

Final disclaimer: this was a complicated set of logic to pull together, and it doesn't help that I'm still getting used to indexing arrays from 1 instead of from 0 (PARI/gp is an outlier in this regard). So, I expect there are still a few bugs. Indeed, I had to re-edit the file several times while writing this post, as I found bugs while making screenshots or validating what I was writing. If I find any more bugs in the coming days, I'll post a small update, perhaps version "6b".

~ Jay Daniel Fox