(33s) OEIS A2516: EARLIEST SEQUENCE WITH A(A(N)) = 2N

[Gerald H]

Takes a real integer from the stack & returns A2516(N) to stack.

https://oeis.org/A002516

Maximum number of recursions is 7, ie any number with a power of 2 greater than 2^7 as a factor will produce the
“XEQ OVERFLOW”
error.

Earliest sequence with
A(A(N)) = 2N.

Code:

1.    LBL C
2.    STO N
3.    CLx
4.    RMDR(N:2)
5.    x=0?
6.    GTO D
7.    CLx
8.    RMDR(N:4)-1
9.    x≠0?
10.    GTO E
11.    CLx
12.    N+2
13.    RTN
1.    LBL E
2.    CLx
3.    (N-2)*2
4.    RTN
1.    LBL D
2.    CLx
3.    N/2
4.    XEQ C
5.    STO N
6.    RCL+N
7.    RTN

C: LN = 62
E: LN = 19
D: LN = 24

[C.Ret]

Hello,

First of all, thank you for sharing this subject, curious by nature, I am glad to discover a type of sequence that I did not know before.

Unfortunately, I am not an HP-33s user, but if I interpret the code correctly, the calculation of the n-th element of the sequence is done by recurrence. Which explains the limitation to 2^7 due to the limited number of nested calls.

I propose below a code for HP-41C, which I hope will be easily adapted for HP-33s, based on an iterative calculation. A loop make it easy to calculate all the elements of the sequence without any limitation as long as their values ​​or index can be represented without rounding.

The principle of this calculation is based on the observation of the structure of the formulas making it possible to calculate each \( u_n \) in the form \( u_n=p\times \left ( a.k+b \right ) \).
The calculation loop runs as long as the indices to be calculated are even. This loop doubles the coefficient \( p \) and halves the indice \( k \).
If the index is odd, the coefficients \( a \) and \( b \) are calculated according to the value of \( ( k \mod 4 ) \).

\( u_n = p\times\left ( a.k+b \right )=\left\{\begin{matrix}
u_{4k }&=&2\times& &u_{2k}&&(m=0)&\textit{even}&:&p&\leftarrow& 2.p&k&\leftarrow& k/2 &\textit{loop} \\
u_{4k+1}&=&& 4.k& + &3&(m=1)&\textit{odd}&:&a&\leftarrow& 4 &b&\leftarrow& 3 &\textit{terminate}\\
u_{4k+2}&=&2\times& &u_{2k+1} & &(m=2)&\textit{even}&:&p&\leftarrow& 2.p&k&\leftarrow& k/2 &\textit{loop} \\
u_{4k+3}&=& & 8.k& + &2&(m=3)&\textit{odd}&:&a&\leftarrow& 8 &b&\leftarrow& 2 &\textit{terminate}\\
\end{matrix}\right. \)

The code only use register R00.

01 LBL"OEIS002516
02   STO 00  1  2  GTO 02       // Initialization
06   LBL 01                  // *** Even loop  
07     RDN  ST* 00  ST* Y       //     p←2p    i←i/2
10     LBL 02                //   loop entrance
11     RCL 00  X=0?  RTN        //     u(0) = 0 STOP  
14     4  ST/ 00  MOD         //     k← i/4    M← i MOD 4
17     X≠Y?  X=0?  GTO 01       //     loop while M in { 0 ; 2 } 
20   X<>Y                   // *** Odd termination  M in { 1 ; 3 }
21   ST* Y  +                // compute a ← 2.M+2 
24   ENTER^  CHS  4  ST/ Y  +    // compute b ← 4-a/4
28   X<>Y  RCL 00  INT  *  +  *  // compute u(n) ← p*(a.k+b)
34 END                     // u(n) so that u(u(n))=2n (OEIS002516)


Usage:
Enter indice n and press [R/S].
The value of u(n) is calculated and left in the display register X:. This makes it easy to verify that u(u(n))=2n; simply press [R/S] again.

Example:
2147483648 [R/S] 6442450944. [R/S] 4294967296.
835584 [R/S] 1605632.000 [R/S] 1671168.000
75861 [R/S] 75863.0000 [R/S] 151722.0000
1024 [R/S] 3072.0000 [R/S] 2048.0000
30 [R/S] 52.0000 [R/S] 60.0000
2 [R/S] 6.0000 [R/S] 4.0000
1 [R/S] 3.0000 [R/S] 2.0000
0 [R/S] 0.0000


Finally, just to satisfy my curiosity, what is a "real integer"?
I know the integers, the real, the rational and the irrational,... as well as the natural integers, the relative integers, positive or negative integers...
Is there any unreal integers ?

Numerous editions to correct broken english and typos.

[Thomas Klemm]

For those who want to run the program in Free42:

Code:

00 { 52-Byte Prgm }
01▸LBL "A2516"
02 STO 00
03 1
04 2
05 GTO 02
06▸LBL 01
07 R↓
08 STO× 00
09 STO× ST Y
10▸LBL 02
11 RCL 00
12 X=0?
13 RTN
14 4
15 STO÷ 00
16 MOD
17 X≠Y?
18 X=0?
19 GTO 01
20 X<>Y
21 STO× ST Y
22 +
23 ENTER
24 +/-
25 4
26 STO÷ ST Y
27 +
28 X<>Y
29 RCL 00
30 IP
31 ×
32 +
33 ×
34 END

Thanks to Gerald H for bringing up such sequences again and again and to C.Ret for posting these excellent programs.
Both are greatly appreciated.

[Gerald H]

Below a direct translation of Thomas' programme to 33s-ese.

Very nicely done, Thomas, congrats.

Not sure why I started using "real integer", perhaps because people started using the phrase "real people" which irritated me enormously. Nice to find members who question my use of English - the expression may have some meaning on 49G, as 77. and 77 are different entities on that calculator.
However, definitely not an innovative type of integer.

Very nice programme, Ret, congrats. My version below also actually gives correct answers until numbers too big for exact representation!

Code:

1.    LBL C
2.    STO N
3.    1
4.    2
5.    GTO N
1.    LBL M
2.    R↓
3.    STO* N
4.    *
5.    LASTx
1.    LBL N
2.    RCL N
3.    x=0?
4.    RTN
5.    4
6.    STO/ N
7.    RMDR
8.    x≠y?
9.    x=0?
10.    GTO M
11.    x<>y
12.    *
13.    LASTx
14.    +
15.    ENTER
16.    +/-
17.    4
18.    /
19.    LASTx
20.    +
21.    x<>y
22.    RCL N
23.    IP
24.    *
25.    +
26.    *
27.    RTN

L: LN = 39
M: LN = 15
N: LN = 105

[John Keith]

(07-31-2022 03:32 PM)Gerald H Wrote:  ... the expression may have some meaning on 49G, as 77. and 77 are different entities on that calculator.

Speaking of which, here is a program for any RPL calculator which returns a list of A002516. Given an integer n on the stack, the length of the list will be n + (1 if n is odd else 2).

Approximate mode should be used on the 49 and 50 because MOD is slow for exact integers.

Code:


\<< \-> n
  \<< 0 3 2 n                            @ Start with a(2)
    FOR k k 2 / PICK 2 *                 @ a(2n) = 2*a(n)
      k DUP 3 + SWAP DUP
      5 - SWAP 2 / 2 MOD * + 2           @ Compute a(2n+1)
    STEP n DUP 2 + SWAP 2 MOD - \->LIST  @ Adjust for n MOD 2
  \>>
\>>

The program is based on the third formula on the OEIS page. A list of 1000 terms takes about 4.5 seconds on my 50g in approximate mode.

[Thomas Klemm]

(07-31-2022 03:32 PM)Gerald H Wrote:  Very nicely done, Thomas, congrats.

All the credits should go to C.Ret.
And then of course to Thomas Okken, who made the marvellous Free42 available to us.
I only reformatted the program so it can be parsed.
And then I thought I might as well post it here.

[Gerald H]

Yes, Thomas, agree with all you write.

I presume you do know that your programme works on a real HP 42S, not only on Free42.

[Thomas Klemm]

(08-01-2022 10:50 AM)Gerald H Wrote:  I presume you do know that your programme works on a real HP 42S, not only on Free42.

Yes of course.
But I assumed that with a real HP-42S you could just as well use the original listing as a source.
The translation from HP-41C code to HP-42S is trivial.
But you can't just copy and paste it into Free42.
So my listing is mainly useful for those who use it.

[C.Ret]

(07-31-2022 03:32 PM)Gerald H Wrote:  Nice to find members who question my use of English - the expression may have some meaning on 49G, as 77. and 77 are different entities on that calculator.
However, definitely not an innovative type of integer.

Thank you very much, I am not a native English speaker and I was afraid of missing something here such as a specific math's object or anything else. Your explanations are clear to me and I think that calling an integer store in the machine as any other real value a real integer makes sense.

I recently read the english version of the HP-35S User manual and I notice they speak about real functions as well. Perhaps is your habit to use real integer link to this ?

(07-31-2022 06:25 PM)John Keith Wrote:  Speaking of which, here is a program for any RPL calculator which returns a list of A002516. Given an integer n on the stack, the length of the list will be n + (1 if n is odd else 2).
[...]
The program is based on the third formula on the OEIS page. A list of 1000 terms takes about 4.5 seconds on my 50g in approximate mode.

I adapt your RPL code to my HP-28S, but I was a bit disappointed by the variable length of the list depending of the argument. As in your RPL code, I used the third formula from the OEIS page. But an HP-28S will not compete again a fast and powerful HP-50g, to compute the list of the first 1000 terms takes about a short minute!

Nevertheless, it doesn't take much more time to also plot or print all these values. Please note the use of flag #1 to alternate a(n) definition on odd k value.


A2516:
« → m « { m } 0 CON 1 SF
   1 m FOR k
      k
         DUP 2 MOD
         « k 2 + IF 1 FC?C THEN k 6 - + 1 SF END »
         « OVER k 2 / GET 2 * »
        IFTE
      PUT
   NEXT » »

Usage: Enter n and execute A2516 to get an array of this length n containing the first n values of the suite from a(1) up to a(n).

7 A2516 returns [ 3 6 2 12 7 4 10 ]

On the HP-28S, this code can with great ease be modified and/or enhanced to allow plotting and printing (here on an HP-82240A IR-printer):

   


I always wonder why HP-50g's users so rarely exploit the amazing plotting capabilities of this great RPL system?

Edited: A typo was in the code, I type an '*' where a '+' is needed in the code.

[John Keith]

I copied and pasted your code to EMU48 and with an input of 7 I got

[ 2 4 3 8 10 6 15 ].

Could there be a typo in your code?

As far as speed is concerned, accumulating numbers on the stack and storing them in a list or array at the end is significantly faster than using GET and PUT on an existing array. Also stack operations such as PICK and ROLL are very fast, sometimes even faster than list processing commands on later calculators.

[C.Ret]

(08-21-2022 01:05 PM)John Keith Wrote:  Could there be a typo in your code?

Sure, I type an '*' where a '+' is needed in the code !!
Sorry for that, corrected code is :

A2516:
« → m « { m } 0 CON 1 SF                       // Build array [ 0 0 ... 0 ] of size { m } and initialze flag n°1
   1 m FOR k                            // For each element of the array
      k                              // k-th element
         DUP 2 MOD                      // parity of k for IF THEN ELSE structure
         « k 2 + IF 1 FC?C THEN k 6 - + 1 SF END »   // k=2n+1: alternate a(k)=a(2n+1)=2n+3=k+2 when n is even   and  a(k)=a(2n+1)=4n-2=2k-4 when n is odd
         « OVER k 2 / GET 2 * »               // k=2n: a(k)=a(2n)=2*a(n) get a(n) from the already built part of the array
        IFTE
      PUT                            // Put a(k) in the array
   NEXT » »

(08-21-2022 01:05 PM)John Keith Wrote:  As far as speed is concerned, accumulating numbers on the stack and storing them in a list or array at the end is significantly faster than using GET and PUT on an existing array. Also stack operations such as PICK and ROLL are very fast, sometimes even faster than list processing commands on later calculators.

I follow your advice and get a faster and shorter code using the exact same algorithm:


A2516:
« → s
   « 1 SF
    1 s FOR k
      IF k 2 MOD
       THEN
         IF 1 FC?C
          THEN 2 k * 4 - 1 SF
          ELSE k 2 + END
       ELSE
         k 2 / PICK 2 * END
    NEXT
    s →LIST » »

The list of the first 544 elements is now obtained in less than 30 seconds: a great improvement.

Thanks

[Gerald H]

Inspired by Ret's etc programme here's my own attempt at derecursivising.

NB Preserves the stack.

Code:

1.    LBL S
2.    x=0?
3.    RTN
4.    STO N
5.    RCL/N
6.    STO C
1.    LBL T
2.    CLx
3.    RMDR(N:2)
4.    x≠0?
5.    GTO U
6.    CLx
7.    2
8.    STO* C
9.    CLx
10.    N/2
11.    STO N
12.    GTO T
1.    LBL U
2.    CLx
3.    C*(N+2+(N-6)*RMDR(IDIV(N:2):2))
4.    RTN

S: LN = 18
T: LN = 60
U: LN = 43

[Gerald H]

A slightly shorter version:

Code:

1.    LBL S
2.    x=0?
3.    RTN
4.    STO N
5.    SGN
6.    STO C
1.    LBL T
2.    CLx
3.    RMDR(N:2)
4.    x≠0?
5.    GTO U
6.    CLx
7.    2
8.    STO* C
9.    STO/ N
10.    GTO T
1.    LBL U
2.    CLx
3.    C*(N+2+(N-6)*RMDR(IDIV(N:2):2))
4.    RTN

S: LN = 18
T: LN = 51
U: LN = 43

[Gerald H]

An improved version of the previous programme.

I fear the section of line T10

SGN(IP(RMDR(N:4)- 0.1))

is awkward - Could someone with superior arithmetic abilities please suggest a more elegant expression?

Code:

1.    LBL S
2.    x=0?
3.    RTN
4.    STO N
5.    SGN
6.    STO C
1.    LBL T
2.    CLx
3.    RMDR(N:-2)+2
4.    STO* C
5.    STO/ N
6.    ACOSH
7.    x≠0?
8.    GTO T
9.    CLx
10.    C*(N+2+(N-6)*SGN(IP(RMDR(N:4)- 0.1))
11.    RTN

S: LN = 18
T: LN = 80