30个整数映射到1个字节，查表法实现小范围内常数时间素性判定-EW帮帮网

想用查表法实现小范围内的常数时间的素性判定，怎样节省存储空间呢？自然是1个bit代表1个整数的素性，同时，偶数不用存，个位为5的整数不用存，只有个位为1、3、7、9的整数才可能是素数，也就是每20个整数中只有8个数可能是素数，正好1个字节。如果我们定义整数的个数除以bit的个数为压缩率，那么此时的压缩率为20/8=2.5.
但本文的标题是30个整数映射到1个字节，即压缩率可以做到30/8=3.75，怎么做呢？
让我们来考虑小于等于30且与30互素的正整数的个数(Euler Phi函数，或称Euler Totient函数)。
先考虑小于等于30且与30不互素的正整数的个数，30的素因子是2、3、5，所以，30以内所有2、3、5倍数都与30不互素。30以内，2的倍数有15个，3的倍数有10个，5的倍数有6个，6(2和3的最小公倍数)的倍数有5个，10(2和5的最小公倍数)的倍数有3个，15(3和5的最小公倍数)的倍数有2个，30(2、3、5的最小公倍数)的倍数有1个，根据容斥原理，与30不互素的正整数的个数为：15+10+6-5-3-2+1=22，那么与30互素的正整数的个数为30-22=8，这8个数依次为1,7,11,13,17,19,23,29. 当然，可以直接用Euler Phi函数的计算公式，Phi(30)=30*(1-1/2)*(1-1/3)*(1-1/5)=8.
于是，不考虑2、3、5这几个小素数，在所有形如[30k+1, 30(k+1)]的区间内(其中k>=0)，只有30k+(1,7,11,13,17,19,23,29)这8个奇数可能是素数，正好用1个字节存储。
2、3、5这几个小素数在编程时可以单独处理，无伤大雅。
很自然地会想到，压缩率能否进一步提高？当然是可以的，考虑n/Phi(n)何时取得最大值即可。将n分解成其素因子的幂的乘积的形式， $n=p_1^{e_1}p_2^{e_2}\cdots p_k^{e_k}$ ，那么
$Phi(n)=n\left(1-\dfrac{1}{p_1}\right)\left(1-\dfrac{1}{p_2}\right)\cdots \left(1-\dfrac{1}{p_k}\right) \\ \dfrac{n}{Phi(n)}=\dfrac{1}{\left(1-\dfrac{1}{p_1}\right)\left(1-\dfrac{1}{p_2}\right)\cdots \left(1-\dfrac{1}{p_k}\right)}$
要让n/Phi(n)取得最大值，那就要让 $\left(1-\dfrac{1}{p_1}\right)\left(1-\dfrac{1}{p_2}\right)\cdots \left(1-\dfrac{1}{p_k}\right)$ 取得最小值，不过，因为素数有无穷多个，这个连乘积并无最小值，根据黎曼zeta函数的无穷乘积表达式(其中的p取遍所有素数)，

$\zeta(s)=\sum\limits_{n=1}^{\infty} \dfrac{1}{n^s}=\prod\limits_{p}\dfrac{1}{1-\frac{1}{p^s}} \\ \zeta(1)=\sum\limits_{n=1}^{\infty} \dfrac{1}{n}=\prod\limits_{p}\dfrac{1}{1-\frac{1}{p}}=\lim\limits_{n\to +\infty}\ln(n)+\gamma=+\infty$

乘积 $\left(1-\dfrac{1}{p_1}\right)\left(1-\dfrac{1}{p_2}\right)\cdots \left(1-\dfrac{1}{p_k}\right)$ 可以随着素因子越来越多而无限接近于0，所以压缩率n/Phi(n)可以趋于无穷大。

比如

n=2*3*5*7=210,
Phi(n)=210*(1-1/2)*(1-1/3)*(1-1/5)*(1-1/7)=48=8*6,
n/Phi(n)=4.375

n=2*3*5*7*11=2310,
Phi(n)=2310*(1-1/2)*(1-1/3)*(1-1/5)*(1-1/7)*(1-1/11)=480=8*60,
n/Phi(n)=4.8125

n=2*3*5*7*11*13=30030,
Phi(n)=30030*(1-1/2)*(1-1/3)*(1-1/5)*(1-1/7)*(1-1/11)*(1-1/13)=5760=8*720,
n/Phi(n)=5.2135……

于是，不考虑2、3、5、7这几个小素数，在所有形如[210k+1, 210(k+1)]的区间内(其中k>=0)，只有
210k+(
     1, 11, 13, 17, 19, 23, 29, 31,
    37, 41, 43, 47, 53, 59, 61, 67,
    71, 73, 79, 83, 89, 97, 101, 103,
    107, 109, 113, 121, 127, 131, 137, 139,
    143, 149, 151, 157, 163, 167, 169, 173,
    179, 181, 187, 191, 193, 197, 199, 209)
这48个奇数可能是素数，可以用6个字节存储。但是，6个字节不是uint32, uint64之类的基本数据类型的长度，只能使用数组存储。而且，必须依赖SIMD指令才能同时对6个字节进行位与运算，不如uint32, uint64等基本类型方便。所以，

从编程方便以及适用于尽可能多的硬件平台的角度来讲，把30个整数映射到1个字节是最佳的。

从提高压缩率的角度来讲，把210个整数映射到6个字节是最佳的。

因为要想压缩率大于210/Phi(210)=4.375，至少要把2310个整数映射到480/8=60个字节，那么需要AVX512指令集才能同时处理至多64字节的数据，然而绝大多数个人电脑的CPU都不具备AVX512指令集。如果不用AVX512指令并行计算，而是逐次计算，那么速度太慢，不值得。

下面是按照“30个整数映射到1个字节”生成整数素性的字节数组的MATLAB代码，其中N=256，就是用256Byte存储了256×30=7680个整数的素性。N可以随意调大。

clc
N=256;
PrimalityTable=zeros(N,1,'uint8');
for k=1:N
    n=(k-1)*30;
    a=uint8(0);
    if isprime(n+1)
        a=bitor(a,uint8(128));
    end
    if isprime(n+7)
        a=bitor(a,uint8(64));
    end
    %-----------------------
    if isprime(n+11)
        a=bitor(a,uint8(32));
    end
    if isprime(n+13)
        a=bitor(a,uint8(16));
    end
    if isprime(n+17)
        a=bitor(a,uint8(8));
    end
    if isprime(n+19)
        a=bitor(a,uint8(4));
    end  
    %-----------------------
    if isprime(n+23)
        a=bitor(a,uint8(2));
    end
    if isprime(n+29)
        a=bitor(a,uint8(1));
    end
    PrimalityTable(k)=a;
    % dec2bin(a)
end

f = fopen('PrimalityTable.txt','w','n','UTF-8');

n_up=floor(N/16);
for n=1:n_up
    for k=1:16
        fprintf(f,"%3d, ",PrimalityTable((n-1)*16+k));
    end
    fprintf(f,"\r\n");
end
for k=1:(N-n_up*16)
    fprintf(f,"%3d, ",PrimalityTable(n_up*16+k));
end

%以下为正确性检验代码
Map0_29To2n = uint8([0, 128, 0, 0, 0, 0, 0, 64, 0, 0, 0, 32, 0, 16, ...
                     0, 0, 0, 8, 0, 4, 0, 0, 0, 2, 0, 0, 0, 0, 0, 1]);
for n=1:(N-1)*30
    k=floor(n/30);
    r=n-k*30;
    if xor(bitand(PrimalityTable(k+1), Map0_29To2n(r+1)), isprime(n))
        n
    end
end

primes(256)
length(primes(N*30))

数组Map0_29To2n的来历如下，本质就是把1,7,11,13,17,19,23,29这些可能为素数的位置挑选出来，

0	0
1	128
2	0
3	0
4	0
5	0
6	0
7	64
8	0
9	0
10	0
11	32
12	0
13	16
14	0
15	0
16	0
17	8
18	0
19	4
20	0
21	0
22	0
23	2
24	0
25	0
26	0
27	0
28	0
29	1

下面是用C#实现的1~7680内全部整数的常数时间素性判定，超过7680的数，就可以使用试除法或者Miller-Rabin算法进行素性判定。

namespace CSharpIsPrime
{
    internal class Program
    {
        public readonly static byte[] Map0_29To2n = [0, 128, 0, 0, 0, 0, 0, 64, 0, 0, 0, 32, 0, 16, 
                                                     0, 0, 0, 8, 0, 4, 0, 0, 0, 2, 0, 0, 0, 0, 0, 1];
        public readonly static byte[] PrimalityTable = [
            127, 251, 247, 126, 109, 219, 188, 159, 171, 242, 120, 207,  87, 101, 183, 121,
            103,  48, 203, 203, 220, 187, 154, 165,  86, 230,  73, 189,  30, 120, 101, 106,
            106, 199, 181, 180, 123,  84,  50, 170, 155, 197,  15, 249, 192,  42, 133,  31,
            116, 191,  34, 151, 102, 111, 200,  92,  29,  50, 212,  92, 162, 136, 253,  42,
             49, 131,  94, 205,  19,  61,  49, 242, 132,  26, 142, 142, 217, 131, 232, 247,
             42, 105,  88,  16, 167, 193,  49,  98,  78, 223, 117, 166,  73, 241,  26, 225,
             75,  73,  27, 129, 166, 100, 199,   5, 136,  28, 227, 100,  60, 129, 215, 153,
            177, 138,  17, 124,  36, 207, 204, 178,  90, 209,  56, 229,  84,  45,  26,  50,
            114, 100, 111, 152,  65,  59, 193, 195,  52, 143,  28,  64, 173, 179, 179,  64,
             77,  82,  41,  48, 234,  50,  94,  12, 194, 208, 143, 211,  34,  54,  36,  31,
            152, 128, 169,  21,  58, 206,  87, 177,  36, 105, 212,  10, 101,  68, 120,  35,
            139,  18,  96,  43,  92, 244,  46,  57, 224,  86, 160,  17, 253,  22, 168, 116,
              6, 170, 199, 237, 138,  25,  16,  40,  97,  90,  85, 162, 178, 146,  14, 228,
             75, 201, 171,  83, 213,  64, 193, 134, 160,  36, 115, 225,  68,  67, 149, 181,
             24,  49, 178,  30, 139, 145, 104,  13, 234, 227,  70,  69,   3,  44,  36,  74,
            117,  90,   2,  76, 177, 132,  16, 194,  44, 109,  75, 109, 155, 152, 135,   6];

        static void Main(string[] args)
        {
            for (int i = 0; i < 256; i++)
            {
                if (IsPrime(i))
                {
                    Console.Write($"{i}, ");
                }
            }
            Console.WriteLine();

            int PrimeCount = 0;
            for (int i = 0; i < (PrimalityTable.Length * 30 + 1); i++)
            {
                if (IsPrime(i))
                {
                    PrimeCount++;
                }
            }
            Console.WriteLine($"PrimeCount: {PrimeCount}");
        }
        public static bool IsPrime(int n)
        {
            ArgumentOutOfRangeException.ThrowIfLessThan(n, 0, "n must be greater than or equal to 0.");
            bool isPrime = false;
            if (n < (PrimalityTable.Length * 30 + 1))
            {
                switch (n)
                {
                    case 0:
                    case 1:
                        return false;
                    case 2:
                    case 3:
                        return true;
                    case 4:
                        return false;
                    case 5:
                        return true;
                    default:
                        if (n % 2 == 0 || n % 3 == 0 || n % 5 == 0)
                        {
                            return false;
                        }
                        else
                        {
                            int quot = n / 30;
                            int rem = n - quot * 30;  // rem = n % 30
                            if ((PrimalityTable[quot] & Map0_29To2n[rem]) != 0)
                            {
                                return true;
                            }
                            else
                            {
                                return false;
                            }
                        }
                }
            }
            else
            {
                // Miller-Rabin primality test
            }
            return isPrime;
        }
    }
}

C#代码的输出结果如下，与MATLAB的验证代码输出结果一致。

2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 
61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 
131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 
197, 199, 211, 223, 227, 229, 233, 239, 241, 251,

PrimeCount: 973

30个整数映射到1个字节，查表法实现小范围内常数时间素性判定

网站公告

今日签到

热门文章

最新发布

0	0
1	128
2	0
3	0
4	0
5	0
6	0
7	64
8	0
9	0
10	0
11	32
12	0
13	16
14	0
15	0
16	0
17	8
18	0
19	4
20	0
21	0
22	0
23	2
24	0
25	0
26	0
27	0
28	0
29	1

0	0
1	128
2	0
3	0
4	0
5	0
6	0
7	64
8	0
9	0
10	0
11	32
12	0
13	16
14	0
15	0
16	0
17	8
18	0
19	4
20	0
21	0
22	0
23	2
24	0
25	0
26	0
27	0
28	0
29	1

0	0
1	128
2	0
3	0
4	0
5	0
6	0
7	64
8	0
9	0
10	0
11	32
12	0
13	16
14	0
15	0
16	0
17	8
18	0
19	4
20	0
21	0
22	0
23	2
24	0
25	0
26	0
27	0
28	0
29	1