Hexdigit product pangrams

from collections import Counter
import itertools
import string
import sys
import time

from sympy import factorint

start_time = time.time()
# Lowercase because hex returns lowercase
HEX_DIGITS = "fedcba9876543210"

factorisations = {}

for num_str in itertools.chain(itertools.permutations(HEX_DIGITS, 3),
                               itertools.permutations(HEX_DIGITS, 5)):
    num_str = ''.join(num_str)
    factorisations[num_str] = factorint(int(num_str, 16))

count = 0
last = None
solutions = set()

numbers_seen = 0
total_factorisations_seen = 0

for digits in itertools.permutations(HEX_DIGITS, 8):
    # Print progress because it's fun to watch progress
    if digits[:5] != last:
        print(digits[:5], last and total_factorisations_seen/numbers_seen, flush=True, file=sys.stderr)
        last = digits[:5]
        numbers_seen = 0
        total_factorisations_seen = 0

    a_str = ''.join(digits[:3])
    b_str = ''.join(digits[3:])
    a = int(a_str, 16)
    b = int(b_str, 16)
    n = a*b

    numbers_seen += 1

    a_factors = factorisations[a_str]
    b_factors = factorisations[b_str]

    combined = Counter(a_factors)
    combined.update(b_factors)

    primes = sorted(combined)
    search_stack = [[1, {p:0 for p in primes}]]
    seen = set()

    while search_stack:
        total_factorisations_seen += 1
        num, num_factorisation = search_stack.pop()
        seen.add(tuple(num_factorisation[p] for p in primes))

        if num > a:
            continue

        c_str = hex(num)[2:].zfill(3)
        d_str = hex(n//num)[2:].zfill(5)

        if len(c_str) == 3 and len(d_str) == 5 and len(set(c_str+d_str+''.join(digits))) == 16:
            t = (a,b,num,n//num)
            t = tuple(sorted(t))

            if t not in solutions:
                print("0x{} * 0x{} == 0x{} * 0x{}".format(a_str, b_str, c_str, d_str), flush=True)
                solutions.add(t)
                count += 1

        for i,p in enumerate(primes):
            if num_factorisation[p] < combined[p]:
                new_num_factorisation = num_factorisation.copy()
                new_num_factorisation[p] += 1

                if tuple(new_num_factorisation[p] for p in primes) in seen:
                    continue

                new_num = num*p

                search_stack.append([new_num, new_num_factorisation])

print(count)
print(time.time() - start_time)