What comes next?

[Thomas Klemm]

This program for the HP-42S will guess the next number of an arbitrary sequence:

Code:

00 { 24-Byte Prgm }
01 1
02 RCL+ 00
03 STO 00
04 X<>Y
05 ENTER
06 GTO 01
07▸LBL 00
08 RCL- IND ST Z
09 STO+ ST Y
10▸LBL 01
11 STO IND ST Z
12 DSE ST Z
13 GTO 00
14 R↓
15 END

Examples

Moser's Circle Problem

A000127

CLRG

1 R/S
1

2 R/S
3

4 R/S
7

8 R/S
15

16 R/S
31

R/S
57

R/S
99

R/S
163

R/S
256

Prime Generating Polynomial

A005846

CLRG

41 R/S
41

43 R/S
45

47 R/S
53

R/S
61

R/S
71

R/S
83

Binomial Coefficient

A000332

CLRG

0 R/S
0

R/S
0

R/S
0

R/S
0

1 R/S
1

R/S
5

R/S
15

R/S
35

R/S
70


The program was inspired by the following video:

[chromos]

It reminded me of my college days, when we were taught that "guessing" the next member of an infinite series is not a problem from mathematics, but from psychology (just put yourself in the mood of the author of the problem and guess which of the popular series he probably had in mind ). Even if we knew the first million members, e.g. 1, 2, 3, ..., 999 998, 999 999, 1 000 000, the next one could be absolutely arbitrary, not necessarily 1 000 001.

[John Keith]

Nice short program! Seems similar to the inverse binomial transform.

[Thomas Klemm]

(10-08-2024 12:14 PM)John Keith Wrote:  Seems similar to the inverse binomial transform.

Instead of the red numbers of Post #7 the following red numbers are stored:

\(
\begin{matrix}
U &: & 0 & & 1 & & 4 & & 10 & & {\color{Red} {20}} & & 35 & & 56 & & \cdots \\
\Delta U &: & & 1 & & 3 & & 6 & & {\color{Red} {10}} & & 15 & & 21 & & \cdots & \\
\Delta^2 U &: & & & 2 & & 3 & & {\color{Red} {4}} & & 5 & & 6 & & \cdots & & \\
\Delta^3 U &: & & & & 1 & & {\color{Red} {1}} & & 1 & & 1 & & \cdots & & & \\
\Delta^4 U &: & & & & & {\color{Red} {0}} & & 0 & & 0 & & \cdots & & & & \\
\cdots &: & & & & & & & \cdots
\end{matrix}
\)

Example

CLRG

0 R/S
0

1 R/S
2

4 R/S
9

10 R/S
20

R/S
35


This leads to the following entries in the registers:

R00: 5
R01: 0
R02: 1
R03: 4
R04: 10
R05: 20
R06: 0
R07: 0
R08: 0
R09: 0

[Thomas Klemm]

(10-08-2024 05:16 AM)chromos Wrote:  It reminded me of my college days, when we were taught that "guessing" the next member of an infinite series is not a problem from mathematics, but from psychology (just put yourself in the mood of the author of the problem and guess which of the popular series he probably had in mind ).

Somewhat related: The Strong Law of Small Numbers

[Gil]

Thanks for the Mathologer video, Thomas.

Here is my version of the Gregory-Newton Interpolation formulae for the HP50G.

Argument:
A non-empty list,
the 1st number being always = f(0),
the 2nd being = f(1),
the 3rd being = f(2), etc.

Example :
{0 1 4 9} NextN

and the result will be:
{ 0 1 4 9 }
f(x): 'x^2'
f(4): 16

Don't forget :
to put a zero for the 1st element in the above example list if you expect to get x². Remember: first element in the list is always f(0), ie in this example 0²=0.

If you put {1 4 9} instead,
then the program supposes that f(0) = 1, f(1) = 4 and f(2) =9, whose solution is obviously not x².

Result for {1 4 9} NextN:
{ 1 4 9 }
f(x): 'x^2+2*x+1'
f(3): 16

Example of Mathologer
{1 2 4 8 16 31} NextN

and the result will be:
{ 1 2 4 8 16 31 }
:f(x): '1/24*x^4-1/12*x^3+11/24*x^2+7/12*x+1'
f(6): 57

Attention again: suppose that the sequence were
{0 1 2 4 8 16 31}

then {0 1 2 4 8 16 31} NextN:
will give the following output:
{ 0 1 2 4 8 16 31 }
f(x): '-1/720*x^6+7/240*x^5-29/144*x^4+37/48*x^3-467/360*x^2+17/10*x'
f(7): 56

From last example, calculate for instance f(2.5):
a) DROP
2.5 'x' STO
EVAL
and you get 2.8291015625
b) or
-105 CF
25 ENTER 10 /
'x' STO
EVAL
and you get '2897/1024'.

Code below slightly changed in comparison to previous versions

Code:


\<< "
Arg
{
 f(0) f(1) f(2)f(n)
}

Ex: 
{ 
0 1 4 9 16 24 
}


Result: f(n+1)
Here: 
30

" DROP DUP DUP DUP SIZE R\->I \-> L s
  \<< s 1 ==
    IF
    THEN DROP L 1 GET DUP FP 0 == { R\->I } IFT "f(x)" \->TAG
    ELSE 1 GET DUP FP 0 == { R\->I } IFT 1 \->LIST 2 s
      FOR i { } 2 s i 2 - -
        FOR j L j GET L j 1 - GET - +
        NEXT DUP 'L' STO 1 GET DUP FP 0 == { R\->I } IFT +
      NEXT -114 CF -105 CF 'x' PURGE 'L' STO 0 0 s 1 -
      FOR i L i 1 + GET DUP 0 ==
        IF
        THEN
        ELSE i R\->I DUP DUP DUP
          CASE 0 ==
            THEN 3 DROPN 1
            END 1 ==
            THEN DROP2 x
            END x 1 ROT 1 -
            FOR i x i R\->I - *
            NEXT SWAP ! /
          END EVAL *
        END +
      NEXT 1 s
      FOR i L i GET TYPE 0 ==
        IF
        THEN s 'i' STO -105 SF
        END
      NEXT EXPAN \->STR 11 14
      FOR i ".E-" i R\->I + "*0" SREPL DROP
      NEXT ".+" "+" SREPL DROP ".-" "-" SREPL DROP ".*" "*" SREPL DROP ".'" "'" SREPL DROP "." "." SREPL 0 ==
      IF
      THEN -105 CF
      END OBJ\-> EVAL EXPAN \->STR "+-" "-" SREPL DROP OBJ\-> "f(x)" \->TAG DUP EVAL s 'x' STO EVAL DUP FP 0 == { R\->I } IFT "f(" s R\->I + ")" + \->TAG 'x' PURGE
    END
  \>>
\>>

[Thomas Klemm]

I'm not sure what your program does, but if you want to transform the polynomial based on binomial coefficients \(\binom{n}{k}\) to the common form based on \(n^k\) the linear transform A131689 can be used.

Example

\(f(n)=\binom{n}{1}+2\binom{n}{2}+\binom{n}{3}\)

[ 0 1 2 1 ]

[[ 1 0 0 0 ]
 [ 0 1 1 1 ]
 [ 0 0 2 6 ]
 [ 0 0 0 6 ]]

/

[ 0 .333333333333 .5 .166666666667 ]

Thus:

\(f(n)=\frac{n}{3}+\frac{n^2}{2}+\frac{n^3}{6}\)

[Gil]

Interesting, but I do not understand where the Matrix in your example comes from.

As mentioned in my 1st post, now with more details and examples, when having
{#1 #2 #3... #n}
with #1=f(0), #2=f(1), #3=f(2),... #n=f(n-1),

then running my program Next_NumberV8 gives 3 outputs:

a) {#1 #2 #3... #n} <again, in stack level 3>;
b) f(x) <the function going through the above points {#1 #2 #3... #n}, in stack level 2>;
c) f(n) <the corresponding number #n+1 following #n, in stack level 1>.

[Thomas Klemm]

(10-10-2024 09:17 PM)Gil Wrote:  Interesting, but I do not understand where the Matrix in your example comes from.

Just write the sequence A131689 downwards as an infinite triangular matrix:

\(
\begin{bmatrix}
1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \cdots \\
& 1 & 1 & 1 & 1 & 1 & 1 & 1 &\cdots \\
& & 2 & 6 & 14 & 30 & 62 & 126 & \cdots \\
& & & 6 & 36 & 150 & 540 & 1806 & \cdots \\
& & & & 24 & 240 & 1560 & 8400 & \cdots \\
& & & & & 120 & 1800 & 16800 & \cdots \\
& & & & & & 720 & 15120 & \cdots \\
& & & & & & & 5040 & \cdots \\
& & & & & & & & \cdots \\
\end{bmatrix}
\)

There is a relative simple formula: \(T(n,k) = k \cdot \left(T(n-1,k-1) + T(n-1,k)\right)\) with \(T(0,0)=1\).

Examples

\(
\begin{align}
30 &= 2 \cdot (1 + 14) \\
1806 & = 3 \cdot (62 + 540) \\
8400 & = 4 \cdot (540 + 1560) \\
\end{align}
\)

[C.Ret]

Bonjour à tous,

A big thank you to Thomas Klemm for this small but very clever program that works wonders to determine arbitrary sequences. Whatever they are!

I think it's great to have had the idea of ​​doing the sum while memorizing the differences. At the end of the program everything is already ready for the next determination. Excellent idea.

I use and abuse this little program that I transcribed for my HP-15C in order to be able to use it in the laboratory.

I had a much longer and much slower program that memorized the elements of the sequence and recalculated the descending and ascending differences at each determination.

And also thank you for the tip 1 RCL+0 STO 0 which is much shorter than the method with an ISG 0 that I used.

It is therefore with great pleasure that I share below a new version for HP-15C directly adapted from the algorithm proposed by Thomas.

Code:

001 {        1 }    1
002 { 45 40  0 }   RCL + 0
003 {    44  0 }   STO 0
004 {    44 25 }   STO I
005 {    43 35 } g CLx
006 {       34 }   X⇄Y
007 {       20 }    ×
008 { 42 21  0 } f LBL 0
009 {    43 36 } g LSTx
010 { 45 30 24 }   RCL - (i)
011 {    44 24 }   STO (i)
012 {       40 }    +
013 { 42  5 25 } f DSE I
014 {    22  0 }   GTO 0

This code is used exactly like Thomas's. So don't forget to do f CLEAR REG before introducing the elements of a new sequence.

[Thomas Klemm]

(10-13-2024 12:35 PM)C.Ret Wrote:  I use and abuse this little program that I transcribed for my HP-15C in order to be able to use it in the laboratory.

Now you made me curious: in what kind of laboratory are you using it?

BTW: Nice adaption of the program for the HP-15C.

[Geoff Quickfall]

Any perfect number algorithms out there?

[John Keith]

(10-13-2024 06:33 PM)Geoff Quickfall Wrote:  Any perfect number algorithms out there?

Perfect numbers are easy (sort of). 2^(p-1)*(2^p - 1) is a perfect number if (2^p - 1) is a Mersenne prime. For example, 2^3-1 = 7 is prime and since 2^(3-1) = 4, 28 = 4*7 is a perfect number.
Of course now you need the ISPRIME? function, which limits you to the HP 49 or later unless you want to roll your own ISPRIME?.

[Thomas Klemm]

The sequence of perfect numbers grows fast: A000396.
Therefore, only the first 6 or 7 can be displayed.

Here's a program for the HP-42S that calculates \(2^{p-1} \cdot \left( 2^p -1 \right)\):

Code:

00 { 13-Byte Prgm }
01 2
02 ENTER
03 1/X
04 X<> ST Z
05 Y↑X
06 ×
07 LASTX
08 1
09 -
10 ×
11 END

Just use it with the first primes but skip 11 since \(2^{11}-1=2047=23\times89\) is not a prime.

Examples

2 R/S
6

3 R/S
28

5 R/S
496

7 R/S
8'128

13 R/S
33'550'336

17 R/S
8'589'869'056

19 R/S
137'438'691'328

[Thomas Klemm]

(10-13-2024 08:49 PM)John Keith Wrote:  Of course now you need the ISPRIME? function, which limits you to the HP 49 or later unless you want to roll your own ISPRIME?.

In case of Mersenne prime numbers we can use the Lucas–Lehmer primality test.

Here's an implementation for the HP-42S:

Code:

00 { 32-Byte Prgm }
01▸LBL "LLT"
02 RCL ST X
03 2
04 STO- ST Z
05 X<>Y
06 Y↑X
07 1
08 -
09 4
10▸LBL 00
11 X↑2
12 2
13 -
14 RCL ST Y
15 MOD
16 DSE ST Z
17 GTO 00
18 END

Examples

3 XEQ "LLT"
0

5 XEQ "LLT"
0

7 XEQ "LLT"
0

11 XEQ "LLT"
1'736

13 XEQ "LLT"
0

17 XEQ "LLT"
0

19 XEQ "LLT"
0

23 XEQ "LLT"
6'107'895

29 XEQ "LLT"
458'738'443

31 XEQ "LLT"
0

[Gilles]

newRPL with LstX library. Enter the first numbers as a list and you get a list with a new guess number. etc.

Code:

«
 DUPDUP RHEAD SWAP
 2 OVER SIZE START
   ΔLIST DUP RHEAD ROT + SWAP 
 NEXT
 + ADD
»

Ex :
{ 1 2 4 8 16 31 } -> { 1 2 4 8 16 31 57 } ->{ 1 2 4 8 16 31 57 99 } ...
A000127 OK
A005846 OK
A000332 OK
...
{ 0 1 4 10 20 } -> { 0 1 4 10 20 35 }

[Thomas Klemm]

Or then use the binomial transform T and add as many 0 as you like to predict the next elements:

{ 1 2 4 8 16 }
T

{ 1 -1 1 -1 1 }

{ 0 0 0 0 0 }
+

{ 1 -1 1 -1 1 0 0 0 0 0 }

T

{ 1 2 4 8 16 31 57 99 163 256 }

---

As a simple program to add what comes next:

Code:

\<< T 0 + T
\>>

[Gil]

I used the method explained in the Mathologer video that you included in your 1st post.

It works fine, but your program is incredibly much shorter.

Could you explain the algorithm you used (as I don't the language of your calculator)?

Thanks for your help.

Regards,
Gil

[Thomas Klemm]

(10-15-2024 10:09 AM)Gil Wrote:  Could you explain the algorithm you used?

From Wikipedia: Binomial transform:

Quote:The binomial transform of a sequence is just the nth forward differences of the sequence, with odd differences carrying a negative sign, namely:

\(
\begin{aligned}
s_{0}&=a_{0}\\
s_{1}&=-(\Delta a)_{0}=-a_{1}+a_{0}\\
s_{2}&=(\Delta ^{2}a)_{0}=-(-a_{2}+a_{1})+(-a_{1}+a_{0})=a_{2}-2a_{1}+a_{0}\\
&\;\;\vdots \\
s_{n}&=(-1)^{n}(\Delta ^{n}a)_{0}
\end{aligned}
\)

where Δ is the forward difference operator.

The only difference from the algorithm in the Mathologer video is the alternating sign due to the factor \((-1)^{n}\).
This is crucial because it makes the transformation an involution, i.e. \(TT=1\).
This allows the same program T to be used in both directions of the transformation.

As for the implementation, I hope the following comments are helpful (right is at the top of the stack):

Code:

\<< { } SWAP                @ S={} A 
  WHILE                     @ S A
    DUP SIZE 1 >            @ S A SIZE(A)>1
  REPEAT                    @ S A
    SWAP OVER HEAD +        @ A S=S+A[0]
    SWAP \GDLIST NEG        @ S A=-ΔA
  END +                     @ S=S+A
\>>                         @ S

In other words, we successively calculate -ΔA and append its first element to the sequence S.
I also recommend stepping through the code in the debugger using a simple example for A.

[Gil]

Thanks for sharing this very compact code to get the coefficients list {c0 x1 c2...}.

Nice ideas yours : the while routine, the use of GDLIST (that I created) and the HEAD command (instead of 1 GET).

I don't use the NEG command in your last code,
building simply the polynomial P(x)=
c0×COMB(x 0) + c1×COMB(x 1) + c2×COMB(x 2) +... , as I understood from the Mathologer video.

My next question is if I may ask: what is the code algorithm you used in your first post to get the binomial coefficient ?
I know of course how to calculate them, but I don't understand where you calculate them in post 1, as the given code there seems to be compact. Did you use a special algorithm to have it coded in such a compact manner?

Last question, what will you get in your calculator or l'émulation for the next number of the following list?
{ 17 72 97 8 32 15 63 97 57 60 83 48 26 12 62 3 49 55 77 97 98 0 89 57 34 92 29 75 13 40 3 2 3 83 69 1 48 87 27 54 92 3 67 28 97 56 63 }

My outputs are the following on my HP50G (EMU48 on the phone):

f(x): '79452528137/73368295464161185998004072384003398572822023372800000000*x^46-62405430595537/34177777390137198446275188998759347161252495360000000000*x^45+6867032884106951/5316543149576897536087251622029231780639277056000000000*x^44-34545108562396661/63042013631346611099848833462797214790189056000000000*x^43+106691370557962873/661179349530767011079125932350358385852416000000000*x^42-923169570689950723/25870912096109496367220670975808100106240000000000*x^41+5116625404961330065763/822442605513880916220376159752884821426176000000000*x^40-311468829373514971297883/352475402363091821237304068465522066325504000000000*x^39+27594326046634577365304​3/2636033992031669603270436409464374427648000000000*x^38-87608265026221671791006869/8324317869573693484011904450940129771520000000000*x^37+3002951632494509448186313​9/32997296059571397594281423048771685580800000000*x^36-4836221201668622333612196203/70708491556224423416317335104510754816000000000*x^35+289882152645079814443270437​89639/6434472731616422530884877494510478688256000000000*x^34-839345543093035778449860512219/3205802337471230543171792686672773120000000000*x^33+5160708636241167083049878373​3509/3823216120984208277412286092994936832000000000*x^32-790243016479486246527834493165571/1274405373661402759137428697664978944000000000*x^31+1068115248911007357425117522​68222691/4193204777853647788129604101994446848000000000*x^30-8431800252037961299841608135968849877/8985438809686388117420580218559528960000000000*x^29+7682535668432642325561622458​09725666423/24725448862516336957591803497967255552000000000*x^28-1963166711585595025861899594374007832651/2119324188215686024936440299825764761600000000*x^27+2042562341192591298725383283​799475884961/82045230011366427310259378122968268800000000*x^26-63148633232678459156494019054764776365418851/104607668264492194820580707106784542720000000000*x^25+21191971223495843964537533​74358695387711789651/160398424672221365391557084230402965504000000000*x^24-17858444764609635577090499051547310644173/68423124970642479286895385772032000000000*x^23+575125831754572115635515365719990​27424387232007/12362724431258168478795901748983627776000000000*x^22-43054411737196706236305010020807767010248926001/575903312015132071993597907561349120000000000*x^21+11146376398461744022298118158​88516792729194246171/1030227035938180706566325145748635648000000000*x^20-572545195255310504916372668049778639167342175967/40666856681770291048670729437446144000000000*x^19+161570098669313148872847271431​80289175362272256437/98278236980944870034287596140494848000000000*x^18-1252209966295262840527254340195582349765162487367/728706650872057859868667297136640000000000*x^17+47399261150573426840777142764422​79198720541863850271/295701871857372358956106502505077145600000000*x^16-168492124329437939303112790767488353520865947410307863/1267293736531595824097599296450330624000000000*x^15+2097376832799371152402216895​03251905819425068750689473/214826146221167953087769681307107328000000000*x^14-924867262101770565596196284887085917299220361455178393/146472372423523604378024782709391360000000000*x^13+21545644845605013092131446491​07615993811498231839781/60276696470585845423055466135552000000000*x^12-111065175127946654779403925649034914832194926918919/631782144683935491623640367104000000000*x^11+30151778145877484822280521461025808​81633138423086659489/4050505360622676774421352242741248000000000*x^10-14144171311503183171149392051950453414030477549299293597/5269790137544809068762473581117440000000000*x^9+12145489715007185763814880131700​4578091817898013556907/14943270644658312752856551242752000000000*x^8-728465154287934961017826027803007940038213284428022953/35863849547179950606855722982604800000000*x^7+8175783823440875967323115918963477​00161404311307203053/19968637217010071263817198845870080000000*x^6-4359514147282944414325206668958078534305004412820119/67563058253041594501637138952192000000*x^5+9251853218738204433802104415183867205​3095671664831/1216909580968584274864143353702400000*x^4-49357091850914086397470724236043481620906093159/790201026602976801859833346560000*x^3+132720057639127275173051837261789837516351​/4202885913130495995387984000*x^2-6869689147511982724488320657561/941958815880242160*x+17'

f(47): -313323742215540

Thanks in advance for your collaboration, Thomas.

Regards,
Gil