Subversion Repositories planix.SVN

Rev

Rev 2 | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
2 - 1 
#include <u.h>
2 
#include <libc.h>
3 
#include <bio.h>
4 
 
5 
double big = 9.007199254740992e15;
6 
 
7 
int pt[] =
8 
{
9 
	  2,  3,  5,  7, 11, 13, 17, 19, 23, 29,
10 
	 31, 37, 41, 43, 47, 53, 59, 61, 67, 71,
11 
	 73, 79, 83, 89, 97,101,103,107,109,113,
12 
	127,131,137,139,149,151,157,163,167,173,
13 
	179,181,191,193,197,199,211,223,227,229,
14 
};
15 
double wheel[] =
16 
{
17 
	10, 2, 4, 2, 4, 6, 2, 6, 4, 2,
18 
	 4, 6, 6, 2, 6, 4, 2, 6, 4, 6,
19 
	 8, 4, 2, 4, 2, 4, 8, 6, 4, 6,
20 
	 2, 4, 6, 2, 6, 6, 4, 2, 4, 6,
21 
	 2, 6, 4, 2, 4, 2,10, 2,
22 
};
23 
uchar table[1000];
24 
uchar bittab[] =
25 
{
26 
	1, 2, 4, 8, 16, 32, 64, 128,
27 
};
28 
 
29 
enum {
30 
	ptsiz	= nelem(pt),
31 
	whsiz	= nelem(wheel),
32 
	tabsiz	= nelem(table),
33 
	tsiz8	= tabsiz*8,
34 
};
35 
 
36 
void	mark(double nn, long k);
37 
 
38 
void
39 
usage(void)
40 
{
41 
	fprint(2, "usage: %s [start [finish]]\n", argv0);
42 
	exits("limits");
43 
}
44 
 
45 
void
46 
ouch(void)
47 
{
48 
	fprint(2, "limits exceeded\n");
49 
	exits("limits");
50 
}
51 
 
52 
void
53 
main(int argc, char *argv[])
54 
{
55 
	int i;
56 
	double k, temp, v, limit, nn;
57 
	char *l;
58 
	Biobuf bin;
59 
 
60 
	ARGBEGIN{
61 
	default:
62 
		usage();
63 
		break;
64 
	}ARGEND;
65 
 
66 
	nn = 0;
67 
	limit = big;
68 
	switch (argc) {
69 
	case 0:
70 
		Binit(&bin, 0, OREAD);
71 
		while ((l = Brdline(&bin, '\n')) != nil) {
72 
			if (*l == '\n')
73 
				continue;
74 
			nn = atof(l);
75 
			if(nn < 0)
76 
				sysfatal("negative start");
77 
			break;
78 
		}
79 
		Bterm(&bin);
80 
		break;
81 
	case 2:
82 
		limit = atof(argv[1]);
83 
		if(limit < nn)
84 
			exits(0);
85 
		if(limit > big)
86 
			ouch();
87 
		/* fallthrough */
88 
	case 1:
89 
		nn = atof(argv[0]);
90 
		break;
91 
	default:
92 
		usage();
93 
		break;
94 
	}
95 
 
96 
	if(nn < 0 || nn > big)
97 
		ouch();
98 
	if(nn == 0)
99 
		nn = 1;
100 
 
101 
	if(nn < 230) {
102 
		for(i=0; i<ptsiz; i++) {
103 
			if(pt[i] < nn)
104 
				continue;
105 
			if(pt[i] > limit)
106 
				exits(0);
107 
			print("%d\n", pt[i]);
108 
//			if(limit >= big)
109 
//				exits(0);
110 
		}
111 
		nn = 230;
112 
	}
113 
 
114 
	modf(nn/2, &temp);
115 
	nn = 2*temp + 1;
116 
/*
117 
 *	clear the sieve table.
118 
 */
119 
	for(;;) {
120 
		for(i = 0; i < tabsiz; i++)
121 
			table[i] = 0;
122 
/*
123 
 *	run the sieve.
124 
 */
125 
		v = sqrt(nn+tsiz8);
126 
		mark(nn, 3);
127 
		mark(nn, 5);
128 
		mark(nn, 7);
129 
		for(i = 0, k = 11; k <= v; k += wheel[i]) {
130 
			mark(nn, k);
131 
			i++;
132 
			if(i >= whsiz)
133 
				i = 0;
134 
		}
135 
/*
136 
 *	now get the primes from the table
137 
 *	and print them.
138 
 */
139 
		for(i = 0; i < tsiz8; i += 2) {
140 
			if(table[i>>3] & bittab[i&07])
141 
				continue;
142 
			temp = nn + i;
143 
			if(temp > limit)
144 
				exits(0);
145 
			print("%.0f\n", temp);
146 
//			if(limit >= big)
147 
//				exits(0);
148 
		}
149 
		nn += tsiz8;
150 
	}
151 
}
152 
 
153 
void
154 
mark(double nn, long k)
155 
{
156 
	double t1;
157 
	long j;
158 
 
159 
	modf(nn/k, &t1);
160 
	j = k*t1 - nn;
161 
	if(j < 0)
162 
		j += k;
163 
	for(; j < tsiz8; j += k)
164 
		table[j>>3] |= bittab[j&07];
165 
}