prime number calculation

Sieve primes upto 2^32

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

class pbuf
{
    public:
	int sizbit;
	unsigned int siz0;
	unsigned char *buf;
	unsigned long long starti;
	pbuf(int bit)
	{
	    sizbit=bit;
	    siz0=1<<(bit-3-1);
	    buf=(unsigned char *)malloc(siz0);
	    printf("# pbuf %.1f alloc\n", siz0*1.0/1e6);
	    clr();
	}

	void save(char *fn)
	{
	    FILE *fp=fopen(fn,"w");
	    fwrite(&sizbit, 1, 4,fp);
	    fwrite(&starti, 1, sizeof(unsigned long long),fp);
	    fwrite(buf, 1, siz0, fp);
	    fclose(fp);
	}

	pbuf(char *fn)
	{
	    FILE *fp=fopen(fn,"r");
	    fread(&sizbit, 1, 4,fp);
	    fread(&starti, 1, sizeof(unsigned long long),fp);
	    siz0=1<<(sizbit-3-1);
	    fread(buf, 1, siz0, fp);
	    fclose(fp);
	}

	void clr()
	{
	    unsigned int i;
	    for(i=0;i<siz0;i++)buf[i]=0xff;
	}

	void set_start(unsigned long long sti)
	{
	    starti=sti;
	    clr();
	}

	int get_flag(unsigned long long ii)
	{
	    unsigned long long i0 = ii -1 - starti;
	    i0/=2;
	    if((buf[i0/8] & (1<<(i0&7)) )!=0) return 1;
	    else return 0;
	}

	void unset_flag(unsigned long long ii)
	{
	    unsigned long long i0 = ii -1 - starti;
	    i0/=2;
	    buf[i0/8] &= 255- (1<<(i0&7));
	}
};

int main(int argc, char **argv)
{
    unsigned long long i,j;
    pbuf p32(32);
    p32.set_start(0);

    int imax=1<<16+1;
    unsigned long long jmax=(1LL<<32)-1;

    for(i=3;i<imax;i+=2)
    {
	if(p32.get_flag(i)==0)continue;

	for(j=i*3LL;j<=jmax;j+=i*2)
	{
	    p32.unset_flag(j);
	}

    }//i

    unsigned int nprime=0;
    unsigned long long pmax=0;
    for(i=3;i<=jmax;i+=2)
    {
	if(p32.get_flag(i)==0)continue;
	nprime++;
	pmax=i;
    }
    printf("Max prime=%lld nprime=%d\n", pmax, nprime+1);
    p32.save("prime32.dat");
}

Sieve primes over 2**32

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
//#define DEB
class pbuf
{
    public:
	int sizbit;
	unsigned int siz0;
	unsigned char *buf;
	unsigned long long starti;
	pbuf(int bit)
	{
	    sizbit=bit;
	    siz0=1<<(bit-3-1);
	    buf=(unsigned char *)malloc(siz0);
	    printf("# pbuf %.1f alloc\n", siz0*1.0/1e6);
	    clr();
	}

	void save(char *fn)
	{
	    FILE *fp=fopen(fn,"w");
	    fwrite(&sizbit, 1, 4,fp);
	    fwrite(&starti, 1, sizeof(unsigned long long),fp);
	    fwrite(buf, 1, siz0, fp);
	    fclose(fp);
	}

	pbuf(char *fn)
	{
	    FILE *fp=fopen(fn,"r");
	    fread(&sizbit, 1, 4,fp);
	    fread(&starti, 1, sizeof(unsigned long long),fp);
	    siz0=1<<(sizbit-3-1);
	    buf=(unsigned char *)malloc(siz0);
	    fread(buf, 1, siz0, fp);
	    fclose(fp);
	}

	void clr()
	{
	    unsigned int i;
	    for(i=0;i<siz0;i++)buf[i]=0xff;
	}

	void set_start(unsigned long long sti)
	{
	    if((sti&1)==1)
	    {
		printf("set start must be even\n");
		exit(1);
	    }
	    starti=sti;
	    clr();
	}

	int get_flag(unsigned long long ii)
	{
	    unsigned long long i0 = ii -1 - starti;
	    i0/=2;

#ifdef DEB
	    if((ii&1)==0)
	    {
		printf("get flag must be odd\n");
		exit(1);
	    }
	    if(ii<starti)
	    {
		printf("get flag range <0\n");
		exit(1);
	    }
	    if(i0/8 >= siz0)
	    {
		printf("get flag range siz0\n");
		exit(1);
	    }
#endif
	    if((buf[i0/8] & (1<<(i0&7)) )!=0) return 1;
	    else return 0;
	}

	void unset_flag(unsigned long long ii)
	{
	    unsigned long long i0 = ii -1 - starti;
	    int i1=i0/2;
#ifdef DEB
	    if((ii&1)==0)
	    {
		printf("unset flag must be odd\n");
		exit(1);
	    }
	    if(ii<starti)
	    {
		printf("unset flag range <0\n");
		exit(1);
	    }
	    if(i0/8 >= siz0)
	    {
		printf("unset flag range >siz0\n");
		exit(1);
	    }
#endif
	    buf[i1/8] &= 255- (1<<(i1&7));
	}

};

#define CBIT 21
#define CSIZE (1LL<<CBIT)
int main(int argc, char **argv)
{
    unsigned long long i,j;
    pbuf p32("prime32.dat");


#if 0
    unsigned long long pmax=0;
    unsigned long long jmax=(1LL<<32) -1;
    printf("%lld\n", jmax);
    for(i=3;i<=jmax;i+=2)
    {
	if(p32.get_flag(i)==0)continue;
	nprime++;
	pmax=i;
    }
    printf("Max prime=%lld nprime=%d\n", pmax, nprime+1);
    exit(1);
#endif

    unsigned long long nprime= 203280221 ;

    unsigned long long istart=1LL<<32;
    unsigned long long iend=istart+CSIZE-1;
    pbuf pblock(CBIT);


    unsigned long long pmax,npmax;
    ////////
    while(1)
    {
	pblock.set_start(istart);

	unsigned long long imax0=1LL<<15;
	while(imax0*imax0<iend) imax0++;

	unsigned long long istart1;
	for(i=3;i<=imax0;i+=2)
	{
	    if(p32.get_flag(i)==0)continue;

	    istart1 = istart/(2*i);
	    istart1 *= 2*i;
	    istart1 -=i;
	    while(istart1<istart)istart1+=i*2;

	    for(j=istart1;j<=iend;j+=i*2)
	    {
		pblock.unset_flag(j);
	    }
	}

	for(i=istart+1;i<=iend;i+=2)
	{
	    if(pblock.get_flag(i))
	    {
		if(nprime>= (1LL<<32))exit(1);
		if((nprime&1)==1 && p32.get_flag(nprime)==1)
		{
		    pmax=i;
		    npmax=nprime;
		}
		nprime++;
	    }
	}
	printf("%.5f  P= %lld Np= %lld\n", istart*100.0/(1LL<<32), pmax, npmax);
	istart += CSIZE;
	iend += CSIZE;
    }
}