Sumas de prefijos de funciones aritméticas (Parte 1)

En este módulo introducimos cómo calcular sumas de prefijos de ciertas funciones aritméticas en tiempo sublineal. Aquí hay algunos ejemplos:

Hecho	Fuente	Nombre	Dificultad	Tags	Solución
	CSES	Sum of Divisors	Fácil		Solución

Hecho	Fuente	Nombre	Dificultad	Tags	Solución
	YS	Sum of Totient Function	Normal		en el módulo

Hecho	Fuente	Nombre	Dificultad	Tags	Solución
	YS	Counting Primes	Normal		en el módulo

Este módulo (parte 1) se enfocará en temas relacionados con los primeros dos problemas foco. El conteo de primos y aplicaciones relacionadas se pospondrán a la parte 2.

Funciones multiplicativas

Las funciones sobre las que quisiéramos calcular sumas de prefijos en los primeros dos problemas foco son ambas multiplicativas.

Definición

Si una función $f: \mathbb{Z}^+ \rightarrow \mathbb{C}$ mapea enteros positivos a números complejos, es una función aritmética.
Si $f(n)$ es una función aritmética, $f(1) = 1$ y $f(p\cdot q) = f(p) \cdot f(q)$ para cualesquiera enteros positivos coprimos $p$ , $q$ , es una función multiplicativa.
Si $f(n)$ es multiplicativa y $f(p \cdot q)$ = $f(p)\cdot f(q)$ para cualesquiera enteros positivos $p$ , $q$ , es una función completamente multiplicativa.

Si $f(n)$ es una función multiplicativa, entonces para un entero positivo $n = \prod_{i=1}^{k} p_i^{k_i}$ , tenemos

f(n) = \prod_{i=1}^{k} f(p_i^{k_i})

Si $f(n)$ es una función completamente multiplicativa, entonces para un entero positivo $n = \prod_{i=1}^{k} p_i^{k_i}$ , tenemos

f(n) = \prod_{i=1}^{k} f(p_i)^{k_i}

Ejemplos

Funciones multiplicativas comunes son

Función suma de divisores: $\sigma_k(n) = \sum_{d|n} d^k$ , que representa la suma de las $k$ -ésimas potencias de los divisores de $n$ . Observar que $\sigma_k(n)$ y $σ^k(n)$ son distintas.
Función conteo de divisores: $\tau(n) = \sigma_0(n) = \sum_{d|n} 1$ , que representa el conteo de divisores de $n$ , también denotada $d(n)$ .
Función suma de divisores: $\sigma(n) = \sigma_1(n) = \sum_{d|n} d$ , que representa la suma de los divisores de $n$ .
Función φ de Euler: $\varphi(n) = \sum_{i=1}^n [(n,i)=1] \cdot 1$ , que representa el conteo de enteros positivos menores o iguales que $n$ y coprimos con $n$ . Además, $\sum_{i=1}^n [(n,i)=1] \cdot i =$ $\frac{n\varphi(n) + [n=1]}{2}$ , $\varphi(n)$ es par.
Función de Möbius: $\mu(n)$ , que sirve como el inverso multiplicativo de la función identidad en la convolución de Dirichlet, $\mu(1) = 1$ , para un número libre de cuadrados $n = \prod_{i=1}^t p_i$ , $\mu(n) = (-1)^t$ , y para un número con factores cuadrados, $\mu(n) = 0$ .
Función unidad: $e(n) = [n=1]$ , que sirve como el elemento identidad en la convolución de Dirichlet, completamente multiplicativa.
Función constante: $I(n) = 1$ , completamente multiplicativa.
Función identidad: $id(n) = n$ , completamente multiplicativa.
Función potencia: $id^k(n) = n^k$ , completamente multiplicativa.

Las dos fórmulas clásicas respecto de la función de Möbius y la función de Euler son:

$[n=1] = \sum_{d|n} \mu(d)$ , interpretando $\mu(d)$ como los coeficientes del principio de inclusión-exclusión lo demuestra.
$n = \sum_{d|n} \varphi(d)$ . Para demostrarlo, podemos contar el número de ocurrencias de $\frac{1}{n}(1 \leq i \leq n)$ en su forma de fracción más simple.

Recursos

Como este módulo pretende servir como una introducción suave a este tema, usualmente describiremos solo la solución más simple que pasa las restricciones dadas. Soluciones con complejidades temporales asintóticamente mejores se pueden encontrar en blogs como los siguientes, aunque hay que tener en cuenta que pueden no ser mucho más rápidas bajo las restricciones dadas.

Recursos
Fuente	Recurso	Notas
CF	[Tutorial] Math note — Möbius inversion
CF	[Tutorial] Math note — Dirichlet convolution
CF	Dirichlet convolution. Part 1: Fast prefix sum computations
CF	Catalog	ver la sección de teoría de números para posts de blog sobre temas relacionados

Calentamiento

Empecemos introduciendo un conjunto de números $Q_N$ con el que a menudo trabajamos al calcular sumas de prefijos de funciones multiplicativas hasta $N$ .

Hecho	Fuente	Nombre	Dificultad	Tags	Solución
	YS	Enumerate Quotients	Fácil		en el módulo

Ejercicio para el lector: ¿qué tan grande puede ser este conjunto?

Respuesta

Sea $s=\lfloor \sqrt N\rfloor$ . Entonces el conjunto puede tener tamaño a lo sumo $2\sqrt N$ ya que cada elemento es o bien a lo sumo $s$ , o de la forma $\lfloor N/i\rfloor$ para algún $i\le s$ .

Una propiedad importante de este conjunto es que si $x\in Q_N$ entonces $Q_x \subseteq Q_N$ .

¿Por qué?

\lfloor \lfloor N/i\rfloor / j\rfloor =\lfloor N/(ij)\rfloor

Implementación

Código

Complejidad temporal: $O(\sqrt N)$

C++


#include <bits/stdc++.h>
using namespace std;
 
template <class T> using V = vector<T>;
#define all(x) begin(x), end(x)
 
using ll = long long;
 
int main() {
	ios::sync_with_stdio(false);
	cin.tie(nullptr);
 
	ll N;
	cin >> N;
 
	auto square = [&](ll s) { return s * s; };
 
	ll s = 1;
	while (square(s + 1) <= N) ++s;
 
	V<ll> a;
	for (ll i = 1; i <= s; ++i) a.push_back(i);
	for (ll i = s; i >= 1; --i) {
		if (i == s && a.back() == N / i) continue;
		a.push_back(N / i);
	}
 
	int k = size(a);
	cout << k << "\n";
	for (int i = 0; i < k; ++i) cout << a[i] << " \n"[i + 1 == k];
}

Ejemplo - Suma de divisores

Hecho	Fuente	Nombre	Dificultad	Tags	Solución
	CSES	Sum of Divisors	Fácil		Solución

Pista: la complejidad temporal es la misma que la del problema anterior.

Solución

Explicación

No es factible calcularlo de forma directa, pero podemos reescribir la suma así:

\sum_{i=1}^{n}\sigma(i)=\sum_{i=1}^{n}\sum_{j=1}^{n}[j|i]\cdot j=\sum_{i=1}^{n}i\cdot\sum_{j=1}^{n}[i|j]=\sum_{i=1}^{n}i\cdot\left\lfloor\frac{n}{i}\right\rfloor

Cuando $i \leq \sqrt{n}$ , hay solo $O(\sqrt{n})$ valores distintos para $\left\lfloor\frac{n}{i}\right\rfloor$ . De forma similar, cuando $i > \sqrt{n}$ , $\left\lfloor\frac{n}{i}\right\rfloor < n$ tiene solo $O(\sqrt{n})$ valores distintos. Para un $\left\lfloor\frac{n}{i}\right\rfloor$ fijo, los valores de $i$ forman un intervalo contiguo, que es $\left[\left\lfloor\frac{n}{\left\lfloor\frac{n}{i}\right\rfloor+1}\right\rfloor+1,\left\lfloor\frac{n}{\left\lfloor\frac{n}{i}\right\rfloor}\right\rfloor\right]$ . Podemos calcular cada una de estas contribuciones en $O(\sqrt n)$ .

Notas adicionales:

La suma del número de divisores de los primeros $n$ enteros positivos se puede calcular de la misma manera.
$\sum_{i=1}^{n}\left\lfloor\frac{n}{i}\right\rfloor\cdot i=\sum_{i=1}^{n}\left\lfloor\frac{n}{i}\right\rfloor\cdot\frac{(\left\lfloor\frac{n}{i}\right\rfloor+1)}{2}$ . Esto es cierto porque ambas cuentan $\sum_{ij\le n}i$ , donde la primera suma recorre $i$ y la segunda recorre $j$ .

Implementación

Código

Complejidad temporal: $\mathcal{O}(\sqrt{n})$

C++

Usamos la clase Modint de AtCoder para simplificar la implementación.


#include <atcoder/modint>
#include <bits/stdc++.h>
 
using namespace std;
using namespace atcoder;
 
using mint = modint1000000007;
 
template <class T> using V = vector<T>;
#define all(x) begin(x), end(x)
 
using ll = long long;
 
int main() {
	ios::sync_with_stdio(false);
	cin.tie(nullptr);
 
	ll n;
	cin >> n;
 
	const mint i2 = mint(1) / 2;
	auto arith = [&](mint l, mint r) { return (l + r) * (r - l + 1) * i2; };
 
	mint ans = 0;
	for (ll l = 1; l <= n;) {
		ll q = n / l;
		ll r = n / q;
		ans += arith(l, r) * mint(q);
		l = r + 1;
	}
 
	cout << ans.val() << "\n";
}

Convolución de Dirichlet

Introducción

La convolución de Dirichlet de funciones aritméticas $f$ y $g$ se define como $(f*g)(n) = \sum_{d|n} f(d) \cdot g(\frac{n}{d})$ . La convolución de Dirichlet cumple conmutatividad, asociatividad y distributividad respecto de la suma. Existe una función identidad $e(n) = [n=1]$ tal que $f*e = f = e*f$ 。Si $f$ y $g$ son funciones multiplicativas, entonces $f*g$ también es multiplicativa.

Una técnica común con la convolución de Dirichlet involucra lidiar con la convolución de una función $f$ y la función identidad $I$ . Por ejemplo, si $n = \prod_{i=1}^{t} p_i^{k_i}$ y $g=f* I$ , entonces $g(n) = \sum_{d|n} f(d)$ . Si $f$ es multiplicativa, entonces tenemos

g(n) = \prod_{i=1}^{t} \sum_{j=0}^{k_i} f(p_i^j).

Si queremos recuperar $f$ a partir de $g$ , podemos escribir $g* \mu = (f* I)* \mu=f* (I* \mu)=f* e=f$ . Es decir, $f(n) = \sum_{d|n} g(d) \cdot \mu\left(\frac{n}{d}\right)$ . Esto se conoce como inversión de Möbius.

Ejemplo - Convolución de Dirichlet y sumas de prefijos

Hecho	Fuente	Nombre	Dificultad	Tags	Solución
	YS	Dirichlet Convolution and Prefix Sums	Normal		en el módulo

Dados $\{f_x | x \in Q_N\}$ y $\{g_x | x \in Q_N\}$ , si definimos $h=f*g$ , podemos calcular $\{h_x | x \in Q_N\}$ en tiempo sublineal.

Solución

Explicación

Para $x\in Q_N$ podemos calcular $h_x$ en $O(\sqrt x)$ ya que si $ij\le x$ , entonces o bien $i\le \sqrt x$ o $j\le \sqrt x$ ( $i$ es chico o $j$ es chico). Podemos contabilizar ambos casos chicos y luego restar su superposición.

Además, $\sum_{x\in Q_N}\sqrt x\le O(\sum_{i=1}^{\sqrt N}\sqrt {N/i})\le O(N^{3/4})$ .

Bonus: es posible reducir esto a $O(N^{2/3})$ (ver este comentario ).

Implementación

Código

Complejidad temporal: $O(N^{3/4})$

C++


#include <atcoder/modint>
#include <bits/stdc++.h>
 
using namespace std;
using namespace atcoder;
 
using mint = modint998244353;
 
template <class T> using V = vector<T>;
#define all(x) begin(x), end(x)
 
using ll = long long;
 
ll sq(ll x) { return x * x; }
 
// BeginCodeSnip{Dirichlet Convolution}
struct Dirichlet {
	ll N;
	int s = 0, k;
	Dirichlet(ll N_) : N(N_) {
		while (sq(s + 1) <= N) ++s;
		k = 2 * s - (s == N / s);
	}
	ll pos_to_val(int i) {
		if (i < s) return i + 1;
		return N / (k - i);
	}
	int val_to_pos(ll l) {
		if (l <= s) return l - 1;
		return k - N / l;
	}
	V<mint> convolve(const V<mint> &f, const V<mint> &g) {
		V<mint> h(k);
		for (int i = 0; i < k; ++i) {
			ll v = pos_to_val(i);
			ll j = 1;
			for (; j * j <= v; ++j) {
				int p = val_to_pos(v / j);
				h[i] += (f[j - 1] - (j == 1 ? 0 : f[j - 2])) * g[p];
				h[i] += (g[j - 1] - (j == 1 ? 0 : g[j - 2])) * f[p];
			}
			--j;
			assert(j > 0);
			h[i] -= f[j - 1] * g[j - 1];
		}
		return h;
	}
};
// EndCodeSnip{}
 
void solve() {
	ll N;
	cin >> N;
	Dirichlet dc(N);
	int k = dc.k;
	V<mint> f(k), g(k);
	for (auto &t : f) {
		int v;
		cin >> v;
		t = v;
	}
	for (auto &t : g) {
		int v;
		cin >> v;
		t = v;
	}
	auto ret = dc.convolve(f, g);
	for (int i = 0; i < k; ++i) { cout << ret.at(i).val() << " \n"[i + 1 == k]; }
}
 
int main() {
	ios::sync_with_stdio(false);
	cin.tie(nullptr);
	int T;
	cin >> T;
	while (T--) solve();
}

Ejemplo - Inverso de Dirichlet y sumas de prefijos

Supongamos que en lugar de hallar $h=f*g$ dados $f$ y $g$ , queremos recuperar $g$ dados $h$ y $f$ .

Hecho	Fuente	Nombre	Dificultad	Tags	Solución
	YS	Dirichlet Inverse and Prefix Sums	Normal		en el módulo

Solución

Explicación

Similar al problema anterior: iteramos sobre $x\in Q_N$ en orden creciente, pero en lugar de derivar $h_x$ de $f_{1\dots x}$ y $g_{1\dots x}$ podemos derivar $g_x$ de $f_{1\dots x}$ , $g_{1\dots x-1}$ y $h_x$ .

Implementación

Código

C++

Complejidad temporal: $O(N^{3/4})$


#include <atcoder/modint>
#include <bits/stdc++.h>
 
using namespace std;
using namespace atcoder;
 
using mint = modint998244353;
 
template <class T> using V = vector<T>;
#define all(x) begin(x), end(x)
 
using ll = long long;
 
ll sq(ll x) { return x * x; }
 
// BeginCodeSnip{Dirichlet Inverse}
struct Dirichlet {
	ll N;
	int s = 0, k;
	Dirichlet(ll N_) : N(N_) {
		while (sq(s + 1) <= N) ++s;
		k = 2 * s - (s == N / s);
	}
	ll pos_to_val(int i) {
		if (i < s) return i + 1;
		return N / (k - i);
	}
	int val_to_pos(ll l) {
		if (l <= s) return l - 1;
		return k - N / l;
	}
	V<mint> invert(const V<mint> &f, const V<mint> &h) {
		// return g s.t. f * g = h
		V<mint> g(k);
		mint inv_f0 = mint(1) / f.front();
		for (int i = 0; i < k; ++i) {
			mint remainder = h.at(i);
			if (i > 0) {
				ll v = pos_to_val(i);
				ll j = 1;
				for (; j * j <= v; ++j) {
					int p = val_to_pos(v / j);
					if (j > 1) remainder -= (f[j - 1] - (j == 1 ? 0 : f[j - 2])) * g[p];
					remainder -= (g[j - 1] - (j == 1 ? 0 : g[j - 2])) * f[p];
				}
				--j;
				assert(j > 0);
				remainder += f[j - 1] * g[j - 1];
			}
			g.at(i) = remainder * inv_f0;
		}
		return g;
	}
};
// EndCodeSnip
 
void solve() {
	ll N;
	cin >> N;
	Dirichlet dc(N);
	int k = dc.k;
	V<mint> f(k);
	for (auto &t : f) {
		int v;
		cin >> v;
		t = v;
	}
	vector<mint> h(k, 1);
	auto ret = dc.invert(f, h);
	for (int i = 0; i < k; ++i) { cout << ret.at(i).val() << " \n"[i + 1 == k]; }
}
 
int main() {
	ios::sync_with_stdio(false);
	cin.tie(nullptr);
	int T;
	cin >> T;
	while (T--) solve();
}

Ejemplo - Suma de totientes

Hecho	Fuente	Nombre	Dificultad	Tags	Solución
	YS	Sum of Totient Function	Normal		en el módulo

Solución

Explicación

Sabemos que $\phi* I=id$ . Las sumas de prefijos de $I$ e $id$ son conocidas, así que todo lo que necesitamos hacer es aplicar la solución de “Inverso de Dirichlet y sumas de prefijos” con $f=I$ y $h=id$ .

Implementación

Código

C++

Complejidad temporal: $O(N^{3/4})$


#include <atcoder/modint>
#include <bits/stdc++.h>
 
using namespace std;
using namespace atcoder;
 
using mint = modint998244353;
 
template <class T> using V = vector<T>;
#define all(x) begin(x), end(x)
 
using ll = long long;
 
ll sq(ll x) { return x * x; }
 
// BeginCodeSnip{Dirichlet Inverse}
struct Dirichlet {
	ll N;
	int s = 0, k;
	Dirichlet(ll N_) : N(N_) {
		while (sq(s + 1) <= N) ++s;
		k = 2 * s - (s == N / s);
	}
	ll pos_to_val(int i) {
		if (i < s) return i + 1;
		return N / (k - i);
	}
	int val_to_pos(ll l) {
		if (l <= s) return l - 1;
		return k - N / l;
	}
	V<mint> invert(const V<mint> &f, const V<mint> &h) {
		// return g s.t. f * g = h
		V<mint> g(k);
		mint inv_f0 = mint(1) / f.front();
		for (int i = 0; i < k; ++i) {
			mint remainder = h.at(i);
			if (i > 0) {
				ll v = pos_to_val(i);
				ll j = 1;
				for (; j * j <= v; ++j) {
					int p = val_to_pos(v / j);
					if (j > 1) remainder -= (f[j - 1] - (j == 1 ? 0 : f[j - 2])) * g[p];
					remainder -= (g[j - 1] - (j == 1 ? 0 : g[j - 2])) * f[p];
				}
				--j;
				assert(j > 0);
				remainder += f[j - 1] * g[j - 1];
			}
			g.at(i) = remainder * inv_f0;
		}
		return g;
	}
};
// EndCodeSnip
 
void solve() {
	ll N;
	cin >> N;
	Dirichlet dc(N);
	int k = dc.k;
	V<mint> f(k), h(k);
	mint i2 = mint(1) / 2;
	for (int i = 0; i < k; ++i) {
		mint v = dc.pos_to_val(i);
		f.at(i) = v;
		h.at(i) = v * (v + 1) * i2;
	}
	auto ret = dc.invert(f, h);
	cout << ret.back().val() << "\n";
}
 
int main() {
	ios::sync_with_stdio(false);
	cin.tie(nullptr);
	solve();
}

Bonus: se puede optimizar esto a $O(N^{2/3})$ aplicando el método introducido en el segundo recurso. Más específicamente, $g$ y sus sumas de prefijos se pueden precomputar hasta $N^{2/3}$ en $O(N^{2/3})$ usando una criba lineal. Luego calcular las sumas de prefijos restantes toma $O(\sum_{i=1}^{N^{1/3}}\sqrt{N/i})=O(N^{2/3})$ de tiempo adicional.

Problemas

Los primeros dos problemas son del primer recurso.

Fuente	Nombre	Dificultad	Solución
HDU	Function	Normal	Solución
SPOJ	Counting Divisors (square)	Difícil	Solución
YS	Counting Square-free Integers	Difícil	Solución