diff --git a/internal/lib/sieve/common.go b/internal/lib/sieve/common.go new file mode 100644 index 0000000..1045e9a --- /dev/null +++ b/internal/lib/sieve/common.go @@ -0,0 +1,32 @@ +/* +Copyright © 2025 filifa + +This program is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with this program. If not, see . +*/ +package sieve + +func updateMultiples(sieve []uint, p uint, n uint) { + for q := p; ; q *= p { + // sieve[a*b] = sieve[a] * sieve[b] if gcd(a,b) = 1 + for i := 2 * q; i < n; i += q { + if i%(p*q) != 0 { + sieve[i] *= sieve[q] + } + } + + if p*q >= n { + break + } + } +} diff --git a/internal/lib/sieve/divisors.go b/internal/lib/sieve/divisors.go index ae4f05d..2a8d9df 100644 --- a/internal/lib/sieve/divisors.go +++ b/internal/lib/sieve/divisors.go @@ -30,24 +30,6 @@ func pow(base uint, exp uint) uint { return result } -func updateMultiples(sieve []uint, x uint, p uint, n uint) { - for q := p; ; q *= p { - // sigma_x(a*b) = sigma_x(a) * sigma_x(b) if gcd(a,b) = 1 - for i := 2 * q; i < n; i += q { - if i%(p*q) != 0 { - sieve[i] *= sieve[q] - } - } - - if p*q >= n { - break - } - - // sigma_x(p^k) = p^(kx) + sigma_x(p^(k-1)) - sieve[p*q] = pow(p*q, x) + sieve[q] - } -} - /* Divisors computes sigma_x(k) for k=1 to n, where sigma_x is the divisor sum function. x sets the power each divisor is raised to. */ @@ -68,7 +50,12 @@ func Divisors(n uint, x uint, buflen uint) chan uint { } sieve[i] = pow(i, x) + 1 - updateMultiples(sieve, x, i, n) + for j := i; i*j < n; j *= i { + // sigma_x(p^k) = p^(kx) + sigma_x(p^(k-1)) + sieve[i*j] = pow(i*j, x) + sieve[j] + } + + updateMultiples(sieve, i, n) ch <- sieve[i] } }()