Generate prime factorisations of factorial numbers

Question

The code was written to generate the factorisation of numbers into primes in ascending order. E.g 3! ='2 * 3', 4!='2^3 * 3'. When I have huge numbers such as 100!, the run time becomes a problem. I want to have a code which will run in less time.

from math import floor,sqrt,factorial
from decimal import Decimal
def prime(x):
   if x==2 or x==3 or x==5:
      return True
   if (x%2==0) or (x%3==0) or (x%5==0) or (x==1):
     return False
   for i in range(6,floor(sqrt(x)),6):
      if x%(i+1)==0:
         return False
      elif i%(i+5)==0:
         return False
   return True
def decomp(n):
   y=factorial(n)
   if y<7:
      pri =list(filter(prime,range(1,y+1)))
   else:
      pri=list(filter(prime,range(1,floor(Decimal(y).sqrt()+1))))
   a={};b=0
   for i in pri:
      while y%i==0:
         y=y//i
         b+=1
         a[i]=b
         b=0
   if len(list(a.keys()))==1:
      if list(a.values())!=[1]: 

         c=str(a.keys())+'^'+str(a.values())
      else:
         c=str(list(a.keys()))
   elif len(list(a.keys()))>1:
      if a[list(a.keys())[0]]>1:
         c=str(list(a.keys())[0])+'^'+str(a[list(a.keys())[0]])
  
      elif a[list(a.keys())[0]]==1:
         c=str(list(a.keys())[0])
      for key in list(a.keys())[1:]:
         if a[key]!=1:
            c+=' * '+str(key)+'^'+str(a[key])
        
         elif a[key]==1:
            c+=' * '+str(key)
        
     
   return c 

Answer 1

Given that you do y=factorial(n) inside your decomp(n) function, I'm gonna assume that this is really about factoring factorials. And my review is gonna be:

• Don't compute the factorial at all.
• Use a better function name.
• Document what the function does.
• Separate the factorization from the presentation. Return something more useful than a string, like pairs of primes and their exponents. That can then be pretty-formatted if desired, or used for computations.
• Test your code, because for example for decomp(10) returns '2 * 3 * 5 * 7', missing the exponents.

Instead of computing the factorial and then factoring that, factor all the numbers from 1 to n and combine their prime factorizations. Or even better... compute, for every prime from 1 to n, its exponent in n!. Let's say we want to find the exponent for 3 in 100!. Every third number from 1 to 100 (i.e., 3, 6, 9, 12, etc) contributes a factor 3, so there are \$\lfloor\frac{100}{3}\rfloor=33\$ numbers that contribute factor 3. Of those, every third number (i.e., 9, 18, 27, etc) contributes another factor 3. Of those, every third number (i.e., 27, 54 and 81) contributes yet another. And 81 even contributes a fourth. So the exponent for 3 in 100! is 33+11+3+1=48.

Here's an O(n) solution doing that. Takes my laptop about 0.33 seconds to factor 1,000,000!. Compare that to just computing math.factorial(10**6), which already takes about 10 seconds.

def factorial_factorization(n):
    """Return the prime factorization of n! as list of (prime, exponent) pairs."""
    
    # Compute primes from 2 to n, see https://cp-algorithms.com/algebra/prime-sieve-linear.html
    lp = [None] * (n+1)
    primes = []
    for i in range(2, n+1):
        if not lp[i]:
            primes.append(i)
            lp[i] = i
        for p in primes:
            if p > lp[i] or i * p > n:
                break
            lp[i * p] = p

    # Find prime exponents for n!
    result = []
    for p in primes:
        e = 0
        m = n // p
        while m:
            e += m
            m //= p
        result.append((p, e))
    return result

Demo:

>>> factorial_factorization(3)
[(2, 1), (3, 1)]

>>> factorial_factorization(4)
[(2, 3), (3, 1)]

>>> factorial_factorization(100)
[(2, 97), (3, 48), (5, 24), (7, 16), (11, 9), (13, 7), (17, 5), (19, 5), (23, 4), (29, 3), (31, 3), (37, 2), (41, 2), (43, 2), (47, 2), (53, 1), (59, 1), (61, 1), (67, 1), (71, 1), (73, 1), (79, 1), (83, 1), (89, 1), (97, 1)]

>>> factorial_factorization(10**6)[:10]
[(2, 999993), (3, 499993), (5, 249998), (7, 166664), (11, 99998), (13, 83332), (17, 62497), (19, 55553), (23, 45453), (29, 35713)]

If you want a less useful but prettier presentation, that's easy to do then:

>>> n = 42
>>> print(f'{n}! =', ' * '.join(f'{p}^{e}' if e > 1 else f'{p}'
                                for p, e in factorial_factorization(n)))
42! = 2^39 * 3^19 * 5^9 * 7^6 * 11^3 * 13^3 * 17^2 * 19^2 * 23 * 29 * 31 * 37 * 41

Checking that:

>>> s = '2^39 * 3^19 * 5^9 * 7^6 * 11^3 * 13^3 * 17^2 * 19^2 * 23 * 29 * 31 * 37 * 41'
>>> (result := eval(s.replace('^', '**')))
1405006117752879898543142606244511569936384000000000

>>> (expect := math.factorial(42))
1405006117752879898543142606244511569936384000000000

>>> result == expect
True

Answer 2

Your prime decomposition could be much (much) more efficient.

Let's look at what you're doing right now : You create a big number (let's say \$100!\$), then get all the primes between between 1 and \$\sqrt (100!)\$ to then check if they are a factor of your number.

One efficient but not "fun" approach would be to cache a predefined table of prime numbers up to a certain high value and use this cache to determine prime factorization.

At the moment, you effectively check if a number is a prime twice, which is costly. Why not simply do this :

def prime_factorization(n):

    primes = []
    for i in range(2,int(sqrt(n))):
        if n == 1:
            return primes

        while n % i == 0:
            n = n // i
            primes.append(i)

Is it efficient? Not so much, but we only look at every number once. And since we start from

How could it be more efficient? Our big problem is the range(2,sqrt(n)). We'd like to be a little bit more specific as to which number we look at. You could create a prime number generator (using the yield keyword) that would return only primes. There are a billion ways to create efficient prime number generators, so I'll skip on this. But then :

def prime_factorization(n):

    primes = []
    for i in prime_generator(from=2, to=sqrt(n)):
        if n == 1:
            return primes

        while n % i == 0:
            n = n / i
            primes.append(i)

Then you'd have something pretty fast and simple.

Lastly, I would recommend wrapping all your code to print your result in another function, as it would make decomp clearer.