(48 42s) prime factorization

[Dirk E.]

Put a natural number 1 < n < 1.E12 on the stack. (Avoid other values or types.)
The program returns it's prime factorization as an algebraic expression - or the number itself, if it's a prime. EVAL brings back n. Display mode STD.

The Argument 'n' (or Remainder 'r' after factors found) and the next factor 'x' to try are kept directly on the stack.
'c' is compared to 'x' as terminating condition or counts how often a factor occurs.
The result is collected in 'f'.

Program variant PC2 (below) needs less storage, but is slower; another PC3 is between them.
Sizes of .hp objects stored by Emu48: PF1 570, PC2 448, PC3 485 Bytes.

PF1

Code:


« DUP √ IP
  "'" → c f
  « 1
    WHILE OVER 1 >
    REPEAT
      IF DUP 7 <
      THEN DUP 3 <
        1 2 IFTE +
      ELSE
        DUP 7 - 30 MOD -
        WHILE DUP2 MOD
          OVER c ≤ AND
        REPEAT 4 +
          IF DUP2 MOD
          THEN 2 +
          END
          IF DUP2 MOD
          THEN 4 +
          END
          IF DUP2 MOD
          THEN 2 +
          END
          IF DUP2 MOD
          THEN 4 +
          END
          IF DUP2 MOD
          THEN 6 +
          END
          IF DUP2 MOD
          THEN 2 +
          END
          IF DUP2 MOD
          THEN 6 +
          END
        END
      END
      IF DUP c >
      THEN DROP DUP
      END
      IF DUP2 MOD NOT
      THEN 0 'c' STO
        WHILE DUP2 MOD NOT
        REPEAT
          SWAP OVER / SWAP
          'c' INCR DROP
        END f
        IF DUP "'" >
        THEN "*" +
        END OVER →STR +
        IF c 1 >
        THEN
          "^" + c →STR +
        END 'f' STO
        OVER √ IP 'c' STO
      END
    END DROP2
    f "'" + OBJ→
  »
»

PC2

Code:


« DUP √ IP
  "'" → c f
  «
    « +
      IF DUP2 MOD NOT
      THEN 0 'c' STO
        WHILE DUP2 MOD NOT
        REPEAT
          SWAP OVER / SWAP
          'c' INCR DROP
        END f
        IF DUP "'" >
        THEN "*" +
        END OVER +
        IF c 1 >
        THEN "^" + c +
        END 'f' STO
        OVER √ IP 'c' STO
      END
    » → w
    « 1 DUP w EVAL 1 w EVAL
      2 w EVAL 2 w EVAL
      DO 4 w EVAL 2 w EVAL
        4 w EVAL 2 w EVAL
        4 w EVAL 6 w EVAL
        2 w EVAL 6 w EVAL
      UNTIL DUP c >
      END
      IF OVER 1 >
      THEN DUP2 - w EVAL
      END DROP2
      f "'" + OBJ→
    »
  »
»

PC3

Code:


« DUP √ IP
  "'" → c f
  «
    « 0 'c' STO
      WHILE DUP2 MOD NOT
      REPEAT
        SWAP OVER / SWAP
        'c' INCR DROP
      END f
      IF DUP "'" >
      THEN "*" +
      END OVER →STR +
      IF c 1 >
      THEN
        "^" + c →STR +
      END 'f' STO
      OVER √ IP 'c' STO
    » → w
    «
      « +
        IF DUP2 MOD NOT
        THEN w EVAL
        END
      » → t
      « 1 DUP t EVAL 1 t EVAL
        2 t EVAL 2 t EVAL
        DO 4 t EVAL 2 t EVAL
          4 t EVAL 2 t EVAL
          4 t EVAL 6 t EVAL
          2 t EVAL 6 t EVAL
        UNTIL DUP c >
        END
        IF OVER 1 >
        THEN DUP2 - t EVAL
        END DROP2
        f "'" + OBJ→
      »
    »
  »
»

For testing how long a program runs put the argument(s) on the stack, then the program name, then run TTST. (Ideally in the same directory.) Note: TTST resets the Display mode to STD.
E.g. 99999989 ENTER 'PC3' ENTER TTST --> "PC3 0.0037774" (i.e. 37.8 sec).

TTST

Code:


« DUP →STR 
    2
    OVER SIZE 1 - SUB
  TIME → p t
  « EVAL
    TIME t HMS-
      7 FIX →STR STD
    p " " + SWAP +
  »
»

Results with Emu48 set to Authentic Calculator speed (h.mmss(.)sss):

Code:


1: 99999989
"PF1 0.0027127"
"PC3 0.0037774"
"PC2 0.0047385"

1: 986837511503
"PF1 0.0002407"
"PC3 0.0001613"
"PC2 0.0001681"
1: '23*29*31*59*61*89*
149'

1: 966215219970
"PF1 0.0002898"
"PC3 0.0002156"
"PC2 0.0002178"
1: '2*3*5*7*11*13*17*19*
23*61*71'

Accordingly PF1 should run about 4:32 min on 9999999967 or 45:13 min on 999999999989.
If not set to Authentic speed, it runs on a mobile Emu48 57.3 sec on 999999999989 (arm64 8x2GHz with powersaver on, 33 s on a desktop).


On the mobile Free42 takes roughly 3s with 999999999989 runing the following 42s prime factorization program in RPN (normally it shows the remainder in Y: and one factor at a time in X:, continues by R/S; it's finished when Y: shows 1 - using X Y Z T and L):

Code:


00 { 112-Byte Prgm }
01▸LBL "PR5"     L   X    Y    Z    T
02 2                 2    n
03▸LBL 00            x    n
04 RCL ST Y          n    x    n
05 X<>Y              x    n    n
06 MOD           x   n%x  n
07 X≠0?
08 GTO 02
09 X<> ST L          x    n
10 STO÷ ST Y
11 STOP              x    n/
12 GTO 00
13▸LBL 01            x    n
14 RCL ST Y          n    x    n
15 SQRT
16 IP                c    x    n
17 X<>Y              x    c    n
18 RTN           L   X    Y    Z    T
19▸LBL 02        x   n%x  n
20 X<> ST L          x    n
21 3                 x+   x    n
22 X=Y?
23 GTO 03
24 X<>Y              x    x+   n
25 R↓                x+   n
26 GTO 00
27▸LBL 03            x    x    n
28 ÷                 1    n
29 XEQ 01            x    c    n
30 STO ST L      x   x    c    n
31▸LBL 04
32 CLX           x   0    c    n
33 RCL ST Z      x   n    c    n
34 4             x   4    n    c    n
35 RCL+ ST L         x+   n    c    n
36▸LBL 05
37 MOD           x   n%x  c    n
38 X=0?
39 GTO 07        L   X    Y    Z    T
40 CLX           x   0    c    n
41 2             x   2    c    n
42 RCL+ ST L         x+   c    n
43▸LBL 06
44 X>Y?
45 GTO 08
46 RCL ST Z          n    x    c    n
47 X<>Y              x    n    c    n
48 MOD           x   n%x  c    n
49 X≠0?
50 GTO 04
51 R↓            x   c    n
52 X<> ST L          x    n
53 STO÷ ST Y
54 STOP              x    n/
55 XEQ 01            x    c    n
56 GTO 06        L   X    Y    Z    T
57▸LBL 07        x   n%x  c    n
58 R↓            x   c    n
59 X<> ST L          x    n
60 STO÷ ST Y
61 STOP              x    n/
62 XEQ 01            x    c    n
63 RCL ST Z          n    x    c    n
64 X<>Y              x    n    c    n
65 GTO 05
66▸LBL 08
67 1                 1    x    c    n
68 RCL ST T          n    1    x    c
69 END

For Free42 the attached .raw.hp may need to be renamed back to .raw.

(After I had learned about a mod 30 loop I had thought about writing one for 42s, but then I decided to remain with mod 6 since everything else would mean much more typing - the disadvantage of my approach here is that the factors need to be noted manually - or the R/S replaced by PRX e.g. Later I took it as learning exercise for User RPL and moved to '48 with my factorization project; quite challenging, because a subroutine is not 'jumped' to, but loaded onto the stack and then 'evaluated'. Without a subroutine - PF1 - indeed is faster on the long run, but it needs more space in memory.)

Since actually Free42 works with 34 digits, it's possible to check higher than 1.E12 values - running time will grow, of course, by factor 10 for each two more places.

On the 48 the binary mode would also allow 64 bits (18 decimal places), but the condition DUP2 MOD NOT in real arithmetic would become DUP2 DUP2 / * - # 0d == much longer and slower...

Many thanks to all who made Emu48 and Free42 possible
-- happy computing! D. :-)