A couple days ago Michael Nielsen posted an image of a one-page paper that gives an algorithm for generating factored random numbers, uniformly distributed from 1 to some designated *N*.

The algorithm does not generate random numbers then factor them. It’s more efficient than that, generating the factorization along with the final result. It does require testing for whether a number is prime, but this is more efficient than factorization.

I thought about trying to code up the algorithm in Python, but then I see that @iconjack beat me to it.

from sympy import isprime from random import random, randint def randfacts(N): while True: n, r, s = N, 1, [] while n > 1: if r > N: break if isprime(n := randint(1,n)): r *= n s.append(n) else: if random() < r/N: return r, s