Writing rational numbers as ratio of factorials of primes

Mathematica, 175 177 169 154 108 bytes

Join@@@Table[x[[1]],{s,{1,-1}},{x,r@#},x[[2]]s]&@*(If[#==1,1,{p,e}=Last@(r=FactorInteger)@#;p^e#0[p!^-e#]]&)

Try it online!

How it works

This is the composition of two functions. The first, which ungolfs to

If[# == 1,
  1,
  {p,e} = Last[FactorInteger[#]];
  p^e * #0[p!^-e * #]
]&

is a recursive function for actually computing the desired factorization. Specifically, given a rational input x, we compute the primes whose factorials should be in the numerator and denominator, and return the fraction with all of those primes multiplied together. (For example, on input 10/9 = 2!*5!/(3!*3!*3!), we return 10/27 = 2*5/(3*3*3).)

We do this by dealing with the largest prime factor at every step: if p^e occurs in the factorization of x, we make sure p!^e occurs in the factorial-factorization, and recurse on x divided by p!^e.

(Earlier, I had a more clever strategy that avoids large numbers by looking at the previous prime number before p, but Mathematica can handle numbers as big as 65521! easily, so there's no point. The old version you can find in the history is much faster: on my computer, it took 0.05 seconds on inputs that this version handles in 1.6 seconds.)

The second function turns the output of the first function into lists of primes.

Join @@@ 
  Table[x[[1]],
    {s,{1,-1}},
    {x,FactorInteger[#]},
    x[[2]]*s
  ]&

For s=1 (positive powers) and s=-1 (negative powers), and for each term {prime,exponent} in the factorization r@#, we repeat the prime number prime exponent*s many times.

Noncompeting version with 109 62 bytes

If[#==1,∇1=1,{p,e}=Last@FactorInteger@#;(∇p)^e#0[p!^-e#]]&

Same as above, but instead of giving output as a list, gives output as an expression, using the ∇ operator (because it has no built-in meaning) as a stand-in for factorials. Thus, an input of 10/9 gives an output of (∇2*∇5)/(∇3)^3 to represent (2!*5!)/(3!)^3.

This is shorter because we skip the second part of the function.

---

+2 bytes: the assignment f=First has to be done in the right place to keep Mathematica from getting upset.

-8 bytes: fixed a bug for integer outputs, which actually made the code shorter.

-15 bytes: FactorInteger returns sorted output, which we can take advantage of.

-46 bytes: we don't actually need to be clever.

05AB1E, 54 53 48 46 40 35 33 32 28 bytes

[D¿÷Z#DÓ€gZD<ØŠQ*DˆR!*]¯øεʒĀ

Try it online! Edit: Saved 2 bytes thanks to @ASCII-only. Saved 1 2 3 4 bytes thanks to @Emigna. (I only need to save one more and I'm down to half my original byte count!) Explanation:

[       Begin an infinite loop
D¿÷     Reduce to lowest terms
Z#      Exit the loop if the (largest) value is 1
DÓ€g    Find the index of the largest prime factor of each value
Z       Take the maximum
D<ØŠ    Convert index back to prime and save for later
Q       Convert to an pair of which value had the largest prime factor
*       Convert to an pair with that prime factor and zero
Dˆ      Save the pair in the global array for later
R!*     Multiply the other input value by the factorial of the prime
]       End of infinite loop
¯ø      Collect all the saved primes
εʒĀ     Forget all the saved 0s

Python 2, 220 202 195 183 bytes

g=lambda a,b:a and g(b%a,a)or b;n,d=input();m=c=()
while n+d>2:
 t=n*d;f=p=2
 while t>p:
	if t%p:p+=1;f*=p
	else:t/=p
 if n%p:c+=p,;n*=f
 else:m+=p,;d*=f
 t=g(n,d);n/=t;d/=t
print m,c

Try it online! Edit: Saved 18 25 bytes thanks to @Mr.Xcoder. Saved 12 bytes thanks to @JonathanFrech.