Numpy切片赋值是否等价于C循环？为何埃氏筛C版难快于Python？

在Numpy已支持SIMD优化的前提下，以下场景仍适合用C改写：

• 无法向量化的复杂逻辑：当算法包含大量分支判断、非规则内存访问，无法用Numpy向量化API表达时，Python循环的解释器开销会成为瓶颈，C实现可彻底消除这部分开销。
• 极致内存控制需求：比如需要位级存储、自定义内存对齐方式，或需精细管理内存分配/释放以避免泄漏时，手动C实现比Numpy的通用数组结构更灵活。
• 系统集成与代码复用：已有成熟C代码库需对接Python，或需与其他C/C++系统（如嵌入式设备、高性能计算框架）集成时。
• 自定义并行策略：当算法并行逻辑无法通过Numpy/Scipy自动优化，需手动控制多线程、多进程或更底层并行调度时。

测试结果 #

在不同输入规模下（如n=1e6、1e7），Numpy版本筛法与手动C实现的耗时差距极小，部分场景下Numpy版本甚至略快。

代码实现 #

C语言实现

#include <stdlib.h>
#include <stdbool.h>
#include <math.h>
#include <stdint.h>

typedef struct
{
    int *primes;
    int size;
} Result;

Result EratosthenesSieveC(int n)
{
    uint8_t *nbs = (uint8_t *)malloc((n + 1) / 8 + 1);
    int i, j, count = 0;
    double limit = sqrt(n);

    // Initialize the array with true everywhere except for 0 and 1
    nbs[0] = nbs[1] = false;
    for (i = 2; i <= n; i++)
    {
        nbs[i / 8] |= 1 << (i % 8);
    }

    // Apply the Sieve of Eratosthenes algorithm
    for (i = 2; i <= limit; i++)
    {
        if (nbs[i / 8] & (1 << (i % 8)))
        {
            for (j = i * i; j <= n; j += i)
            {
                if (nbs[j / 8] & (1 << (j % 8)))
                {
                    nbs[j / 8] &= ~(1 << (j % 8));
                    count++; // how many numbers are not prime
                }
            }
        }
    }

    // Allocate memory for the array of primes
    int *primes = (int *)malloc((n + 1 - count - 2) * sizeof(int));

    // Store the prime numbers in the 'primes' array
    int index = 0;
    for (i = 2; i <= n; i++)
    {
        if (nbs[i / 8] & (1 << (i % 8)))
        {
            primes[index++] = i;
        }
    }

    // Free the memory used by the 'nbs' array
    free(nbs);

    Result result;
    result.primes = primes;
    result.size = n - count - 1;

    return result;  // Note: I never really deallocate the memory of result which might be causing memory leaks??
}

Numpy高效实现

import numpy as np

def EratosthenesSieveFullNumpy(n):
    nbs = np.ones(n+1, dtype=bool)
    nbs[:2] = 0
    for i in range(2, int(n**0.5)+1):
        if nbs[i]:
            nbs[i*i::i] = 0
    return np.where(nbs)[0]

Numpy低效实现（反例）

import numpy as np

def EratosthenesSieveNumpy(n):
    nbs = np.ones(n+1, dtype=bool)
    nbs[:2] = 0
    for i in range(2, int(n**0.5)+1):
        if nbs[i]:
            for j in range(i*i, n+1, i):
                nbs[j] = 0
    return np.where(nbs)[0]

Python调用C代码的方式

import ctypes

class Result(ctypes.Structure):
    _fields_ = [('primes', ctypes.POINTER(ctypes.c_int)), ('size', ctypes.c_int)]

lib = ctypes.cdll.LoadLibrary(my_path_of_compiled_shared_library)
lib.EratosthenesSieveC.restype = Result

def EratosthenesSieveC(n, lib):
    result = lib.EratosthenesSieveC(n)
    return result.primes[:result.size]

性能差距的原因分析 #

1. SIMD利用效率差异：Numpy的nbs[i*i::i] = 0是向量化操作，底层用SIMD批量清零连续内存块；而你的C实现采用位存储，每次操作需计算字节索引和位偏移，位运算无法被SIMD有效批量处理，抵消了内存紧凑的优势。
2. 循环内额外开销：C代码中count的统计需要在循环内判断位是否已清零，增加了指令数；Numpy赋值操作直接覆盖，无需额外判断。
3. 编译器优化不足：若编译C代码时未启用最高级优化（如-O3 -mavx2），编译器不会自动生成SIMD指令，进一步拉大差距。

C代码优化建议 #

• 改用uint64_t数组存储，利用SIMD批量清零64位块，减少位操作开销。
• 移除循环内的count统计，最后遍历数组统计质数数量，简化循环逻辑。
• 编译时启用编译器优化：添加-O3 -mavx2参数，让编译器自动生成SIMD指令。
• 确保数组内存对齐，提升缓存命中率。

内容的提问来源于stack exchange，提问作者FluidMechanics Potential Flows