sieve.go

package primegen

import (
	"runtime"
	"sync"
)

var two = [32]uint32{
	0x00000001, 0x00000002, 0x00000004, 0x00000008,
	0x00000010, 0x00000020, 0x00000040, 0x00000080,
	0x00000100, 0x00000200, 0x00000400, 0x00000800,
	0x00001000, 0x00002000, 0x00004000, 0x00008000,
	0x00010000, 0x00020000, 0x00040000, 0x00080000,
	0x00100000, 0x00200000, 0x00400000, 0x00800000,
	0x01000000, 0x02000000, 0x04000000, 0x08000000,
	0x10000000, 0x20000000, 0x40000000, 0x80000000,
}

func doit4(a []uint32, x, y int, start int64) {
	var i, i0, y0 int

	x += x
	x += 15
	y += 15

	start += 1000000000
	for start < 0 {
		start += int64(x)
		x += 30
	}
	start -= 1000000000
	i = int(start)

	for i < _B {
		i += x
		x += 30
	}

	for {
		x -= 30
		if x <= 15 {
			return
		}
		i -= x

		for i < 0 {
			i += y
			y += 30
		}

		i0 = i
		y0 = y
		for i < _B {
			pos := uint32(i) >> 5
			data := uint32(i) & 31
			i += y
			y += 30
			bits := a[pos]
			data = two[data]
			bits ^= data
			a[pos] = bits
		}
		i = i0
		y = y0
	}
}

func doit6(a []uint32, x, y int, start int64) {
	var i, i0, y0 int

	x += 5
	y += 15

	start += 1000000000
	for start < 0 {
		start += int64(x)
		x += 10
	}
	start -= 1000000000
	i = int(start)
	for i < _B {
		i += x
		x += 10
	}

	for {
		x -= 10
		if x <= 5 {
			return
		}
		i -= x

		for i < 0 {
			i += y
			y += 30
		}

		i0 = i
		y0 = y
		for i < _B {
			pos := uint32(i) >> 5
			data := uint32(i) & 31
			i += y
			y += 30
			bits := a[pos]
			data = two[data]
			bits ^= data
			a[pos] = bits
		}
		i = i0
		y = y0
	}
}

func doit12(a []uint32, x, y int, start int64) {
	var i, i0, y0 int

	x += 5

	start += 1000000000
	for start < 0 {
		start += int64(x)
		x += 10
	}
	start -= 1000000000
	i = int(start)
	for i < 0 {
		i += x
		x += 10
	}

	y += 15
	x += 10

	for {
		for i >= _B {
			if x <= y {
				return
			}
			i -= y
			y += 30
		}
		i0 = i
		y0 = y
		for (i >= 0) && (y < x) {
			pos := uint32(i) >> 5
			data := uint32(i) & 31
			i -= y
			y += 30
			bits := a[pos]
			data = two[data]
			bits ^= data
			a[pos] = bits
		}
		i = i0
		y = y0
		i += x - 10
		x += 10
	}
}

var deltainverse = [60]int{
	-1, _B32 * 0, -1, -1, -1, -1, -1, _B32 * 1, -1, -1, -1, _B32 * 2, -1, _B32 * 3, -1,
	-1, -1, _B32 * 4, -1, _B32 * 5, -1, -1, -1, _B32 * 6, -1, -1, -1, -1, -1, _B32 * 7,
	-1, _B32 * 8, -1, -1, -1, -1, -1, _B32 * 9, -1, -1, -1, _B32 * 10, -1, _B32 * 11, -1,
	-1, -1, _B32 * 12, -1, _B32 * 13, -1, -1, -1, _B32 * 14, -1, -1, -1, -1, -1, _B32 * 15,
}

func squarefree1big(buf [16][]uint32, base uint64, q uint32, qq uint64) {
	var i uint64
	bound := base + 60*_B

	for qq < bound {
		if bound < 2000000000 {
			i = qq - uint64(uint32(base)%uint32(qq))
		} else {
			i = qq - (base % qq)
		}
		if (i & 1) == 0 {
			i += qq
		}

		if i < _B*60 {
			pos := uint32(i)
			n := deltainverse[pos%60]
			if n >= 0 {
				pos /= 60
				idx := n + int(pos>>5)
				buf[idx/_B32][idx%_B32] |= two[pos&31]
			}
		}

		qq += uint64(q)
		q += 1800
	}
}

func squarefree1(buf [16][]uint32, L uint64, q uint32) {
	var i uint32

	base := 60 * L
	qq := q * q
	q = 60*q + 900

	for qq < _B*60 {
		if base < 2000000000 {
			i = qq - (uint32(base) % qq)
		} else {
			i = qq - uint32(base%uint64(qq))
		}
		if (i & 1) == 0 {
			i += qq
		}

		if i < _B*60 {
			qqhigh := qq / 60
			ilow := i % 60
			ihigh := i / 60

			qqhigh += qqhigh
			for ihigh < _B {
				n := deltainverse[ilow]
				if n >= 0 {
					idx := n + int(ihigh>>5)
					buf[idx/_B32][idx%_B32] |= two[ihigh&31]
				}

				ilow += 2
				ihigh += qqhigh
				if ilow >= 60 {
					ilow -= 60
					ihigh += 1
				}
			}
		}

		qq += q
		q += 1800
	}

	squarefree1big(buf, base, q, uint64(qq))
}

func squarefree49big(buf [16][]uint32, base uint64, q uint32, qq uint64) {
	var i uint64
	bound := base + 60*_B

	for qq < bound {
		if bound < 2000000000 {
			i = qq - uint64(uint32(base)%uint32(qq))
		} else {
			i = qq - (base % qq)
		}
		if (i & 1) == 0 {
			i += qq
		}

		if i < _B*60 {
			pos := uint32(i)
			n := deltainverse[pos%60]
			if n >= 0 {
				pos /= 60
				idx := n + int(pos>>5)
				buf[idx/_B32][idx%_B32] |= two[pos&31]
			}
		}

		qq += uint64(q)
		q += 1800
	}
}

func squarefree49(buf [16][]uint32, L uint64, q uint32) {
	var i uint32

	base := 60 * L
	qq := q * q
	q = 60*q + 900

	for qq < _B*60 {
		if base < 2000000000 {
			i = qq - (uint32(base) % qq)
		} else {
			i = qq - uint32(base%uint64(qq))
		}
		if (i & 1) == 0 {
			i += qq
		}

		if i < _B*60 {
			qqhigh := qq / 60
			ilow := i % 60
			ihigh := i / 60

			qqhigh += qqhigh
			qqhigh += 1
			for ihigh < _B {
				n := deltainverse[ilow]
				if n >= 0 {
					idx := n + int(ihigh>>5)
					buf[idx/_B32][idx%_B32] |= two[ihigh&31]
				}

				ilow += 38
				ihigh += qqhigh
				if ilow >= 60 {
					ilow -= 60
					ihigh += 1
				}
			}
		}

		qq += q
		q += 1800
	}

	squarefree49big(buf, base, q, uint64(qq))
}

/* squares of primes >= 7, < 240 */
var qqtab = [49]uint32{
	49, 121, 169, 289, 361, 529, 841, 961, 1369, 1681,
	1849, 2209, 2809, 3481, 3721, 4489, 5041, 5329, 6241, 6889,
	7921, 9409, 10201, 10609, 11449, 11881, 12769, 16129, 17161, 18769,
	19321, 22201, 22801, 24649, 26569, 27889, 29929, 32041, 32761, 36481,
	37249, 38809, 39601, 44521, 49729, 51529, 52441, 54289, 57121,
}

/* (qq * 11 + 1) / 60 or (qq * 59 + 1) / 60 */
var qq60tab = [49]uint32{
	9, 119, 31, 53, 355, 97, 827, 945, 251, 1653,
	339, 405, 515, 3423, 3659, 823, 4957, 977, 6137,
	1263, 7789, 1725, 10031, 1945, 2099, 11683, 2341, 2957,
	16875, 3441, 18999, 21831, 22421, 4519, 4871, 5113, 5487,
	31507, 32215, 35873, 6829, 7115, 38941, 43779, 9117, 9447, 51567, 9953, 56169,
}

func squarefreetiny(a []uint32, Lmodqq []uint32, d uint32) {
	for j := 0; j < 49; j++ {
		qq := qqtab[j]
		k := qq - 1 - ((Lmodqq[j] + qq60tab[j]*d - 1) % qq)
		for k < _B {
			pos := k >> 5
			data := k & 31
			k += qq
			bits := a[pos]
			data = two[data]
			bits |= data
			a[pos] = bits
		}
	}
}

type todo struct {
	index, f, g, k int
}

var for4 = []todo{
	{0, 2, 15, 4}, {0, 3, 5, 1}, {0, 3, 25, 11}, {0, 5, 9, 3},
	{0, 5, 21, 9}, {0, 7, 15, 7}, {0, 8, 15, 8}, {0, 10, 9, 8},
	{0, 10, 21, 14}, {0, 12, 5, 10}, {0, 12, 25, 20}, {0, 13, 15, 15},
	{0, 15, 1, 15}, {0, 15, 11, 17}, {0, 15, 19, 21}, {0, 15, 29, 29},
	{3, 1, 3, 0}, {3, 1, 27, 12}, {3, 4, 3, 1}, {3, 4, 27, 13},
	{3, 6, 7, 3}, {3, 6, 13, 5}, {3, 6, 17, 7}, {3, 6, 23, 11},
	{3, 9, 7, 6}, {3, 9, 13, 8}, {3, 9, 17, 10}, {3, 9, 23, 14},
	{3, 11, 3, 8}, {3, 11, 27, 20}, {3, 14, 3, 13}, {3, 14, 27, 25},
	{4, 2, 1, 0}, {4, 2, 11, 2}, {4, 2, 19, 6}, {4, 2, 29, 14},
	{4, 7, 1, 3}, {4, 7, 11, 5}, {4, 7, 19, 9}, {4, 7, 29, 17},
	{4, 8, 1, 4}, {4, 8, 11, 6}, {4, 8, 19, 10}, {4, 8, 29, 18},
	{4, 13, 1, 11}, {4, 13, 11, 13}, {4, 13, 19, 17}, {4, 13, 29, 25},
	{7, 1, 5, 0}, {7, 1, 25, 10}, {7, 4, 5, 1}, {7, 4, 25, 11},
	{7, 5, 7, 2}, {7, 5, 13, 4}, {7, 5, 17, 6}, {7, 5, 23, 10},
	{7, 10, 7, 7}, {7, 10, 13, 9}, {7, 10, 17, 11}, {7, 10, 23, 15},
	{7, 11, 5, 8}, {7, 11, 25, 18}, {7, 14, 5, 13}, {7, 14, 25, 23},
	{9, 2, 9, 1}, {9, 2, 21, 7}, {9, 3, 1, 0}, {9, 3, 11, 2},
	{9, 3, 19, 6}, {9, 3, 29, 14}, {9, 7, 9, 4}, {9, 7, 21, 10},
	{9, 8, 9, 5}, {9, 8, 21, 11}, {9, 12, 1, 9}, {9, 12, 11, 11},
	{9, 12, 19, 15}, {9, 12, 29, 23}, {9, 13, 9, 12}, {9, 13, 21, 18},
	{10, 2, 5, 0}, {10, 2, 25, 10}, {10, 5, 1, 1}, {10, 5, 11, 3},
	{10, 5, 19, 7}, {10, 5, 29, 15}, {10, 7, 5, 3}, {10, 7, 25, 13},
	{10, 8, 5, 4}, {10, 8, 25, 14}, {10, 10, 1, 6}, {10, 10, 11, 8},
	{10, 10, 19, 12}, {10, 10, 29, 20}, {10, 13, 5, 11}, {10, 13, 25, 21},
	{13, 1, 15, 3}, {13, 4, 15, 4}, {13, 5, 3, 1}, {13, 5, 27, 13},
	{13, 6, 5, 2}, {13, 6, 25, 12}, {13, 9, 5, 5}, {13, 9, 25, 15},
	{13, 10, 3, 6}, {13, 10, 27, 18}, {13, 11, 15, 11}, {13, 14, 15, 16},
	{13, 15, 7, 15}, {13, 15, 13, 17}, {13, 15, 17, 19}, {13, 15, 23, 23},
	{14, 1, 7, 0}, {14, 1, 13, 2}, {14, 1, 17, 4}, {14, 1, 23, 8},
	{14, 4, 7, 1}, {14, 4, 13, 3}, {14, 4, 17, 5}, {14, 4, 23, 9},
	{14, 11, 7, 8}, {14, 11, 13, 10}, {14, 11, 17, 12}, {14, 11, 23, 16},
	{14, 14, 7, 13}, {14, 14, 13, 15}, {14, 14, 17, 17}, {14, 14, 23, 21},
}

var for6 = []todo{
	{1, 1, 2, 0}, {1, 1, 8, 1}, {1, 1, 22, 8}, {1, 1, 28, 13},
	{1, 3, 10, 2}, {1, 3, 20, 7}, {1, 7, 10, 4}, {1, 7, 20, 9},
	{1, 9, 2, 4}, {1, 9, 8, 5}, {1, 9, 22, 12}, {1, 9, 28, 17},
	{5, 1, 4, 0}, {5, 1, 14, 3}, {5, 1, 16, 4}, {5, 1, 26, 11},
	{5, 5, 2, 1}, {5, 5, 8, 2}, {5, 5, 22, 9}, {5, 5, 28, 14},
	{5, 9, 4, 4}, {5, 9, 14, 7}, {5, 9, 16, 8}, {5, 9, 26, 15},
	{8, 3, 2, 0}, {8, 3, 8, 1}, {8, 3, 22, 8}, {8, 3, 28, 13},
	{8, 5, 4, 1}, {8, 5, 14, 4}, {8, 5, 16, 5}, {8, 5, 26, 12},
	{8, 7, 2, 2}, {8, 7, 8, 3}, {8, 7, 22, 10}, {8, 7, 28, 15},
	{11, 1, 10, 1}, {11, 1, 20, 6}, {11, 3, 4, 0}, {11, 3, 14, 3},
	{11, 3, 16, 4}, {11, 3, 26, 11}, {11, 7, 4, 2}, {11, 7, 14, 5},
	{11, 7, 16, 6}, {11, 7, 26, 13}, {11, 9, 10, 5}, {11, 9, 20, 10},
}

var for12 = []todo{
	{2, 2, 1, 0}, {2, 2, 11, -2}, {2, 2, 19, -6}, {2, 2, 29, -14},
	{2, 3, 4, 0}, {2, 3, 14, -3}, {2, 3, 16, -4}, {2, 3, 26, -11},
	{2, 5, 2, 1}, {2, 5, 8, 0}, {2, 5, 22, -7}, {2, 5, 28, -12},
	{2, 7, 4, 2}, {2, 7, 14, -1}, {2, 7, 16, -2}, {2, 7, 26, -9},
	{2, 8, 1, 3}, {2, 8, 11, 1}, {2, 8, 19, -3}, {2, 8, 29, -11},
	{2, 10, 7, 4}, {2, 10, 13, 2}, {2, 10, 17, 0}, {2, 10, 23, -4},
	{6, 1, 10, -2}, {6, 1, 20, -7}, {6, 2, 7, -1}, {6, 2, 13, -3},
	{6, 2, 17, -5}, {6, 2, 23, -9}, {6, 3, 2, 0}, {6, 3, 8, -1},
	{6, 3, 22, -8}, {6, 3, 28, -13}, {6, 4, 5, 0}, {6, 4, 25, -10},
	{6, 6, 5, 1}, {6, 6, 25, -9}, {6, 7, 2, 2}, {6, 7, 8, 1},
	{6, 7, 22, -6}, {6, 7, 28, -11}, {6, 8, 7, 2}, {6, 8, 13, 0},
	{6, 8, 17, -2}, {6, 8, 23, -6}, {6, 9, 10, 2}, {6, 9, 20, -3},
	{12, 1, 4, -1}, {12, 1, 14, -4}, {12, 1, 16, -5}, {12, 1, 26, -12},
	{12, 2, 5, -1}, {12, 2, 25, -11}, {12, 3, 10, -2}, {12, 3, 20, -7},
	{12, 4, 1, 0}, {12, 4, 11, -2}, {12, 4, 19, -6}, {12, 4, 29, -14},
	{12, 6, 1, 1}, {12, 6, 11, -1}, {12, 6, 19, -5}, {12, 6, 29, -13},
	{12, 7, 10, 0}, {12, 7, 20, -5}, {12, 8, 5, 2}, {12, 8, 25, -8},
	{12, 9, 4, 3}, {12, 9, 14, 0}, {12, 9, 16, -1}, {12, 9, 26, -8},
	{15, 1, 2, -1}, {15, 1, 8, -2}, {15, 1, 22, -9}, {15, 1, 28, -14},
	{15, 4, 7, -1}, {15, 4, 13, -3}, {15, 4, 17, -5}, {15, 4, 23, -9},
	{15, 5, 4, 0}, {15, 5, 14, -3}, {15, 5, 16, -4}, {15, 5, 26, -11},
	{15, 6, 7, 0}, {15, 6, 13, -2}, {15, 6, 17, -4}, {15, 6, 23, -8},
	{15, 9, 2, 3}, {15, 9, 8, 2}, {15, 9, 22, -5}, {15, 9, 28, -10},
	{15, 10, 1, 4}, {15, 10, 11, 2}, {15, 10, 19, -2}, {15, 10, 29, -10},
}

type doitParams struct {
	bufIdx     int
	f          func([]uint32, int, int, int64)
	start, end int
	ffor       []todo
	q          uint32
}

var sieveParams = []doitParams{
	{0, doit4, 0, 16, for4, 1},
	{3, doit4, 16, 32, for4, 13},
	{4, doit4, 32, 48, for4, 17},
	{7, doit4, 48, 64, for4, 29},
	{9, doit4, 64, 80, for4, 37},
	{10, doit4, 80, 96, for4, 41},
	{13, doit4, 96, 112, for4, 49},
	{14, doit4, 112, 128, for4, 53},

	{1, doit6, 0, 12, for6, 7},
	{5, doit6, 12, 24, for6, 19},
	{8, doit6, 24, 36, for6, 31},
	{11, doit6, 36, 48, for6, 43},

	{2, doit12, 0, 24, for12, 11},
	{6, doit12, 24, 48, for12, 23},
	{12, doit12, 48, 72, for12, 47},
	{15, doit12, 72, 96, for12, 59},
}

func doit(pg *Primegen, Lmodqq []uint32, paramsIdx chan int, wg *sync.WaitGroup) {
	for idx := range paramsIdx {
		if idx == -1 {
			break
		}
		p := sieveParams[idx]
		buf := pg.buf[p.bufIdx]
		for i := p.start; i < p.end; i++ {
			p.f(buf, p.ffor[i].f, p.ffor[i].g, int64(uint64(p.ffor[i].k)-pg._L))
		}
		squarefreetiny(buf, Lmodqq, p.q)
	}
	wg.Done()
}

func (pg *Primegen) sieve() {
	var i, j int
	Lmodqq := make([]uint32, 49)

	buf := pg.buf
	L := pg._L

	if L > 2000000000 {
		for i = 0; i < 49; i++ {
			Lmodqq[i] = uint32(L % uint64(qqtab[i]))
		}

	} else {
		for i = 0; i < 49; i++ {
			Lmodqq[i] = uint32(L % uint64(qqtab[i]))
		}
	}

	for j = 0; j < 16; j++ {
		for i = 0; i < _B32; i++ {
			buf[j][i] = uint32(0xffffffff)
		}
	}

	// Create GOMAXPROCS doit threads.
	maxprocs := runtime.GOMAXPROCS(-1)
	cp := make(chan int)
	var wg sync.WaitGroup
	wg.Add(maxprocs)
	for i = 0; i < maxprocs; i++ {
		go doit(pg, Lmodqq, cp, &wg)
	}

	// Send params index to threads.
	for i := 0; i < len(sieveParams); i++ {
		cp <- i
	}
	// Notify threads that there's no more work.
	for i = 0; i < maxprocs; i++ {
		cp <- -1
	}
	wg.Wait()

	squarefree49(buf, L, 247)
	squarefree49(buf, L, 253)
	squarefree49(buf, L, 257)
	squarefree49(buf, L, 263)
	squarefree1(buf, L, 241)
	squarefree1(buf, L, 251)
	squarefree1(buf, L, 259)
	squarefree1(buf, L, 269)
}