add sundaram algorithm

bitset.go

@@ -51,10 +51,15 @@ func (s BitSet) byteFor(i uint64) (b *byte, mask byte) {
 	return
 }
 
-// NewBitSet returns a new BitSet big enough to hold at least l values
-func NewBitSet(l uint64) BitSet {
-	numBlocks := uint64(math.Floor(float64(l)/Byte)) + 1
-	return make(BitSet, int(numBlocks))
+// BitSetSize returns the number of bytes required for the BitSet to contain at
+// least n values
+func BitSetSize(n uint64) uint64 {
+	return uint64(math.Floor(float64(n)/Byte)) + 1
+}
+
+// NewBitSet returns a new BitSet big enough to hold at least n values
+func NewBitSet(n uint64) BitSet {
+	return make(BitSet, int(BitSetSize(n)))
 }
 
 // Max returns the highest value that can be set

cmd/primes/main.go

@@ -33,7 +33,7 @@ func init() {
 		},
 		cli.StringFlag{
 			Name:  "algorithm, a",
-			Usage: "which algorithm to use [eratosthenes]",
+			Usage: "which algorithm to use [eratosthenes, sundaram]",
 		},
 		cli.StringFlag{
 			Name:  "profile",
@@ -61,6 +61,8 @@ func before(c *cli.Context) error {
 	switch algopt {
 	case "eratosthenes", "":
 		algo = primes.EratosthenesAlgo
+	case "sundaram":
+		algo = primes.SundaramAlgo
 	default:
 		cli.ShowAppHelp(c)
 		log.Fatalf("unknown algorithm: %s\n", algopt)

eratosthenes.go

@@ -1,17 +1,17 @@
 package primes
 
-// Eratosthenes is a bitset that sets bits at each index to 1 if the index
-// number itself is a prime
+// Eratosthenes is a Sieve that is calculated using the sieve of eratosthenes
+// algorithm
 type Eratosthenes struct {
 	BitSet
 }
 
 // NewEratosthenes returns a new Eratosthenes calculated for all values from 0
-// to the nearest byte greater than l
-func NewEratosthenes(l uint64) Sieve {
+// to the nearest byte greater than n
+func NewEratosthenes(n uint64) Sieve {
 	// initialize the sieve with all bits set but unset 0 and 1 (as they are
 	// not-prime)
-	s := Eratosthenes{NewBitSet(l).SetAll().Unset(0).Unset(1)}
+	s := Eratosthenes{NewBitSet(n).SetAll().Unset(0).Unset(1)}
 
 	for i := uint64(2); i <= sqrt(s.Max()); i++ {
 		if !s.IsPrime(i) {

eratosthenes_test.go

@@ -1,7 +1,6 @@
 package primes
 
 import (
-	"math"
 	"testing"
 
 	. "github.com/smartystreets/goconvey/convey"
@@ -29,22 +28,3 @@ func TestEratosthenes(t *testing.T) {
 		}
 	})
 }
-
-// a very naïve approach to testing for primes
-func isPrime(val uint64) bool {
-	if val == 0 || val == 1 {
-		return false
-	}
-
-	if val == 2 {
-		return true
-	}
-
-	for i := uint64(2); i <= uint64(math.Sqrt(float64(val))); i++ {
-		if val%i == 0 {
-			return false
-		}
-	}
-
-	return true
-}

primes.go

@@ -20,8 +20,10 @@ func Between(a, b uint64, algo SieveAlgo) []uint64 {
 	switch algo {
 	case EratosthenesAlgo:
 		s = NewEratosthenes(b)
+	case SundaramAlgo:
+		s = NewSundaram(b)
 	default:
-		fmt.Fprintf(os.Stderr, "unknown sieve algorithm: %v", algo)
+		fmt.Fprintf(os.Stderr, "unknown sieve algorithm: %v\n", algo)
 		return nil
 	}
 

primes_test.go

@@ -1,6 +1,7 @@
 package primes
 
 import (
+	"math"
 	"testing"
 
 	. "github.com/smartystreets/goconvey/convey"
@@ -18,7 +19,7 @@ func TestPrimes(t *testing.T) {
 	Loop:
 		for _, algo := range SieveAlgos {
 			switch algo {
-			case SundaramAlgo, AtkinAlgo:
+			case AtkinAlgo:
 				Print("TODO(jrubin) enable prime test for ", algo, "\n")
 				continue Loop
 			}
@@ -31,3 +32,22 @@ func TestPrimes(t *testing.T) {
 		}
 	})
 }
+
+// a very naïve approach to testing for primes
+func isPrime(val uint64) bool {
+	if val == 0 || val == 1 {
+		return false
+	}
+
+	if val == 2 {
+		return true
+	}
+
+	for i := uint64(2); i <= uint64(math.Sqrt(float64(val))); i++ {
+		if val%i == 0 {
+			return false
+		}
+	}
+
+	return true
+}

sundaram.go

@@ -0,0 +1,46 @@
+package primes
+
+// Sundaram is a Sieve that is calculated using the sieve of sundaram algorithm
+type Sundaram struct {
+	BitSet
+}
+
+// NewSundaram returns a new Sundaram calculated for all values from 0 to the
+// nearest byte greater than n
+func NewSundaram(n uint64) Sieve {
+	s := NewBitSet(n)
+
+	imax := sqrt(s.Max())
+	for i := uint64(1); i <= imax; i++ {
+		jmax := (n - i) / (1 + 2*i)
+		for j := i; j <= jmax; j++ {
+			s.Set(i + j + 2*i*j)
+		}
+	}
+
+	ret := Sundaram{NewBitSet(n).Set(2).Set(3)}
+	imax = (s.Max() - 1) / 2
+	for i := uint64(2); i <= imax; i++ {
+		if !s.IsSet(i) {
+			ret.Set(2*i + 1)
+		}
+	}
+
+	return ret
+}
+
+// IsPrime returns if value is a prime
+func (s Sundaram) IsPrime(i uint64) bool {
+	return s.IsSet(i)
+}
+
+// ListPrimes returns the set of all primes in the sieve
+func (s Sundaram) ListPrimes() []uint64 {
+	return s.ListSet()
+}
+
+// Len returns the size of the BitSet in bytes
+func (s Sundaram) Len() uint64 {
+	// twice the bitset size since 2 were required
+	return s.BitSet.Len() * 2
+}

sundaram_test.go

@@ -0,0 +1,30 @@
+package primes
+
+import (
+	"testing"
+
+	. "github.com/smartystreets/goconvey/convey"
+)
+
+func TestSundaram(t *testing.T) {
+	Convey("Sundaram algorithm should work", t, func() {
+		max := uint64(10000)
+		s := NewSundaram(max)
+		So(s.Len(), ShouldEqual, 2502)
+		numPrimes := 0
+		for i := uint64(0); i <= max; i++ {
+			if s.IsPrime(i) {
+				numPrimes++
+				So(isPrime(i), ShouldBeTrue)
+			} else {
+				So(isPrime(i), ShouldBeFalse)
+			}
+		}
+		So(numPrimes, ShouldEqual, 1229)
+		ps := s.ListPrimes()
+		So(len(ps), ShouldEqual, 1230)
+		for _, p := range ps {
+			So(isPrime(p), ShouldBeTrue)
+		}
+	})
+}