基于素数指数的大数压缩算法性能优化求助

我希望通过素数指数实现数字压缩，采用基于连续素数序列、仅输出指数值来重构原数的方案。编写的Python函数在处理大数时效率极低，耗时可达数小时；扩充素数序列也未改善效果。需要优化现有代码的效率与结构，或提供适配任意大小数字的替代方案。

原代码 #

def get_primes(n):
    primes = []
    num = 2

    while len(primes) < n:
        is_prime = all(num % i != 0 for i in range(2, int(num**0.5) + 1))
        if is_prime:
            primes.append(num)
        num += 1

    return primes

def get_max_prime_amount(num):
    result = num // 3  # Use integer division for max_prime
    return result

def factorize_with_errors(original_number, primes, max_prime):
    error_count = 0
    soriginal_number = original_number
    tCount = 0
    PrimeArray = {}
    ExponentArray = {}

    for tPrimeCount in primes:
        tCount += 1
        PrimeArray[tCount] = tPrimeCount
        ExponentArray[tCount] = 0

    max_prime_number = PrimeArray[tCount]
    tCurrentPrimeCount = 1
    tFactorisation = False

    while not tFactorisation:
            if error_count == 99999999:
                breakpoint()

            tQuotient = gmpy2.mpfr(gmpy2.div(original_number,PrimeArray[tCurrentPrimeCount]))

        if tQuotient.is_integer():
            ExponentArray[tCurrentPrimeCount] += 1

            if tQuotient > max_prime_number:
                if tCurrentPrimeCount == max_prime:
                    tCurrentPrimeCount = 1
                    error_count += 1
                    original_number = int(tQuotient) - 1
                else:
                    original_number = int(tQuotient)
            elif tQuotient <= max_prime_number:
                if tQuotient == 1:
                    tFactorisation = True
                elif tQuotient in primes:
                    for tKey in PrimeArray:
                        if PrimeArray[tKey] == tQuotient:
                            ExponentArray[tKey] += 1
                            tFactorisation = True
                            break
                elif tQuotient not in primes:
                    tCurrentPrimeCount = 1
                    FinalFactorisation = False
                    while not FinalFactorisation:
                        tQuotientNew = tQuotient / PrimeArray[tCurrentPrimeCount]

                        if tQuotientNew == 1:
                            ExponentArray[tCurrentPrimeCount] += 1
                            tQuotient = tQuotientNew
                            tFactorisation = True
                            break
                        elif tQuotientNew.is_integer():
                            ExponentArray[tCurrentPrimeCount] += 1
                            tQuotient = tQuotientNew
                            break
                        else:
                            tCurrentPrimeCount += 1
                break

        elif not tQuotient.is_integer():
            if tCurrentPrimeCount == max_prime:
                tCurrentPrimeCount = 1
                error_count += 1
                soriginal_number -= 1
                original_number = soriginal_number
                for sKey in ExponentArray:
                    ExponentArray[sKey] = 0
            elif tCurrentPrimeCount < max_prime:
                tCurrentPrimeCount += 1

    ExponentArray["error_count"] = error_count
    return ExponentArray


original_number = 288684097887703
input_amount_len = len(str(original_number))
max_prime = get_max_prime_amount(input_amount_len)
primes = get_primes(max_prime)
result = factorize_with_errors(original_number, primes, max_prime)
print(result)

1. 素数生成优化 #

原素数生成函数采用暴力逐个判断，效率极低。替换为埃氏筛法，结合素数定理估算上限，大幅提升素数生成速度：

import math

def get_primes(n):
    if n == 0:
        return []
    # 估算第n个素数的上限，留有余量
    upper_bound = int(n * (math.log(n) + math.log(math.log(n)))) if n > 1 else 3
    sieve = [True] * (upper_bound + 1)
    sieve[0] = sieve[1] = False
    for i in range(2, int(math.isqrt(upper_bound)) + 1):
        if sieve[i]:
            sieve[i*i : upper_bound+1 : i] = [False]*len(sieve[i*i : upper_bound+1 : i])
    primes = [i for i, is_p in enumerate(sieve) if is_p]
    # 若估算上限不足，动态扩展
    while len(primes) < n:
        upper_bound *= 2
        sieve = [True] * (upper_bound + 1)
        sieve[0] = sieve[1] = False
        for p in primes:
            start = max(p*p, ((primes[-1] // p) + 1)*p)
            sieve[start : upper_bound+1 : p] = [False]*len(sieve[start : upper_bound+1 : p])
        primes.extend([i for i in range(primes[-1]+1, upper_bound+1) if sieve[i]])
    return primes[:n]

2. 因式分解逻辑重构 #

原函数逻辑混乱，存在浮点数精度问题与低效查找。以下优化实现基于整数运算，避免修改原数：

import gmpy2

def factorize(num, primes):
    num = gmpy2.mpz(num)
    # 构建素数到索引的映射，实现O(1)查找
    prime_index = {p: idx for idx, p in enumerate(primes, 1)}
    max_idx = len(primes)
    exponents = [0] * (max_idx + 1)  # 索引从1开始对应素数序列

    for idx, p in enumerate(primes, 1):
        if p*p > num:
            break
        while num % p == 0:
            exponents[idx] += 1
            num = num // p
    # 处理剩余的大素数
    if num > 1:
        if num in prime_index:
            exponents[prime_index[num]] += 1
        else:
            raise ValueError(f"剩余数 {num} 不在提供的素数序列中，无法分解")
    # 转换为原输出格式的字典
    exponent_dict = {i: exponents[i] for i in range(1, max_idx+1)}
    exponent_dict["error_count"] = 0
    return exponent_dict

3. 整体流程优化 #

原get_max_prime_amount逻辑缺乏依据，改用素数定理估算所需素数数量，确保覆盖原数所有素因数：

def estimate_required_primes(num):
    num = gmpy2.mpz(num)
    if num == 1:
        return 1
    max_p = num
    # 素数定理反向估算所需素数数量，留有余量
    pi = int(max_p / math.log(max_p)) if max_p > 2 else 1
    return pi + int(math.sqrt(pi))

# 使用示例
original_number = 288684097887703
required_primes_count = estimate_required_primes(original_number)
primes = get_primes(required_primes_count)
result = factorize(original_number, primes)
print(result)

4. 替代方案：适配超大数的因式分解 #

对于超大数，试除法效率不足，可采用Pollard's Rho算法进行快速因式分解：

import random
import gmpy2

def pollards_rho(n):
    if n % 2 == 0:
        return 2
    if n % 3 == 0:
        return 3
    if n % 5 == 0:
        return 5
    while True:
        c = gmpy2.mpz(random.randint(1, n-1))
        f = lambda x: (gmpy2.powmod(x, 2, n) + c) % n
        x, y, d = 2, 2, 1
        while d == 1:
            x = f(x)
            y = f(f(y))
            d = gmpy2.gcd(abs(x-y), n)
        if d != n:
            return d

def factorize_large(n):
    factors = {}
    def _factorize(n):
        if n == 1:
            return
        if gmpy2.is_prime(n):
            factors[n] = factors.get(n, 0) + 1
            return
        d = pollards_rho(n)
        _factorize(d)
        _factorize(n//d)
    _factorize(gmpy2.mpz(n))
    return factors

分解得到因数后，再映射到连续素数的指数序列即可。

内容的提问来源于stack exchange，提问作者Tone