%%HP: T(3)A(R)F(.);
@ by Mark Adler
@ FACNUM (BYTES = 277.5, #74FFh)
@ Given an integer, return its prime factorization as a list.
@ For example, 16353 FACNUM returns { 3 3 23 79 }.
\<<
  @ initialize factor list
  { } SWAP

  @ For this part of the program the stack is: n factorlist, where the two are
  @ kept so that the product of the list times n is the original integer.

  @ factor out 2's
  WHILE DUP 2 MOD NOT REPEAT
    2 / SWAP 2 + SWAP END

  @ factor out 3's
  WHILE DUP 3 MOD NOT REPEAT
    3 / SWAP 3 + SWAP END

  @ start factor search at 5
  5

  @ The stack is now: k n factorlist, where the second two are maintained as
  @ before, and k is the largest 5 mod 6 integer less than or equal to the
  @ last factor found (execpt initially when it is set to 5).

  @ search from k to sqrt(n) for factors---k must be 5 mod 6
  WHILE OVER 1 \=/ REPEAT  @ do while n is not one
    OVER \v/ FLOOR         @ go up to the floor of the square root
    IFERR             @ (divide by zero used as a loop breaker)
      FOR i           @ look at factors that are 1 and 5 mod 6
  IF DUP i MOD NOT THEN
    i 0 / END         @ if 5 mod 6 divides n, then cause error
  IF DUP i 2 + MOD NOT THEN
    i 2 + 0 / END          @ if 1 mod 6 divides n, then cause error
      6 STEP
    THEN DROP              @ got an error---trash the zero
    ELSE DUP               @ end of loop---n divides n
    END                    @ after this, stack is: factor n factorlist
    ROT OVER + ROT ROT          @ add factor to list (factor n factorlist')
    SWAP OVER / SWAP       @ divide out divisor (factor n' factorlist')
    DUP 6 MOD 5 - 2 / +         @ set k' to largest 5 mod 6 <= divisor
  END                 @ find next factor (k' n' factorlist')

  @ return list, dropping n and k
  DROP2
\>>