1: package de.luschny.math.primes;
2:
3: import de.luschny.math.arithmetic.Xint;
4: import de.luschny.math.util.PositiveRange;
5:
6: import java.io.BufferedWriter;
7: import java.io.FileWriter;
8: import java.io.IOException;
9: import java.io.PrintWriter;
10: import java.util.Iterator;
11: import java.util.concurrent.atomic.AtomicInteger;
12: import java.util.logging.Logger;
13:
14: /**
15: * Sieve.
16: *
17: * @author Peter Luschny
18: * @version 2004-10-20
19: */
20:
21: public class PrimeSieve implements IPrimeSieve
22: {
23: int[] primes;
24: PositiveRange sieveRange;
25: PositiveRange primeRange;
26: Logger logger;
27:
28: /**
29: * Constructs a prime sieve for the integer range [1,n].
30: *
31: * @param n
32: * The upper bound of the range to be sieved.
33: */
34:
35: public PrimeSieve(int n)
36: {
37: primes = new int[getPiHighBound(n)];
38: // Note: This forces n>=1
39: sieveRange = new PositiveRange(1, n);
40:
41: int numberOfPrimes = makePrimeList(n);
42: primeRange = new PositiveRange(1, numberOfPrimes);
43: }
44:
45: /**
46: * Enables logging of this component if the logger object is not null.
47: *
48: * @param log
49: * The Logger used to log messages for this component.
50: */
51:
52: public void enableLogging(Logger log)
53: {
54: this.logger = log;
55: }
56:
57: /**
58: * Prime number sieve, Eratosthenes (276-194 b. t.) This implementation
59: * considers only multiples of primes greater than 3, so the smallest value
60: * has to be mapped to 5.
61: * <p/>
62: * Note: There is no multiplication operation in this function.
63: *
64: * @param composite
65: * After execution of the function this boolean array includes
66: * all composite numbers in [5,n] disregarding multiples of 2 and 3.
67: */
68:
69: private static void sieveOfEratosthenes(final boolean[] composite)
70: {
71: int d1 = 8;
72: int d2 = 8;
73: int p1 = 3;
74: int p2 = 7;
75: int s1 = 7;
76: int s2 = 3;
77: int n = 0;
78: int len = composite.length;
79: boolean toggle = false;
80:
81: while (s1 < len) // -- scan sieve
82: {
83: if (!composite[n++]) // -- if a prime is found
84: { // -- cancel its multiples
85: int inc = p1 + p2;
86:
87: for (int k = s1; k < len; k += inc)
88: {
89: composite[k] = true;
90: }
91:
92: for (int k = s1 + s2; k < len; k += inc)
93: {
94: composite[k] = true;
95: }
96: }
97:
98: if (toggle = !toggle) // Never mind, it's ok.
99: {
100: s1 += d2;
101: d1 += 16;
102: p1 += 2;
103: p2 += 2;
104: s2 = p2;
105: }
106: else
107: {
108: s1 += d1;
109: d2 += 8;
110: p1 += 2;
111: p2 += 6;
112: s2 = p1;
113: }
114: }
115: }
116:
117: /**
118: * Transforms the sieve of Eratosthenes into the sequence of prime numbers.
119: *
120: * @param n
121: * Upper bound of the sieve.
122: * @return Number of primes found.
123: */
124:
125: private int makePrimeList(int n)
126: {
127: boolean[] composite = new boolean[n / 3];
128:
129: sieveOfEratosthenes(composite);
130:
131: int[] prime = this.primes; // -- on stack for eff.
132: boolean toggle = false;
133: int p = 5, i = 0, j = 2;
134:
135: prime[0] = 2;
136: prime[1] = 3;
137:
138: while (p <= n)
139: {
140: if (!composite[i++])
141: {
142: prime[j++] = p;
143: }
144: // -- never mind, it's ok.
145: p += (toggle = !toggle) ? 2 : 4;
146: }
147:
148: return j; // number of primes
149: }
150:
151: /**
152: * Get a high bound for pi(n), the number of primes less or equal n.
153: *
154: * @param n
155: * Bound of the primes.
156: * @return A simple estimate of the number of primes <= n.
157: */
158:
159: private static int getPiHighBound(long n)
160: {
161: if (n < 17)
162: {
163: return 6;
164: }
165: return (int) Math.floor(n / (Math.log(n) - 1.5));
166: }
167:
168: /**
169: * Get the n-th prime number.
170: *
171: * @param n
172: * The index of the prime number.
173: * @return The n-th prime number.
174: */
175:
176: public int getNthPrime(int n)
177: {
178: primeRange.containsOrFail(n);
179: return primes[n - 1];
180: }
181:
182: /**
183: * Checks if a given number is prime.
184: *
185: * @param cand
186: * The number to be checked.
187: * @return True if and only if the given number is prime.
188: */
189:
190: public boolean isPrime(int cand)
191: {
192: // The candidate is interpreted as an one point interval!
193: int primeMax = primeRange.getMax();
194: int start = PrimeIteration.indexOf(primes, cand - 1, 0, primeMax);
195: int end = primes[start] == cand ? start + 1 : start;
196: return start < end;
197: }
198:
199: /**
200: * The default iteration of the full sieve.
201: *
202: * @return An iteration over all prime numbers found by the sieve.
203: */
204:
205: public IPrimeIteration getIteration()
206: {
207: return new PrimeIteration(this);
208: }
209:
210: /**
211: * Gives the iteration of the prime numbers in the given interval.
212: *
213: * @param low
214: * Lower bound of the iteration.
215: * @param high
216: * Upper bound of the iteration.
217: * @return An iteration of the prime numbers which are in the interval [low,
218: * high].
219: */
220:
221: public IPrimeIteration getIteration(int low, int high)
222: {
223: return new PrimeIteration(this, new PositiveRange(low, high));
224: }
225:
226: /**
227: * Gives the iteration of the prime numbers in the given range.
228: *
229: * @param range
230: * The range of the iteration.
231: * @return The prime number iteration.
232: */
233:
234: public IPrimeIteration getIteration(PositiveRange range)
235: {
236: return new PrimeIteration(this, range);
237: }
238:
239: /**
240: * Gives the product of the prime numbers in the given interval.
241: *
242: * @param low
243: * Lower bound of the iteration.
244: * @param high
245: * Upper bound of the iteration.
246: * @return The product of the prime numbers which are in the interval [low,
247: * high].
248: */
249:
250: public Xint getPrimorial(int low, int high)
251: {
252: return new PrimeIteration(this, new PositiveRange(low, high)).primorial();
253: }
254:
255: /**
256: * Gives the product of the prime numbers in the given range.
257: *
258: * @param range
259: * The range of the iteration.
260: * @return The product of the prime numbers in the range.
261: */
262:
263: public Xint getPrimorial(PositiveRange range)
264: {
265: return new PrimeIteration(this, range).primorial();
266: }
267:
268: // ----------------------- inner class --------------------------
269:
270: /**
271: * PrimeIteration.
272: *
273: * @author Peter Luschny
274: * @version 2004-09-12
275: */
276:
277: private static class PrimeIteration implements IPrimeIteration
278: {
279: private final PrimeSieve sieve;
280: private final PositiveRange sieveRange;
281: private final PositiveRange primeRange;
282: private final int start;
283: private final int end;
284: private int current;
285: private AtomicInteger state;
286:
287: /**
288: * Constructs the iteration for the passed sieve.
289: *
290: * @param sieve
291: * The sieve which is to be enumerated.
292: */
293:
294: PrimeIteration(final PrimeSieve sieve)
295: {
296: this.sieve = sieve;
297: current = start = 0;
298: end = sieve.primeRange.getMax();
299: sieveRange = sieve.sieveRange;
300: primeRange = sieve.primeRange;
301: state = new AtomicInteger(0);
302: }
303:
304: /**
305: * Constructs an iteration over a subrange of the range of the sieve.
306: *
307: * @param sieve
308: * Prime number sieve to be used.
309: * @param sieveRange
310: * Range of iteration.
311: */
312:
313: PrimeIteration(final PrimeSieve sieve, final PositiveRange sieveRange)
314: {
315: this.sieve = sieve;
316: this.sieveRange = sieveRange;
317:
318: sieve.sieveRange.containsOrFail(sieveRange);
319:
320: int sieveMax = sieve.primeRange.getMax();
321: int primeMin = indexOf(sieve.primes, sieveRange.getMin() - 1, 0, sieveMax - 1);
322: int primeMax = indexOf(sieve.primes, sieveRange.getMax(), primeMin, sieveMax);
323:
324: if (primeMin == primeMax) // there are no primes in this range
325: {
326: start = end = 0;
327: primeRange = new PositiveRange(0, 0);
328: }
329: else
330: {
331: start = primeMin;
332: end = primeMax;
333: primeRange = new PositiveRange(primeMin + 1, primeMax);
334: }
335:
336: current = primeMin;
337: state = new AtomicInteger(0);
338: }
339:
340: /**
341: * Provides an iterator over the current prime number list. This
342: * iterator is thread save.
343: *
344: * @return An iterator over the current prime number list.
345: */
346:
347: public Iterator<Integer> iterator()
348: {
349: PrimeIteration result = this;
350:
351: if (0 != state.getAndIncrement())
352: {
353: result = new PrimeIteration(sieve, sieveRange);
354: }
355:
356: result.current = result.start;
357:
358: return result;
359: }
360:
361: /**
362: * Returns The next prime number in the iteration.
363: *
364: * @return The next prime number in the iteration.
365: */
366:
367: public Integer next()
368: {
369: return sieve.primes[current++];
370: }
371:
372: /**
373: * Checks the current status of the finite iteration.
374: *
375: * @return True iff there are more prime numbers to be enumerated.
376: */
377:
378: public boolean hasNext()
379: {
380: return current < end;
381: }
382:
383: /**
384: * The (optional operation) to remove from the underlying collection the
385: * last element returned by the iterator is not supported.
386: *
387: * @throws UnsupportedOperationException
388: */
389:
390: public void remove()
391: {
392: throw new UnsupportedOperationException();
393: }
394:
395: /**
396: * Identifies the index of a prime number. Uses a (modified!) binary
397: * search.
398: *
399: * @param data
400: * List of prime numbers.
401: * @param value
402: * The prime number given.
403: * @param low
404: * Lower bound for the index.
405: * @param high
406: * Upper bound for the index.
407: * @return The index of the prime number.
408: */
409:
410: static int indexOf(final int[] data, int value, int low, int high)
411: {
412: while (low < high)
413: {
414: // int mid = low + ((high - low) / 2);
415: // Probably faster, and arguably as clear is:
416: // int mid = (low + high) >>> 1;
417: // In C and C++ (where you don't have
418: // the >>> operator), you can do this:
419: // mid = ((unsigned) (low + high)) >> 1;
420:
421: int mid = (low + high) >>> 1;
422:
423: if (data[mid] < value)
424: {
425: low = mid + 1;
426: }
427: else
428: {
429: high = mid;
430: }
431: }
432:
433: if (low >= data.length)
434: {
435: return low;
436: }
437:
438: if (data[low] == value)
439: {
440: low++;
441: }
442:
443: return low;
444: }
445:
446: /**
447: * Computes the number of primes in the iteration range.
448: *
449: * @return Cardinality of primes in iteration range.
450: */
451:
452: public int getNumberOfPrimes()
453: {
454: if (0 == primeRange.getMax()) // If primeRange is empty...
455: {
456: return 0;
457: }
458:
459: return primeRange.size();
460: }
461:
462: /**
463: * Computes the density of primes in the iteration.
464: *
465: * @return Ratio of the number of primes relative to the number of
466: * integers in this iteration.
467: */
468:
469: public double getPrimeDensity()
470: {
471: // Note: By construction sieveRange.size is always != 0.
472: return (double) primeRange.size() / (double) sieveRange.size();
473: }
474:
475: /**
476: * Gives the interval [a,b] of the sieve.
477: *
478: * @return sieved interval.
479: */
480:
481: public PositiveRange getSieveRange()
482: {
483: return (PositiveRange) sieveRange.clone();
484: }
485:
486: /**
487: * Gives the range of the indices of the prime numbers in the
488: * enumeration.
489: *
490: * @return Range of indices.
491: */
492:
493: public PositiveRange getPrimeRange()
494: {
495: return (PositiveRange) primeRange.clone();
496: }
497:
498: /**
499: * Gives the prime numbers in the iteration as an array.
500: *
501: * @return An array of prime numbers representing the iteration.
502: */
503:
504: public int[] toArray()
505: {
506: int primeCard = primeRange.size();
507: int[] primeList = new int[primeCard];
508:
509: System.arraycopy(sieve.primes, start, primeList, 0, primeCard);
510:
511: return primeList;
512: }
513:
514: /**
515: * Computes the product of all primes in the range of this iteration.
516: *
517: * @return The product of all primes in the range of the iteration.
518: */
519:
520: public Xint primorial()
521: {
522: if (0 == primeRange.getMax()) // if primeRange is empty...
523: {
524: return Xint.ONE;
525: }
526: return Xint.product(sieve.primes, start, primeRange.size());
527: }
528:
529: /**
530: * Writes the primes iteration to a file.
531: *
532: * @param fileName
533: * Name of the file to write to.
534: * @param format
535: * The print-format requested.
536: */
537:
538: public void toFile(String fileName, PrintOption format)
539: {
540: try
541: {
542: PrintWriter primeReport = new PrintWriter(new BufferedWriter(new FileWriter(fileName)));
543:
544: switch (format)
545: {
546: case CommaSeparated:
547: writeComma(primeReport);
548: break;
549: case FormattedText:
550: writeFormatted(primeReport);
551: break;
552: case Xml:
553: writeXml(primeReport);
554: break;
555: default:
556: writeComma(primeReport);
557: break;
558: }
559:
560: primeReport.close();
561: }
562: catch (IOException e)
563: {
564: if (null != sieve.logger)
565: {
566: sieve.logger.severe(e.toString());
567: }
568: }
569: }
570:
571: /**
572: * Prints the list of prime numbers in a formatted way to a file.
573: *
574: * @param file
575: * The file to be printed to.
576: */
577:
578: private void writeFormatted(PrintWriter file)
579: {
580: file.println("SieveRange " + sieveRange.toString() + " : " + sieveRange.size());
581: file.println("PrimeRange " + primeRange.toString() + " : " + primeRange.size());
582: file.println("PrimeDensity " + getPrimeDensity());
583:
584: int primeOrdinal = start;
585: int primeLow = sieve.primes[start];
586: int lim = (primeLow / 100) * 100;
587:
588: for (int prime : this)
589: {
590: primeOrdinal++;
591:
592: if (prime >= lim)
593: {
594: lim += 100;
595:
596: file.println();
597: file.print('<');
598: file.print(primeOrdinal);
599: file.print(".> ");
600: }
601:
602: file.print(prime);
603: file.print(' ');
604: }
605:
606: file.println();
607: }
608:
609: /**
610: * Prints the list of primes in a comma separated format to a file.
611: *
612: * @param file
613: * The file to be printed to.
614: */
615:
616: private void writeComma(PrintWriter file)
617: {
618: for (int prime : this)
619: {
620: file.print(prime);
621: file.print(',');
622: }
623: }
624:
625: /**
626: * Prints the list of prime numbers in XML format to a file.
627: *
628: * @param file
629: * The file to be printed to.
630: */
631:
632: private void writeXml(PrintWriter file)
633: {
634: int primeOrdinal = start;
635:
636: file.println("<?xml version=\"1.0\" encoding=\"ISO-8859-1\"?>");
637: file.println("<primes>");
638:
639: for (int prime : this)
640: {
641: primeOrdinal++;
642: file.print("<ol>(");
643: file.print(primeOrdinal);
644: file.print(',');
645: file.print(prime);
646: file.println(")</ol>");
647: }
648:
649: file.println("</primes>");
650: }
651:
652: } // endOfPrimeIteration
653: } // endOfPrimeSieve