+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Calculators (and very old HP Computers) (/forum-3.html)
+--- Forum: General Forum (/forum-4.html)
+--- Thread: OEIS A228297: Need Help to Find Algorithm (/thread-9092.html)
OEIS A228297: Need Help to Find Algorithm - Gerald H - 09-16-2017 02:41 PM
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.
RE: OEIS A228297: Need Help to Find Algorithm - Thomas Okken - 09-16-2017 03:20 PM
http://www.poslarchive.com/math/math/papers/11-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.
RE: OEIS A228297: Need Help to Find Algorithm - Gerald H - 09-16-2017 04:41 PM
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
;
;RE: OEIS A228297: Need Help to Find Algorithm - Thomas Okken - 09-16-2017 06:27 PM
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.
RE: OEIS A228297: Need Help to Find Algorithm - Gilles59 - 09-16-2017 06:46 PM
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' STOWith 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 ?
RE: OEIS A228297: Need Help to Find Algorithm - Thomas Okken - 09-16-2017 11:24 PM
(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. :-)
RE: OEIS A228297: Need Help to Find Algorithm - Gerald H - 09-17-2017 05:15 AM
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.
RE: OEIS A228297: Need Help to Find Algorithm - Gerald H - 09-17-2017 11:55 AM
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.
RE: OEIS A228297: Need Help to Find Algorithm - Thomas Okken - 09-17-2017 03:30 PM
(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' STOI 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.
RE: OEIS A228297: Need Help to Find Algorithm - Gerald H - 09-17-2017 04:08 PM
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.
RE: OEIS A228297: Need Help to Find Algorithm - Thomas Okken - 09-17-2017 04:51 PM
(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 ENDThis 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. :-)
RE: OEIS A228297: Need Help to Find Algorithm - Gerald H - 09-18-2017 10:44 AM
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.
RE: OEIS A228297: Need Help to Find Algorithm - Paul Dale - 09-18-2017 11:39 AM
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
endConversion to a calculator is left as an exercise
The % operation is mod or remainder.
Pauli
RE: OEIS A228297: Need Help to Find Algorithm - Gerald H - 09-18-2017 04:34 PM
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. ENDRE: OEIS A228297: Need Help to Find Algorithm - Gerald H - 09-19-2017 05:30 AM
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
»