9base

revived minimalist port of Plan 9 userland to Unix
git clone git://git.suckless.org/9base
Log | Files | Refs | README | LICENSE

primes.c (2126B)


      1 #include	<u.h>
      2 #include	<libc.h>
      3 
      4 #define	ptsiz	(sizeof(pt)/sizeof(pt[0]))
      5 #define	whsiz	(sizeof(wheel)/sizeof(wheel[0]))
      6 #define	tabsiz	(sizeof(table)/sizeof(table[0]))
      7 #define	tsiz8	(tabsiz*8)
      8 
      9 double	big = 9.007199254740992e15;
     10 
     11 int	pt[] =
     12 {
     13 	  2,  3,  5,  7, 11, 13, 17, 19, 23, 29,
     14 	 31, 37, 41, 43, 47, 53, 59, 61, 67, 71,
     15 	 73, 79, 83, 89, 97,101,103,107,109,113,
     16 	127,131,137,139,149,151,157,163,167,173,
     17 	179,181,191,193,197,199,211,223,227,229,
     18 };
     19 double	wheel[] =
     20 {
     21 	10, 2, 4, 2, 4, 6, 2, 6, 4, 2,
     22 	 4, 6, 6, 2, 6, 4, 2, 6, 4, 6,
     23 	 8, 4, 2, 4, 2, 4, 8, 6, 4, 6,
     24 	 2, 4, 6, 2, 6, 6, 4, 2, 4, 6,
     25 	 2, 6, 4, 2, 4, 2,10, 2,
     26 };
     27 uchar	table[1000];
     28 uchar	bittab[] =
     29 {
     30 	1, 2, 4, 8, 16, 32, 64, 128,
     31 };
     32 
     33 void	mark(double nn, long k);
     34 void	ouch(void);
     35 
     36 void
     37 main(int argc, char *argp[])
     38 {
     39 	int i;
     40 	double k, temp, v, limit, nn;
     41 
     42 	if(argc <= 1) {
     43 		fprint(2, "usage: primes starting [ending]\n");
     44 		exits("usage");
     45 	}
     46 	nn = atof(argp[1]);
     47 	limit = big;
     48 	if(argc > 2) {
     49 		limit = atof(argp[2]);
     50 		if(limit < nn)
     51 			exits(0);
     52 		if(limit > big)
     53 			ouch();
     54 	}
     55 	if(nn < 0 || nn > big)
     56 		ouch();
     57 	if(nn == 0)
     58 		nn = 1;
     59 
     60 	if(nn < 230) {
     61 		for(i=0; i<ptsiz; i++) {
     62 			if(pt[i] < nn)
     63 				continue;
     64 			if(pt[i] > limit)
     65 				exits(0);
     66 			print("%d\n", pt[i]);
     67 			if(limit >= big)
     68 				exits(0);
     69 		}
     70 		nn = 230;
     71 	}
     72 
     73 	modf(nn/2, &temp);
     74 	nn = 2.*temp + 1;
     75 /*
     76  *	clear the sieve table.
     77  */
     78 	for(;;) {
     79 		for(i=0; i<tabsiz; i++)
     80 			table[i] = 0;
     81 /*
     82  *	run the sieve.
     83  */
     84 		v = sqrt(nn+tsiz8);
     85 		mark(nn, 3);
     86 		mark(nn, 5);
     87 		mark(nn, 7);
     88 		for(i=0,k=11; k<=v; k+=wheel[i]) {
     89 			mark(nn, k);
     90 			i++;
     91 			if(i >= whsiz)
     92 				i = 0;
     93 		}
     94 /*
     95  *	now get the primes from the table
     96  *	and print them.
     97  */
     98 		for(i=0; i<tsiz8; i+=2) {
     99 			if(table[i>>3] & bittab[i&07])
    100 				continue;
    101 			temp = nn + i;
    102 			if(temp > limit)
    103 				exits(0);
    104 			print("%.0f\n", temp);
    105 			if(limit >= big)
    106 				exits(0);
    107 		}
    108 		nn += tsiz8;
    109 	}
    110 }
    111 
    112 void
    113 mark(double nn, long k)
    114 {
    115 	double t1;
    116 	long j;
    117 
    118 	modf(nn/k, &t1);
    119 	j = k*t1 - nn;
    120 	if(j < 0)
    121 		j += k;
    122 	for(; j<tsiz8; j+=k)
    123 		table[j>>3] |= bittab[j&07];
    124 }
    125 
    126 void
    127 ouch(void)
    128 {
    129 	fprint(2, "limits exceeded\n");
    130 	exits("limits");
    131 }