OEIS A244953: Ridiculous Formula

[Gerald H]

The entry in OEIS

https://oeis.org/A244953

contains the formula

a(n) = 1 + n + ( 2*(1 + n) - (1 + (-1)^n)*(1 + 2*i^(n*(n+1))) )/4, where i = sqrt(-1)

which I have great difficulty to believe is anywhere near the simplest solution to the problem.

Can someone please suggest a simpler calculation?

[Thomas Klemm]

How about:

\(
\begin{align}
a(n) = \frac{3n + 4 - (n \mod 4 - 2)^2}{2}
\end{align}
\)

Code:

00 { 20-Byte Prgm }
01 3
02 RCL× ST Y
03 4
04 +
05 X<>Y
06 4
07 MOD
08 2
09 -
10 X↑2
11 -
12 2
13 ÷
14 END

[Gerald H]

Very nice, Thomas, thank you.

LN for formula on 33s reduced from 45 to 25, & the two formulae give the same answers!

[C.Ret]

(08-20-2022 11:34 AM)Thomas Klemm Wrote:  How about:\( \begin{align} a(n) = \frac{3n + 4 - (n \mod 4 - 2)^2}{2} \end{align} \)

Very nice formula that will be efficient for large \( n \), but actually lead to a longer program than the simple original definition : \( a(n)=\begin{align}\sum_{k=0}^{n}mod(-k,4)\end{align} \)

A001 LBL A
A002 INPUT N
A003 0
A004 RCL ST Y  CHS  4  Rmdr  +  DSE ST Y  GTO A004
A011 RTN

Of course, this shorter version of the code need more time to compute a(n) for large n.

PS: Step A002 is optional, it is just the equivalent of a "N=?" PROMPT operation. Just to remember the user to enter a real integer.

[Thomas Klemm]

We can use the following program to generate the whole sequence:

Code:

00 { 6-Byte Prgm }
01 STO- ST L
02 R↑
03 RCL- ST L
04 END

Initialisation

0 +
3 ENTER
5 ENTER
6 ENTER

Run

Just hit repeatedly:

R/S

[John Keith]

(08-20-2022 10:35 PM)Thomas Klemm Wrote:  We can use the following program to generate the whole sequence:

Code:

00 { 6-Byte Prgm }
01 STO- ST L
02 R↑
03 RCL- ST L
04 END

That's amazing, only 4 instructions!

Your formula in post #2 comes out to 14 objects in RPL, which is 3 shorter than any formula that I could come up with. You may want to add it to the OEIS page since it is shorter and simpler than the formulas there.

[Thomas Klemm]

(08-21-2022 02:21 PM)John Keith Wrote:  That's amazing, only 4 instructions!

It's just this recursive formula: \(a(n) = a(n-1) + a(n-4) - a(n-5)\)
Or then rather: \(a(n) = a(n-4) - \left( a(n-5) - a(n-1) \right) \)

Quote:Your formula in post #2 comes out to 14 objects in RPL …

Is this what you refer to?

Code:

\<< 3 OVER * 4 + SWAP
4 MOD 2 - SQ - 2 /
\>>

Quote:You may want to add it to the OEIS page since it is shorter and simpler than the formulas there.

Meanwhile it has been updated.

[John Keith]

(08-22-2022 03:06 PM)Thomas Klemm Wrote:  Is this what you refer to?

Code:

\<< 3 OVER * 4 + SWAP
4 MOD 2 - SQ - 2 /
\>>

Meanwhile it has been updated.

Yes, exactly that. Thanks for updating the formula.

[DavidM]

I know the goal here was to find the formula for a(n) instead of computing the entire sequence, but I thought I'd include this as a possible solution to creating the list a(0)..a(n) when given n as the argument.

Target platform is 50g with ListExt installed:

Code:

%%HP: T(3)A(R)F(.);
\<<
  3 0 LSEQR                @ creates { 3 2 1 0 }
  OVER 4 / CEIL LMRPT      @ repeat the above enough for at least n elements
  SWAP LTAKE               @ trim the increments list to n elements
  0 ::+ LSCAN              @ start with 0 and build based on increments provided
\>>
@ Size: 58.0 bytes
@ CRC: 3319h

a(5000) ~10.8 seconds in exact mode, ~7.6 seconds in approximate mode.

[Thomas Klemm]

In the same spirit we can use the following in Clojure:

Code:

(->> [0 3 2 1]
     cycle
     (reductions +)
     (take 50))

This leads to:

Code:

(0 3 5 6 6 9 11 12 12 15 17 18 18 21 23 24 24 27 29 30 30 33 35 36 36 39 41 42 42 45 47 48 48 51 53 54 54 57 59 60 60 63 65 66 66 69 71 72 72 75)

It takes 0.032141 msecs for 5000 on a Macbook Pro.

[John Keith]

(08-24-2022 08:23 PM)DavidM Wrote:  

Code:

%%HP: T(3)A(R)F(.);
\<<
  3 0 LSEQR                @ creates { 3 2 1 0 }
  OVER 4 / CEIL LMRPT      @ repeat the above enough for at least n elements
  SWAP LTAKE               @ trim the increments list to n elements
  0 ::+ LSCAN              @ start with 0 and build based on increments provided
\>>
@ Size: 58.0 bytes
@ CRC: 3319h

An elegantly simple program and a great demonstration of ListExt. One would expect the Clojure version to be similar since Clojure and ListExt are strongly influenced by Haskell, where it would be something like
scanl (+) 0 cycle [3, 2, 1, 0].

Thanks also because I never realized that FLOOR and CEIL work with rational numbers. IP and FP do also, returning approximate numbers, whereas FLOOR and CEIL return exact integers.

[DavidM]

(08-25-2022 03:32 PM)John Keith Wrote:  ...I never realized that FLOOR and CEIL work with rational numbers.

I will confess that I didn't, either, until I tested this for exact/approximate mode compatibility. I initially worked on this in approximate mode and fully expected a problem at that step. I was pleasantly surprised.

[Allen]

Good morning and thank you for an interesting discussion.

Thomas, nice work on the formula!! From a different approach, I arrived at a surprisingly similar routine (perhaps the exact same formula but factored differently?) that saves two steps (5 bytes) in RPL and 2 bytes in RPN by reusing the intermediate value \( \begin{align} 3n \mod 4 \end{align} \) .

Originally, I spent some time trying to implement a different formula, which still might be possible,
\(
\begin{align}
a(n) = 6 \lfloor{\frac n 4}\rfloor + \frac{7(n \mod 4) - (n \mod 4)^2}{2}
\end{align}
\)

but arrived a something very similar to Thomas' elegant formula from last week:

\(
\begin{align}
a(n) = \frac{3n + 4(3n \mod 4) - (3n \mod 4)^2}{2}
\end{align}
\)

The formula looks uglier, but is slightly more efficient to calculate with only 6 mathematical operations and 2 less commands in RPL.

Code:


<< 3 * DUP 4 MOD DUP  4 - * - 2 / >>

(40 bytes, checksum #913Bh)


Somehow I feel the stack lifts could be optimized here to replicate the DUP function from RPL, but I've gotten rusty during covid.

Code:


00 { 18-Byte Prgm }
01 3
02 ×
03 RCL ST X
04 4
05 MOD
06 RCL ST X
07 4
08 -
09 ×
10 -
11 2
12 ÷
13 .END.

or python

Code:


example=[0, 3, 5, 6, 6, 9, 11, 12, 12, 15, 17, 18, 18, 21, 23, 24, 24, 27, 29, 30, 30, 33, 35, 36, 36, 39, 41, 42, 42, 45, 47, 48, 48, 51, 53, 
54, 54, 57, 59, 60, 60, 63, 65, 66, 66, 69, 71, 72, 72, 75, 77, 78, 78, 81, 83, 84, 84, 87, 89, 90, 90, 93, 95, 96, 96, 99]
def f(n):
    i=n%4
    return (3*n + 4*i - i**2)//2
assert [f(x) for x in range(len(example))]==example

[Thomas Klemm]

(08-27-2022 12:22 PM)Allen Wrote:  Somehow I feel the stack lifts could be optimized here to replicate the DUP function from RPL

This is what I came up with:

Code:

00 { 17-Byte Prgm }
01 3
02 ×
03 RCL ST X
04 4
05 MOD
06 STO- ST L
07 RCL× ST L
08 +
09 2
10 ÷
11 END

Thanks for your contribution.

[Thomas Klemm]

Or then longer, but saving one byte:

Code:

00 { 16-Byte Prgm }
01 3
02 ×
03 RCL ST X
04 4
05 MOD
06 STO- ST L
07 LASTX
08 ×
09 +
10 2
11 ÷
12 END

[Allen]

(08-28-2022 02:31 PM)Thomas Klemm Wrote:  Or then longer, but saving one byte:

Code:

00 { 16-Byte Prgm }

Nice work! I love how combining 2 ideas saved 20% from the original program size!!!