diff --git a/junk/faulhaber/faulhaber.js b/junk/faulhaber/faulhaber.js new file mode 100644 index 0000000..26b78d1 --- /dev/null +++ b/junk/faulhaber/faulhaber.js @@ -0,0 +1,111 @@ +/* + * Copyright (C) 2026 filifa + * + * License: 0BSD (https://spdx.org/licenses/0BSD.html) + */ +let exp = null; +let s = null; + +function stirling2(p) { + if (p === 0) { + return [1n]; + } + + let s = [0n]; + for (let i = 1; i < p + 1; i++) { + s.push(1n); + for (let j = i - 1; j > 0; j--) { + s[j] = BigInt(j) * s[j] + s[j-1]; + } + } + + return s; +} + +function faulhaber(n, s) { + const limit = BigInt(s.length); + if (limit === 1n) { + return n; + } + + let q = n + 1n; + let total = 0n; + for (let k = 1n; k < limit; k++) { + q *= n + 1n - k; + total += s[k] * q / (k + 1n); + } + + return total; +} + +function calculate() { + const limit = document.querySelector("#limit"); + const exponent = document.querySelector("#exponent"); + + const e = parseInt(exponent.value) + if (exp === null || exp !== e) { + exp = e; + s = stirling2(exp); + } + + const n = BigInt(limit.value); + const r = faulhaber(n, s); + const result = document.querySelector("#result"); + result.value = r; +} + +function calculateStirling() { + const limit = document.querySelector("#stirling-limit"); + const p = parseInt(limit.value) + s = stirling2(p); + + const table = buildTable(s); + const result = document.querySelector("#stirling-result"); + while (result.firstChild) { + result.removeChild(result.firstChild); + } + result.appendChild(table); +} + +function buildTable(values) { + const table = document.createElement("table"); + const caption = document.createElement("caption"); + caption.textContent = "Stirling numbers S(p,k) from k=0 to k=p"; + table.appendChild(caption); + + const tbody = document.createElement("tbody"); + table.appendChild(tbody); + + const n = values.length; + const ks = Array(n); + for (let i = 0; i < n; i++) { + ks[i] = i; + } + + const nrow = buildRow("k", ks); + tbody.appendChild(nrow); + + const vrow = buildRow("S(p,k)", values); + tbody.appendChild(vrow); + + return table; +} + +function buildRow(htext, values) { + const row = document.createElement("tr"); + const head = document.createElement("th"); + head.textContent = htext; + row.appendChild(head); + + const n = values.length; + for (let i = 0; i < n; i++) { + const elem = document.createElement("td"); + elem.textContent = values[i]; + row.appendChild(elem); + } + + return row; +} + +document.querySelector("#enter").addEventListener("click", () => calculate()) +document.querySelector("#stirling-enter").addEventListener("click", () => calculateStirling()) diff --git a/junk/faulhaber/index.html b/junk/faulhaber/index.html new file mode 100644 index 0000000..fc5ab78 --- /dev/null +++ b/junk/faulhaber/index.html @@ -0,0 +1,1079 @@ + + + + + + + + + + Sums of Powers + + + +
+ dairydemon.net +
+ +
+
+
+

Sums of Powers

+

Posted

+
+ +

This page is all about how we can efficiently compute large sums of + powers:

+
+ + + + 1 + p + + + + + 2 + p + + + + + 3 + p + + + + ⋯ + + + + n + p + + + +
+ +
+

Demo

+

Note that large exponents will take longer to calculate, but your + upper bound can be quite large - even thousands of digits!

+
+ + + + + + + +
+ + +
+ +
+

How it works

+

A straightforward approach to computing a sum like this would be + something like

+ +
+ 
function sumOfPowers(n, p) {
+	let total = 0;
+	for (let k = 1; k <= n; k++) {
+		total += Math.pow(k, p);
+	}
+	return total;
+}
+
A simple algorithm for computing sums of powers
+
+ +

This works fine for small values of n, but + this algorithm requires we perform n + exponentiations and n additions. In other + words, the time to run grows as n increases, + which makes computing sums with a high upper bound take a long time. + This is not the algorithm the demo uses.

+ +

In fact, the algorithm used in the demo looks quite different. + It's not that many lines of code, but it's relying on some nifty (and + slightly tricky) mathematics. If all you do is read the code, it's + hard to believe it works, so rather than simply presenting it right + now, I want to introduce it piece by piece and explain it as we + go.

+ +
+

Special cases

+

If we want to compute a sum where each term's exponent is 1, it + turns out there's a simple formula we can use + (read more about why it works here):

+
+ + + 1 + + + 2 + + + 3 + + + ⋯ + + + n + = + + + n + ( + n + + + 1 + ) + + 2 + + + +
+ +

The nice thing about this formula is the number of operations we + perform is no longer dependent on our upper bound - we'll always + perform one addition, one multiplication, and one division. That + means we can easily determine the sum of the numbers from 1 to some + extreme limit like 1,000,000,000,000,000 much faster than with the + naive algorithm.

+ +

There are similar formulas for computing sums when our exponent + is 2 or 3:

+
+ + + + 1 + 2 + + + + + 2 + 2 + + + + + 3 + 2 + + + + ⋯ + + + + n + 2 + + = + + + n + ( + n + + + 1 + ) + ( + 2 + n + + + 1 + ) + + 6 + + + + + + + 1 + 3 + + + + + 2 + 3 + + + + + 3 + 3 + + + + ⋯ + + + + n + 3 + + = + + + ( + + + n + ( + n + + + 1 + ) + + 2 + + ) + + 2 + + + +
+ +

It turns out there's a pattern that lets us construct + polynomials like these for any exponent, but the pattern + isn't super obvious.

+
+ +
+

Faulhaber's formula

+

+ Faulhaber's formula + gives the general pattern these polynomials follow. The most + commonly given version of the formula involves a special sequence + called the + Bernoulli numbers. +

+ +

Of course, there are algorithms to compute the Bernoulli + numbers, but since they're rational numbers, there are tradeoffs to + consider. If we use floats to store the sequence, we need to be + mindful of rounding errors; additionally, after the 258th term, the + Bernoulli numbers become too large to store in a 64 bit float. This + means to compute sums of large powers precisely, we'll want a + rational data type, which not every language has handy.

+ +

I suppose if you need the coefficients of the polynomial for + some reason, this is the formula for you, but if all you want is + the sum of powers for some n and + p, we can sidestep computing the Bernoulli + numbers by taking another approach.

+
+ +
+

Stirling numbers

+

To begin explaining this alternate approach, we need to + introduce the + Stirling numbers of the second kind + (as you might expect, there are Stirling numbers of the + first kind too, but they aren't super relevant here). + These numbers give the number of ways to partition a set of + n elements into k + non-empty subsets.

+ +

For instance, there's 6 ways to partition the letters a, b, c, + and d into 3 non-empty subsets:

+
    +
  • {{a,b},{c},{d}}
    • +
  • {{a,c},{b},{d}}
    • +
  • {{a},{b,c},{d}}
    • +
  • {{a,d},{b},{c}}
    • +
  • {{a},{b,d},{c}}
    • +
  • {{a},{b},{c,d}}
    • +
+ +

The Stirling numbers of the second kind (or just "Stirling + numbers" from here on) happen to have a nice recursive formula that + we can apply to compute them.

+
+ + + + { + + n + k + + } + + = + k + + { + + + n + − + 1 + + k + + } + + + + + { + + + n + − + 1 + + + k + − + 1 + + + } + + + +
+ +

To compute sums of powers for a given exponent + p, we'll need to know

+
+ + + + { + + p + 0 + + } + + , + + { + + p + 1 + + } + + , + + { + + p + 2 + + } + + , + … + , + + { + + p + p + + } + + + +
+ +

(We'll learn why in a bit.) Here's a function to get all those + values.

+ +
+ 
function stirling2(p) {
+	if (p === 0) {
+		return [1n];
+	}
+
+	let s = [0n];
+	for (let i = 1; i < p + 1; i++) {
+		s.push(1n);
+		for (let j = i - 1; j > 0; j--) {
+			s[j] = BigInt(j) * s[j] + s[j-1];
+		}
+	}
+
+	return s;
+}
+
An algorithm for computing Stirling numbers of the second kind
+
+ +

This algorithm is quadratic in the value of + p, which isn't great, but it's simple to + implement (and I don't know of any faster way).

+ +
+

Stirling numbers calculator

+

Be warned that these numbers can grow quite large, even for + small values of p.

+ +
+ + + + +
+ + +
+
+ +
+

Connection to exponents

+

To understand what the Stirling numbers have to do with our + problem, it helps to first introduce the notation for + falling factorials.

+
+ + + + x + + n + + + = + x + ( + x + − + 1 + ) + ( + x + − + 2 + ) + ⋯ + ( + x + − + n + + + 1 + ) + + +
+ +

It turns out there is a very elegant identity relating + exponents, Stirling numbers, and falling + factorials:1

+
+ + + + x + n + + = + + { + + n + 0 + + } + + + x + + 0 + + + + + + { + + n + 1 + + } + + + x + + 1 + + + + + + { + + n + 2 + + } + + + x + + 2 + + + + + ⋯ + + + + { + + n + n + + } + + + x + + n + + + + +
+ +

We can apply this formula to every term of our power sum and + rearrange to get

+
+ + + + + + + 1 + p + + + + + 2 + p + + + + + 3 + p + + + + ⋯ + + + + n + p + + = + + + + + + { + + p + 0 + + } + + ( + + 1 + + 0 + + + + + + 2 + + 0 + + + + + + 3 + + 0 + + + + + ⋯ + + + + n + + 0 + + + ) + + + + + + + + + + + + + { + + p + 1 + + } + + ( + + 1 + + 1 + + + + + + 2 + + 1 + + + + + + 3 + + 1 + + + + + ⋯ + + + + n + + 1 + + + ) + + + + + + + + + + + + + { + + p + 2 + + } + + ( + + 1 + + 2 + + + + + + 2 + + 2 + + + + + + 3 + + 2 + + + + + ⋯ + + + + n + + 2 + + + ) + + + + + + + + + + + + ⋯ + + + + + + + + + + + + + { + + p + p + + } + + ( + + 1 + + p + + + + + + 2 + + p + + + + + + 3 + + p + + + + + ⋯ + + + + n + + p + + + ) + + + + + +
+ +

This might seem more complicated at first, but despite the + similar notation to normal exponentiation, it's much easier to + compute sums of falling factorials. For + p>0,

+
+ + + + 1 + + p + + + + + + 2 + + p + + + + + + 3 + + p + + + + + ⋯ + + + + n + + p + + + = + + + ( + n + + + 1 + + ) + + + p + + + 1 + + + + + + p + + + 1 + + + + +
+

(For p=0, the sum is + just n.)2 This lets us + simplify.3

+ +
+ + + + 1 + p + + + + + 2 + p + + + + + 3 + p + + + + ⋯ + + + + n + p + + = + + { + + p + 0 + + } + + n + + + + { + + p + 1 + + } + + + + ( + n + + + 1 + + ) + + 2 + + + + 2 + + + + + { + + p + 2 + + } + + + + ( + n + + + 1 + + ) + + 3 + + + + 3 + + + + + { + + p + 3 + + } + + + + ( + n + + + 1 + + ) + + 4 + + + + 4 + + + + ⋯ + + + + { + + p + p + + } + + + + ( + n + + + 1 + + ) + + + p + + + 1 + + + + + + p + + + 1 + + + + +
+ +

The first term is 0 unless + p=0, which is a trivial + case, so we really only need to consider the remaining terms.

+
+ +
+

Implementation

+

This formula might still seem a little complicated, but you may + be surprised by how simple it is to implement. We already have a + method to compute the Stirling numbers we need. With those values + available, we can write a simple loop to compute each term and add + to our total.

+ +
+ 
function faulhaber(n, s) {
+	const limit = BigInt(s.length);
+	if (limit === 1n) {
+		return n;
+	}
+
+	let q = n + 1n;
+	let total = 0n;
+	for (let k = 1n; k < limit; k++) {
+		q *= n + 1n - k;
+		total += s[k] * q / (k + 1n);
+	}
+
+	return total;
+}
+
A function that takes an upper bound and Stirling + numbers to compute a sum of powers
+
+
+ +

Note that the loop is dependent on p and not + n, which allows us to compute sums with + extremely large upper bounds very quickly.

+ +
+

Final thoughts

+

This approach spends most of its time computing the Stirling + numbers. Additionally, it seems like the performance is largely + hampered by how large the Stirling numbers get - computing them + with a modulus is considerably quicker. The demo is a little + clever, in that it only recomputes the Stirling numbers if the + value of the exponent changes.

+ +

If you would rather use the more commonly known Faulhaber + formula, you'll need to compute the Bernoulli numbers. + Brent and Harvey + give a simple algorithm that is as efficient as how we're + calculating the Stirling numbers. The needed binomial coefficients + can also be calculated easily.

+ +

Furthermore, Brent and Harvey give another algorithm for the + Bernoulli numbers that is even more optimized (although it is + considerably more complicated). I believe using that algorithm in + conjunction with the more traditional formula would be more + efficient than the Stirling number approach.

+ +

However, in other ways, I like this formula better. For one + thing, whichever algorithm you use, Bernoulli numbers are a bit + harder to calculate, and require more lines of code. You'll have to + compute the binomial coefficients, as well.

+ +

One nice thing about the Stirling approach is the terms are + integers throughout. I think that makes this algorithm a little + better for modular calculations.

+
+
+ +
+

Footnotes

+
    +
  1. 1 For proofs of this identity, see chapter 1.9 of + Enumerative Combinatorics by Stanley and chapter 6.1 + of Concrete Mathematics by Graham, Knuth, and + Patashnik.
    2. +
  2. 2 See chapter 2.6 of Concrete + Mathematics for an explanation of this identity.
    3. +
  3. 3 A version of this formula is also mentioned in + chapter 6.5 of Concrete Mathematics.
    4. +
+
+
+
+ + + + + + diff --git a/junk/sieves/index.html b/junk/sieves/index.html new file mode 100644 index 0000000..361e2c5 --- /dev/null +++ b/junk/sieves/index.html @@ -0,0 +1,1166 @@ + + + + + + + + + + Sieving Multiplicative Functions + + + +
+ dairydemon.net +
+ +
+
+
+

Sieving Multiplicative Functions

+

Posted

+
+ +

The + sieve of Eratosthenes + is a well-known algorithm for computing all the prime numbers up to a + given bound. There are also well-known variants on the sieve. In + particular, a segmented version can be implemented to limit memory + usage.

+ + Furthermore, a similar approach can be used to "sieve" the values of + multiplicative functions. + For instance, + this page + discusses computing values of + Euler's totient function + with an approach that results in the same time complexity as the + standard sieve.

+ +

+ This Codeforces page + also discusses an approach for sieving any multiplicative + function. However, the provided code has some downsides: it's not + segmented, and the modifications needed for a given multiplicative + function are not always obvious.

+ +

Here, I'll discuss a variant of the sieve that +

+

+ +
+

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 + + ( + a + b + ) + + = + f + + ( + a + ) + + f + + ( + b + ) + + + . + This means if we can factor a number n into + coprime integers x and y, + we can find f(n) by + instead finding f(x) + and f(y) 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(n) 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:

+ +

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

+ +
+

An example

+

Suppose we want to know how many divisors 120 has. We + could try counting every divisor by hand, but that would + be pretty tedious and error-prone. We could also write a simple + loop to have a computer find the answer, but that approach wouldn't + scale well for larger numbers. Instead, we'll take advantage of the + fact that the + number of divisors function + is multiplicative to find the answer.

+ +

We'll use d(n) + to denote the function. It has a convenient formula for prime + powers:1

+ +
+ + + d + + ( + + p + k + + ) + + = + k + + + 1 + + +
+ +

This lets us answer our question by factoring + 120:2

+ +
+ + + d + + ( + 120 + ) + + = + d + + ( + + 2 + 3 + + × + 3 + × + 5 + ) + + = + d + + ( + + 2 + 3 + + ) + + d + + ( + 3 + ) + + d + + ( + 5 + ) + + = + 4 + × + 2 + × + 2 + = + 16 + + +
+ +

Sure enough, the 16 divisors are 1, 2, 3, 4, 5, 6, 8, 10, 12, + 15, 20, 24, 30, 40, 60, and 120.

+ +

Suppose, however, that we want to find the number of divisors + for each number up to 120. + There's no general algorithm to efficiently factor a number, + so it isn't a great idea to factor every number to get the + solution. We'll take a different approach, but it will still take + advantage of the multiplicative property.

+
+ +
+

What about additive functions?

+

There's another class of functions called + additive functions + with similar properties to multiplicative functions, namely that + + + f + + ( + a + b + ) + + = + f + + ( + a + ) + + + + f + + ( + b + ) + + + + whenever a and b are + coprime. We can sieve these functions in nearly the exact same + manner as multiplicative functions (e.g. the + prime omega functions + in the demo are additive). When we get to the code, you'll see the + provided functions take a boolean indicating whether the function + is additive to apply the proper logic.

+
+
+ +
+

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;
+}
+
This function performs a segmented sieve on any + multiplicative (or additive) function.
+
+ +

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];
+}
+
The initArrays function initializes arrays + that keep track of the value of the function we're sieving and + the running quotients.
+
+ +

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:

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+ Initial values of sieve and quo +
n303132333435363738394041424344454647484950515253545556575859
sieve111111111111111111111111111111
quo303132333435363738394041424344454647484950515253545556575859
+ +

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
+	}
+}
+
These functions are where most of the computation + happens.
+
+ +

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:

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+ Values of sieve and quo after sieving 2 +
n303132333435363738394041424344454647484950515253545556575859
sieve216121312141213121512131214121
quo1531133173593719395412143114523473492551135327557572959
+ +

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

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+ Values of sieve and quo after sieving 3 +
n303132333435363738394041424344454647484950515253545556575859
sieve4162219122414133211012231814221
quo531111173513719135417431152347149251713531557192959
+ +

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

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+ Values of sieve and quo just before + calling unusuals +
n303132333435363738394041424344454647484950515253545556575859
sieve8162249122818136211036231828221
quo131111171137191314114311123471111713531111192959
+

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.3 + 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);
+		}
+	}
+}
+
The last step is to finalize the values of any unusual + numbers (numbers n with a prime factor greater than + sqrt(n)).
+
+ +
+

Analysis

+

It's pretty easy to see that initArrays and + unusuals are + O(r-l).4 + Analyzing sievePowers and its impact is a bit more + involved.

+ +

Looking at sieveMultiples, its main loop performs a + fixed number of arithmetic operations + + + O + + ( + + + r + − + l + + q + + ) + + + + times. sieveMultiples, in turn, is called by + sievePowers at values of + q=p, + q=p2, + q=p3, + and so on until pk + exceeds r. This means the total time + complexity of sievePowers is4

+ +
+ + + O + + ( + + + r + − + l + + p + + + + + + r + − + l + + + p + 2 + + + + + + + r + − + l + + + p + 3 + + + + + … + ) + + = + O + + ( + + + r + − + l + + p + + ) + + + +
+ +

With this in mind, we can evaluate sieveSegment, + which calls sievePowers for every prime less than the + square root of r. Thanks to Euler, we know that + the sum of the reciprocal of primes up to k + is asymptotic to + + + + log + ⁡ + + + + log + ⁡ + + + k + + . + Therefore, the time complexity of sieveSegment is

+ +
+ + + O + + ( + + + r + − + l + + 2 + + + + + + r + − + l + + 3 + + + + + + r + − + l + + 5 + + + + … + ) + + = + O + + ( + + ( + r + − + l + ) + + + log + ⁡ + + + + log + ⁡ + + + + + r + + + + ) + + = + O + + ( + + ( + r + − + l + ) + + + log + ⁡ + + + + log + ⁡ + + + r + ) + + + +
+ +

This is the same complexity as the method given + here + for sieving the totient function.

+
+
+ +
+

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;
+}
+
Functions that can be sieved in the demo.
+
+ +

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. 1 This is because the only numbers that divide + pk are + 1, p, + p2, ..., up to + pk, which is + k+1 total numbers.
    2. + +
  2. 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. + +
  3. 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. + +
  4. 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.
    5. +
+
+
+
+ + + + + + diff --git a/junk/sieves/sieves.js b/junk/sieves/sieves.js new file mode 100644 index 0000000..d787c0d --- /dev/null +++ b/junk/sieves/sieves.js @@ -0,0 +1,222 @@ +/* + * Copyright (C) 2026 filifa + * + * License: 0BSD (https://spdx.org/licenses/0BSD.html) + */ +function sievePrimes(n) { + const isPrime = Array(n); + + isPrime[0] = false; + isPrime[1] = false; + for (let i = 2; i < n; i++) { + isPrime[i] = true; + } + + for (let i = 2; i < Math.sqrt(n); i++) { + if (!isPrime[i]) { + continue; + } + + for (let j = i * i; j < n; j += i) { + isPrime[j] = false; + } + } + + return isPrime; +} + +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; +} + +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]; +} + +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; +} + +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); + } + } +} + +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 + } +} + +function getFunc(func) { + switch (func.value) { + case "num-divisors": + return [numDivisors, false]; + case "sum-divisors": + return [sumDivisors, false]; + case "totient": + return [totient, false]; + case "mobius": + return [mobius, false]; + case "radical": + return [radical, false]; + case "little-omega": + return [littleOmega, true]; + case "big-omega": + return [bigOmega, true]; + default: + throw new Error("function not implemented!"); + } +} + +function buildTable(l, r, values) { + const table = document.createElement("table"); + const caption = document.createElement("caption"); + caption.id = "result-caption"; + caption.textContent = "Results"; + table.appendChild(caption); + + const tbody = document.createElement("tbody"); + table.appendChild(tbody); + + const ns = Array(r-l); + for (let i = l; i < r; i++) { + ns[i - l] = i; + } + + const nrow = buildRow("n", l, r, ns); + tbody.appendChild(nrow); + + const vrow = buildRow("f(n)", l, r, values); + tbody.appendChild(vrow); + + return table; +} + +function buildRow(htext, l, r, values) { + const row = document.createElement("tr"); + const head = document.createElement("th"); + head.textContent = htext; + row.appendChild(head); + + for (let i = l; i < r; i++) { + if (i === 0) { + continue; + } + + const elem = document.createElement("td"); + elem.textContent = values[i - l]; + row.appendChild(elem); + } + + return row; +} + +function calculate() { + const lower = document.querySelector("#lower"); + const upper = document.querySelector("#upper"); + const func = document.querySelector("#function"); + + const [f, additive] = getFunc(func); + + const l = parseInt(lower.value); + const r = parseInt(upper.value); + const s = Math.floor(Math.sqrt(r)); + const primes = sievePrimes(s + 1); + + const sieve = sieveSegment(l, r, f, primes, additive); + + const table = buildTable(l, r, sieve); + const result = document.querySelector("#result"); + while (result.firstChild) { + result.removeChild(result.firstChild); + } + result.appendChild(table); +} + +document.querySelector("#enter").addEventListener("click", () => calculate()) diff --git a/styles/faulhaber.css b/styles/faulhaber.css new file mode 100644 index 0000000..ddd5e48 --- /dev/null +++ b/styles/faulhaber.css @@ -0,0 +1,33 @@ +#demo-inputs { + display: grid; + grid-template-columns: 1fr 1fr; +} + +#limit-label, +#exponent-label { + text-align: end; +} + +#enter { + grid-column: 2; +} + +#result { + word-break: break-all; + max-height: 24lh; + overflow: auto; + display: inline-block; +} + +#stirling-result { + margin: 1vh; +} + +/* table styling */ +caption { + text-align: left; +} + +td { + text-align: right; +} diff --git a/styles/junk.css b/styles/junk.css new file mode 100644 index 0000000..8263636 --- /dev/null +++ b/styles/junk.css @@ -0,0 +1,47 @@ +.math-block { + overflow-x: auto; + overflow-y: hidden; +} + +pre { + overflow: auto; + background-color: black; + color: lime; + padding: 1ch; + border: 1vh solid gray; +} + +hgroup p { + font-style: italic; +} + +figcaption { + text-align: center; +} + +#footnote-list { + list-style-type: none; +} + +/* table styling */ +table { + display: block; + overflow: auto; + border-collapse: collapse; +} + +caption { + font-weight: bold; +} + +tbody { + border: 2px solid black; +} + +th { + background-color: mediumturquoise; +} + +th, td { + border: 1px solid black; +} diff --git a/styles/sieves.css b/styles/sieves.css new file mode 100644 index 0000000..4b4e744 --- /dev/null +++ b/styles/sieves.css @@ -0,0 +1,26 @@ +/* demo styling */ +#demo-inputs { + display: grid; + grid-template-columns: 1fr 1fr; +} + +label { + text-align: end; +} + +#enter { + grid-column: 2; +} + +#result { + margin: 1vh; +} + +/* table styling */ +#result-caption { + text-align: left; +} + +td { + text-align: right; +}