针对超大n（10^18及以上）的特定约束特殊半素数高效迭代方案优化问询

我正在做一个和素数计数相关的项目，需要遍历一类特殊的半素数。给定一个n，我要遍历满足以下所有条件的半素数k=pq：

• k ≤ n^(2/3)
• p、q都是素数，且p < q
• q ≤ √(n/p)
• p ≤ ∛n

关键是要能在常数时间内重构k、p、q。我得强调一下，我不是要统计这类数的数量，而是要遍历/生成它们来做后续处理。

我能想到的最佳方法是埃拉托斯特尼筛法的变种。因为目标n非常大，大概在10^18甚至更大（之后会转成C++实现），所以必须用分段筛法，这样空间占用也更友好。

用13900K处理器跑Java代码的话，n=10^15时耗时57秒，找到了144611615个符合条件的k。代码如下：

import java.util.*;
public class Result {
    static final int SEGMENT_SIZE = 1 << 16;

    public static void main(String[] args) {
        long start = System.nanoTime();
        long n = (long) 1e15;
        System.out.println(S(n));
        long end = System.nanoTime();
        System.out.println("Iteration complete in: " + (end - start) / (1000000L) + " ms");
    }

    static long S(long n) {
        long limit = n;
        long[] basePrimes = SoE((long) Math.cbrt(limit));
        long segstart = 0;
        long seglim = (long) Math.pow(n, 2.0/3.0);
        long[][] sieve = new long[2][SEGMENT_SIZE];
        long count = 0;
        long low = segstart;

        while (low <= seglim) {
            long high = Math.min(low + SEGMENT_SIZE - 1, seglim);
            long segLen = (int) (high - low + 1);
            Arrays.fill(sieve[0], 1L);
            Arrays.fill(sieve[1], 0);
            long localsqrt = (long) Math.sqrt(high);

            for (long p : basePrimes) {
                if (p > localsqrt) break;
                // mark multiples of p
                long start = Math.max(((low + p - 1) / p) * p, p * p);
                for (long m = start; m <= high; m += p) {
                    int idx = (int) (m - low);
                    sieve[0][idx] *= p;
                    sieve[1][idx] += 1;
                }
                // mark multiples of p^2 as non-squarefree
                long p2 = p * p;
                long start2 = ((low + p2 - 1) / p2) * p2;
                for (long m = start2; m <= high; m += p2) {
                    int idx = (int) (m - low);
                    sieve[0][idx] = 0;
                }
            }

            // scan segment for semiprimes
            for (int i = 0; i < segLen; i++) {
                long k = low + i;
                if (sieve[1][i] != 1 || sieve[0][i] <= 1) continue;
                // squarefree and exactly one small prime
                long p = sieve[0][i];
                long q = k / p;
                if (q > Math.sqrt(limit / p)) continue;
                count++;
            }
            low = high + 1;
        }
        return count;
    }

    public static long[] SoE(long n) {
        if (n < 2) return new long[0];
        boolean[] isPrime = new boolean[(int) n + 1];
        Arrays.fill(isPrime, true);
        isPrime[0] = isPrime[1] = false;
        for (int i = 2; i * i <= n; i++) {
            if (!isPrime[i]) continue;
            for (int j = i * i; j <= n; j += i) {
                isPrime[j] = false;
            }
        }
        int count = 0;
        for (int i = 2; i <= n; i++) if (isPrime[i]) count++;
        long[] primes = new long[count];
        int idx = 0;
        for (int i = 2; i <= n; i++) {
            if (isPrime[i]) primes[idx++] = i;
        }
        return primes;
    }
}

代码说明 #

我用了一个两行的二维数组来实现分段筛：

• 第一行初始化为1，用来记录对应整数的素因子乘积
• 第二行初始化为0，用来记录对应整数的素因子个数

对于每个≤∛n的素数p：

1. 从p²开始标记p的倍数，更新筛子的两行数据（乘积乘p，计数加1）
2. 标记p²的倍数，把第一行设为0，因为这类数包含平方因子，不符合我们的要求

处理完所有基础素数后，我们筛选符合条件的k：

• 如果第一行是0：直接跳过（有平方因子）
• 如果第一行是1：说明这个数的最小素因子大于∛n，不符合要求
• 如果第二行计数>2：说明至少有3个素因子，不是半素数，跳过

剩下的情况是第二行计数为1或2：

• 计数为2时，通过数学推导可以知道，此时这个数必然还有额外的素因子，不可能是符合要求的半素数，所以不用考虑
• 最终只需要看计数为1的情况：此时第一行的数值就是素因子p，k/p得到q，再检查q是否满足q ≤ √(n/p)的约束，符合的话就是我们要找的半素数

另外，素数本身不会被标记，因为我们是从p²开始标记倍数的，所以素数的第二行计数为0，会被过滤掉，这正好符合我们的需求。

内容来源于stack exchange