Mini-Challenge: Rudin-Shapiro Sequence

[John Keith]

A simple programming challenge to stretch your brains before Valentin's next round: Just write the most efficient program you can for your favorite HP calculator to return at least 1000 terms of the Rudin-Shapiro sequence as an array or list. To be clear, I am referring to what the Wikipedia article calls r(n) but you may also submit a program for u(n) if you like for "extra credit"

The winner will be determined by the lowest product of (program size in bytes) * (execution time in seconds).

For the RPL series, speed is to be timed on an HP 50g in approximate mode. For other calculators including the HP-71B and the Prime, execution time can only be compared with other programs on the same model of calculator.

Older calculators can play along too, just compute as many terms as memory allows and scale the execution time accordingly.

[Gerald H]

Edit 2022-12-11: Supplied missing "→". CKSUM & SIZE remain correct.

For 48G, 49G & 50g.

Code:

CKSUM # 9B57h

SIZE 162.

«
  « DUP 2. >
    IF
    THEN DUP 4. / IP SWAP 4. MOD DUP 1. >
      IF
      THEN SWAP DUP + 1. + ←P EVAL SWAP 2. MOD +
      ELSE DROP ←P EVAL
      END
    ELSE DROP 0.
    END
  » → ←P
  « ←P EVAL
  »
»

[C.Ret]

For HP-28S,one algorithm, three submissions.

Version with extra \(u_k\) computations:

Uk:
« 0 SWAP BIN R→B →STR OVER
  WHILE 1 + 999 SUB DUP "11" POS DUP
  REPEAT ROT 1 + ROT ROT END
  DROP2 »

Usage:375775 Uk returns 10.

Rk:
« -1 SWAP Uk ^ »

LST1000:
« 0 999 FOR k k Rk NEXT
  1E3 →LIST
  960 .1 BEEP »

The list of the first 1000 terms of the Rudin-Shapiro sequence obtained in about 2'46".


Version without extra \(u_k\):
Rk:
« 1 SWAP BIN R→B →STR 2
  WHILE 1 + 999 SUB DUP "11" POS DUP
  REPEAT ROT NEG ROT ROT END
  DROP2 »

Usage:375775 Rk returns 1.

LST1000:
« 0 999 FOR k k Rk NEXT
  1E3 →LIST
  960 .1 BEEP »

The list of the first 1000 terms of the Rudin-Shapiro sequence obtained in about 2'26".


Stand alone version:
LST1000:
« 0 999 FOR k
    1 k BIN R→B →STR 2
    WHILE 1 + 999 SUB DUP "11" POS DUP
    REPEAT ROT NEG ROT ROT END
    DROP2 NEXT
  1E3 →LIST
  960 .1 BEEP »

The list of the first 1000 terms of the Rudin-Shapiro sequence obtained in about 2'18".


Next challenge, trying to improuve speed efficiency of the list build-up by using the recursive algorithm.

P.S.: All system settings expected to be default values, especially binary words size at 64 bits.

[John Keith]

(12-10-2022 05:54 AM)Gerald H Wrote:  «
« DUP 2. >
IF
THEN DUP 4. / IP SWAP 4. MOD DUP 1. >
IF
THEN SWAP DUP + 1. + ←P EVAL SWAP 2. MOD +
ELSE DROP ←P EVAL
END
ELSE DROP 0.
END
» → ←P
« ←P EVAL
»
»[/code]

There seems to be a missing right arrow before the second from last ←P as indicated above. It now runs and returns u(n).

Neat program, I'm still trying to wrap my head around compiled local variables.

[Gerald H]

Correct, will re-edit.

Yes, compiled local variables are useful & the programme is meant as a demonstration of their utility, sadly I can't think of a previous example in programmes available on the forum.

[Gerald H]

I've found a previous reference to "compiled local variables".

https://www.hpmuseum.org/forum/thread-14...#pid125758

[Allen]

Not quite to 1000, but I found a LCG with multiplier 14, increment 2 and prime 3449 that to my surprise will generate the first 33 bit entries of the Zero-one version of Golay-Rudin-Shapiro sequence in only 44 bytes, including the "RSS" label.

Code:


00 { 44-Byte Prgm }
01>LBL "RSS"
02 2
03 STO 00
04 33
05 STO 01
06>LBL 00
07 RCL 00
08 14
09 ×
10 2
11 +
12 3449
13 MOD
14 STO 00
15 64
16 ÷
17 1
18 AND
19 PRX
20 DSE 01
21 GTO 00
22 .END.

Yields:

Code:


           0.0000    ***
           0.0000    ***
           0.0000    ***
           1.0000    ***
           0.0000    ***
           0.0000    ***
           1.0000    ***
           0.0000    ***
           0.0000    ***
           0.0000    ***
           0.0000    ***
           1.0000    ***
           1.0000    ***
           1.0000    ***
           0.0000    ***
           1.0000    ***
           0.0000    ***
           0.0000    ***
           0.0000    ***
           1.0000    ***
           0.0000    ***
           0.0000    ***
           1.0000    ***
           0.0000    ***
           1.0000    ***
           1.0000    ***
           1.0000    ***
           0.0000    ***
           0.0000    ***
           0.0000    ***
           1.0000    ***
           0.0000    ***
           0.0000    ***

[Allen]

Output with the LCG internal state and the bit mask for a_n = bit7 .. except for a few entries, it's correct up through a_50

Code:


n  LCG       LCG bits    correct?
-2    0 00000 0 000000 N/A
-1    2 00000 0 000010 N/A   <-- initial state from program 
 0   30 00000 0 011110 True  <-- a_0 = 0
 1  422 00011 0 100110 True
 2 2461 10011 0 011101 True
 3 3415 11010 1 010111 True
 4 2975 10111 0 011111 True
 5  264 00010 0 001000 True
 6  249 00001 1 111001 True
 7   39 00000 0 100111 True
 8  548 00100 0 100100 True
 9  776 00110 0 001000 True
10  519 00100 0 000111 True
11  370 00010 1 110010 True
12 1733 01101 1 000101 True
13  121 00000 1 111001 True
14 1696 01101 0 100000 True
15 3052 10111 1 101100 True
16 1342 01010 0 111110 True
17 1545 01100 0 001001 True
18  938 00111 0 101010 True
19 2787 10101 1 100011 True
20 1081 01000 0 111001 True
21 1340 01010 0 111100 True
22 1517 01011 1 101101 True
23  546 00100 0 100010 True
24  748 00101 1 101100 True
25  127 00000 1 111111 True
26 1780 01101 1 110100 True
27  779 00110 0 001011 True
28  561 00100 0 110001 True
29  958 00111 0 111110 True
30 3067 10111 1 111011 True
31 1552 01100 0 010000 True
32 1036 01000 0 001100 True
33  710 00101 1 000110 False  <-- starts going off the rails here
34 3044 10111 1 100100 False
35 1230 01001 1 001110 True
36 3426 11010 1 100010 False
37 3129 11000 0 111001 True
38 2420 10010 1 110100 True
39 2841 10110 0 011001 True
40 1837 01110 0 101101 True
41 1577 01100 0 101001 True
42 1386 01010 1 101010 False
43 2161 10000 1 110001 True
44 2664 10100 1 101000 True
45 2808 10101 1 111000 True
46 1375 01010 1 011111 False
47 2007 01111 1 010111 True
48  508 00011 1 111100 True
49  216 00001 1 011000 True
50 3026 10111 1 010010 True
51  978 00111 1 010010 False
52 3347 11010 0 010011 False
53 2023 01111 1 100111 True
54  732 00101 1 011100 False
55 3352 11010 0 011000 False
56 2093 10000 0 101101 True
57 1712 01101 0 110000 True
58 3276 11001 1 001100 False
59 1029 01000 0 000101 False
60  612 00100 1 100100 True
61 1672 01101 0 001000 False

[John Keith]

Very cool program, Allen! I like those off-the-wall solutions.

[John Keith]

Since interest in this challenge seems to have waned, I will post my solutions tomorrow afternoon. I have a 97 byte program which runs on any RPL calculator, and a faster 55 byte program which requires the 48G or newer. Let's see if anyone can do better, you have 24 hours!

[John Keith]

As promised, here are my solutions to this challenge.

First, a 97 byte program for any RPL calculator. This is a simple implementation of the recurrence from the Wikipedia article. Given a number n on the stack, returns 2n+1 terms. An input of 500 returns 1001 terms in 3.25 seconds on my 50g.

Code:


\<< \-> n
  \<< 1 1 1 n
    FOR k k PICK DUP k 2 MOD
     \<< NEG
     \>> IFT
    NEXT n 2 * 2 + \->LIST
  \>>
\>>

Next, a larger but faster version which avoids the 2 MOD test but requires extra code to avoid problems with odd number inputs. 121 bytes. An input of 500 returns 1000 terms in 1.94 seconds on my 50g.

Code:


\<< DUP DUP 2 MOD + \-> n m
  \<< 1 1 1 m
    FOR k k PICK DUP NEG k 1. + PICK DUP 2
    STEP m 2 * 2 + \->LIST 1 n 2 * SUB
  \>>
\>>

Next, a completely different sort of program which requires the 48G or newer because it uses the listability of functions (called "automatic list processing" in the manuals). Given a number n on the stack, returns 2^n terms. For example an input of 10 will return 1024 terms in 1.22 seconds on my 50g. 55 bytes.

Code:


\<< { 1 } DUP ROT 2 SWAP   @ Do n-1 iterations
  START DUP2 + ROT ROT     @ Concatenate and move to top
    0 SWAP - +             @ Negate second part and concatenate
  NEXT +                   @ Concatenate last two lists
\>>

Finally, a similar program for u(n) aka A014081. Also returns 2^n terms, 55.5 bytes. An input of 10 returns 1024 terms in 1.19 seconds.

Code:


\<< { 0 } DUP ROT 2 SWAP
  START DUP2 + ROT ROT 1 ADD +
  NEXT +
\>>

[ttw]

The Thue-Morse sequence isn't that hard (it's very fast on the Hp50g as one can compute en masse with X = {0 1} and follow up with X = X + (1 - X). The Nth element is the parity of the number of bits set in the binary representation of N. The anti-Thue sequence (I don't know an accepted name) is the sequence that counts the number of 0's in the number N.

[Allen]

@ John,
Thank you for an interesting problem!! These challenges are always a good cause to learn more about other areas of math I knew little about. While I was trying to squeeze more juice out of this, I learned some interesting things about the Shapiro Polynomials and some other items.

One of the potential avenues I explored was using run length encoding for A020987. There are some really interesting properties of this function and these polynomials they represent!!

Regarding the run length, for n up to 2**20 there are no more than 4 0's or 1' in a row. I think it's possible to prove that there are never more than 4 in a row, and Brillhard and Morton have made a really heavy dent in proving some related theorems- generating 16 pages of math, but without a grand unified proof.

(humorous aside: the authors of this paper say on p.3

"We find ourselves being distracted from the original question by our conjecture, which is interesting in it's own right. Let's indulge ourselves and pursue this pattern question, ignoring the original question for the time being."

The irony of this statement, both in the context of this series, this original forum post, and their paper seems appropriately recursive! In fact, this may be one of the best quotes I've ever seen with respect to applied math research.

Perhaps the second-best math quote I've seen is on p. 15: "What strange numbers!" - Seems like the kind of thing Ramanujan or Erdős, or Hawking might put on their gravestones.)


For the first 1024*1024, here are the counts of the run lengths:

Code:


 0    262144 -> 2**18
 2    131072 -> 2**17
 3    65536   -> 2**16
 1    65536   -> 2**16

This seems a beautiful coincidence!

In the original sequence count for 2**20 is not quite as nice, but close:

Code:


0  524800   -> 2**19 + 2**9
1  523776   -> 2**19 - 2**9

For various reasons, I encoded a run length of 1 as 0 and 2 as 1 and so on, so the first 20 entries of A020987 would be run-length encoded as:

Code:


[0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 1, 0, 1, 0, 0, 0, 1] -> [2, 0, 1, 0, 3, 2, 0, 0, 2, 0]

when packed as binary encoded, the sequence for the first 2*20 A020978 entries fed into zlib.compress yields only 2791 bytes. Due to the semi-cyclical nature of even the 1X compressed data, feeding those bytes back into zlib again returns only 274 bytes!!

A 500-fold reduction in size (when 2X compressed) can sometimes be a bellwether of deeper mathematical mysteries- or at least some very long repeated sub-sequences.

Behold, 2**20 entries of A020987 are here:

Code:


789cab98f3f6cac1ac24461196c079d57fff9ff775611174693aa158b2562db9272aecbaaddff19f​
16b17ec7dfef9ffef9ebf3dbfff59fbf4ede7e5cbebebef6cfbf3dc79feffff5fbf7e7f9f1efafbe​
fefff5effcf8ff779fff7efdfcbfd9ffdde7ffbffcfefdedfbf8fdff97ffb3fffbefebbf2fbfff97​
ff9f15ff7dfdf7f53fafffcdfef7ebfcfe2febdfdfb7da7f6fffb3ebff5edbff3ffdcf74ffbefa59​
f1bfd7ffcf8529fcbcff6dfce7d7bf97ffdcfaff2b501aa2b4fef6fe67ebfffff1feffeb5fe6bf52​
3c66625788d54cecd653e64ea2cd24c19d44fb7d80c39316eea44578d222de4772fa1c2ae13954d2​
e750c9ef43257d0e95f01c4d9f23337d8e969fd44e9ffff97f2ffab205001f2570d6

Or I could examine 1 line of python and stop procrastinating my Saturday morning house cleaning:

Code:


A020987_gen=[[x for x in zip(bin(x),bin(x)[1:])].count(('1','1'))%2 for x in range(limit)]

[John Keith]

(12-17-2022 02:30 AM)ttw Wrote:  The Thue-Morse sequence isn't that hard (it's very fast on the Hp50g as one can compute en masse with X = {0 1} and follow up with X = X + (1 - X). The Nth element is the parity of the number of bits set in the binary representation of N. The anti-Thue sequence (I don't know an accepted name) is the sequence that counts the number of 0's in the number N.

The sequence that counts the number of 0's in binary n is A023416. The parity of the number of 0's is A059448. I don't know if there are special names for either sequence. Both can be computed by similar small, fast programs.

This was actually the reason behind my posting of this challenge. I have observed that nearly all sequences describing properties of binary numbers have a fractal pattern, that is they are self-similar at different scales. This self-similarity leads to fast and simple algorithms as you noted for the Thue-Morse sequence.

[John Keith]

(12-17-2022 04:03 PM)Allen Wrote:  @ John,
Thank you for an interesting problem!! These challenges are always a good cause to learn more about other areas of math I knew little about. While I was trying to squeeze more juice out of this, I learned some interesting things about the Shapiro Polynomials and some other items.

One of the potential avenues I explored was using run length encoding for A020987. There are some really interesting properties of this function and these polynomials they represent!!

Regarding the run length, for n up to 2**20 there are no more than 4 0's or 1' in a row. I think it's possible to prove that there are never more than 4 in a row, and Brillhard and Morton have made a really heavy dent in proving some related theorems- generating 16 pages of math, but without a grand unified proof.

...snip...

For the first 1024*1024, here are the counts of the run lengths:

Code:


 0    262144 -> 2**18
 2    131072 -> 2**17
 3    65536   -> 2**16
 1    65536   -> 2**16

This seems a beautiful coincidence!

In the original sequence count for 2**20 is not quite as nice, but close:

Code:


0  524800   -> 2**19 + 2**9
1  523776   -> 2**19 - 2**9

For various reasons, I encoded a run length of 1 as 0 and 2 as 1 and so on, so the first 20 entries of A020987 would be run-length encoded as:

Code:


[0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 1, 0, 1, 0, 0, 0, 1] -> [2, 0, 1, 0, 3, 2, 0, 0, 2, 0]

Thanks for your interesting observations! I don't think that the patterns you found in the run length are coincidences at all, see my post above in response to ttw. Of course, I don't have a grand unified theory to explain the patterns either.

As it happens, the way I discovered the pattern that led to my HP48G program is by looking at A020987. The Rudin-Shapiro sequence (A020985) is a confusing wall of ones and minus signs but seeing it as a sequence of ones and zeros made the pattern much clearer visually.

One additional note on the subject of run-length encoding: the ListExt commands LRPCT AND LRPC→ perform run-length encoding and decoding respectively. They are very fast.

[Allen]

A 42S program to print (n...1) entries: ( for n up to 1023)

Code:


00 { 48-Byte Prgm }
01 STO 01
02>LBL 00
03 RCL 01
04 PRX
05 STO 00
06 2
07 ÷
08 RCL 00
09 AND
10 ENTER
11 STO 00
12 8
13 ÷
14 RCL 00
15 64
16 ÷
17 XOR
18 RCL 00
19 XOR
20 7
21 AND
22 150
23 X<>Y
24 ROTXY
25 1
26 AND
27 PRX
28 DSE 01
29 GTO 00
30 .END.

The bit twiddling was adapted from here.

Code:


    # find consecutive ones
    x =  x & (x>>1)
  
    #then get parity
    x ^= (x >> 3)
    x ^= (x >> 6)
    x &= 0x7
    return (0x96 >> x) % 2

[Allen]

(12-16-2022 10:41 PM)John Keith Wrote:  

Code:


\<< { 1 } DUP ROT 2 SWAP   @ Do n-1 iterations
  START DUP2 + ROT ROT     @ Concatenate and move to top
    0 SWAP - +             @ Negate second part and concatenate
  NEXT +                   @ Concatenate last two lists
\>>

Some nice programming! I just noticed you can save 10% (5 bytes) of the program size by using the NEG function rather than "0 SWAP -"

(50 bytes)

Code:


\<< { 1 } DUP ROT 2 SWAP   @ Do n-1 iterations
  START DUP2 + ROT ROT     @ Concatenate and move to top
    NEG +                       @ NEG-ate second part and concatenate
  NEXT +                    @ Concatenate last two lists
\>>

[John Keith]

(12-18-2022 08:30 PM)Allen Wrote:  I just noticed you can save 10% (5 bytes) of the program size by using the NEG function rather than "0 SWAP -"

(50 bytes)

Code:


\<< { 1 } DUP ROT 2 SWAP   @ Do n-1 iterations
  START DUP2 + ROT ROT     @ Concatenate and move to top
    NEG +                       @ NEG-ate second part and concatenate
  NEXT +                    @ Concatenate last two lists
\>>

As strange as it may seem, 0 SWAP - is faster than NEG for lists of numbers. I considered the speed gain to be worth the extra 5 bytes since the operation is inside the loop. Also, on the 49/50 you can replace ROT ROT with UNROT and save another 2.5 bytes.

[Gilles]

newRPL HP50g:

Code:

«
  { 1 } DUP ROT
  2 SWAP START
    DUP2 ADD UNROT NEG ADD 
  NEXT
  ADD 
»

Nota : + and ADD are switched in NewRPL, wich is more logical but loose compatibility

Size: 68 bytes

1024 items, Time : 0.022 sec (usb connected)
8192 items, Time : 0.161 sec

PC newRPL : 65536 items, Time : 0 sec (!!)

[John Keith]

Dang that's fast!