C语言求一千亿以前的所有质数(二)

风吹无痕

上一篇文章, 我们讲解了欧拉筛的算法和原理, 那么这篇文章我们就来讲解具体如何用C语言实现欧拉筛来求 1000 亿以前的所有质数.
首先呢, 我们得创建一个数组, 用来保存每个数是否被划去. 我们这里为了节省内存, 用每一位是 0 或 1 来表示是否被划去. 那理论上我们需要
\begin{align} \frac{100000000000}{8}=12500000000=11.6415\,\text{GB} \end{align}
的内存空间, 但是很不幸, 我的编译器不允许申请超过 int 最大值的连续的内存空间. 所以我们用一下二维数组.


定义一个包含 12 个 8 字节数的数组, 该数组的每一个数都代表一个数组的内存地址.


让这地址指向一个 1\,\text{GB} 大小的 long long unsigned int 数组,


先把 long long unsigned int 数组, 的每一位都赋值为 1 , 代表最开始都没有被划去.
一个 long long unsigned int  是 8 字节一共 64 位, 因此每个 long long unsigned int 可以保存一共 64 个数的状态


对于 0\sim63 我们特殊处理一下.
创建了数组之后, 下一步就是开始划掉每一个合数了.


对于每一个在 [2,5\times10^{10}] 的每一个数我们都要遍历, 分别处理其是质数(没有划去)和其是合数(被划去)的情况. 而对于每一个数, 我们都计算其在哪个数组的第几个数的第几位, 存在结构体 locateStruct 中, 保存在 location 中.


getLocation 函数用来计算每一个数的位置.


对于每一个数, 我们都要判断其是否为质数,


isPrime 函数接收对应数的二进制信息存在的 long long unsigned int 和其位置信息, 返回是否为质数.


如果为质数, 我们就划掉所有该质数和其之前质数(包括本身)的乘积.


getNextPrime 接收上一个质数的位置信息, 返回下一个质数




通过位运算求出这个质数, 并且记录下这次的质数的位置信息.


其实应该把这个函数改为接收储存上一个质数的结构体和内存地址数组, 而不是直接接收这么多数字, 可以让代码精简点.
但是因为这个函数是我从求 21 亿质数一直用过来的, 懒得大改了.
如果为合数, 逻辑就要稍微复杂一点, 需要判断是否被质因子整除, 如果整除, 要先划去其和该质因子的乘积再退出.


我们也单独判断一下是否越界, 如果越界就直接 break 了.


OK, 经过上面的大循环, 我们就划去了所有合数, 下一步就是把质数读出来.


用位运算计算出质数


其实这个函数也可以直接接收结构体而不是单独的数来简化代码, 但是我从 21 亿一直用过来的函数就懒得改了.
返回的是一个结构体


最后, 再输出一下最后一个质数和质数的个数


<hr/>





<hr/><hr/>文章封面


蛍ちゃん
<hr/>//  primesEurl
//  Created by Pandora on 2022/11/20.

#include <stdio.h>
#include <stdbool.h>
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <unistd.h>
#include <string.h>

struct previousPrimeStruct{
    long long unsigned int thisPrime;
    long long unsigned int indexLLU;
    unsigned indexBit;
    unsigned indexArray;
};

struct readFromBitsStruct{
    long long unsigned int primeInnerCount;
    bool maxiumReached;
    long long unsigned int theLastPrime;
};

struct locateStruct{
    long long unsigned int indexArray;
    long long unsigned int indexLLU;
    long long unsigned int indexBit;
    long long unsigned int number;
};

typedef struct previousPrimeStruct PPS;
typedef struct readFromBitsStruct RFBS;
typedef struct locateStruct LS;

LS *getLocation(long long unsigned int number,LS *locationS){
    locationS->indexArray=number>>33;
    locationS->indexLLU=(number&((1ULL<<33)-1))>>6;
    locationS->indexBit=number&63ULL;
    locationS->number=number;
    return locationS;
}

PPS *getNexePrime(long long unsigned int *primArray,long long unsigned int indexArray,long long unsigned int indexLLU,unsigned indexBit,PPS *primeStruct){
    indexBit++;
    while(!(*(&primArray[indexArray]+indexLLU)&(1ULL<<(64-indexBit-1)))){
        indexBit++;
        if(indexBit>63){
            indexBit=1;
            indexLLU++;
        }
        if(indexLLU>=(1<<27)){
            indexBit=1;
            indexLLU=0;
            indexArray++;
        }
    }
    long long unsigned int thisPrime=((indexLLU<<6)+indexBit)^(indexArray<<33ULL);
    primeStruct->thisPrime=thisPrime;
    primeStruct->indexLLU=indexLLU;
    primeStruct->indexBit=indexBit;
    return primeStruct;
}

bool isPrime(LS *location,long long unsigned int primeLLU){
    if(primeLLU&(1ULL<<(64-location->indexBit-1))){
        return true;
    }
    return false;
}

RFBS *readFromBits(long long unsigned int bitData,long long unsigned int indexArray,long long unsigned int index,RFBS *readFromBitsS){
    long long unsigned int prime=index<<6;
    long long unsigned int toAnd=1ULL<<63;
    unsigned primeInnerCount=0;
    readFromBitsS->maxiumReached=false;
    for(int innerIndex=0;innerIndex<64;innerIndex++){
        if(((prime+innerIndex)^(indexArray<<33))>100000000000ULL){
            readFromBitsS->maxiumReached=true;
            readFromBitsS->primeInnerCount=primeInnerCount;
            return readFromBitsS;
        }
        if(bitData&toAnd){
            printf("%llu\n",(prime+innerIndex)^(indexArray<<33));
            readFromBitsS->theLastPrime=(prime+innerIndex)^(indexArray<<33);
            primeInnerCount++;
        }
        toAnd=toAnd>>1;
    }
    readFromBitsS->primeInnerCount=primeInnerCount;
    return readFromBitsS;
}

int main(int argc, const char * argv[]) {
    long long unsigned int *primArray[12]={0};
    for(int indexPoint=0;indexPoint<12;indexPoint++){
        primArray[indexPoint]=(long long unsigned int*)malloc(1<<30);
        memset(primArray[indexPoint],255,1<<30);
    }
    *(primArray[0]+0)=3824771065533498388ULL;
    long long unsigned int number=2;
    long long unsigned int prime=2;
    long long unsigned int indexLLU=1;
    long long unsigned int toGetOut=0;
    long long unsigned int primeLLU=0;
    unsigned indexBit=1;
    unsigned indexArray=0;
    PPS *primeStruct=malloc(sizeof(PPS));
    LS *location=malloc(sizeof(LS));
    while(number<50000000001ULL){
        indexArray=0;
        indexLLU=0;
        indexBit=2;
        prime=2;
        location=getLocation(number,location);
        primeLLU=*(primArray[location->indexArray]+location->indexLLU);
        if(isPrime(location,primeLLU)){
            while(prime<=number && number*prime<100000000001ULL){
                toGetOut=number*prime;
                location=getLocation(toGetOut,location);
                *(primArray[location->indexArray]+location->indexLLU)=*(primArray[location->indexArray]+location->indexLLU)&(~(1ULL<<(64-location->indexBit-1)));
               
                primeStruct=getNexePrime(*primArray,indexArray,indexLLU,indexBit,primeStruct);
                prime=primeStruct->thisPrime;
                indexArray=primeStruct->indexArray;
                indexLLU=primeStruct->indexLLU;
                indexBit=primeStruct->indexBit;
            }
        }
        else{
            toGetOut=number*prime;
            location=getLocation(toGetOut,location);
            *(primArray[location->indexArray]+location->indexLLU)=*(primArray[location->indexArray]+location->indexLLU)&(~(1ULL<<(64-location->indexBit-1)));
            while(number%prime && number*prime<100000000001ULL){
                primeStruct=getNexePrime(*primArray,indexArray,indexLLU,indexBit,primeStruct);
                prime=primeStruct->thisPrime;
                indexArray=primeStruct->indexArray;
                indexLLU=primeStruct->indexLLU;
                indexBit=primeStruct->indexBit;
               
                toGetOut=number*prime;
                location=getLocation(toGetOut,location);
                if(location->indexArray>11){
                    break;
                }
                *(primArray[location->indexArray]+location->indexLLU)=*(primArray[location->indexArray]+location->indexLLU)&(~(1ULL<<(64-location->indexBit-1)));
            }
        }
        number++;
    }
    RFBS *readFromBitsS=malloc(sizeof(RFBS));
    readFromBitsS->maxiumReached=false;
    indexLLU=0;
    indexArray=0;
    long long unsigned int primeCount=0;
    while(!(readFromBitsS->maxiumReached)){
        primeLLU=*(primArray[indexArray]+indexLLU);
        readFromBitsS=readFromBits(primeLLU,indexArray,indexLLU,readFromBitsS);
        indexLLU++;
        primeCount=primeCount+readFromBitsS->primeInnerCount;
        if(indexLLU>=(1<<27)){
            indexLLU=0;
            indexArray++;
        }
    }
    printf("THE LAST PRIME IS: %llu\n",readFromBitsS->theLastPrime);
    printf("THE NUMBER OF PRIMES UNDER 100,000,000,000 IS %llu.\n",primeCount);
    return 0;
}