a new challenge



#2

A new year deserves a new challenge!

Write a program on whatever calculator you like to answer the following question: Are there 4 consecutive 4-digit integers that are each evenly divisible by the sum of their digits? If so, list them (there MAY be more than one set, and I probably wouldn't pose this challenge if there were no sets).

For example, if 1001 is evenly divisible by 2 and 1002 is evenly divisible by 3 and 1003 is evenly divisible by 4 and 1004 is evenly divisible by 5, you've got a set. Of course, only one of those 4 conditions is true in this example, so 1001-1004 don't work.

So post your program and tell how long it took. This challenge has no practical value that I can think of other than it's just fun to develop the algorithm and get it to work, and that's reason enough for me.


#3

How about this for the 12c:

01-       9    9
02- 9 9
03- 9 9
04- 9 9
05- 44 1 STO 1
06- 35 clx
07- 44 0 STO 0
08- 9 9
09- 44 5 STO 5
10- 9 9
11- 44 4 STO 4
12- 9 9
13- 44 3 STO 3
14- 9 9
15- 44 2 STO 2
16- 45 1 RCL 1
17- 45 2 RCL 2
18- 45 3 RCL 3
19- 40 +
20- 45 4 RCL 4
21- 40 +
22- 45 5 RCL 5
23- 40 +
24- 10 /
25- 43 24 FRAC
26- 43 35 x=0?
27-43,33 55 GTO 55
28- 35 Clx
29- 44 0 STO 0
30- 1 1
31- 44 30 1 STO- 1
32- 44 30 2 STO- 2
33- 16 CHS
34- 45 2 RCL 2
35- 43 34 x<=y?
36-43,33 38 GTO 38
37-43,33 16 GTO 16
38- 44 40 3 STO+ 3
39- 45 3 RCL 3
40- 43 34 x<=y?
41-43,33 43 GTO 43
42-43,33 14 GTO 14
43- 44 40 4 STO+ 4
44- 45 4 RCL 4
45- 43 34 x<=y?
46-43,33 48 GTO 48
47-43,33 12 GTO 12
48- 45 5 RCL 5
49- 1 1
50- 30 -
51- 43 35 x=0?
52-43,33 00 GTO 00
53- 44 5 STO 5
54-43,33 10 GTO 10
55- 1 1
56- 44 40 0 STO+ 0
57- 45 0 RCL 0
58- 4 4
59- 43 34 x<=y?
60-43,33 62 GTO 62
61-43,33 30 GTO 30
62- 45 1 RCL 1
63- 43 31 PSE
64- 3 3
65- 44 0 STO 0
66-43,33 30 GTO 30

It finds three solutions in about 80 seconds on my current model (i.e. fast) 12c. I'll leave the solutions as an exercise for the reader :-) However, I've verified them as being correct. Only the smallest of the four values is displayed.

It would be easy to speed this up to about thirty seconds by changing line 01 and line 08 to 3. If you are only interested in a single solution, change the pause to GTO 00 and watch it speed up even more.

- Pauli

[edit: wrote largest instead of smallest]


Edited: 6 Jan 2011, 12:30 a.m.


#4

Good work Paul. I love the 12c, and especially the super-fast new one.

I solved this first on the 32sii, then on the 30b. I'll post the code in a day or so. The speed of the 30b revealed some 5-digit solutions that would have taken much too long on the 32sii (probably require a new set of batteries!).

Hint: one of the 4-digit solutions is an IBM mainframe! That can't be a coincidence.

Don


#5

Quote:
The speed of the 30b revealed some 5-digit solutions that would have taken much too long on the 32sii (probably require a new set of batteries!).

Seven of them??

Thirty nine six digit ones. And 234 seven digit ones :-)

I didn't do this on a calculator though.


- Pauli


#6

Quote:
Seven of them??

Actually, I stopped it at 45,639 and it had found 4 5-digit series solutions by then. I was hoping to find more than 4 consecutive numbers somewhere, have you found that yet?

Did you use Mathematica or something?

Don


#7

Hi Don,

Quote:
I was hoping to find more than 4 consecutive numbers somewhere, have you found that yet?

In the range from 10 to 1E6, there're indeed only two 5-consecutive sets: starting at 131052 and 491424.
There're 22 more, if you go to 1E7.

In the range from 10 to 1E8, I found only 9 6-consecutive sets. (First one at 10000095.) Two of which are actually 7-consecutive (10000095 and 41441423).

After 3 minutes I found one 9-consecutive set at 12432xxxx. The last digits have vanished. Maybe to give your hp 30b something to do. ;-)

Oliver


Edited: 6 Jan 2011, 1:08 p.m.


#8

Quote:
After 3 minutes I found one 8-consecutive set at 12432xxxx. The last digits have vanished. Maybe to give your hp 30b something to do. ;-)

Well, my 30b was indeed up to the task! There are actually 9 consecutive numbers, beginning at 124324220 and going to 124324228.

Amazing!


#9

Don,

Yes, I noticed the jump from 7 to 9 consecutive numbers and found it intriguing enough to hunt a little more.

Feed your 30b this number 62162100000000 to make him feel great! (This one is a starting point, not a solution.)

:-)

Edited: 6 Jan 2011, 4:23 p.m.


#10

Wow, another niner'.

Amazing.

#11

Paul,

I'm curious.

Would you care to post your "non-calculator" program and timings?

On my Mac using CalcPad, the MorphEngine implementation reported for ND1 finds the 7-digit solutions in 10.8s, and the 1383 solutions for 8-digit numbers in 120s.

Cheers.

[EDIT: For completeness: With the improved algorithm, the timings are 0.6s and 6.1s, respectively.]


Edited: 8 Jan 2011, 3:13 p.m. after one or more responses were posted


#12

Quote:
Would you care to post your "non-calculator" program and timings?

Sure. The code is not in the least bit elegant and is trading copious amounts of memory for speed. It is in c++ for the bit vector library support.

Here is my initial version that does the problem as stated in essentially zero time on my old slow laptop:

#include <vector>
#include <stdio.h>

#define MAX 10000

int main() {
std::vector<bool> v(MAX, false);
for (int n = 1000; n<MAX; n++) {
int d1 = n%10;
int d2 = (n/10)%10;
int d3 = (n/100)%10;
int d4 = (n/1000)%10;
int ds = d1+d2+d3+d4;
if (n % ds == 0)
v[n] = true;
}
for (int n = 1003; n<MAX; n++)
if (v[n-3] && v[n-2] && v[n-1] && v[n])
printf("%d : %d\n", n-3, n);
return 0;
}


Here is a version that checks from 10,000,000 up to 100,000,000 in under 4 seconds on my rather old mac laptop:

#include <vector>
#include <stdio.h>

#define MAX 100000000

int main() {
std::vector<bool> v(MAX, false);
for (int n = 10000000; n<MAX; n++) {
int d1 = n%10;
int d2 = (n/10)%10;
int d3 = (n/100)%10;
int d4 = (n/1000)%10;
int d5 = (n/10000)%10;
int d6 = (n/100000)%10;
int d7 = (n/1000000)%10;
int d8 = (n/10000000)%10;
int ds = d1+d2+d3+d4+d5+d6+d7+d8;
if (n % ds == 0)
v[n] = true;
}
for (int n = 1006; n<MAX; n++)
if (v[n-6] && v[n-5] && v[n-4] && v[n-3] && v[n-2] && v[n-1] && v[n])
printf("%d : %d\n", n-3, n);
return 0;
}

The extension to search for the nine sequences should be obvious.
That program takes about forty five seconds.


- Pauli


#13

Nice, Paul! Thanks. I was hoping for Mathematica (as I'd like to see how it's done there) but this, too, is clean, compact and easy to read.


#14

Here is another version that much more closely mimics the 12c program above (counting up not down though):

#include <stdio.h>

int main() {
unsigned int a, b, c, d, e, f, g, h, i;
unsigned int n = 100000000;
unsigned int count = 0;
for (a=1; a<10; a++)
for (b=0; b<10; b++)
for (c=0; c<10; c++)
for (d=0; d<10; d++)
for (e=0; e<10; e++)
for (f=0; f<10; f++)
for (g=0; g<10; g++)
for (h=0; h<10; h++)
for (i=0; i<10; i++) {
const unsigned int s = a+b+c+d+e+f+g+h+i;
if (n % s == 0) {
count++;
if (count >= 9)
printf("%u - %u\n", n-count+1, n);
} else
count = 0;
n++;
}
return 0;
}


Seven seconds to find both 9 sequences.


- Pauli

#15

For reference my program does the full iteration from 9999 down to 1000. The timings I've given are for the entire range.


- Pauli

#16

Here's an implementation and timings for ND1.

This function finds the next solution vector, given a starting number:

function (x) { // arg is starting number
  function crossSumOfDecDigits(x) {
var v = String(x);
      var sum = 0;
      for (var i=v.length-1; i>=0; i--)
        sum += v[i] * 1.0;
      return sum;
  }

  function nextNCrossSumDivisibles(n, x) {
   var solutionVec = [];
   do {
     if (!(x % crossSumOfDecDigits(x)))
       solutionVec.push(x); // add candidate num to solution
     else
       solutionVec = []; // reset solution
     ++x;
   }
   while (solutionVec.length < n);

   return solutionVec;
  }

  return nextNCrossSumDivisibles(4, x);
}

Given the starting point 1E5, the first solution vector, [10307 10308 10309 10310], comes up after 5ms.

Given the starting point 1E6, the first solution vector, [1010022 1010023 1010024 1010025], comes up after 0.17s.

The following function uses the function above (named "crossSumDivisibles" in the same user folder) and returns a vector of solution vectors, given a search range, and reports timings:

function (start, end) {    
  Bench.start();
      
  var solutions = [];
  for (var i=start; i<end; i++) {
    var solution = crossSumDivisibles(i);
    solutions.push(solution);
    i = solution[0];
  }
  solutions.pop(); // we're computing one solution too many; drop it

  Bench.elapsedMsg();

  return solutions;
}

The 3 solutions for 4-digit numbers, as per challenge, are found in 0.13s, but somehow, mysteriously, are kept secret.

The 7 solutions for 5-digit numbers are found in 1.3s.

The 39 solutions for 6-digit numbers are found in 14.2s.

The 234 solutions for 7-digit numbers are found in 149s.

That's on my iPhone 3GS.

RPL on ND1 would run significantly slower on this task. (But I shall try running it, if someone posts a UserRPL solution.)

Cheers.


Edited: 6 Jan 2011, 12:34 p.m.


#17

Quote:
The 3 solutions for 4-digit numbers, as per challenge, are found in 0.13s, but somehow, mysteriously, are kept secret.

Maybe you'll have to pay Steve Jobs a couple of million $ to unlock the secret!

Very nice work. It really shows the fundamental differences between the new breed of "app" calculators running on i-devices and Androids versus the traditional dedicated RPN calculators.

Thanks.

Don

#18

I was able to speed up this implementation by a factor 20, by

- maintaining a running count of the digit sum instead of deriving it from x

- no longer using an array to hold candidates

The new code is

function (x) {   // arg is starting number
    function digitSumOf(x) {
        var s = String(x);
        var sum = 0;
        for (var i=s.length-1; i>=0; i--)
            sum += s[i] * 1.0;
        return sum;
    }
    
    function nextNivenSequence(n, x) {
        var digitSum = digitSumOf(x);
        var sequenceLength = 0;
        do {
            if (!(x % digitSum))
                ++sequenceLength;
            else
                sequenceLength = 0;
            
            if (++x % 10)
                ++digitSum;
            else {
                if (x % 100)
                    digitSum -= 8;
                else if (x % 1000)
                    digitSum -= 17;
                else if (x % 10000)
                    digitSum -= 26;
                else
                    digitSum = digitSumOf(x);
            }
        }
        while (sequenceLength < n);
        
        return x-sequenceLength;
    }

    return nextNivenSequence(4, x);
}

The algorithm retains the ability to start the search from any given number, is not confined to a certain number of digits, and will look for an adjustable count of consecutive numbers.

New timings:

The 3 4-digit solutions: 8ms

The 234 7-digit solutions: 6.5s


Edited: 8 Jan 2011, 5:23 a.m.

#19

Hello Don,

Interesting problem!

I was glad my first HP-33s program found the solution in less than five seconds... until I read Paul Dale's post and realized there were three sets of them. Not so nice, too many labels, not so optimized, anyway this second try finds the three sets in about 260 seconds (after 3 s, 132 s and 260 s, respectively).

C0001 LBL C
C0002 CLx
C0003 STO S
C0004 STO T
C0005 1000
D0001 LBL D
D0002 CF 1
D0003 RCL T
D0004 4
D0005 -
D0006 x=0?
D0007 GTO E
F0001 LBL F
F0002 x<>y
F0003 ENTER
F0004 ENTER
F0005 ENTER
F0006 10
F0007 INT/
F0008 STO+ S
F0009 LASTx
F0010 INT/
F0011 STO+S
F0012 LASTx
F0013 INT/
F0014 STO+S
F0015 Rv
F0016 RCS S
F0017 9
F0018 *
F0019 -
F0020 RMDR
F0021 0
F0022 x#y?
F0023 STO T
F0024 STO S
F0025 Rv
F0026 x=0?
F0027 SF 1
F0028 CLx
F0029 1
F0029 FS? 1
F0031 STO+ T
F0032 +
F0033 GTO D
E0001 LBL E
E0002 x<>y
E0003 LASTx
E0004 -
E0005 STOP
E0006 LASTx
E0007 +
E0008 x<>y
E0009 GTO F

Reference:

http://oeis.org/A007953

http://oeis.org/A054899

#20

It's a HP 50g solution but slow a bit. The 3 solutions: 52 sec, and the full loop from 1001 to 9999 is 3 min 50 sec.


\<< 0. {} 'CLST' STO 
1001. 9999.
FOR M
M
DUP
1000. /
DUP IP SWAP FP
1. 2. FOR N
10. *
DUP IP SWAP FP
NEXT
10. * + + + MOD
IF NOT THEN
1. + DUP
IF 4. == THEN
M 3. - 'CLST' STO+
DUP -
END
ELSE
DUP -
END
NEXT
DROP CLST REVLIST
\>>

#21

A restructured 50g version. A bit faster.
The solution: from 52 sec reduced to 44 sec
and the full loop from 3' 50" to 3' 14".

\<< 
{} 'CLST' STO 0.
1. 9. FOR A
0. 9. FOR B
0. 9. FOR C
0. 9. FOR D
A 1000. * B 100 *
C 10 * D + + + DUP
A B C D + + + MOD
IF NOT THEN
SWAP
1. + DUP
IF 4. == THEN
SWAP 3. - 'CLST' STO+
DUP -
ELSE
SWAP DROP
END
ELSE
SWAP DROP DUP -
END
NEXT
NEXT
NEXT
NEXT
DROP CLST REVLIST
\>>

#22

George: I bet if you find a way to write this without STO+, this'll run ~10 faster.


#23

Thanks Oliver. But in this case this STO+ playing three times only till the program is running.


#24

George,

Sorry about the bad lead! I didn't read your code properly and was just parsing it for slow instructions to explain the unexpected slow speed. (Surely a 50g must run this faster than a 12c!)

What about changing "B 100" to "B 100." and "C 10" to "C 10."?

Oliver

Edited: 8 Jan 2011, 9:43 a.m.


#25

The "B 100 and C 10" it's a text mistake only, I'm sorry for it. My 50g is in Approx mode and the program ran with "B 100. C 10." values.
The low speed related the level of my knowledge I think.


#26

This is my third HP 50g version. I just found an interesting topic on the HP48 forum about the local variables.
With this method the three solutions was found in 31" and the full loop's time was 2' 15".

\<< 1000. 1. \-> up dn 
\<< {} 'CLST' STO 0.
1. 9. FOR A
0. 9. FOR B
0. 9. FOR C
0. 9. FOR D
1. up + 'up' STO
1. dn + 'dn' STO
up dn MOD
IF NOT THEN
1. + DUP
IF 4. == THEN
up 3. - 'CLST' STO+
DUP -
END
ELSE
DUP -
END
NEXT
-9. dn + 'dn' STO
NEXT
-9. dn + 'dn' STO
NEXT
-9. dn + 'dn' STO
NEXT
DROP CLST REVLIST
\>>
\>>

Thanks Don for this challange. I'm waiting for the next one.

#27

It seems I missed the link.

https://groups.google.com/d/topic/comp.sys.hp48/Ujngf2y134A/discussion

With this link and with Raymund's method made my last version

#28

Programs that used only the stack used to be faster on the HP-48G/GX and the HP-49G. Since the original Saturn microprocessor is emulated on the HP-49G+/50g, I don't think this applies to the HP-50g. The 50g program below is base on the Turbo Pascal code somewhere else in this thread and takes roughly one minute to find all three solutions, starting at 1000. On my HP-48GX it took about 150 seconds, that is, the 50g is only 2.5 times faster in this case.

Regards,

Gerson.

%%HP: T(3)A(D)F(,);
\<< 0, ROT CLLCD
DO DUP DUP 0, OVER
DO 10, / IP SWAP OVER + SWAP DUP
UNTIL 10, <
END DROP 9, * - MOD
IF NOT
THEN SWAP 1, + DUP
IF 4, \>=
THEN OVER 2, DISP 1, WAIT
END SWAP
ELSE NIP 0, SWAP
END 1, + 3, PICK OVER
UNTIL <
END 3, DROPN
\>>

On the HP-48G/GX, replace NIP with SWAP DROP.

<< << 1000. 3035. DC >> TEVAL >> -> s: 60.1381


#29

Thanks for the nice program, Gerson. It takes one minute for me too. In this moment I don't know what is it doing but I will visit it's run with my debugger. :-) It seems my third version with the local variables is faster. If you modify it's middle section:


       IF 4. == THEN
up 3. - 'CLST' STO+
DUP -
CLST SIZE IF 3. == THEN KILL END <---- Here.
END

In this case the program takes 31 seconds and the three solutions are in the list named CLST.


#30

George,

At half a minute, it looks like you're reaching the limit of the HP-50g for the particular algorithm you've been using. Well done!

Quote:
In this moment I don't know what is it doing but I will visit it's run with my debugger. :-)

I haven't written any remarks, instead I am providing a Turbo Pascal program which has virtually the same structure, so it can be more easily followed. The bold-face text in the RPL listing is borrowed code from hpsolo I found at comp.sys.hp48 ( "sum of digits" contest).

Gerson.

%%HP: T(3)A(D)F(,);
\<< 0, ROT CLLCD
DO DUP DUP MANT 0,
DO OVER IP + SWAP FP MANT SWAP
UNTIL OVER NOT
END +
MOD
IF NOT
THEN SWAP 1, + DUP
IF 4, \>=
THEN OVER 2, DISP 1, WAIT
END SWAP
ELSE NIP 0, SWAP
END 1, + 3, PICK OVER
UNTIL <
END 3, DROPN
\>>

<< << 1000. 3035. DC >> TEVAL >> -> s: 55.3735


Program Niven;
var n1, n2, n: real;
c, i, m: byte;
k: char;
function Sod(x: Real): byte;
var s, t: real;
begin
s:=0;
t:=x;
repeat
t:=Int(t/10);
s:=s+t
until t<10;
Sod:=Trunc(x-9*s)
end;
begin
ClrScr;
Write('n1, n2, m: ');
ReadLn(n1, n2, m);
WriteLn;
n:=n1;
c:=0;
repeat
if Frac(n/Sod(n)) = 0
then
begin
c:=c+1;
if c>=m then WriteLn (n-i+1:18:0)
end
else
c:=0;
n:=n+1;
until n=n2
end.


#31

Gerson,

Thanks for the lists. I visited the first one. It was a good occasion to learn the stack manipulation. On trial i change the older DUP DUP and 3. PICK commands with the subsequent DUPDUP and PICK3 ones. The program's run time reduced with three hard seconds. :-)

Edited: 11 Jan 2011, 5:25 p.m.


#32

A small improvement and a version for the HP-28S:

HP-50g:

%%HP: T(3)A(D)F(,);
\<< { } 0, 1000,
DO DUP DUPDUP 10, / IP DUP 10, / IP DUP 10, / IP + + 9, * - MOD
IF NOT
THEN SWAP 1, + DUP
IF 4, \>=
THEN UNROT + LASTARG NIP ROT
END SWAP
ELSE NIP 0, SWAP
END 1, + DUP 10000,
UNTIL >
END DROP2
\>>

<< DC >> TEVAL -> {1017. 2025. 3033.}
s: 171.2803

( 38.3768 seconds when up to 3033)

HP-28S:

\<< { } 0 1000
DO DUP DUP DUP 10
/ IP DUP 10 / IP DUP
10 / IP + + 9 * -
MOD
IF NOT
THEN SWAP 1 +
DUP
IF 4 \>=
THEN OVER 4
ROLL + ROT ROT
END SWAP
ELSE SWAP DROP 0
SWAP
END 1 + DUP 3033
UNTIL >
END DROP2
\>>

DC -> { 3033 2025 1017 } (after 3.00 minutes)


#33

Gerson,


These lot of program lists are very interesting and instructive me.
Thanks, George


#34

Since positive multiples of 9, except 0, are Niven numbers, they can be skipped in the sum of digits routine. However, for 4-digit numbers it is faster to simply do the sum. This might be useful when testing larger numbers though. The contents THEN and the ELSE statements have been exchanged and the common SWAP instruction is now duely out of the structure. I hope there are no more obvious mistakes.

\<< { } 0. 1000.
DO DUP DUPDUP 10. / IP DUP 10,
/ IP DUP 10. / IP + + 9. * - MOD
IF
THEN NIP 0.
ELSE SWAP 1. + DUP
IF 4. \>=
THEN UNROT + LASTARG NIP
ROT
END
END SWAP 1. + DUP 3033.
UNTIL >
END DROP2
\>>

37.6573 sec


\<< { } 0. 1000,
DO DUPDUP 9. MOD
IF
THEN
DUPDUP 10. / IP DUP 10.
/ IP DUP 10. / IP + + 9. * - MOD
ELSE DROP 0.
END
IF
THEN NIP 0.
ELSE SWAP 1. + DUP
IF 4. \>=
THEN UNROT + LASTARG NIP
ROT
END
END SWAP 1. + DUP 3033.
UNTIL >
END DROP2
\>>

38.5524 sec


#35

Quote:
I hope there are no more obvious mistakes.

No, the lists are correct. (37.8970 and 39.1434. Your HP is faster a bit. :-)

Gerson, in this case (4 digit numbers) maybe isn't the fastest way: Simple to count the numbers one by one? (Oliver's method for example.)


#36

George,

I know that is not the best algorithm, I was just trying to optimize for speed the one I had already started. But you're right, there is no need to use a formula in this case. As we can see from the table, the sums of digits follow a quite predictable pattern. In the program below, the sum of digits is computed with that pattern in mind, but then again it is equivalent to a program that simply would sum the digits one by one in the innermost loop. So, no gain this time. The RPL program cycles from 1000 to 4000 and leaves the answers in the stack.

%%HP: T(3)A(D)F(.);
\<< 1000. 0. \-> n c
\<< 0. 2.
FOR i 0. 9.
FOR j 0. 9.
FOR k 1. 10.
FOR l n i j k l + + + / FP
IF
THEN 0. 'c' STO
ELSE 1. 'c' STO+ c 4. \>= { n } IFT
END 1. 'n' STO+
NEXT
NEXT
NEXT
NEXT
\>>
\>>


<< DC >> TEVAL -> 1017.
2025.
3033.
s:59.4524

------------------

n SOD(n) n SOD(n) n SOD(n)

1000 1 2090 11 3001 4
1001 2 2091 12 3002 5
1002 3 2092 13 3003 6
1003 4 2093 14 3004 7
1004 5 2094 15 3005 8
1005 6 2095 16 3006 9
1006 7 2096 17 3007 10
1007 8 2097 18 3008 11
1008 9 2098 19 3009 12
1009 10 2099 20 3010 4
1010 2 2100 3 3011 5
1011 3 2101 4 3012 6
1012 4 2102 5 3013 7
1013 5 2103 6 3014 8
1014 6 2104 7 3015 9
1015 7 2105 8 3016 10
1016 8 2106 9 3017 11
1017 9 2107 10 3018 12
1018 10 2108 11 3019 13
1019 11 2109 12 3020 5
1020 3 2110 4 3021 6
1021 4 2111 5 3022 7


------------------

Program Niven;

var c, i, j, k, l, sod: byte;
n: integer;

begin
ClrScr;
n:=1000;
c:=0;
for i:=0 to 9 do
for j:=0 to 9 do
for k:=0 to 9 do
for l:=1 to 10 do
begin
sod:=i+j+k+l;
if Frac(n/sod)=0 then
begin
c:=c+1;
if c>=4 then WriteLn(n:4)
end
else
c:=0;
n:=n+1
end;
end.


1017
2025
3033

>

------------------


#37

Quote:
I was just trying to optimize for speed the one I had already started.

Gerson,

I visited your programs and I tried to modify the first one from these. It was a good idea from you: to restructure the IF-END part. I changed the DO UNTIL structure to a FOR NEXT one and use the loop variable to count the dividend and so the stack commands decreased a bit.

\<< { } 0. 
1000. 3033. FOR n
n DUPDUP
10. / IP DUP 10. / IP DUP
10. / IP + + 9. * -
MOD
IF THEN
DROP 0.
ELSE
1. + DUP
IF 4. \>= THEN
DROP n 3. - + 0.
END
END
NEXT
DROP
\>>

34.6757 sec

Edited: 13 Jan 2011, 1:44 p.m.


#38

Great!

A couple of seconds less:

\<< { } 0. 1000. 3033. 
FOR n n DUPDUP MANT IP OVER .1
* IP DUP .1 * IP + + 9. * - MOD
IF
THEN DROP 0.
ELSE 1. + DUP
IF 4. \>=
THEN DROP n 3. - + 0.
END
END
NEXT DROP
\>>

32.5687 sec (average of three runs)

145.1378 sec (the whole cycle)


#39

Another couple of seconds:

\<< { } 0. 1000. 3033.
FOR n n DUP MANT DUP IP SWAP
FP MANT DUP IP SWAP FP MANT DUP
IP SWAP FP MANT + + + MOD
IF
THEN DROP 0.
ELSE 1. + DUP
IF 4. \>=
THEN DROP n 3. - + 0.
END
END
NEXT DROP
\>>

30.6459 sec (average of three runs)

141.292 sec (the whole cycle)


#40

Gerson, lenity me please! :-)

I just tried to understand your previous solution. I'll print the whole conversation between us and I'll read it a lot of days. And this link too.

Quote:
Another couple of seconds:

I understand: Isn't it your real goal, but however it was a nice solution from you.


George
#41

Here's another RPL version:

\<< {} 0 \-> results len \<<
   27
   1000 9999 FOR x

     IF x 10 MOD THEN
       1 +
     ELSE
       IF x 100 MOD THEN 8 -
       ELSE IF x 1000 MOD THEN 17 -
            ELSE 26 -
            END
       END
     END
     
     DUP x SWAP
     IF MOD THEN
       0 'len' STO
     ELSE
       IF len 3 == THEN
         'results' x len - STO+
       ELSE
         'len' 1 STO+
       END
     END
   NEXT
   results
\>> \>>

This one allows you to adjust the range to search but the initial digit sum (maintained on stack) needs to be in sync. Also, if you want to search into 5 digits, you need to add a "x 10000 MOD" case.

Run-time:

25.8s Emu48-based calculator on iPhone

10.8s ND1 (v1.3) on iPhone

105.3s 50g (thanks, George, for reporting)

Here's the equivalent JavaScript program:

function () { 
    var results = [], len = 0;
    var digitSum = 27;

    for (var x=1000; x<9999; x++) {
        if (x % 10)
            digitSum += 1;
        else {
            if (x % 100) digitSum -= 8;
            else if (x % 1000) digitSum -= 17;
            else if (x % 10000) digitSum -= 26;
        }
        
        if (x % digitSum)
            len = 0;
        else {
            if (len == 3)
                results.push(x-len);
            else
                len += 1;
        }
    }
    
    return results;
}

ND1 runs this one in 11ms. That is, 1000x faster than the RPL program. (2500x faster than Emulator on iPhone, ~200x faster than Emu48-emulated 49G on a Mac. ~10,000x faster than 50g!)

I'm working on technology that will transform ("code-morph") UserRPL into JavaScript for the obvious benefit.


Edited: 11 Jan 2011, 4:11 p.m. after one or more responses were posted


#42

The full loop is from 1000 to 9999: 105.3351 sec.
Found the three solutions is: 23.7689 sec. Congratulations Oliver!

Edited: 11 Jan 2011, 3:56 p.m.

#43

Quote:
Surely a 50g must run this faster than a 12c!

Oliver, maybe isn't it a 12c with an ARM processor?
http://en.wikipedia.org/wiki/HP-10C_series#HP-12C

#44

You're citing Wikipedia on this instead of this Forum, outrageous!
I think that everyone who regularly posts here should be deeply insulted, I know I am :)

Edited: 9 Jan 2011, 6:18 p.m.


#45

I'm sorry, Katie. I'm a new member and I don't want to insult enybody. I don't do it once more.

George


#46

Quote:

You're citing Wikipedia on this instead of this Forum, outrageous! I think that everyone who regularly posts here should be deeply insulted, I know I am :)

Quote:
I'm sorry, Katie. I'm a new member and I don't want to insult enybody. I don't do it once more.

George --

I'm quite sure that Katie's statement was "tongue-in-cheek" humor, not meant to be taken seriously.

I honestly didn't know about that fairly-new Wikipedia page you cited, and was surprised not to recognize any of the authors.

-- Karl


Edited: 10 Jan 2011, 2:34 a.m.


#47

Ok, Karl. But I'm really a greenhorn here and caution me please if I muss anything. :-)


#48

You're doing fine George, not to worry. Karl is right, Katie was just making some fun, hence her smiley face.

Thanks for your great contributions here (and everyone, really). We all benefit from the insights and knowledge of everyone. It is very interesting to see all of these different approaches to the same problem. Thanks for joining in.

Don


#49

George,

Don and Karl have that right, I was just kidding around. Anyway whoever wrote the Wikipedia article got a least some of their information from this site and from members' own sites (Eric Smith's in particular) if not this Forum.

Welcome, and thanks for your contributions!

-Katie


#50

Katie,

If i had known this! The next time preferably i'll citing the Bible. :-)

George

#51

In last december 5227 articles were removed from the hungarian Wikipedia site. All of them were copied from a literary history cyclopaedia.

#52

Quote:
Oliver, maybe isn't it a 12c with an ARM processor?

George,

Right. The 30b (also ARM) is faster still. Both surprising to me.

I'd love to see a SysRPL implementation...

Oliver

#53

I thought I would offer another solution for the hp-49g+ or the hp-50g using User RPL.  I was not concerned with speed but wanted to learn more about programing using lists and this turned out to be a great learning tool for me.  I only evaluated the 4-digit case.  

I have generated 2 sets of times, one for the emulator and one for the calc. The emulator is running on a 32bit Linux laptop under Wine. Times for full range for 4-digit numbers is 43 min 34 sec. A real battery eater! For the emulator 2 min 57 sec. The calc is in exact mode.

Listing

«
8996 { 1000 1001 1002 1003} {} {} {}
n a d1 d2 b
« TIME NEG 1 SF
1 n FOR k
a 'd1' STO
1 4 FOR j
d1 j GET 1000 / DUP @ CHANGE
IP SWAP FP
10 * DUP
IP SWAP FP
10 * DUP
IP SWAP FP
10 *
+ + + d1 j ROT PUT 'd1' STO @ CHANGE
NEXT
a d1 / FP 'd2' STO
CASE
d2 4 GET 0 > THEN k 2 + 'k' STO 2 SF END @ CHANGE
d2 3 GET 0 > THEN k 1 + 'k' STO 2 SF END @ CHANGE
d2 2 GET 0 > THEN 2 SF END @ CHANGE
d2 1 GET 0 > THEN 2 SF END @ CHANGE
END
IF 2 FC?C THEN
IF 1 FS?C THEN a 'b' STO
ELSE a b + 'b' STO END
END
{ 1000 1001 1002 1003} k ADD 'a' STO
NEXT
TIME HMS+ 4 FIX "TIME H.MMSS" + 2 FIX
b SORT
»
»


#54

Ronald, try set the calculator to approx mode. Maybe your program will run faster I think.


#55

George,
Thanks for the advice. For some strange reason I can't figure out, when in approximate mode the time for the full set is 1 hr 1 min and 3 sec. Go figure! Again I'm not concerned with speed. I just had a thought that because most of the work is being done on integers that this may cause exact mode to be faster. I don't really know but maybe someone as knowledge with hp calculators as John Myers may have an answer.


#56

I was imprecise I'm sorry. The goal is: Convert every numbers to floating point form in your program. From this: 100 to this: 100.

The easiest way is: First set the HP to Approx, then push your program list to the stack, call it to the edit line and press Enter to push it to the stack again. The floating point number is really faster, try it.


#57

George,
Your suggestion proved to be correct. By changing most integers to reals the time for the full 4-digit range is now 12 min 42 sec. I use reals and approximate mode in the future if I am concerned with speed. Thank you.

#58

I noticed that the \-> in the second line did not come through. Very important for defining local variables. Sorry I did not catch it before I posted.

#59

Many thanks to Paul, Oliver, Gerson, and George for posting your solutions to this problem. Very clever work.

I developed solutions for the HP-32sii, HP-30b, and the HP-12c+. Below I show the code for the 32sii and the 30b; the code for the 12c+ is just the same algorithm adjusted slightly for the 12c+ instruction set.

The three solution sets are: (1014, 1015, 1016, 1017), (2022, 2023, 2024, 2025) and (3030, 3031, 3032, 3033).

The 32sii took 9.5 minutes to find the third solution. The 30b took 18 seconds, and the 12c+ took 48 seconds.

Among 4-digit numbers, the highest number of consecutive numbers is 4 (the three solutions shown above). When you open it up to larger numbers, there are some examples of 9 consecutive numbers that meet the criteria. I find that amazing.

Here is the code for the 32sii and 30b:

For the 32sii, enter starting number to check then XEQ F

A = sum of digits accumulator
B = number being tested
C = number of consecutive numbers divisible by SOD

lbl F entry point
clrvars so A and C are initialized to 0
sto B current number to test
lbl A begin loop to get SOD
10
/
fp
sto+A will x 10 later
lastx
ip
x=/=0
goto A loop until all digits summed
rcl B
10
rcl*A gives SOD
/ see if divisible by SOD
fp
x=/=0
goto B number not divisible by SOD
1
sto+C is divisible, count
4
rcl C
x<y if fewer than 4, don't display
goto C
view B the 4th consecutive number (or 5th or 6th...!) number
view C the frequency (found at least 4 in a row)
lbl C get ready to test a new number
0
sto A SOD accumulator
1
sto+B next number to test
rcl B so it will be in x
goto A try next number (program never ends until you R/S)
lbl B this number was not divisible by SOD, so clear C
0
sto C
goto C


30b solution:

This is a direct translation of the 32sii program to the 30b.

Reg. 0 = 32sii reg. A
Reg. 1 = 32sii reg. B
Reg. 2 = 32sii reg. C

Label 00 = 32sii label A
Label 01 = 32sii label B
Label 02 = 32sii label C

Press Amort to run after keying in number to start with.

(0 Amort)
sto 1
0
sto 0
sto 2
swap
lbl 00
10
/
math up = (yes, all that for fp)
sto+0
ans
math up up = (ip)
0
swap
?=/=
gt 00
rcl 1
10
rcl*0
/
math up =
gt 01
1
sto+2
rcl 2
4
?<
gt 02
rcl 1
r/s
rcl 2
r/s
lbl 02
0
sto 0
1
sto+1
rcl 1
goto 00
lbl 01
0
sto 2
goto 02



#60

Don,

Thanks for literally giving us some joy (see link below). Thanks also for the didactly explained algorithm.

http://en.wikipedia.org/wiki/Harshad_number

http://oeis.org/A141769

Too bad (or too good) no joy beyond 20 consecutive Harshad numbers :-)

Regards,

Gerson.


#61

Thanks for those links Gerson, I wasn't even aware that there was such a thing as a Harshad number. I think my little 30b would have run out of juice looking for that 20th consecutive Harshad number!

Don


#62

Quote:
I wasn't even aware that there was such a thing as a Harshad number.

Quote:
This challenge has no practical value that I can think of other than it's just fun to develop the algorithm and get it to work, and that's reason enough for me.

I wasn't aware of them either. You were quite right about their ludic purpose. Thanks again.

Gerson.


#63

Using Don's HP32Sii program on my HP42S, it took a dismal 59 minutes to find the third solution. I did realize how slow the old technology is.


#64

The following is a straightfoward conversion of Don's HP-32SII program to the HP-42S, except for the sum of digits algorithm. All three solutions to his original problem come up instantaneaously on Free42. If line 28 is changed to 5, it will take 4.6 s for the first 5-number sequence solution, starting with 10 (Free42 @ 1.86 GHz)

00 { 64-Byte Prgm }
01>LBL "DC"
02 STO 02
03 CLX
04 STO 01
05 STO 03
06>LBL 01
07 RCL 02
08 ENTER
09>LBL 04
10 10
11 /
12 IP
13 STO+ 01
14 X#0?
15 GTO 04
16 X<>Y
17 ENTER
18 ENTER
19 RCL 01
20 9
21 *
22 -
23 MOD
24 X#0?
25 GTO 02
26 1
27 STO+ 03
28 4
29 RCL 03
30 X<Y?
31 GTO 03
32 RCL 02
33 RCL 03
34 STOP
35>LBL 03
36 CLX
37 STO 01
38 1
39 STO+ 02
40 GTO 01
41>LBL 02
42 CLX
43 STO 03
44 GTO 03
45 .END.

Usage:

Starting number XEQ DC for the first solution then R/S for the next ones.

#65

Quote:
Using Don's HP32Sii program on my HP42S, it took a dismal 59 minutes to find the third solution.

Still based upon Don's program (except for the sum of digits formula), the following does it in one-third of the time:

00 { 63-Byte Prgm }
01>LBL "DC"
02 STO 02
03 CLX
04 STO 03
05>LBL 01
06 RCL 02
07 ENTER
08>LBL 00
09 10
10 BASE/
11 STO+ ST Z
12 X#0?
13 GTO 00
14 CLX
15 RCL ST Y
16 9
17 RCL* ST T
18 -
19 MOD
20 X#0?
21 GTO 02
22 1
23 STO+ 03
24 4
25 RCL 03
26 X<Y?
27 GTO 03
28 RCL 02
29 RCL 03
30 STOP
31 GTO 03
32>LBL 02
33 CLX
34 STO 03
35>LBL 03
36 1
37 STO+ 02
38 CLX
39 STO 01
40 GTO 01
41 .END.

Keystrokes Display Elapsed time

1000 XEQ DC y: 1,017 x: 4 0 min 11.0 s
R/S y: 2,025 x: 4 9 min 44.4 s
R/S y: 3,033 x: 4 19 min 21.2 s


#66

About 15% faster:

   HP-42S                        			 HP-12C+

00 { 65-Byte Prgm }
01>LBL "DC" 01 STO 2
02 STO 02 02 1
03 CLX 03 0
04 STO 03 04 STO 0
05>LBL 01 05 CLx
06 RCL 02 06 STO 1
07 ENTER 07 STO 3
08 ENTER 08 RCL 2
09 10 09 RCL 0
10>LBL 00 10 /
11 BASE/ 11 INTG
12 STO+ ST Z 12 STO+ 1
13 LASTX 13 RCL 0
14 X<=Y? 14 x<=y
15 GTO 00 15 GTO 10
16 Rv 16 RCL 2
17 Rv 17 RCL 2
18 ENTER 18 RCL 1
19 ENTER 19 9
20 9 20 *
21 RCL* ST T 21 -
22 - 22 /
23 MOD 23 FRAC
24 X#0? 24 x=0
25 GTO 02 25 GTO 27
26 1 26 GTO 36
27 STO+ 03 27 1
28 4 28 STO+ 3
29 RCL 03 29 3
30 X<Y? 30 RCL 3
31 GTO 03 31 x<=y
32 RCL 02 32 GTO 38
33 RCL 03 33 RCL 2
34 STOP 34 R/S
35 GTO 03 35 GTO 38
36>LBL 02 36 CLx
37 CLX 37 STO 3
38 STO 03 38 1
39>LBL 03 39 STO+ 2
40 1 40 CLx
41 STO+ 02 41 STO 1
42 CLX 42 GTO 08
43 GTO 01 43 GTO 00
44 .END.


Keystrokes Display Elapsed time Keystrokes Display Elapsed time

1000 XEQ DC y: 1,017 x: 4 0 min 09.4 s 1000 R/S 1,017.000000 0 min 01.0 s
R/S y: 2,025 x: 4 8 min 11.5 s R/S 2,025.000000 0 min 37.6 s
R/S y: 3,033 x: 4 16 min 14.8 s R/S 3.033.000000 1 min 14.1 s

#67

Agreed! This was a pretty cool challenge, and I loved seeing the different implementations and timings.

Thanks Don!

Bruce

#68

Yes, thanks for those links! Nice to read up on this after getting one's hands dirty.

I was only able to give Don a 14-digit number yesterday, because I tripped over an obvious pattern after finding the first 9-consecutive numbers:

  X{0}n2[0-8]

where

     X: special numbers { 1243242, 6216210, 21135114, 31081050, ... }
{0}n: 0 to n "0" digits
[0-8]: the digits "0" through "8", making the set

Example: 310810500000020, 310810500000021, 310810500000022, ...

That is, there're infinitely many 9-consecutive numbers following this pattern. (The Wikipedia page confirms this in general.)

All 9-consecutive Niven numbers I found are Nivenmorphic. (They also shared the same digit sum, 20.)

Anyone with a mathematical mind who cares to prove, or disprove, that this applies to all 9- (and higher?) consecutive Niven numbers?

The Wikipedia page is silent on that...


#69

Quote:
I was only able to give Don a 14-digit number yesterday, because I tripped over an obvious pattern after finding the first 9-consecutive numbers:

  X{0}n2[0-8]

Nice finding! I think my HP notebook (T2350 @ 1.86GHz) would have to be at least 1,000 times faster to check all through the 1e14 to 1e15 range for these sequences in a day's time, using the TurboBCD program below. Perhaps I little less if I switched to C :-)

Program Niven;

var n1, n2, n: real;
c, i, m: byte;
k: char;

function Sod(x: Real): byte;

var s, t: real;

begin
s:=0;
t:=x;
repeat
t:=Int(t/10);
s:=s+t
until t<10;
Sod:=Trunc(x-9*s)
end;

begin
ClrScr;
Write('n1, n2, m: ');
ReadLn(n1, n2, m);
WriteLn;
n:=n1;
c:=0;
repeat
if Frac(n/Sod(n)) = 0
then
begin
c:=c+1;
if c>=m then
begin
for i:=m downto 1 do
WriteLn (n-i+1:18:0);
repeat
read(kbd,k)
until not KeyPressed;
WriteLn
end;
end
else
c:=0;
n:=n+1;
until n=n2
end.


----------------------------------------

n1, n2, m: 3.108105e14 3.108106e14 9

310810500000020
310810500000021
310810500000022
310810500000023
310810500000024
310810500000025
310810500000026
310810500000027
310810500000028

#70

Good old HP41C:

Somewhat different approach: No digit counting.
The divisor is increased and decreased in line with the increasing dividend. Correction with -9 at each 10/100/1000 cycle.

I watched until solution 2 (202x .. 202y), then my battery died
So I did not test it to the end!

01LBL "D0111"
02 TIME
03 CLA
04 FIX 4
05 ATIME
06 AVIEW
07 FIX 0
08 CLA
09 1.01
10 STO 01
11 STO 02
12 STO 03
13 STO 04
14 STO 05
15 1.004
16 STO 06
17 STO 07
18 STO 08
19 0
20 STO 09
21 999
22 STO 10
23LBL 01
24 RCL 05
25 STO 06
26 1
27 ST+ 09
28 ST+ 10
29 RCL 10
30 RCL 09
31 /
32 FRC
33 X=0?
34 GTO 03
35 RCL 08
36 STO 07
37 GTO 02
38LBL 03
39 ISG 07
40 GTO 02
41 CLA
42 FIX 4
43 TIME
44 ATIME
45 FIX 0
46 RCL 10
47 3
48 -
49 ": " (Remark: Append)
50 ARCL X
51 " .. " (Remark: Append)
52 ARCL 10
53 AVIEW
54 1
55 ST- 07
56LBL 02
57 ISG IND 06
58 GTO 01
59 9
60 ST- 09
61 RCL 05
62 STO IND 06
63 ISG 06
64 GTO 02
65 .END.

BR. Raymund


#71

The same approach as my previous post but a solution in C.
Looks as a candidate for a recursive approach with a variabele amount of digits....

void func2()
{
int Dividend = 999;
int Divisor = 0;
int Hit = 0;
for(int i=0; i<=9; i++){
for(int j=0; j<=9; j++){
for(int k=0; k<=9; k++){
for (int l=0; l<=9; l++) {
Dividend++;
Divisor++;
int r = Dividend%Divisor;
if(!r)
{
Hit ++;
if(Hit == 4)
{
printf("%d .. %d\n",Dividend-3,Dividend);
Hit--;
}
}
else
{
Hit = 0;
}
}
Divisor -= 9;
}
Divisor -= 9;
}
Divisor -= 9;
}
}

BR. Raymund


#72

Quote:
for(int i=0; i<=9; i++){

????
shouldn't it be replaced with

for(int i=1; i<=9; i++){

?

it is a 10% optimizing

Patrice


#73

Hello Patrice

Since we have 10 digits we need 10 steps from 0 up to and including 9. Therefore from 0 .. <=9. Variant: from 1 .. <=10 but this doesn't cover the digits - although it has the same functionality.

And .. I tried your proposal to ckeck if my considerations are OK (was not 100% sure). Result, starting with 1 doesn't lead to the solution (with the program as presented).

Thanks for the review comment, it forced me to think about my assumption!

BR. Ray


#74

Hi Raymound,

add

printf("%d\n",Dividend);

at the end of your program, I think you will have a surprise!

Patrice


#75

Hi Patrice,

Yes! a surprise. I'll look at it next week. Thanx !!

I'm pleased to know that this is not the first bug (TheFirstBug)

BR

Raymund

#76

Hi Patrice

of course, starting at 1000 .... Yes, your remark

Quote:
for(int i=1; i<=9; i++){

is fully correct!

BR and good night

Ray

#77

Hmm, among the 3-digit numbers I find one instance of 4 consecutive integers that satisfy this. It includes a number that any computer scientist would recognize.


#78

Thanks, Don, for a new mini-challenge -- in trying to understand George's second 50g program, I tried modifying it to search the three-digit numbers, and it got the one solution you mentioned. Just
eliminate the D loop and the other D, change "A 1000." etc. to "A 100. * B 10. * C" and eliminate one plus sign in each group of three, and also one NEXT. It took only 15 seconds to find it.


#79

Glenn: Could you post here this modified program list please?


#80

Ok George, now with your new program I am sure to lose some sleep!
Here is your "trimmed down" program for 3-digit numbers:
<< {} 'CLST' STO 0.
1. 9. FOR A
0. 9. FOR B
0. 9. FOR C
A 100. * B 10. * C + +
DUP A B C + + MOD
IF NOT THEN
SWAP 1. + DUP
IF 4. == THEN
SWAP 3. - 'CLST' STO+ DUP -
ELSE SWAP DROP END
ELSE SWAP DROP DUP -
END
NEXT
NEXT
NEXT
DROP CLST REVLIST >>
Thanks for continuing this challenge! Cheers, Glenn


#81

Thanks Glenn. It's a good challenge, not too easy to discontinue it. :-)

#82

Hello, maybe I am writing a bit late but for those interested the link http://oeis.org/ offers you the chance to write some numbers and the page tries to find the sequence for them. If you write 1014, 2022, 3033 it finds the answer. Best wishes!


#83

Quote:
1014, 2022, 3033

That didn't work for me, so I tried 1014, 2022, and 3030, that worked. Very interesting site, thanks much.

Don

#84

Here is a solution to the challenge in User RPL that takes a completely different approach. On a 50g, it runs in 123 seconds. I
haven't attempted any coding optimization - the time is due entirely
to the algorithm - so I'm sure it can be improved.

Instead of checking each possible 4-digit number to see whether it's
divisible by the sum of its digits, this program starts with a
possible sum of digits for four consecutive numbers. It then checks
the much smaller set of numbers that can be divided by this set.

For example, the sum of the digits in the first solution (beginning at 1014) are 6,7,8, and 9. It turns out that any set of numbers N, N+1, N+2, N+3 that are divisible by 6,7,8 and 9, must be of the form N = 6 + 504k where k is an integer. So for the sums {6 7 8 9} we only have to check every 504'th number.

The program is quite long and this will be a long explanation. I hope you'll find it interesting. I implemented the program is a directory for simplicity (and debugging!) and I'll refer to the various programs in the directory.

Throughout this explanation, I'll use N, N+1, N+2 and N+3 (or N0, N1, N2 and N3) to refer to the sequential 4-digit numbers, and I'll use s0, s1, s2 and s3 to refer to the sum of the digits in each of these numbers.

Enumerating the sums of digits

To begin, lets figure the possible values for the sum of the digits of four sequential numbers between 1000 and 9999. Clearly the smallest sum is 1 (for 1,000) and the largest sum is 36 (for 9999). So any sum must between 1 and 36.

For most sequences of N, the sum of their digits is also a
sequence. E.G, for N = {1001, 1002, 1003, 1004 }, s = {2 3 4 5}.
So one possible set of sums is:

    {s s+1 s+2 s+3} for s= 1 through 33
Now consider N = {1207 1208 1209 1210}. The digit sums are s = {10 11 12 4}. In general, when the one's digit rolls over from 9 to 0, you lose 9 in the digit sum, but you gain 1 because the 10's digit increments. So another possible set of sums is:
    {s s+1 s+2 s-6} for s=8 to 33

The roll-over from 9 to 0 could occur anywhere in the sequence, so
other possibilities are:

   (s s+1 s-7 s-6) for s=9 to 34
(s s-8 s-7 s-6) for s=10 to 35
The numbers in N can also roll over from 99 to 00 or 999 to 000, these give the possible sums:
   {s s-17 s-16 s-15} for s=19 to 35
{s s+1 s-16 s-15} for s= 18 to 34
{s s+1 s+2 s-15 } for s= 17 to 33

{s s-26 s-25 s-24} for s=27 to 36
{s s+1 s-25 s-24} for s=26 to 35
{s s+1 s+2 s-24} for s=25 to 34

In my solution, all of these lists are calculated by the SUMS program. The result of SUMS is a list of lists where each sublist is 4 numbers that represent the sums of the digits of 4 consecutive numbers.

Restricting the start and end values

The challenge asks us to find a 4-digit number - one that's between
1000 and 9999, but if we start with a possible sum of the digits in N, we can use that information to restrict the upper and lower bounds of the search. For example, suppose we're considering the sums 10, 11, 12, and 13. The smallest 4-digit number whose digits add up to 10 is 1009 and the largest is 9001, so there's no point in checking numbers less than 1009 or larger than 9001.

The subprogram BIGandSMALL calculates the largest and smallest 4-digit numbers whose digits sum to each of the values between 1 and 36. It leaves the results on the stack in two arrays. The main program calls BIGandSMALL once and stores the results in global variables BIGGEST and SMALLEST. For example, BIGGEST[12] is the biggest 4 digit number whose digits sum to 12. Later when we check a particular set of sums, we use BIGGEST and SMALLEST to get the upper and lower bounds of the loop.

Mathematics

Suppose we're given a list {s0 s1 s2 s3} that is a solution for the
sequence of numbers {N N+1 N+2 N+3}. That means that N is evenly
divisible by s0, N+1 is evenly divisible by s1 etc. Mathmatically:

N   == 0 mod s0
N+1 == 0 mod s1
N+2 == 0 mod s2
N+3 == 0 mod s3

or equivalently:
N == 0 mod s0
N == -1 mod s1
N == -2 mod s2
N == -3 mod s3

(I'm using "==" to mean "congruent to," which is usually indicated by a symbol with three horizontal lines).

What are the values of N that satisfy these equations? The first
equation means that N = k0*s0 for any integer k0. Plugging that into the second equation gives k0*s0 = -1 mod s1. Can we simplify this into the form k0 = r mod s? The answer is yes.

Solving modulo equations - there's an app for that!

In a more general sense, how do you solve Ax == B mod C? This is done with the Linear Congruence Theorem. The procedure goes like this:

  1. If (B / GCD(A,C)) not an integer,then there is no solution
  2. Let D = GCD(A,C).
  3. Find r and s such that rA + sC = D
  4. The solution is x == (rB/D) mod (C/D)

How do you find r and s? It's built-in to the 50G! The IEGCD command does it for you.

Here is the function that solves a single modulo equation:

@ This solves an equation of the form ax = b mod n
@ Where "=" really means "is congruent to".
@ It returns r s such that x = r mod s, or
@ equivalently x = s*k + r
@ arguments are a b n (must be integers)
@ Results are r s 1. if there's a solution, or
@ 0. if no solution exists.
@ See http://en.wikipedia.org/wiki/Linear_congruence_theorem
@ and the IEGCD command for details.
SLVMOD
\<< \-> a b n \<<
@ If GCD(a,n) doesn't divide b evenly then no solution
IF b a n GCD / DUP IP \=/ THEN
0.
ELSE
@ It's good!
a n IEGCD DROP
b * OVER /
n ROT /
1. @ Indicates success
END
\>>
\>>

Solving simultaneous modulo equations

Let's return to the original problem. We want to find N such that it
simultaneously solves:

N ==  0 mod s0
N == -1 mod s1
N == -2 mod s2
N == -3 mod s3
More generally, suppose we want to solve two simultaneous equations:
X == b mod a, and
X == r mod s
we do it as follows:
X = a*k0 + b	     (for any integer k0)
a*k0 + b == r mod s (substituting for X into the second equation
k0*a == (r-b) mod s (subtract b)
We can use SLVMOD on this equation to find:
   k0 == u mod v, or equivalently
k0 = k1*v + u for any integer k1
and substituting for k0 into the equation above:
X = a*(k1*v+u) + b
= (a*v)k1 + (a*u + b)
X == (a*u+b) mod (a*v)
The program SSLVMOD solves these simultaneous modulo equations:
@ Given two equations:
@ x = r mod s, and
@ x = b0 mod a0
@ This program returns b1 and a1 such that
@ x = b1 mod a1
@ Arguments: b0 a0 r s
@ results: b1 a1 1. (if there is a solution), or
@ 0. if not
SSLVMOD
\<< \-> b a r s
\<< a r b - s
IF SLVMOD
THEN a * SWAP a * b + SWAP 1.
ELSE 0.
END
\>>
\>>

Putting it all together

The somewhat poorly named FACTS program takes a list of possible sums of digits {s0 s1 s2 s3} and computes a and b such that N == a mod b and

N ==  0 mod s0
N == -1 mod s1
N == -2 mod s2
N == -3 mod s3
Remember, this means that if {s0 s1 s2 s3) are a possible sum-of-digits, then the vales N == a mod b are the only possible numbers where N, N+1, N+2 and N+3 are evenly divisible by s0, s1, s2, and s3 respectively.

For example, given the list of possible sums-of-digits {4 5 6 7}, FACTS determines that N == 4 mod 420 are the only values that will work.

The CHKSUM prorgam takes a list of possible digit sums and checks for a solution to the challenge that matches that list. In other words, it checks to see if there is a sequence of four 4-digit numbers that are (1) evenly divisible by the numbers in the list, AND (2) whose sum of digits equals the values in the list.

CHKSUM uses FACTS to determine the possible values that satisfy the
first criteria. Then it uses the BIGGEST and SMALLEST arrays arrays
to look up the largest and smallest numbers that satisfy the
sum-of-digits criteria. Finally, it uses a FOR loop to check each of the possible values to see if it matches the sum-of-digits criteria.

The DIGITS program takes an integer and computes the sum of its digits.

The MAINP program is the main entry point to the algorithm. It
computes and stores the biggest and smallest 4-digit numbers for each possible sum. Then it calls SUMS to create the list of all possible sums of digits of sequences of 4-digit numbers. Finally, it uses DOLIST to call CHKSUM on each of these possible sums.

Conclusions

I think this is a good example of how you can squeeze performance out of a program by exploiting the properties of the problem that you're trying to solve. By starting with the potential sum of digits and then working forward to the possible 4-digit number (rather than the other way around), this program cuts down on the work needed:

  1. Computing the biggest and smallest numbers with the given sum cuts the range of numbers that you need to check.
  2. Computing the possible numbers that satisfy the divisibility test cuts the numbers down even more.

The program also demonstrates that it's worth stepping back from a
problem and asking "is there a different approach that might work."


Program Listing


%%HP: T(3)A(R)F(.);
DIR
MAINP
\<< BIGandSMALL 'SMALLEST' STO 'BIGGEST' STO SUMS 1.
\<< CHKSUM
\>> DOLIST
\>>
DIGITS
\<< 0. SWAP
WHILE DUP
REPEAT DUP 10. MOD ROT + SWAP 10. / IP
END DROP
\>>
FACTS
\<< 0. \-> L n
\<< 0. 1. 1. L 1.
\<<
IF SWAP
THEN n 'n' 1. STO- SWAP SSLVMOD
ELSE DROP 0.
END
\>> DOLIST
\>>
\>>
CHKSUM
\<< \-> L
\<< L
IF FACTS
THEN L HEAD { } \-> a b N RES
\<< SMALLEST N GET a - b / CEIL b * a + BIGGEST N GET
IF DUP2 \<=
THEN
FOR I
IF I DIGITS N ==
THEN 1. 1. 3.
FOR K I K + DIGITS L K 1. + GET == AND
NEXT
IF
THEN I 'RES' STO+
END
END b
STEP
IF RES SIZE
THEN RES
END
ELSE DROP DROP
END
\>>
END
\>>
\>>
SUMS
\<< 1. 33.
FOR i i DUP 1. + DUP 1. + DUP 1. + 4 \->LIST
NEXT 8. 33.
FOR i i DUP 1. + DUP 1. + i 6. - 4 \->LIST
NEXT 9. 34.
FOR i i DUP 1. + i 7. - i 6. - 4 \->LIST
NEXT 10. 35.
FOR i i DUP 8. - DUP 1. + DUP 1. + 4 \->LIST
NEXT 19. 35.
FOR i i DUP 17. - DUP 1. + DUP 1. + 4 \->LIST
NEXT 18. 34.
FOR i i DUP 1. + DUP 17. - DUP 1. + 4 \->LIST
NEXT 17. 33.
FOR i i DUP 1. + DUP 1. + DUP 17. - 4 \->LIST
NEXT 28. 35.
FOR i i DUP 26. - DUP 1. + DUP 1. + 4 \->LIST
NEXT 27. 34.
FOR i i DUP 1. + DUP 26. - DUP 1. + 4 \->LIST
NEXT 26. 33.
FOR i i DUP 1. + DUP 1. + DUP 26. - 4 \->LIST
NEXT 186. \->LIST
\>>
BIGandSMALL
\<< 0. 36. NDUPN \->ARRY DUP \-> bigres smallres
\<< 1. 36.
FOR i 0. 1000. i \-> val mult cur
\<<
WHILE cur
REPEAT cur 9. MIN 'cur' OVER STO- mult * 'val' STO+ 'mult' 10. STO/
END 'bigres' i val PUT 1. 'mult' STO i 1. - 'cur' STO 1000. 'val' STO
WHILE cur
REPEAT cur 9. MIN 'cur' OVER STO- mult * 'val' STO+ 10. 'mult' STO*
END 'smallres' i val PUT
\>>
NEXT bigres smallres
\>>
\>>
SLVMOD
\<< \-> a b n
\<<
IF b a n GCD / DUP IP \=/
THEN 0.
ELSE a n IEGCD DROP b * OVER / n ROT / 1.
END
\>>
\>>
SSLVMOD
\<< \-> b a r s
\<< a r b - s
IF SLVMOD
THEN a * SWAP a * b + SWAP 1.
ELSE 0.
END
\>>
\>>
END

#85

David, I'm speechless, and that doesn't happen that often. I'm going to print your post and read it carefully, I really like what you did here.

This thread will fade into archive oblivion one day, but your post cannot. This needs to be immediately available for all to see and benefit from. Please enter this post as an article.

This is truly great work and I look forward to reading it in detail, even though I'm not an RPL guy. You (and Allen) may force me to buy and learn the 50g!

Don


#86

Thanks Don. I really enjoyed this challenge.

To me, the best part of my solution came at the deepest part of the math when I needed to solve an equation in modulus math. I spent days trying to work out a solution by myself but when I finally turned to wikipedia to see how it's done, it turns out that the 50g has the critical function (IEGCD) built in!

I've posted the challenge and my solution as an article. I noted the dates of the original posts so any curious future readers can find the thread in the archives.

Dave

#87

Beautiful. Thank you for this and the nice presentation.


Possibly Related Threads...
Thread Author Replies Views Last Post
  RE: 35s sorting routine challenge - Gene's Challenge Miguel Toro 4 588 08-01-2007, 08:36 AM
Last Post: Miguel Toro

Forum Jump: