OEIS A228297: Need Help to Find Algorithm

[Gerald H]

I can't manage to find a way of representing this series

https://oeis.org/A228297

through some formula or recursion.

Would be very grateful for suggestions.

[Thomas Okken]

http://www.poslarchive.com/math/math/pap...-42139.pdf

The abstract appears to indicate that it describes the solution for the case at hand, as well as similar recurrence relations with an odd number of terms.

[Gerald H]

Thank you, Thomas for your suggestion, resulting in the programme below, which does function, but for large input, eg 20, it becomes extremely slow (31.5s).

Suggestions welcome.

Code:

Size: 180.

CkSum: # 2911h

::
  CK1&Dispatch
  BINT1
  ::
    %ABSCOERCE
    '
    ::
      DUP
      BINT6
      #<
      ITE
      ::
        BINT0
        #MAX
      ;
      ::
        DUP
        BINT6
        #=
        casedrop
        BINT5
        DUP
        DUP#1-
        1GETLAM
        EVAL
        #-
        1GETLAM
        EVAL
        OVER
        #1-
        DUP#1-
        1GETLAM
        EVAL
        #-
        1GETLAM
        EVAL
        #+
        OVER
        #2-
        DUP#1-
        1GETLAM
        EVAL
        #-
        1GETLAM
        EVAL
        #+
        OVER
        #3-
        DUP#1-
        1GETLAM
        EVAL
        #-
        1GETLAM
        EVAL
        #+
        SWAP
        #4-
        DUP#1-
        1GETLAM
        EVAL
        #-
        1GETLAM
        EVAL
        #+
      ;
    ;
    DUP1LAMBIND
    EVAL
    ABND
    UNCOERCE
  ;
;

[Thomas Okken]

That paper is beyond me, and I'm not familiar with SysRPL, but here's a memoized implementation in C which is quite fast:

Code:

#include <stdio.h>
#include <stdlib.h>

int *m = NULL;
int sz = 0;

int a(int n) {
    int r;
    if (n < 0)
        return 0;
    if (n < sz) {
        r = m[n];
        if (r != -1)
            return r;
    } else {
        int i;
        m = realloc(m, (n + 1) * sizeof(int));
        for (i = sz; i <= n; i++)
            m[i] = -1;
        sz = n + 1;
    }
    if (n <= 5)
        r = n;
    else
        r = a(n-a(n-1)) + a(n-1-a(n-2)) + a(n-2-a(n-3)) + a(n-3-a(n-4)) + a(n-4-a(n-5));
    m[n] = r;
    return r;
}

int main(int argc, char *argv[]) {
    int i, n;
    for (i = 1; i < argc; i++)
        if (sscanf(argv[i], "%d", &n) != 1)
            printf("can't parse \"%s\"\n", argv[i]);
        else
            printf("a(%d) = %d\n", n, a(n));
    return 0;
}

If you throw very large arguments at it, you may get a segmentation fault caused by the process running out of stack space. ulimit -s unlimited prevents that, although the program does need O(n) memory, for the stack and the memoization array, so you're going to hit some limit eventually.

[Gilles59]

Code:

« -> n
  «
    CASE 
      'n<=s+1' THEN 1 END 
      'n<=s+k' THEN n s - END 
      `∑(j=1,k,a(n-j-(s-1)-a(n-j)))`
    END
  »
»
'a' STO

With 0 's' STO , 5 'k' STO
I get correct result for a(1) to a(5)

but a(6) returns "THEN error :undefined name"

What is the problem ?

[Thomas Okken]

(09-16-2017 06:46 PM)Gilles59 Wrote:  

Code:

« -> n
  «
    CASE 
      'n<=s+1' THEN 1 END 
      'n<=s+k' THEN n s - END 
      `∑(j=1,k,a(n-j-(s-1)-a(n-j)))`
    END
  »
»
'a' STO

With 0 's' STO , 5 'k' STO
I get correct result for a(1) to a(5)

but a(6) returns "THEN error :undefined name"

What is the problem ?

I'm not familiar with the backquote notation you're using in the default case of your CASE statement (but I don't know RPL very well so that's probably just me). I replaced the backquotes with straight single quotes, and added EVAL to force the expression to be evaluated, and I changed the first case to

'n<=s' THEN 0 END

Then the program works, although, like Gerald's program, it becomes too slow very quickly because of the exponentially increasing run time. You really need to memoize the calculation to keep its run time fast enough, but because of my lack of familiarity with RPL, adapting my C code is left as an exercise to the reader. :-)

[Gerald H]

I should say that for input N I only want the Nth element returned, not the whole series up to N, which may lead to speedier algorithms.

[Gerald H]

Thank you, Thomas, for your suggestion, but I fear any algorithm that uses a complicated recursion is doomed to get bogged down & be inefficient.

Similarly Gilles, a nicely small programme, but the recursions will prohibit use for large numbers.

I don't think any recursive programme will be successful.

[Thomas Okken]

(09-17-2017 11:55 AM)Gerald H Wrote:  Thank you, Thomas, for your suggestion, but I fear any algorithm that uses a complicated recursion is doomed to get bogged down & be inefficient.

Similarly Gilles, a nicely small programme, but the recursions will prohibit use for large numbers.

I don't think any recursive programme will be successful.

I disagree:

Code:

<< -> N
  <<
    IF N 0 <=
    THEN 0
    ELSE
      IF 'M' DUP EVAL SAME
      THEN [ -1 ] 'M' STO
      END
      IF N M SIZE LIST-> DROP DUP ROT <
      THEN N 1 ->LIST M SWAP RDM SWAP 1 + N
        FOR I I -1 PUT
        NEXT 'M' STO
      ELSE DROP
      END M N GET DUP
      IF -1 ==
      THEN DROP
        IF N 5 <=
        THEN N
        ELSE 'A(N-A(N-1))+A(N-1-A(N-2))+A(N-2-A(N-3))+A(N-3-A(N-4))+A(N-4-A(N-5))' EVAL
        END DUP M SWAP N SWAP PUT 'M' STO
      END
    END
  >>
>>

'A' STO

I don't have a real calculator to test this on, but in m48 on my iPhone, set to emulate a 48GX at "authentic speed," this evaluates a(100) in 137 seconds.

It looks like evaluating A(N) always requires evaluating all of A(1) through A(N-1), so the algorithm could also be written as an iteration, which would be equivalent in terms of the calculations it performs, but probably much faster.

[Gerald H]

Thomas, I installed your programme on a real 50g & have checksum A401 & size 560.5.

For input 100 the programme returned 81 in 334s.

While this is the best published programme so far (we don't know what the others have & not published), I would like the range of input on the real calculator to be at least in the millions.

[Thomas Okken]

(09-17-2017 04:08 PM)Gerald H Wrote:  I would like the range of input on the real calculator to be at least in the millions.

That would require an algorithm that doesn't store all of A(1)...A(N). I don't know if that's possible.

Iterative implementation for HP-42S:

Code:

00 { 134-Byte Prgm }
01▸LBL "AA"
02 X>0?
03 GTO 00
04 0
05 RTN
06▸LBL 00
07 5
08 X<>Y
09 X≤Y?
10 RTN
11 STO "N"
12 1
13 +
14 1
15 NEWMAT
16 STO "REGS"
17 5
18▸LBL 01
19 STO IND ST X
20 DSE ST X
21 GTO 01
22 6
23 STO "I"
24▸LBL 02
25 RCL "I"
26 STO "J"
27 DSE "J"
28 RCL- IND "J"
29 RCL IND ST X
30 RCL "J"
31 DSE "J"
32 RCL- IND "J"
33 X<>Y
34 RCL+ IND ST Y
35 RCL "J"
36 DSE "J"
37 RCL- IND "J"
38 X<>Y
39 RCL+ IND ST Y
40 RCL "J"
41 DSE "J"
42 RCL- IND "J"
43 X<>Y
44 RCL+ IND ST Y
45 RCL "J"
46 DSE "J"
47 RCL- IND "J"
48 X<>Y
49 RCL+ IND ST Y
50 STO IND "I"
51 1
52 STO+ "I"
53 RCL "I"
54 RCL "N"
55 X≥Y?
56 GTO 02
57 R↑
58 END

This calculates a(100) is 97 seconds, on a real HP-42S.

UPDATE: The program doesn't avoid tight-memory situations; it really should do CLV "REGS" before NEWMAT, and do CLX after STO "REGS", to make sure nothing is competing with REGS for the calculator's RAM.
Still, even it its current form, running in Free42 on an iPhone 5s, it calculates a(2,000,000) in under a minute, so Gerald's requirement can be met, as long as you use the right equipment. :-)

[Gerald H]

OK, I have now stored your programme on an otherwise empty 42S (due to memory problems on the unit I customarily use) & confirm your timing for input 100 & correctness of answer.

Bravo!

Given enough computing power & sufficient energy the problem may be solved, although you could add longevity to your requirements.

However, I'll stick to the real 42S or 50g.

As for specifications, different for different models:

For 50g: Given any natural number input the correct answer is returned in a reasonable time, ie time characteristic of algorithm should be less than linear.

For 42S: Given any natural number < 5*10^11input the correct answer is returned in < 20s.

[Paul Dale]

This Lua program seems to produce the correct sequence:

Code:

for i = 1, 10000 do
    k = 1
    repeat
        print(i)
        k = k * 5
    until i % k ~= 0
end

Conversion to a calculator is left as an exercise
The % operation is mod or remainder.


Pauli

[Gerald H]

I have run up against a new problem.

If you use

N 5 / IP

to find the integer quotient of N divided by 5 on a 42S an erroneous answer is returned for large numbers, eg for

999999999999

the answer is

200000000000

whereas it should be

199999999999

I have a programme to do the calculation correctly without disturbing the stack but believe it is inefficient.

Could someone suggest a more efficient programme that leaves the stack intact?

Code:

0.    { 29-Byte Prgm }
1.    LBL “IQT5”
2.    R↑
3.    STO 01
4.    R↓
5.    STO 02
6.    5
7.    MOD
8.    +/-
9.    RCL+ 02
10.    5
11.    /
12.    R↑
13.    X<> 01
14.    R↓
15.    END

[Gerald H]

Here my programme for 42S

Store 5 in variable "5".
Unreliable for input > 5*10^11.

Code:

0.    { 50-Byte Prgm }
1.    LBL “Qs0k5”
2.    ENTER
3.    ENTER
4.    RCL/ "5"
5.    IP
6.    - 
7.    LBL 00
8.    ENTER
9.    ENTER
10.    X<> 00
11.    CLX
12.    X<> 00
13.    LBL 02
14.    RCL/ "5"
15.    IP
16.    STO+ 00
17.    X≠0
18.    GTO 02
19.    X<> 00
20.    +
21.    RCL- ST Z
22.    X≥0?
23.    GTO 01
24.    SIGN
25.    -
26.    GTO 00
27.    LBL 01
28.    R↓
29.    END

& for 49G

Code:

Qs0k5

« DUPDUP 5 IQUOT -
  WHILE DUP 0 OVER
    DO 5 IQUOT DUP
ROT + SWAP DUP NOT
    UNTIL
    END DROP +
PICK3 <
  REPEAT 1 +
  END NIP
»