Primes

main.go

package main

import (
	"container/heap"
	"fmt"
	"math"
	"math/big"
	"sort"
	"time"
)

func main() {
	// azimuths := []int{1, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 49, 53, 59, 61, 67, 71, 73, 77, 79, 83, 89, 91, 97, 101, 103, 107, 109, 113, 119, 121, 127, 131, 133, 137, 139, 143, 149, 151, 157, 161, 163, 167, 169, 173, 179, 181, 187, 191, 193, 197, 199, 203, 209, 211, 217, 221, 223, 227, 229, 233, 239, 241, 247, 251, 253, 257, 259, 263, 269, 271, 277, 281, 283, 287, 289, 293, 299, 301, 307, 311, 313, 317, 319, 323, 329, 331, 337, 341, 343, 347, 349, 353, 359}
	azimuths := generateAzimuths()
	// smallPrimes := []int{7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359}
	var limit int = 1e5
	smallPrimes := serialPrimes(azimuths, int(math.Sqrt(float64(limit))))
	fmt.Printf("Small Primes: %d\n", len(smallPrimes))

	// sp := 23
	start := time.Now()
	// spirals := bruteForceSpirals(sp, azimuths)
	// azimuthMap := make(map[int]PriorityQueue)
	azimuthSpiralQueue := make(PriorityQueue, 0, len(smallPrimes)*len(azimuths))
	for _, p := range smallPrimes {
		spirals := predictedSpirals(p, azimuths)
		azimuthSpiralQueue = append(azimuthSpiralQueue, spirals...)

		// for _, s := range spirals {
		// 	if _, ok := azimuthMap[s.Azm]; !ok {
		// 		azimuthMap[s.Azm] = make(PriorityQueue, 0, len(smallPrimes))
		// 	}
		// 	azimuthMap[s.Azm] = append(azimuthMap[s.Azm], s)
		// }
	}
	fmt.Printf("Generating spirals: %s\n", time.Since(start))

	heap.Init(&azimuthSpiralQueue)
	// for _, v := range azimuthMap {
	// heap.Init(&v)
	// }

	var primeCount int
	start = time.Now()
	pq := azimuthSpiralQueue
	for i := 23; i < limit; i += 360 {
		for i > pq[0].Value {
			// fmt.Printf("Incrementing spiral: %d, [%d] -> %d\n", pq[0].Value, pq[0].Prime, pq[0].Next)
			pq.MultIncrStep(pq[0], i)
		}

		if i == pq[0].Value {
			// fmt.Printf("Overlapping spiral: %d, [%d]\n", i, pq[0].Prime)

			// Found a non-prime along this azimuth
			pq.IncrStep(pq[0])
			continue
		}

		primeCount++
		// fmt.Printf("Found Prime: %d\n", i)
	}
	fmt.Printf("Single Azimuth Elapsed: %s\n", time.Since(start))
	fmt.Printf("Primes: %d\n", primeCount)

	azimuthSpiralQueue.Reset()

	missingPrimes := map[int]int{8641: 0, 1801: 0, 6481: 0, 6841: 0, 2521: 0, 6121: 0, 9001: 0, 2161: 0, 7561: 0, 9721: 0}
	tmp := big.NewInt(0)

	start = time.Now()
	filterTable := make(map[int]int)
	filterTable[1] = 0
	for _, item := range azimuthSpiralQueue {
		if _, ok := missingPrimes[item.Value]; ok {
			fmt.Printf("Found missing prime! %d [%d] on spiral %d\n", item.Value, item.Prime, item.Azm)
		}

		for item.Value < limit {
			if tmp.SetInt64(int64(item.Value)).ProbablyPrime(0) {
				fmt.Printf("Found prime on spiral: %d [%d] on azm %d\n", item.Value, item.Prime, item.Azm)
			}

			filterTable[item.Value] = 0
			item.Value += item.StepSize
			item.Next += item.StepSize
		}
	}
	fmt.Printf("Generating Filter: %s\n", time.Since(start))
	fmt.Printf("Filter Size: %d\n", len(filterTable))

	start = time.Now()
	primes := []int{2, 3, 5}
	for index := 0; index < len(azimuths); index++ {
		for candidate := azimuths[index]; candidate < limit; candidate += 360 {
			if _, ok := filterTable[candidate]; ok {
				if tmp.SetInt64(int64(candidate)).ProbablyPrime(0) {
					fmt.Printf("Found a prime in the filter!!! %d\n", candidate)
				}
				continue
			}
			if candidate%360 != 23 {
				continue
			}
			primes = append(primes, candidate)
		}
	}
	fmt.Printf("Generating Primes: %s\n", time.Since(start))

	sort.Ints(primes)
	fmt.Printf("Primes: %v\n", primes)
	fmt.Printf("Found %d Primes\n", len(primes))

	// pq := azimuthMap[23]
	// for pq.Len() > 0 {
	// 	s := pq[0]
	// 	if s.Value%360 == 0 {
	// 		heap.Pop(&pq)
	// 		fmt.Printf("%d,%d,%d,%d,%d,%d,%d\n", s.Value, s.Prime, s.Value/s.Prime, s.Value/s.Prime+360, s.Prime*360, s.Next, s.Azm)
	// 	}
	// 	pq.IncrStep(s)
	// }

	// sort.Slice(azimuthMap[23], func(i, j int) bool {
	// 	return azimuthMap[23][i].Value < azimuthMap[23][j].Value
	// })

	// fmt.Println("value,prime,div,stepsize,next,azm")
	// for _, s := range azimuthSpiralQueue {
	// 	if s.Value%360 == 0 {
	// 		fmt.Printf("%d,%d,%d,%d,%d,%d\n", s.Value, s.Prime, s.Value/s.Prime, s.Prime*360, s.Next, s.Azm)
	// 	}
	// }
}

func generateAzimuths() []int {
	azm := []int{}
	lastDigits := []int{1, 3, 7, 9}
	for digit := 0; digit < len(lastDigits); digit++ {
		for a := lastDigits[digit]; a < 360; a += 10 {
			if a%3 == 0 {
				continue
			} else if a == 1 {
				// continue
			}
			azm = append(azm, a)
		}
	}
	return azm
}

func serialPrimes(azimuths []int, sqrtN int) []int {
	primes := []int{2, 3, 5}
	tmp := big.NewInt(0)
	for index := 0; index < len(azimuths); index++ {
		if azimuths[index] > sqrtN {
			continue
		}
		for candidate := azimuths[index]; candidate < sqrtN; candidate += 360 {
			if tmp.SetInt64(int64(candidate)).ProbablyPrime(0) {
				primes = append(primes, candidate)
			}
		}
	}
	sort.Ints(primes)
	return primes
}

func predictedSpirals(sp int, azimuths []int) []*AzimuthSpiral {
	var spirals []*AzimuthSpiral
	for _, azm := range azimuths {
		if azm == 1 {
			continue
		}
		i := sp * azm
		next := sp * (i/sp + 360)
		spirals = append(spirals, &AzimuthSpiral{Index: 0, Value: i, Prime: sp, StepSize: sp * 360, Next: next, Azm: i % 360, Steps: 0})
	}
	return spirals
}

// AzimuthSpiral iterates each prime along an azimuth
type AzimuthSpiral struct {
	Index    int // Index in the heap
	Value    int // Value of the spiral
	Prime    int // Prime for this spiral
	StepSize int // Step size. This value is prime * 360.
	Next     int // Value + Step size
	Azm      int // Azimuth of each step
	Steps    int // Number of steps taken
}

func (a *AzimuthSpiral) Step() {
	a.Value += a.StepSize
	a.Next += a.StepSize
	a.Steps++
}

func (a *AzimuthSpiral) Reset() {
	a.Value -= a.StepSize * a.Steps
	a.Next = a.StepSize + a.Value
	a.Steps = 0
}

// A PriorityQueue implements heap.Interface and holds Items.
type PriorityQueue []*AzimuthSpiral

func (pq PriorityQueue) Len() int { return len(pq) }

func (pq PriorityQueue) Less(i, j int) bool {
	return pq[i].Value < pq[j].Value
}

func (pq PriorityQueue) Swap(i, j int) {
	pq[i], pq[j] = pq[j], pq[i]
	pq[i].Index = i
	pq[j].Index = j
}

func (pq *PriorityQueue) Push(x interface{}) {
	n := len(*pq)
	item := x.(*AzimuthSpiral)
	item.Index = n
	*pq = append(*pq, item)
}

func (pq *PriorityQueue) Pop() interface{} {
	old := *pq
	n := len(old)
	item := old[n-1]
	old[n-1] = nil  // avoid memory leak
	item.Index = -1 // for safety
	*pq = old[0 : n-1]
	return item
}

// update modifies the priority and value of an Item in the queue.
func (pq *PriorityQueue) IncrStep(item *AzimuthSpiral) {
	item.Step()
	heap.Fix(pq, item.Index)
}

// update modifies the priority and value of an Item in the queue.
func (pq *PriorityQueue) MultIncrStep(item *AzimuthSpiral, watermark int) {
	for item.Value < watermark {
		item.Step()
	}
	heap.Fix(pq, item.Index)
}

func (pq *PriorityQueue) Reset() {
	for _, item := range *pq {
		item.Reset()
	}
	heap.Init(pq)
}