Sieving Multiplicative Functions

Demo

Note that for large ranges, outputting the full result will chew up a lot of memory.

Multiplicative functions

Multiplicative functions are functions with a special property: if two numbers $a$ and $b$ are coprime (i.e. their greatest common divisor is 1), then $f\left(ab\right)=f\left(a\right)f\left(b\right)$. This means if we can factor a number $n$ into coprime integers $x$ and $y$, we can find $f\left(n\right)$ by instead finding $f\left(x\right)$ and $f\left(y\right)$ and multiplying them.

Going a step further, if we have a prime factorization of a number $n$, along with a formula for $f$ at prime powers, we can easily compute $f\left(n\right)$ by evaluating $f$ at each prime power and multiplying.

It turns out there are lots of multiplicative functions, including many important functions in number theory, such as:

• the divisor function
• Euler's totient function
• the Mobius function

Furthermore, many of these functions have very simple formulas for prime power inputs, making them easy to evaluate if you have a prime factorization.

The algorithm

Here's the routine for sieving multiplicative functions:

function sieveSegment(l, r, f, primes, additive) {
	const [sieve, quo] = initArrays(l, r, additive);

	for (let i = 0; i < primes.length; i++) {
		if (i*i > r) {
			break;
		}

		if (!primes[i]) {
			continue;
		}

		sievePowers(i, l, r, f, sieve, quo, additive);
	}

	unusuals(f, sieve, quo, additive);

	return sieve;
}

l and r are the lower and upper bounds for the sieve, respectively. f is the function we are sieving. primes is an array indicating all the prime numbers up to the square root of r, which we can precompute with the standard sieve. Finally, additive is a boolean indicating if the function is additive or multiplicative.

The first step of the function is to initialize two arrays of size r - l called sieve and quo. sieve is what ultimately gets returned and stores the values of the function we're sieving. quo is an internal array we use to finalize the values of sieve at the very end of the function.

function initArrays(l, r, additive) {
	const sieve = Array(r - l);
	const quo = Array(r - l);
	for (let i = 0; i < r-l; i++) {
		quo[i] = l + i;
		if (additive) {
			sieve[i] = 0;
		} else {
			sieve[i] = 1;
		}
	}

	return [sieve, quo];
}

Say we want to calculate the number of divisors function from 30 to 60. After initArrays runs, here's what the sieve and quo arrays look like:

Next, we'll call the following functions for every prime number in the primes array. This is where the function spends most of its time.

function sievePowers(p, l, r, f, sieve, quo, additive) {
	let i = 1;
	for (let q = p; q < r; q *= p) {
		const y = f(p, i);
		sieveMultiples(p, q, l, r, y, sieve, quo, additive);
		i++;
	}
}

function sieveMultiples(p, q, l, r, y, sieve, quo, additive) {
	const low = Math.ceil(l / q);
	for (let i = low * q; i < r; i += q) {
		if (i % (p*q) === 0) {
			continue;
		}

		if (additive) {
			sieve[i - l] += y;
		} else {
			sieve[i - l] *= y;
		}

		quo[i - l] /= q
	}
}

In short, these functions evaluate our function at prime powers only, then update the values in the sieve array for non-prime powers using the multiplicative property. We also update the quo array by dividing out the factors we find along the way.

Here's what the arrays will look like after sieving with 2:

And here's what the arrays will look like after subsequently sieving with 3:

After we've gone through the entire primes array, the two arrays will look like this:

At this point, there's one last step. Numbers can have up to one prime factor greater than their own square root. When that happens, they're called unusual numbers.³ But since we only require the primes array to hold primes up to the square root of r, any unusual number will have one last prime factor that we never evaluated. Fortunately, we can identify those numbers easily using the quo array, so we just need to do one more pass through that array to finalize everything.

function unusuals(f, sieve, quo, additive) {
	for (let i = 0; i < quo.length; i++) {
		const q = quo[i];
		if (q === 1) {
			continue;
		}

		if (additive) {
			sieve[i] += f(q, 1);
		} else {
			sieve[i] *= f(q, 1);
		}
	}
}

Final thoughts

This all ends up being more lines of code than writing a single sieve for, say, totient or the number of divisors. However, if you want to sieve multiple functions, it can be tricky to make the right modifications to an existing sieve to work with a different function.

In contrast, with this approach, you can use the same routine and simply pass a different function as an argument. For example, here's all the functions available in the above demo:

function numDivisors(p, k) {
	return k + 1;
}

function sumDivisors(p, k) {
	return (p**(k+1) - 1) / (p - 1)
}

function totient(p, k) {
	return p**(k-1) * (p - 1);
}

function radical(p, k) {
	return p;
}

function littleOmega(p, k) {
	return 1;
}

function bigOmega(p, k) {
	return k;
}

function mobius(p, k) {
	if (k === 1) {
		return -1;
	}

	return 0;
}

Many multiplicative and additive functions have simple formulas for prime powers, which makes adding a new function to sieve very easy with this approach.

Footnotes

1. ¹ This is because the only numbers that divide $p^k$ are $1$, $p$, $p^2$, ..., up to $p^k$, which is $k+1$ total numbers.
2. ² You might be thinking at this point that factoring is also not very efficient. Sit tight; this is only an example to get you comfortable with multiplicative functions. We're going to take advantage of this property later.
3. ³ "Unusual" is actually a bit of a misnomer - most numbers have a prime factor greater than their square root, so unusual numbers are actually more common than so-called "usual" numbers!
4. ⁴ Well, technically unusuals and sievePowers are also dependent on the time complexity of the function we're sieving, but generally that complexity isn't too bad, so I'm going to gloss over that bit.