This commit is contained in:
filifa 2025-12-05 22:29:56 -05:00
parent 0dd446fe74
commit 592fadb6d7
1 changed files with 18 additions and 15 deletions

View File

@ -18,21 +18,21 @@ package sieve
// NOTE: these formulas come from https://web.archive.org/web/20151224013638/http://www.maa.org/sites/default/files/Walter_D22068._Stangl.pdf
func updatePowersOfTwoCoprime(sieve []uint, n uint) {
func qrPowersOfTwoCoprime(sieve []uint, n uint) {
for q := uint(8); 2*q < n; q *= 2 {
// q(2^n) = 2^(n-3) for n >= 3
sieve[2*q] = 2 * sieve[q]
}
}
func updatePowersOfOddPrimesCoprime(sieve []uint, p uint, n uint) {
func qrPowersOfOddPrimesCoprime(sieve []uint, p uint, n uint) {
for q := p; p*q < n; q *= p {
// q(p^n) = (p^n - p^(n-1)) / 2
sieve[p*q] = (p*q - q) / 2
}
}
func updatePowersOfTwo(sieve []uint, n uint) {
func qrPowersOfTwo(sieve []uint, n uint) {
k := 1
for q := uint(1); 2*q < n; q *= 2 {
if k%2 == 0 {
@ -47,13 +47,10 @@ func updatePowersOfTwo(sieve []uint, n uint) {
}
}
func updatePowersOfOddPrimes(sieve []uint, p uint, n uint) {
k := 2
for q := p; p*q < n; q *= p {
if p == q {
// s(p^2) = (p^2 - p + 2) / 2
sieve[p*q] = (p*p - p + 2) / 2
} else if k%2 == 0 {
func qrPowersOfOddPrimes(sieve []uint, p uint, n uint) {
k := 3
for q := p * p; p*q < n; q *= p {
if k%2 == 0 {
// s(p^n) = (p^(n+1) + p + 2) / (2*(p+1)) for even n
sieve[p*q] = (p*p*q + p + 2) / (2 * (p + 1))
} else {
@ -68,7 +65,7 @@ func updatePowersOfOddPrimes(sieve []uint, p uint, n uint) {
/*
QuadraticResidues computes the number of quadratic residues modulo k for k=1 to n.
see https://oeis.org/A046073
see https://oeis.org/A000224 and https://oeis.org/A046073
*/
func QuadraticResidues(n uint, coprime bool, buflen uint) chan uint {
sieve := make([]uint, n)
@ -97,10 +94,11 @@ func sieveQRCoprime(sieve []uint, n uint, ch chan uint) {
}
if i == 2 {
updatePowersOfTwoCoprime(sieve, n)
qrPowersOfTwoCoprime(sieve, n)
} else {
// q(p) = (p - 1) / 2
sieve[i] = (i - 1) / 2
updatePowersOfOddPrimesCoprime(sieve, i, n)
qrPowersOfOddPrimesCoprime(sieve, i, n)
}
updateMultiples(sieve, i, n, false)
@ -116,10 +114,15 @@ func sieveQR(sieve []uint, n uint, ch chan uint) {
}
if i == 2 {
updatePowersOfTwo(sieve, n)
qrPowersOfTwo(sieve, n)
} else {
// s(p) = (p + 1) / 2
sieve[i] = (i + 1) / 2
updatePowersOfOddPrimes(sieve, i, n)
if i*i < n {
// s(p^2) = (p^2 - p + 2) / 2
sieve[i*i] = (i*i - i + 2) / 2
}
qrPowersOfOddPrimes(sieve, i, n)
}
updateMultiples(sieve, i, n, false)