OEIS featured in The New York Times

[Valentin Albillo]

.
Hi, all,

Today (Sunday, May 21th, 2023) The New York Times features an article commemorating OEIS (On-Line Encyclopedia of Integer Sequences) 50th anniversary:

What Number Comes Next? The Encyclopedia of Integer Sequences Knows

It makes for a good read and it's nice to see a math website being featured in a main USA newspaper. Have a look at a sequence mentioned in the article, namely:

The Sisyphus sequence: start the sequence S with a(1) = 1 and extend S with a(n)/2 when a(n) is even, otherwise with a(n) + the smallest prime not yet added: 1, 3, 6, 3, 8, 4, 2, 1, 8, 4, 2, 1, 12, 6, 3, 16, 8, 4, 2, 1, 18, ...

Generating the sequence is trivial, here's for instance code for the HP-71B:

SISYPHUS 76 bytes

10 DESTROY ALL @ STD @ N=1 @ P=1
20 LOOP @ DISP N; @ IF NOT MOD(N,2) THEN N=N DIV 2 ELSE P=FPRIM(P+1) @ N=N+P
30 END LOOP

>RUN
1 3 6 3 8 4 2 1 8 4 2 1 12 6 3 16 8 4 2 1 18 9 28 14 7 30 15 44 22 ...

and there are a number of interesting things to explore in the sequence.
V.

[Joe Horn]

(05-21-2023 11:44 PM)Valentin Albillo Wrote:  20 LOOP @ DISP N; @ IF NOT MOD(N,2) THEN N=N DIV 2 ELSE P=FPRIM(P+1) @ N=N+P

Cool! Unfortunately, my memories of the 71B are rapidly fading with age. Is FPRIM built-in, or is it in the Math ROM, or in the JPC ROM, or something else? Thanks for refreshing my memory.

[Valentin Albillo]

(05-22-2023 12:21 AM)Joe Horn Wrote:  

(05-21-2023 11:44 PM)Valentin Albillo Wrote:  20 LOOP @ DISP N; @ IF NOT MOD(N,2) THEN N=N DIV 2 ELSE P=FPRIM(P+1) @ N=N+P

Cool! Unfortunately, my memories of the 71B are rapidly fading with age. Is FPRIM built-in, or is it in the Math ROM, or in the JPC ROM, or something else? Thanks for refreshing my memory.

You're welcome. FPRIM is in the fantastic JPC ROM, as are LOOP and END LOOP.

Best regards.
V.

[Fernando del Rey]

This recent post from Valentín seems to have gone unfairly unnoticed. It's not one of his great challenges, but it's still an interesting topic to the calculator-minded or mathematically-minded folk that hang out in this forum, I believe.

While I was not able to read the NYT article, which is blocked to non-subscribers, I've been exploring the OEIS site which was new to me.

It's quite a source of interesting and fun stuff!

It would have been fantastic to have OEIS during the heyday of the HP calcs, when we were looking for any idea we could transform into a program in our beloved little toys.

Thanks, Valentín, for posting this entertaining and informative thread!

[jonmoore]

I always enjoy Valentin's posts, so not I'm not certain how I missed this one.

I've been aware of OEIS for some time, mainly because of Mathematica and Maple listings, but it was nice to see Valentin's approach to that particular number sequence using the 71B. There's always something new to learn with the manner in which Valentin approaches things. Not that I've kept his program listing for adding them the missing trigonometric functions on any of my 12C's, I was totally transfixed whilst studying the listing.

I hadn't seen the NYT piece either. I used to get access via my Apple One subscription, but the NYT pulled out a few years back. Living in the UK and being a lifetime reader (and subscriber for most of it) of The Guardian, I wasn't going to pay the NYT directly for the privilege! The Washington Post has always been more my US newspaper taste anyway.

On a separate note, concerning the JPC ROM. I've been tempted to download the compendiums of Paris Chapter 48 programs that Eric publishes on hpcalc.org. But there's no handy index (like the one for Joe's Goodies), so it would involve an awful lot of manual sifting. Does anyone know if there's anything worth investigating within the various discs that isn't already published on the rest of hpcalc.org?

No worries if not, it's just that I've read great things about the Paris chapter over the years.

[J-F Garnier]

(05-30-2023 06:05 PM)jonmoore Wrote:  On a separate note, concerning the JPC ROM. I've been tempted to download the compendiums of Paris Chapter 48 programs that Eric publishes on hpcalc.org.

What are you referring to, exactly?
I can't find anything about the JPC ROM or the PPC Paris Chapter on hpcalc.org.

J-F

[jonmoore]

(05-30-2023 06:25 PM)J-F Garnier Wrote:  

(05-30-2023 06:05 PM)jonmoore Wrote:  On a separate note, concerning the JPC ROM. I've been tempted to download the compendiums of Paris Chapter 48 programs that Eric publishes on hpcalc.org.

What are you referring to, exactly?
I can't find anything about the JPC ROM or the PPC Paris Chapter on hpcalc.org.

J-F

The JPC ROM is for the 71B and I first encountered it here:

http://www.jeffcalc.hp41.eu/emu71/jpcrom.html

When I saw it mentioned in this thread, it acted as an aide-mémoire ref the PPC-Paris User Club.

For the HP historians, this is a great read - it's a link within the above page, but it's easy to miss:

http://www.jeffcalc.hp41.eu/emu71/files/jpcromstory.pdf

[rprosperi]

(05-30-2023 07:19 PM)jonmoore Wrote:  

(05-30-2023 06:25 PM)J-F Garnier Wrote:  What are you referring to, exactly?
I can't find anything about the JPC ROM or the PPC Paris Chapter on hpcalc.org.

J-F

The JPC ROM is for the 71B and I first encountered it here:

http://www.jeffcalc.hp41.eu/emu71/jpcrom.html

When I saw it mentioned in this thread, it acted as an aide-mémoire ref the PPC-Paris User Club.

For the HP historians, this is a great read - it's a link within the above page, but it's easy to miss:

http://www.jeffcalc.hp41.eu/emu71/files/jpcromstory.pdf

So, only fair to warn you that the person who created that website is the person you are answering...

[jonmoore]

(05-30-2023 07:48 PM)rprosperi Wrote:  So, only fair to warn you that the person who created that website is the person you are answering...

Thanks Bob, I had a suspicion by the barbed language and Gallic surname that the French Connection was a distinct possibility. Tell him, I adore his emulator.

[Valentin Albillo]

      
Hi again,

Thanks a lot to Fernando del Rey and jonmoore for their interest and very kind words, much appreciated. Re the OEIS (On-Line Encyclopedia of Integer Sequences), I heartily recommend reading these PDF documents to all people new to it:

On-Line Encyclopedia of Integer Sequences (25-page PDF)
The On-Line Encyclopedia of Integer Sequences - An illustrated guide with many unsolved problems (37-page PDF)
A Handbook of Integer Sequences - Fifty Years Later (23-page PDF)

Now for a little self-quoting:

I Wrote:[...] Generating the sequence is trivial [...]

1 3 6 3 8 4 2 1 8 4 2 1 12 6 3 16 8 4 2 1 18 9 28 14 7 30 15 44 22 ...

and there's a number of interesting things to explore in the sequence.

Indeed there is, and presently I'll mention a couple. First of all, the sequence is conjectured to include all positive integers 1, 2, 3, ..., as elements, and we can investigate this alleged fact by using this 5-liner for the HP-71B, which takes as input K, the maximum number of elements to generate, and will output the pairs (number, index) for up to the first 60 numbers (1..60), where index is the lowest one where number first appears in the sequence.

It will also output how many numbers weren't found within the first K elements (plus timing) and finally it will explicitly list all numbers not found:

10 DESTROY ALL @ INPUT K @ INTEGER D(K),M(50) @ N=1 @ P=1 @ SETTIME 0
20 FOR C=1 TO K @ IF N>K THEN STD ELSE IF NOT D(N) THEN D(N)=C
30 IF MOD(N,2) THEN P=FPRIM(P+1) @ N=N+P ELSE N=N DIV 2
40 NEXT C @ C=0 @ FOR I=1 TO MIN(K,60) @ IF D(I) THEN DISP (I,D(I)); ELSE C=C+1 @ M(C)=I
50 NEXT I @ DISP @ DISP "Not found:";C;TIME$ @ FOR I=1 TO C @ DISP M(I); @ NEXT I @ DISP

Let's make a sample run specifying a maximum of 100 elements:

>RUN ? 100


(1,1) (2,7) (3,2) (4,6) (5,71) .. (48,42) (51,58) (53,73) (55,88) (58,33)
Not found: 21 (timing)
13 17 19 25 26 27 32 33 34 36 38 41 47 49 50 52 54 56 57 59 60

This means that 1 first appears at index 1, 2 at index 7, 3 at index 2, 4 at index 6, 5 at index 71, and so on until 55 appears at index 88 and 58 at index 33. We also learn that as many as 21 numbers from 1..60 didn't appear among the first 100 elements of the sequence, the first of them being 13, 17, ... and the last ones being ..., 59, 60.

Of course, increasing the maximum number of elements to explore will probably result in finally locating the first occurrence of some or all the numbers not found within the first 100 elements of the sequence. Indeed, running the above program for 100, 200, ...., 2000 elements we get these results:

#Elem  #  Not found
--------------------------------------------------------------------------
 100  21  13 17 19 25 26 27 32 33 34 36 38 41 47 49 50 52 54 56 57 59 60
 200  12  13 25 26 27 32 36 41 49 50 52 54 60
 500   8  25 27 32 36 49 50 54 60
1000   6  25 27 36 50 54 60
2000   4  27 36 54 60

and in fact only 4 numbers from 1..60 failed to make an appearance within the first 2,000 elements, corroborating the conjecture that all positive integers will eventually appear in the sequence.

Now we can try to locate the first occurrence of the 4 missing numbers (27, 36, 54, 60) in a faster way and using minimal memory by running this 4-line HP-71B program, which will accept the number to locate and will search for it among the first 1,000,000 elements (if not using a very fast emulator, you might want to reduce this upper limit to save time; generating a million elements requires primes up to 4,761,697):

10 DESTROY ALL @ STD @ INPUT K @ N=1 @ P=1 @ SETTIME 0
20 FOR C=1 TO 1000000 @ IF N=K THEN DISP N;C;TIME$ @ END
30 IF MOD(N,2) THEN P=FPRIM(P+1) @ N=N+P ELSE N=N DIV 2
40 NEXT C @ DISP "Not found: ";TIME$

Let's try a sample run with an easily found number, 25:

>RUN
? 25 ->  25  1154 (timing)

which quickly tells us that 25 appears in the sequence at index 1,154. Now let's consider the four numbers which weren't located earlier within the first 2,000 elements, namely 27, 36, 54 and 60, to try and locate their first occurrence in the sequence, if possible:

>RUN
? 27 ->  27  161336  (Emu71/Win @972x:  28", go71b @128x:  3'36", Physical: 7h 41')
? 60 ->  60  614667  (Emu71/Win @972x: 3'3", go71b @128x: 23'14", Physical: ~ 50h)

so 27 appears for the first time in the sequence at index 161,336 and 60 at index 614,667. As for 54, being twice 27, we don't need to search for it, its index will be one less than the one for 27, namely 161,335.

Unfortunately, the search up to index 1,000,000 fails for 36 because its first occurrence in the sequence happens to be at index 77,534,485,877 !!. Even worse, if considering numbers up to 100 instead of up to 60, the number 97 first appears at index 17,282,073,747,557 !!!

Additionally, a second conjecture is that every positive integer appears in the sequence not just once but an infinite number of times. We can check it out by finding multiple indexes for any given input number, simply delete the END statement at line 20 and reduce to 100,000 the maximum index to search up to. Line 20 will then look like this:

20 FOR C=1 TO 100000 @ IF N=K THEN DISP N;C;TIME$

Let's try a few selected numbers (you may break execution before reaching index 100,000, to save time):

>RUN
?  1  ->  1  8  12  20  742 ...             (the very next appearance is at index 513,152,128)
?  2  ->     7  11  19  741 ...             (ditto at index 513,152,127)
?  3  ->  2  4  15  46  95  6355 ...
?  5  ->  71  4849 ...                            (no others up to index 1,000,000)
?  7  ->  25  114  123  446  7104 ...
? 13  ->  345  418  4621 ...                  (no others up to index 1,000,000)
? 25  ->  1154  1519  10359  13330 ... (no others up to index 1,000,000)
? 85  ->  140  3161  72349...

so it seems that indeed there might be infinite appearances of any positive integer but the distance from one appearance to the next increases at an extremely fast rate.

Additional Trivia:

● A few other nice appearances are 31416 at indexes 6,768 and 6,923, 11111 at 12,497, 55555 at 56,551, 100000 at 26,488 and last but not least, 2023 at 2,165. On the other hand, 1992 does not appear within the first 1,000,000 elements.

● There's a number of solutions of A(n) = n, i.e. numbers whose index in the sequence equals the number itself. They can be found very easily with a trivial modification of my second program above, and the first ones are n = 1, 16, 787, 427447 and no others up to index 1,000,000.

Regards.
V.

[jason32]

New York Times featuring the OEIS is a testament to the profound impact of mathematical sequences on various fields. It highlights how the exploration of patterns and numbers can inspire new insights and innovations.

[Gil]

Just a possible, simple adaptation for the HP50G calculator:

Code:

\<< DUP 1 1 \-> n l p
  \<< 1 SWAP 1 \=/
    IF
    THEN 1 n 1 -
      START l 2 MOD NOT
        IF
        THEN l 2 /
        ELSE l p NEXTPRIME DUP 'p' STO +
        END DUP 'l' STO
      NEXT
    END n \->LIST
  \>>
\>>

As argument,
enter the number n
for the the number of elements to be calculated.

But in fact, no necessity to save the last element calculated, we can just play on its duplication. The above program becomes then:

Code:

\<< DUP 1 \-> n p
  \<< 1 SWAP 1 \=/
    IF
    THEN 1 n 1 -
      START DUP DUP 2 MOD NOT
        IF
        THEN 2 /
        ELSE p NEXTPRIME DUP 'p' STO +
        END
      NEXT
    END n \->LIST
  \>>
\>>

[John Keith]

(07-24-2024 09:08 AM)Gil Wrote:  But in fact, no necessity to save the last element calculated, we can just play on its duplication. The above program becomes then:

Code:

\<< DUP 1 \-> n p
  \<< 1 SWAP 1 \=/
    IF
    THEN 1 n 1 -
      START DUP DUP 2 MOD NOT
        IF
        THEN 2 /
        ELSE p NEXTPRIME DUP 'p' STO +
        END
      NEXT
    END n \->LIST
  \>>
\>>

Nice program, simple and fast. Here is a slightly modified version to illustrate a couple of memory saving tips. This is not meant as criticism of your program, just as information that may be of help to you and others.

Code:


\<< DUP 1 \-> n p
  \<< 1 SWAP 1
    IF \=/
    THEN 2 n
      START DUPDUP 2 MOD
        IF NOT
        THEN 2 /
        ELSE p NEXTPRIME DUP 'p' STO +
        END
      NEXT
    END n \->LIST
  \>>
\>>

My first tip turns out not to save any memory in this case but deserves mention. I have always heard that in IF..THEN..(ELSE)..END structures, having exactly one object between IF and ELSE minimizes memory use. It seems that your method of having nothing between IF and ELSE uses no more memory. I do think that my version is more readable but that may not be the case with other programs.

My next tip is to use 2 n START rather than 1 n 1 - START which does the same number of iterations but saves 5 bytes. Note that this may not work with FOR loops where the value of the looping variable would be different, but is always safe to use with START loops.

A final, simple tip is to use the extended stack operations available on the HP 49 and 50. In this case I replaced DUP DUP with DUPDUP, saving 2.5 bytes. Others include UNROT, PICK3 and NIP. These extended stack operations are faster than their equivalents, as well as saving memory.

[Gil]

Thanks for the tips, John, above all regarding the cleaner & more logical 2 n START... NEXT, that, besides memory saving, for 1 million loops/start seems — logically — to reduce the execution time by a factor of 1.001 in comparison to 1 n 1 - START... NEXT.

IF "condition" THEN... ELSE...
is, yes, clearer than "condition" IF THEN... ELSE: an habit from RPN logic.

DUP DUP is just safer — to include potential users of
the old HP48 that can't use the new, faster command DUPDUP.

[Gil]

By the way, it seems on HP50G/EMU48 that putting the condition/result before the IF THEN END is faster for the following program:

Code:

\<<
  \<< 1 100000
    START 1
      IF
      THEN
      END
    NEXT
  \>> TEVAL
\>>

than putting the condition/result just after the IF as here:

Code:

\<<
  \<< 1 100000
    START
      IF 1
      THEN
      END
    NEXT
  \>> TEVAL
\>>

[John Keith]

(07-25-2024 11:59 PM)Gil Wrote:  By the way, it seems on HP50G/EMU48 that putting the condition/result before the IF THEN END is faster ...

Testing on my physical 50g and EMU48, the difference is slight, less than 0.5%. However,

Code:

\<<
  \<< 1. 100000.
    START 1.
      IF
      THEN
      END
    NEXT
  \>> TEVAL
\>>

using approximate numbers, is 8 times as fast. Of course, not all programs can use approximate numbers but if the numbers are < 10^12 and exact integers are not required, the time savings can be significant. By the way, the decimal points following integers do not cause problems with the HP-48, they are simply ignored.

Finally, apologies to Valentin for derailing his thread on the OEIS (of which I am a huge fan) with discussions of his least favorite programming language.

[Gil]

My apologises too to Valentin.

Regards,
Gil

[C.Ret]

Bonjour à toutes et à tous.

I would like to intervene on behalf of Gil and John. Even though their discussion regarding the use this really esoteric language has made a small digression; their interventions allowed me to notice this very interesting subject that I had not seen last year.

Determining The Sisyphus sequence seems easy with a machine equipped with a JPC-ROM that, like other specific pocket computers, allows you to easily have the prime numbers in sequence. But what about a machine without this very opportunistic feature?

This is the question I asked myself and to try to explore the possible methods, I took my SHARP PC-1211 out of its hard case.

The first program proposed by Valentin led me to this adaptation:

 1: D=4.2424626 , W=10-Ǣ-7 , N=1 , P=1
 2: PRINT N : N=N/2 : IF N>INT N GOSUB 5 : N=2N+P
 3: GOTO 2 

 5: IF P<7 LET P=1+P+(P>2 : RETURN 
 6: P=P+INT D , D=10D-W*INT D , F=5
 7: F=F+2 , Q=P/F : IF Q<F RETURN 
 8: GOTO 6+(Q>INT Q

where Ǣ stand for the Exp key (10-exponentiation key)

As the SHARP PC-1211 is not fast, I tried to optimize things by determining the next prime number from a rudimentary generator of quasi-prime numbers. In reality, integers that are not multiples of 2, 3 or 5. Which facilitates primality tests based on divisibility by odd factors, thus starting only from 7.

My code seems to work and list the integers of the sequence quite quickly (for this SHARP).

But I have trouble finding an efficient method to adapt the codes given by Valentin later. In particular, how to find with a PC-1211 the results obtained by an HP-71B without having to wait a century?

Any ideas or suggestions are of course welcome