怎么样求N以下所有素数

fly

GMP有没有自带这么一个函数：给定一个整数N，输出N以下所有的素数N[x]。
matlab里面就有类似的函数
如果GMP没有的话，这个函数要怎么写？

其实我的目的是想随机得到N以下的一个素数，如果有上面的函数的话，我只要从N[x]里面选一个就行了。

关于随机得到N以下的一个素数，我自己是这么写的
        gmp_randinit_default(state);                   设置随机初始状态
        do
        {
                mpz_urandomm(r, state, N);            得到N以下随机一个数，赋给r
                mpz_nextprime(r, r);                        求得r的下一次素数，赋给r
        }
        while(mpz_cmp(r, N)>=0);                        如果r>N，返回循环重新生成

如果只要得到一个随机素数的话，感觉还可以。不过我是用到很多个随机素数，所以感觉每次都要这么生成，效率不高。

求牛人帮我解答问题！

tprime

//计算2^31 以内素数个数, 打印每个结果需自己写点代码
//
#include <math.h>
#include <stdio.h>
#include <memory.h>
#include <time.h>

#define ST             7
#define MOVE           4
#define BLOCKSIZE      (510510 * 4) // 2 * 3 * 5 * 7 * 11 * 13 * 17 = 510510
#define MAXM           ((BLOCKSIZE >> MOVE) + 1)
#define SET_BIT(a, n)  a[n >> 3] |= (1 << (n & 7))

typedef unsigned int uint;
typedef unsigned char uchar;
static uint Prime[6800];
static uchar NumBit0[1 << 8];
static uchar BaseTpl[MAXM * 2];
static const uint MAXN = (-1u / 2);

static int sieve( )
{
        int primes = 1;
        const uint maxp = (1 << 16) + 1000;
        for (uint p = 3; p < maxp; p += 2) {
                if (0 == (BaseTpl[p >> MOVE] & (1 << (p / 2 & 7)))) {
                        Prime[primes++] = p;
                        for (uint i = p * p / 2; i < maxp / 2; i += p)
                                SET_BIT(BaseTpl, i);
                }
        }
        return primes;
}

static void initBit( )
{
        memset(BaseTpl, 0, sizeof(BaseTpl));
        for (int k = 1; k < ST; k++)
        for (int p = Prime[k] / 2; p < sizeof(BaseTpl) * 8; p += Prime[k])
                SET_BIT(BaseTpl, p);

        for (int i = 1; i < sizeof(NumBit0) / sizeof(NumBit0[0]); i++)
                NumBit0[i] = NumBit0[i >> 1] + (i & 1);
        for (int j = 0; j < sizeof(NumBit0) / sizeof(NumBit0[0]); j++)
                NumBit0[j] = 8 - NumBit0[j];
}

static void inline
crossOutFactorMod6(uchar bitarray[], const uint start,
                const uint leng, uint factor)
{
        uint s2, s1 = factor - start % factor;
        s1 += factor - factor * (s1 % 2);
        if (start <= factor)
                s1 = factor * factor;

        const char skip[] = {0, 2, 1, 1, 0, 1};
        const char mrid = ((start + s1) / factor) % 6 - 1;
        s1 = s1 / 2 + skip[mrid] * factor;
        s2 = s1 + skip[mrid + 1] * factor;

        for (factor += 2 * factor; s2 <= leng / 2; ) {
                SET_BIT(bitarray, s1); s1 += factor; //6k + 1
                SET_BIT(bitarray, s2); s2 += factor; //6k + 5
        }
        if (s1 <= leng / 2)
                SET_BIT(bitarray, s1);;
}

static int inline setSieveTpl(uchar bitarray[], uint start, int leng)
{
        const int offset = start % BLOCKSIZE;
        memcpy(bitarray, BaseTpl + (offset >> MOVE), (leng >> MOVE) + 2);
        if (start < 32)
                *((short*)bitarray) = 0x3461; //0x 0011 0100 0110 0001

        bitarray[0] |= (1 << (offset % 16 / 2)) - 1;
        leng += offset % 16;
        *((uint*)bitarray + (leng >> 6)) |= ~((1u << (leng / 2 & 31)) - 1);
        return offset % 16;
}

static int segmentedSieve(uint start, int leng)
{
        uchar bitarray[MAXM + 64];
        const int offset = setSieveTpl(bitarray, start, leng);
        start -= offset, leng += offset;

        const int maxp = (int)sqrt((double)start + leng) + 1;
        for (int i = ST, p = Prime[i]; p < maxp; p = Prime[++i])
                crossOutFactorMod6(bitarray, start, leng, p);

        int primes = 0;
        for (int k = 0; k <= leng >> MOVE; k++)
                primes += NumBit0[bitarray[k]];
        return primes;
}

int countPrimes(uint start, uint stop)
{
        int primes = 0;
        if (start <= 2 && stop >= 2)
                primes = 1;

        const int blocksize = 63 << 14;
        if (start + blocksize > stop)
                return primes + segmentedSieve(start, stop - start + 1);

        primes += segmentedSieve(start, blocksize - start % blocksize);
        for (uint i = start / blocksize + 1; i < stop / blocksize; i++)
                primes += segmentedSieve(i * blocksize, blocksize);

        primes += segmentedSieve(stop - stop % blocksize, stop % blocksize + 1);
        return primes;
}

int main(int argc, char* argv[])
{
        uint start = 0, stop = MAXN;
        clock_t tstart = clock();

        sieve();
        initBit();

        int primeCnt = countPrimes(0, MAXN);
        const char *ouput = "PI(%u, %u) = %d, time use %u ms\n";
        printf(ouput, start, stop, primeCnt, clock() - tstart);

        while (scanf("%u %u", &start, &stop) == 2 && stop >= start) {
                tstart = clock();
                primeCnt = countPrimes(start, stop);
                printf(ouput, start, stop, primeCnt, clock() - tstart);
        }

        return 0;
}
//PI[1, 2147483647] = 105097565, time use 1592 ms
//

fly

这个程序能生成2^2048内的素数吗？

我发现你没有用到GMP大数库啊！

g99

LZ头像好热阿

liangbch

题目是求N以下的所有素数。3#又问"这个程序能生成$2^2048$内的素数吗?", 如果我没有理解错的话，你是想要$2^2048$以内的所有素数，这实在是一个"天才"的的想法，让我来算一下，需要多少硬盘才能存下这么多素数。
1. 以最朴素的存储方法，存储1个2048bit的素数需要2048/8=256B
2. 根据素数定理可知，$2^2048$以内素数个数为 $2^2048/ln(2^2048),$约等于$2^2038$个素数。需要$2^2038$×256B=$2^2046$B=$2^2006$TB.
3. 假设硬盘的容量为2TB(当前主流硬盘容量),并假设技术足够进步，使得这样的硬盘的质量仅为1个电子的质量(约为$10^-30$千克),则所需的硬盘总质量为 $2^2006$TB/2TB×$10^-30$千克=$3.67×10^573$千克。

4. $3.67×10^573$千克 是个什么概念呢？我们知道地球的质量约为$6×10^24$千克。太阳的质量约为$2×10^30$千克，宇宙$1.53×10^54$千克(来源http://wenwen.soso.com/z/q136056252.htm)，那个$3.67×10^573$千克比宇宙的质量实在是大的太多太多太多了。

liangbch

退一步，即使求$2^64$以内的所有素数，所需的存储容量也相当可观。
$2^64$以内的所有素数大约有$4×10^17$个，按每个素数需要8字节，需要$4×10^17×8=3.2×10^18B=3.2×10^6TB$,以每块硬盘2T，每块硬盘重量0.5公斤计算，硬盘的总重量仍达$3.2×10^6"TB"//(2"TB")xx0.5KG=800$吨。

tprime

3# fly

  在我看来2^60 已经是天文数字， 楼主对数字大小还没有概念， 还是认为
现在的计算机是以光速运行。 即便是能生成2^60以内素数，
你要这么大范围数以内素数do what？循序扫描一遍什么也不干也，要等个十年八年的。

fly

唉呀，我本来就是要找很大的素数，达到密码上要求的安全。
2^60从密码上来说太小了，微不足道，不够安全，而且如果是2^60随便找个软件就行了，要不然还用GMP干嘛！
5楼6楼的东西我看得头晕，本来我的N是2^2048的，怎么到你那么变成了22048了，下次写的东西都不知你写什么，宇宙有1.53×1054千克（＝1吨）那么轻吗？你忽悠人啊！

还有啊，麻烦看一下一楼的内容，不要只是看标题嘛。

fly

如果只要得到一个随机素数的话，感觉还可以。不过我是用到很多个随机素数，所以感觉每次都要这么生成，效率不高


大哥，这才是我真正想要滴，不能被理解，好郁闷啊

liangbch

8# fly

严重怀疑你的浏览器不能显示Latext公式，请double check.