50g Mini-Challenge: Number of positive divisors of x!

+- HP Forums (http://www.hpmuseum.org/forum)
+-- Forum: HP Calculators (and very old HP Computers) (/forum-3.html)
+--- Forum: General Forum (/forum-4.html)
+--- Thread: 50g Mini-Challenge: Number of positive divisors of x! (/thread-9195.html)

50g Mini-Challenge: Number of positive divisors of x! - Joe Horn - 09-29-2017 06:49 PM

The HP 50g program << ! DIVIS SIZE >> returns the number of positive divisors of the input factorial. However, this brute-force program becomes impossibly slow for medium-sized inputs, and runs out of memory for any input that's interestingly large. Just look how this slows down:

<< 12 ! DIVIS SIZE >> returns 792 in 6.1 seconds
<< 13 ! DIVIS SIZE >> returns 1584 in 17.5 seconds
<< 14 ! DIVIS SIZE >> returns 2592 in 39.3 seconds
<< 15 ! DIVIS SIZE >> returns 4032 in 86.2 seconds!

(Unrelated note: Prime's CAS returns 4032 for size(idivis(15!)) in 0.7 seconds)

The winner of this challenge will be the 49G/49g+/50g RPL program that returns the exact number of positive divisors of X! the fastest. Testing will focus on large inputs.

Although half the fun of this challenge will consist in thinking up different ways of doing it (obviously avoiding the factorial and DIVIS functions), some algorithmic hints can be found here if you totally get stuck: https://oeis.org/A027423

Many thanks to Gerald H for his many OEIS-related postings and programs which inspired this mini-challenge.

RE: 50g Mini-Challenge: Number of positive divisors of x! - Thomas Okken - 09-29-2017 08:48 PM

I realize that this doesn’t qualify since it’s a 42S program, but I couldn’t resist:

Code:
00 { 67-Byte Prgm }
01▸LBL "PF"
02 STO 00
03 1
04 STO 01
05 STO 04
06 2
07 GTO 02
08▸LBL 00
09 2
10 STO+ 01
11 1
12 STO 02
13▸LBL 01
14 2
15 STO+ 02
16 RCL 02
17 X↑2
18 RCL 01
19 X<Y?
20 GTO 02
21 LASTX
22 MOD
23 X=0?
24 GTO 00
25 GTO 01
26▸LBL 02
27 STO 03
28 RCL 00
29 X<Y?
30 GTO 04
31 1
32 X<>Y
33▸LBL 03
34 RCL÷ 03
35 IP
36 STO+ ST Y
37 X≠0?
38 GTO 03
39 R↓
40 STO× 04
41 GTO 00
42▸LBL 04
43 RCL 04
44 END

Constant memory use; run time something like O(n*sqrt(n)), I think, but that analysis is harder than writing the program. :-)

RE: 50g Mini-Challenge: Number of positive divisors of x! - Thomas Okken - 09-30-2017 09:01 AM

HP-67: 52 → 8087040000 in 69 s.

(I really should get that card reader restored. It’s been 40 years since the last time I used a calculator that would forget everything when you turned it off!)

UPDATE:

HP-25: 52 → 8087040000 in 54 s. :-)

Code:
01     23 00  STO 0
02        01  1
03     23 01  STO 1
04     23 04  STO 4
05        02  2
06     13 24  GTO 24
07        02  2
08  23 51 01  STO + 1
09        01  1
10     23 02  STO 2
11        02  2
12  23 51 02  STO + 2
13     24 02  RCL 2
14     15 02  x^2
15     24 01  RCL 1
16     14 41  x<y
17     13 24  GTO 24
18     14 73  LASTx
19        71  ÷
20     15 01  FRAC
21     15 71  x=0
22     13 07  GTO 07
23     13 11  GTO 11
24     23 03  STO 3
25     24 00  RCL 0
26     14 41  x<y
27     13 40  GTO 40
28        01  1
29     23 05  STO 5
30        22  R↓
31     24 03  RCL 3
32        71  ÷
33     14 01  INT
34  23 51 05  STO + 5
35     15 61  x≠0
36     13 31  GTO 31
37     24 05  RCL 5
38  23 61 04  STO × 4
39     13 07  GTO 07
40     24 04  RCL 4
41     13 00  GTO 00


RE: 50g Mini-Challenge: Number of positive divisors of x! - Thomas Okken - 09-30-2017 10:33 AM

My algorithm is based on the prime factorization of n!. When a number has prime factorization Π(i=1,j,p[i]^k[i]), i.e. it is the product of j distinct primes and each prime p[i] occurs k[i] times in the factorization, then the number of divisors is Π(i=1,j,k[i]+1).

To find this prime factorization of n!, realize that it can only contain primes that are less than or equal to n, so the search space is pleasantly small. And how often does a prime p occur in n!? Answer: floor(n/p) of the factors 1, 2, 3, ... , n are divisible by p; floor(n/p^2) are divisible by p^2, etc. Keep dividing n by p and rounding toward zero, until you reach zero; the sum of all the intermediate results equals the exponent k[i] of the prime p[i], and the number of possible contributions of p[i] to divisors of n! is thus k[i]+1. Multiply all those k[i]+1 factors together, and you end up with the total number of divisors.

The first part of the program finds primes using a simple but reasonably efficient algorithm, only checking divisors <= sqrt(p). Hence the O(n*sqrt(n)) part of my time estimate. Once a prime has been found, the repeated divisions take O(log(n)). Working this out more accurately would require dusting off my rusty calculus skills... Maybe later. :-)

RE: 50g Mini-Challenge: Number of positive divisors of x! - Gerald H - 09-30-2017 11:29 AM

Thank you, Thomas.

Here a SysRPL programme following Thomas' suggestion which correctly processes 100 in

0.5604s

a definite improvement on my first attempt, and bestfit gives the input time relationship as linear.

Size: 124.

CkSum: # 6E3Ch

Code:
::
  CK1&Dispatch
  # FF
  ::
    ZINT 1
    SWAP
    ZINT 0
    BEGIN
    FPTR2 ^Prime+
    2DUP
    Z>=
    WHILE
    ::
      2DUP
      ZINT 0
      SWAPROT
      BEGIN
      OVER
      FPTR2 ^ZQUOText
      ROTOVER
      FPTR2 ^RADDext
      3UNROLL
      FPTR2 ^DupQIsZero?
      UNTIL
      2DROP
      ZINT 1
      FPTR2 ^RADDext
      4ROLL
      FPTR2 ^RMULText
      3UNROLL
    ;
    REPEAT
    2DROP
  ;
;


RE: 50g Mini-Challenge: Number of positive divisors of x! - Joe Horn - 09-30-2017 03:17 PM

Wow, very impressive! Two quick questions:

(1) FPTR 6 119 (^RADDext) and FPTR 6 118 (^QAdd) both resolve to address 6:5253B, so they are really identical functions. Is there a reason that you prefer using ^RADDext instead of ^QAdd in this program? Same question for the identical functions ^RMULText and ^QMul. (It might help if I knew what the "Q" and the "R" stand for in those function names.)

(2) Your second System RPL program evaluates an input of 100 in half a second, so when I keyed up the same program in User RPL (exactly the same, step-by-step, just in URPL), I expected it to run slower, but it ran MUCH slower than I anticipated. URPL is always slower than SysRPL, but this difference seems crazy:

320 SysRPL --> 1.5 seconds
320 URPL --> 18.6 seconds!

Why is this SysRPL program so much faster? Is it simply because SysRPL is always faster, especially when looping is involved?

Always eager to learn! Thanks in advance for your insights!

RE: 50g Mini-Challenge: Number of positive divisors of x! - Gerald H - 09-30-2017 04:31 PM

The User programme below took

14.9141s

to process 320.

To question 1, ^RMULText etc is just habit, to q 2 I have surely less idea than you.

How great a speed improvement occurs is variable, loops do seem to go much faster in Sys.

Code:
« 1 SWAP 0
  WHILE NEXTPRIME
DUP2 ≥
  REPEAT DUP2 0
OVER 4 ROLL 4 ROLL
    WHILE DUP2
IQUOT DUP 0 >
    REPEAT 5 ROLL +
4 ROLLD PICK3 *
    END 4 DROPN 1 +
4 ROLL * UNROT
  END DROP2
»