Programming puzzle: Longest list of regular numbers?

[pier4r]

Description:
A regular number is a number that divides evenly a power of 60 (so at least one of the following numbers 60, 3600, 216000, \( 60^4 \) , \( 60^5 \) and so on)

Challenge #1:
Print on the screen all the regular numbers found in 120 seconds without gaps (so, without focusing on a subset of the regular numbers. One has to return consecutive regular numbers). If the program fails to produce any regular number after 120 seconds (unlikely, since the first is 2), one can just report the first regular number found and how much time was needed.
Note: No quick start allowed. One has to start from 2.

For example by hand using my el506w as support (or using the properties of regular numbers, those are faster to apply mentally) I can do this in 120 seconds:

2,3,4,5,6,8,9,10,12,15,16,18,20

Challenge #2 (longer):
Find the 1429th regular number, return the time needed for the computation (I'm not sure if this would stretch thin some old calculators).
Note: No quick start allowed. One has to start finding numbers from 2.

Challenge #2a (for very fast hw/solutions):
Like #2, but one has to find the 9733rd regular number. (for computers, one has the original challenge that raises the bar a bit)

Source of the idea.

[pier4r]

So first attempt for challenge 1.

126 seconds to get the first 60 numbers on the stack (then I collected them in a list manually). Hp 50g.

Code:


DIR
  @ Problem: www.hpmuseum.org/forum/thread-8178.html
  
  @Description: A regular number is a number that divides evenly a power of 
  @60 (so at least one of the following numbers 60, 3600, 216000, 604 , 605 
  @and so on) 
  @
  @Challenge #1: Print on the screen all the regular numbers found in 120 
  @seconds without gaps (so, without focusing on a subset of the regular 
  @numbers. One has to return consecutive regular numbers). If the program 
  @fails to produce any regular number after 120 seconds (unlikely, since 
  @the first is 2), one can just report the first regular number found and 
  @how much time was needed. Note: No quick start allowed. One has to start 
  @from 2. 
  @
  @For example by hand using my el506w as support (or using the properties 
  @of regular numbers, those are faster to apply mentally) I can do this in 
  @120 seconds: 
  @
  @2,3,4,5,6,8,9,10,12,15,16,18,20 
  @
  @Challenge #2 (longer): Find the 1429th regular number, return the time 
  @needed for the computation (I'm not sure if this would stretch thin some 
  @old calculators). Note: No quick start allowed. One has to start finding 
  @numbers from 2. 
  @
  @Challenge #2a (for very fast hw/solutions): Like #2, but one has to find 
  @the 9733rd regular number. (for computers, one has the original 
  @challenge that raises the bar a bit) 


  
  @in short: smart ways to compute regular numbers, that divides evenly at least
  @ one power of 60.
  
  mainTimed
  \<<
    @assuming the input on the stack
    'regNumC1v1' TEVAL
  \>>
  
  regNumC1v1
  @basic approach, dividing power of 60.
  @idea is to create the sequence without gaps of the regular numbers
  
  @ so far 21 50 16.04.2017 it works, pretty slow though. 126 seconds to get
  @ the list up to 400. 60 regular numbers
  \<<
    @the idea using factors seems smarter than basic division, but it is not easy
    @to formalize, so for the moment, basic divisions or checking factors every time.
    
    2 @counterV
    0 @numElV
    0 @factorV
    10 @uFbreak
    \->
    @input
    upperLimit
    @local var
    counterV
    numElV
    factorV
    uFbreak
    \<<
      -3 CF
        @this for FACTORS
    
      2 upperLimit
      FOR counter
        counter FACTORS
          @as result list of factors on the stack
          
        OBJ\->
          @exploded list, number of elements on the stack L1
        'numElV' STO
        
        @we check the factors (in even positions)
        @preparation and then loop
        2 'counterV' STO
        uFbreak CF
        WHILE
          uFbreak FC?
          counterV numElV \<=
          AND @ we don't quit the while and we did not finish going through the list.
        REPEAT
          counterV PICK 'factorV' STO
          IF
            factorV 2 \=/
            factorV 3 \=/
            AND
            factorV 5 \=/
            AND
          THEN
            @not a regular number
            uFbreak SF
          END
          2 'counterV' STO+
        END
        
        numElV DROPN
          @clear stack
        
        uFbreak FC? counter IFT
          @if there was no break in the factor check loop
          @leave the number on the stack
      NEXT
    \>>
  \>>
END

the second challenge fails to produce the 1429th number. I'm waiting since 8h at least, so the code is obviously slow. I do suppose it is the exact mode with FACTORS. Let's try to take it away.

[e_emil]

Greetings from a new user on the forum

105 seconds on a 50g was enough to find 580 first "regular" numbers, put them in order on the stack and tag each number by its consecutive number.

580'th regular number is 2073600

RPL code:
<< 20 2. LN 3. LN / 2. LN 5. LN / 0. 0. -> smax a b m n
<< 1. smax
FOR s 0. s
FOR x 0. s 1. + x - a * FLOOR
FOR y s x - y a / - b * DUP CEIL 0. MAX SWAP b + FLOOR DUP ROT ==
IF THEN 5. SWAP ^ 2. x ^ 3. y ^ * * 1. 'm' STO+
ELSE DROP
END
NEXT
NEXT
m ->LIST SORT 'k+n' 'k' 1. m 1. SEQ 2. << ->TAG >> DOLIST
OBJ-> DROP m 'n' STO+ 0. 'm' STO
NEXT
>>
>>

Description of the code: the outer loop (FOR s) finds all regular numbers between 2^s and 2^(s+1). The list is then sorted and output on the stack. The code proceeds to the next value of s, until smax (20 in the example above).

In other words, the code above finds all regular numbers smaler than 2^21, and it seems that there are 580 such numbers.

The number immediately after the last one found, will be always a power of two and it can be added manually:

580th number = 2073600 (last line of the program)
581st number = 2^21 = 2097152

I was thinking geometrically when writing the code: in a 3-d coordinate system, every point with non-negative integer coordinates corresponds to a regular number, i.e. (x,y,z) corresponds to 2^x*3^y*5^z. The inequality s <= x + y*(ln 3/ln2) + z(ln 5/ln2) < s+1 represents a "slab" in the coordinate system which contains all regular numbers greater than 2^s but smaller than 2^(s+1). (You can say that "slabs" are perpendicular to the direction of increasing regular numbers.) My code searches each such slab for points of integral coordinates, then generates a sorted list of 2^x3^y5^z and proceeds to the nest "slab" a bit further from the origin. "smax" in the beginning is the total number of slabs to be searched.

[pier4r]

neat! Can you explain (editing the previous post) your approach? Otherwise with all those stack operations I have to go through it.

[e_emil]

Challenge 2 - solved in 334 seconds on a 50g

According to my program, the 1429th number is 604661760

(1427: 600000000
1428: 603979776
1429: 604661760
1430: 607500000)

Execution time: 5 minutes and 34 seconds (50g). The first number put on the stack in the program (smax) has to be changed to 29. The last number found was 1062882000 (number 1544). The 1545th number would be then 2^30 = 1073741824

[Dave Britten]

Tip: These are just the "Hamming numbers".

A quick-and-dirty TI-86 program for a simplistic algorithm got me H(395)=295,245.

A similar program for my Casio fx-9860g Slim failed at about 24 seconds, as it had already found 999 and hit the maximum length of a list. Would have to work around this by garbage collecting unneeded entries from the beginning of the list to keep it short.

[pier4r]

E_emil neat idea the one of coordinate system. I do believe it is quite efficient.

I will try to replace the call to factors in my code to see if there is a speed up


Anyone else with hp 41, 42 , 48 family, prime, 17 b solver, other HPs, various casio-s, various TIs ?
Someone with the newRPL ?

[e_emil]

Thanks for hint, Dave! Good to hear that well-known algorithms exist

By the way, my nSpire CX CAS solves challenge 2 in just 15 seconds! (1543 first regular numbers, which means all of them up to 2^30). It's exactly the same algorithm as in the RPL function in one of the previous posts. I'm not quite familiar with the nSpire language and I couldn't find a command to sort a list within a program, so the function returns a partly-sorted answer and the final sorting has to be done outside of the program.

   

regular(29) returns in 15 seconds (compare it with more than 5 minutes on a 50g!!!). SortA is done in a fraction of a second.

The list returned by the function is "partly-sorted" in a sense that the order of powers of two is preserved, but the numbers between any two consecutive powers of two are not in order.

And here's the code:

Define regular(smax)=
(C) finds all regular numbers up to (but not including) 2^(smax+1)
(C) output list has to be sorted after the function returns
Func
:Local a,b,x,y,zlow,z,s,m,n,lpart,lall
:a:=((ln(2.))/(ln(3.)))
:b:=((ln(2.))/(ln(5.)))
:lpart:={}
:lall:={}
:For s,1,smax
:For x,0,s
: For y,0,floor((s+1-x)*a)
: zlow:=(s-x-((y)/(a)))*b
: z:=floor(zlow+b)
: If ceiling(zlow)=z Then
: lpart:=augment(lpart,{2^(x)*3^(y)*5^(z)})
: EndIf
: EndFor
:EndFor
(C) The list "lpart" should have to be sorted at this point but for some reason nSpire's SortA doesn't work in a function
(C) so the output list will be sorted after the function returns
:lall:=augment(lall,lpart)
:lpart:={}
:EndFor
:Return lall
:EndFunc

Update to challenge #1: the nSpire CX CAS with the above code finds 6909 first numbers (all up to 2^51) in about 90 seconds. Sorting of the output list using "SortA" is still around a second or two. "Resource exhaustion" happens shortly after this point, so to solve challenge 3, the function should not store intermediate numbers in a list.

[pier4r]

7000 numbers do exhaust the 64 mb of the Nspire?

Anyway thanks for the results. I tried to change the usage of factors with another function and I got even slower. I like the 3d approach, it is more elegant than the algorithms mentioned by Wikipedia, although maybe a bit slower but really nice.

[Dave Britten]

Still trying to get the Casio to find terms beyond H(998), but Casio BASIC is terrible and throws a syntax error if you try to do something simple like Seq(List 26[X],X,L,N) to extract a sublist.

It's irritating how many things on a Casio just randomly can't be used in certain contexts. This will probably work fine on the TI, but it's way slower.

EDIT: Got the Casio playing ball. Even with disgusting manual list copying in a For loop, it's still fast.

H(1429)=604,661,760 (37 sec.)

Still wrestling with H(9733). I know the program will slow down as it runs. List garbage collection becomes more frequent as the look-back range widens, but it's coming up on the 8-minute mark, so either I've underestimated how much it slows down, or I've coded some clever infinite loop.

[e_emil]

(04-17-2017 11:24 AM)pier4r Wrote:  7000 numbers do exhaust the 64 mb of the Nspire?

The innermost loop appends every new number at the end of the existing list which may not be an efficient way of managing memory on the nSpire? I'll try with a pre-allocated list of sufficient length and see if there's any difference. Anyway, speed difference compared with the 50g is remarkable.

I read the wikipedia page on Hamming numbers and it appears that they use a very similar geometric approach to estimate the number of regular numbers not greater than N: the inequality (ln 2)x + (ln 3)y + (ln 5)z <= ln N is satisfied by all regular numbers that are less or equal than N (assumming that a number of the form 2^x3^y5^z corresponds to the point (x,y,z)), which is the same idea that I used in the algorithm

EDIT: Challenge 3 solved!

the 9733rd number is 198135565516800000 (exactly - zeros at the end are _not_ due to rounding). Can be also written as 2^30*3^10*5^5.

9732: 197912092999680000
9733: 198135565516800000
9734: 200385994162176000

Time: 4 minutes 20 seconds on the nSpire CX CAS, using a modification of my previous program (calculating but not storing intermediate results). Had to switch to the "exact" mode, because the numbers already have 17 digits. The greatest power of two smaller than H(9733) is 2^57, so "smax" has to be 57 or more.

[pier4r]

(04-17-2017 11:32 AM)Dave Britten Wrote:  Still wrestling with H(9733). I know the program will slow down as it runs. List garbage collection becomes more frequent as the look-back range widens, but it's coming up on the 8-minute mark, so either I've underestimated how much it slows down, or I've coded some clever infinite loop.

Cool. I got interested by the regular/humming numbers because after a while they need a platform to handle big numbers, and handling big numbers by itself is non trivial (surely they cannot fit the CPU registers after a while). So maybe you got slowed down also due to this problem. I do believe in general it could be used as benchmark, I may add it to the wiki4hp.

Besides, which algorithm are you applying? The one mentioned on wikipedia or some elaborated by you?

Edit: could you elaborate on extracting sublists. Is there no built in function?

[pier4r]

(04-17-2017 12:28 PM)e_emil Wrote:  I read the wikipedia page on Hamming numbers and it appears that they use a very similar geometric approach to estimate the number of regular numbers not greater than N: the inequality (ln 2)x + (ln 3)y + (ln 5)z <= ln N is satisfied by all regular numbers that are less or equal than N (assumming that a number of the form 2^x3^y5^z corresponds to the point (x,y,z)), which is the same idea that I used in the algorithm

Yes the fastest algorithm known (but I do like yours more, what about putting it in the talk page? So maybe it will be included one day) is quite smart but it requires a good memory management, indeed on /r/dailyprogrammer the problem for many was to use a proper data structure.

Quote:Time: 4 minutes 20 seconds on the nSpire CX CAS, using a modification of my previous program (calculating but not storing intermediate results). Had to switch to the "exact" mode, because the numbers already have 17 digits. The greatest power of two smaller than H(9733) is 2^57, so "smax" has to be 57 or more.

Well done and yes, I chose a "big enough" number (well I was unsure at first) to push the mathematical environment to handle integers bigger than the canonical 32 or 64 bits or 10-15 digits.

I would be really interested in RPN programmable calculators, prime and so on.

[Han]

Only just now read the thread, but this should be useful: https://oeis.org/A051037

Additional resources: https://rosettacode.org/wiki/Hamming_numbers

[pier4r]

@Han, nice. But I do not see RPL or RPN keystroke programmign on the list! (neither casio basic, ti basic, ti lua ), time to fix it! Sure one can adapt algorithms.

(04-16-2017 09:05 PM)e_emil Wrote:  Challenge 2 - solved in 334 seconds on a 50g

According to my program, the 1429th number is 604661760

You algorithm is very interesting, I adapted one using the known algorithm how to solve it and still it is slower than yours (I'm not sure how much does it cost to put all the numbers in a list though).

378 seconds for 1429 numbers.

For the ti nspire I saw that you used Ti basic, have you tried with Lua? Should be a bit faster and more capable.

Edit_note: I wanted to put the list of 1429 numbers in a GROB and then print it, but the 50g does not have enough memory.

code for 50g.

Code:


  regNumC2v3
  @ using known algorithm
  @ not so original but goot to understand refined solutions if one does not
  @ find them by themselves.
  @ See Dijkstra
  @ www.reddit.com/r/dailyprogrammer/comments/60b63i/weekly_27_mini_challenges/df6hd31/
  @ the idea is that a "tree" is built, with branches emanating from 2,3 and 5.
  @ and since it is a tree, the branch emanates from the last useful value used by a
  @ certain factor. For more info, see the clues above.
  
  @ 141 seconds for 600 numbers
  @ 378 seconds for 1429 numbers.
  \<<
    @the idea using factors seems smarter than basic division, but it is not easy
    @to formalize, so for the moment, basic divisions or checking factors every time.
    
    1 @twoStackLvl
    1 @threeStackLvl
    1 @fiveStackLvl
    0 @newValue
    0 @numFound
    10 @uFbreak
    11 @uFbool
    \->
    @input
    upperLimit
    @local var
    twoStackLvl
    threeStackLvl
    fiveStackLvl
    newValue
    numFound
    uFbreak
    \<-uFbool
    \<<
      @TODO: save and restore flags.
    
      -3 SF
      -105 SF
        @faster with reals instead of exact mode.
    
      1
        @starting value on the stack
    
      WHILE
        numFound upperLimit <
      REPEAT
        
        @building the next regular number
        twoStackLvl PICK 2 *
        threeStackLvl 1 + PICK 3 *
          @ we need to compensate for the number just added on the stack
        MIN
        fiveStackLvl 1 + PICK 5 *
          @ we need to compensate for the number just added on the stack
        MIN
        
        DUP 'newValue' STO
        
        @ now the next regular number is on the stack
        @ so we need to update the values, especially
        @ stack pointers because everything is lifted by one.
        1 'numFound' STO+
        1 'twoStackLvl' STO+
        1 'threeStackLvl' STO+
        1 'fiveStackLvl' STO+
        
        @now we update the last usable value for the branches
        @relative to every factor compared to the last value
        twoStackLvl PICK 2 * newValue \<= 
        \<< 'twoStackLvl' 1 STO- \>>
        IFT
        
        threeStackLvl PICK 3 * newValue \<= 
        \<< 'threeStackLvl' 1 STO- \>>
        IFT
        
        fiveStackLvl PICK 5 * newValue \<= 
        \<< 'fiveStackLvl' 1 STO- \>>
        IFT
      END
      
      @create a list with the result.
      numFound \->LIST
    \>>
  \>>

[Han]

(04-17-2017 03:30 PM)pier4r Wrote:  You algorithm is very interesting, I adapted one using the known algorithm how to solve it and still it is slower than yours (I'm not sure how much does it cost to put all the numbers in a list though).

378 seconds for 1429 numbers.

Creating a list is expensive due to the memory allocations and possible garbage collection. Also, I wonder if your use of IFT is also compounding the time. I wonder if IF THEN END is faster than IFT.

[pier4r]

Good question. I did not write IF THEN ELSE due to more verbosity (and since we are in RPL, IFT is understandable too for small code.

Someone should launch a loop to test a condition with IF THEN ELSE and IFT.

Did a test.

Programs:




just the FOR is stretched up to 1000 to be sure. Then TEVAL.

IF then else : 40.90 seconds
IFT: 41.80 seconds

[e_emil]

Just 2% difference, you should probably repeat it several times to make sure it's not random

At last, I generated a complete list of (almost) 10,000 first regular numbers in one variable! I got around the memory problem on the nspire by allocating a table of zeros beforehand and making the program fill in the table without changing the structure. The entire calculation took about 13 minutes (just creating a table of zeros was half a minute, a total surprise to me). Sorting the output table is no more than 2 seconds using the built-in command.

   

the code:

Define hamming(smax)=
(C) A global table 'h' has to be defined with enough space for all generated numbers.
(C) The table is not sorted at output. The calculator has to be in exact mode if smax>=39
Prgm
:Local a,b,x,y,zlow,z,s,m
:a:=((ln(2.))/(ln(3.)))
:b:=((ln(2.))/(ln(5.)))
:m:=0
:For s,1,smax
:For x,0,s
: For y,0,floor((s+1-x)*a)
: zlow:=(s-x-((y)/(a)))*b
: z:=floor(zlow+b)
: If ceiling(zlow)=z Then
: m:=m+1
: 2^(x)*3^(y)*5^(z)→h[m]
: EndIf
: EndFor
: EndFor
: EndFor
:EndPrgm

Learning lua has never been my priority so far...

[Han]

I used the HP48 and solved #2 in 3 minutes 24 seconds. I think it could probably be improved a slight bit with pure stack storage.

Also, the 1429th number is actually 603979776 and not 604661760 (this is the 1430th number). I think e_emil is off by one. Place index on the stack and run the program below:

EDIT: Fixed a few typos (I typed this out by hand and did some sloppy copy/paste/edit).

Code:


\<<
  2 2 1 2 3 5 \-> n i j k x2 x3 x5
  \<<
    1 2 n
    IF DUP2 > THEN DROP2 ELSE

      FOR l
        x2 x3 MIN x5 MIN DUP l 1 + ROLLD
        IF x2 OVER == THEN
          1 'i' STO+ i PICK 2 * 'x2' STO
        END
        IF x3 OVER == THEN
          1 'j' STO+ j PICK 3 * 'x3' STO
        END
        IF x5 == THEN
          1 'k' STO+ k PICK 5 * 'x5' STO
        END
      NEXT

    END
  \>>
\>>

Output is basically all the terms in the sequence, with the n-th term placed on the n-th level of the stack.

[pier4r]

You just solved 2 in 204 seconds on a Hp 48? Impressive.

I got the same number of Emil because I did not count 1

Ah clever you saved already the multiplication results instead of doing the over and over