在Python中处理大计算的内存使用

Question

我正在用Python做一些计算，但内存不够用了。所以我想读写一个文件来释放内存。我需要一个像非常大的列表对象，因此我想每个对象写一行到文件里，然后就可以读写这些行，而不是直接在内存中操作。行的顺序对我很重要，因为我会用行号作为索引。所以我在想，怎么在Python中替换行，而不影响其他行的位置（其实，移动行也没关系，只要它们能回到我预期的位置就行）。

编辑

我在帮一个朋友，他的Python水平和我差不多。这段代码是用来找出能整除给定非质数的最大质数。这个代码在处理到一百万以下的数字时能正常工作，但一旦超过这个范围，我的内存就不够用了，因为在尝试生成数字列表时。

# a comes from a user input
primes_upper_limit = (a+1) / 2
counter = 3L
numbers = list()
while counter <= primes_upper_limit:
    numbers.append(counter)
    counter += 2L

counter=3
i=0
half = (primes_upper_limit + 1) / 2 - 1
root = primes_upper_limit ** 0.5
while counter < root:
    if numbers[i]:
        j = int((counter*counter - 3) / 2)
        numbers[j] = 0
        while j < half:
            numbers[j] = 0
            j += counter
    i += 1
    counter = 2*i + 3
primes = [2] + [num for num in numbers if num]
for numb in reversed(primes):
    if a % numb == 0:
        print numb
        break

另一个编辑

那如果为每个索引写不同的文件呢？比如创建十亿个文件，文件名是长整型数字，文件里只放一个数字？

Answer 1

对于一个只有十二位数字的数，比如在Project Euler #3中，不需要复杂的整数分解方法，也不需要把中间结果存储到硬盘上。可以用下面这个算法来找出n的因数：

1. 先设定f = 2。
2. 如果n等于1，就停止。
3. 如果f乘以f大于n，打印n并停止。
4. 用f去除n，记下商q和余数r。
5. 如果r等于0，打印q，然后用q去除n，回到第2步。
6. 如果r不等于0，就把f加1，然后回到第3步。

这个方法就是通过每个整数进行试除，直到达到平方根为止，这样就能知道剩下的因数是质数。每找到一个因数就会打印出来。

Answer 2

pytables 是一个非常适合处理和存储大量数据的工具。不过，首先可以参考smci的回答中的建议，尽量减少需要存储的数字数量。

Answer 3

你想找出一个数 a 的最大质因数。（Project Euler 问题 3）你现在的算法和实现步骤如下：

1. 生成一个列表 numbers，包含所有在范围内的候选质数（3 <= n <= sqrt(a)，或者像你现在做的那样用 (a+1)/2）
2. 对 numbers 列表进行筛选，得到一个小于等于 sqrt(a) 的质数列表 {p}
3. 试除法：用每个质数 p 测试 a 是否能被整除。把所有的质因数 {q} 存储下来。
4. 打印所有的因数 {q}；我们只需要最大的一个。

关于这个算法，我有以下几点评论。筛选和试除法在处理大数时并不高效，正如 Owen 和我所提到的。对于很大的 a（比如十亿或万亿），你真的应该使用 NumPy。无论如何，关于实现这个算法的一些建议：

1. 你知道吗，你只需要测试到 √a，也就是 int(math.sqrt(a))，而不是像你现在那样用 (a+1)/2？
2. 其实不需要先生成一个巨大的候选列表 numbers，然后再筛选出质数——这个列表不适合扩展。你可以直接构建质数列表 primes。可以使用 while/for 循环和 xrange(3,sqrt(a)+2,2)（这会给你一个迭代器）。虽然你提到的 xrange() 在 2**31L 时会溢出，但结合平方根的观察，你仍然可以成功分解到 2**62。
3. 一般来说，这种方法不如获取 a 的质因数分解有效，也就是说，每次找到一个质因数 p | a 后，你只需要继续筛选剩下的因数 a/p 或 a/p² 或 a/p³ 等等。除了非常大的质数（或伪质数）这种特殊情况外，这样做会大大减少你需要处理的数字的大小。
4. 另外，你只需要生成一次质数列表 {p}；之后只需存储它并进行查找，而不是每次都重新生成。所以我建议把 generate_primes(a) 和 find_largest_prime_divisor(a) 分开。分解会大有帮助。

这是我对你代码的重写，但在处理十亿以上的数字（a > 10**11 +1）时性能仍然会下降，因为要保持筛选后的列表。我们可以用 collections.deque 来替代列表，以获得更快的 O(1) 的 append() 操作，但这只是一个小优化。

# Prime Factorization by trial division

from math import ceil,sqrt
from collections import deque

# Global list of primes (strictly we should use a class variable not a global)
#primes = deque()
primes = []

def is_prime(n):
    """Test whether n is divisible by any prime known so far"""
    global primes
    for p in primes:
         if n%p == 0:
             return False #  n was divisible by p
    return True # either n is prime, or divisible by some p larger than our list    
def generate_primes(a):
    """Generate sieved list of primes (up to sqrt(a)) as we go"""
    global primes
    primes_upper_limit = int(sqrt(a))
    # We get huge speedup by using xrange() instead of range(), so we have to seed the list with 2
    primes.append(2)
    print "Generating sieved list of primes up to", primes_upper_limit, "...",
    # Consider prime candidates 2,3,5,7... in increasing increments of 2
    #for number in [2] + range(3,primes_upper_limit+2,2):
    for number in xrange(3,primes_upper_limit+2,2):
        if is_prime(number): # use global 'primes'
            #print "Found new prime", number
            primes.append(number) # Found a new prime larger than our list
    print "done"    
def find_largest_prime_factor(x, debug=False):
    """Find all prime factors of x, and return the largest."""
    global primes
    # First we need the list of all primes <= sqrt(x)    
    generate_primes(x)
    to_factor = x # running value of the remaining quantity we need to factor
    largest_prime_factor = None
    for p in primes:
        if debug: print "Testing divisibility by", p
        if to_factor%p != 0:
            continue
        if debug: print "...yes it is"
        largest_prime_factor = p
        # Divide out all factors of p in x (may have multiplicity)
        while to_factor%p == 0:
            to_factor /= p
        # Stop when all factors have been found
        if to_factor==1:
            break
    else:
        print "Tested all primes up to sqrt(a), remaining factor must be a single prime > sqrt(a) :", to_factor
    print "\nLargest prime factor of x is", largest_prime_factor
    return largest_prime_factor