Puntos de retículo dentro de un polígono que no es de retículo

Para polígonos de retículo existe la fórmula de Pick para enumerar los puntos de retículo dentro del polígono. ¿Qué hay de los polígonos con vértices arbitrarios?

Procesemos cada una de las aristas del polígono de forma individual, y después de eso podemos sumar las cantidades de puntos de retículo bajo cada arista considerando sus orientaciones para elegir un signo (como al calcular el área de un polígono usando trapecios).

Primero de todo debemos notar que si la arista actual tiene extremos en $A=(x_1;y_1)$ y $B=(x_2;y_2)$ entonces se puede representar como una función lineal:

$$y=y_1+(y_2-y_1) \cdot \dfrac{x-x_1}{x_2-x_1}=\left(\dfrac{y_2-y_1}{x_2-x_1}\right)\cdot x + \left(\dfrac{y_1x_2-x_1y_2}{x_2-x_1}\right)$$

$$y = k \cdot x + b,~k = \dfrac{y_2-y_1}{x_2-x_1},~b = \dfrac{y_1x_2-x_1y_2}{x_2-x_1}$$

Ahora realizaremos una sustitución $x=x’+\lceil x_1 \rceil$ de modo que $b’ = b + k \cdot \lceil x_1 \rceil$. Esto nos permite trabajar con $x_1’=0$ y $x_2’=x_2 - \lceil x_1 \rceil$. Denotemos $n = \lfloor x_2’ \rfloor$.

No sumaremos puntos en $x = n$ ni en $y = 0$ por la integridad del algoritmo. Se pueden agregar manualmente después. Así tenemos que sumar $\sum\limits_{x’=0}^{n - 1} \lfloor k’ \cdot x’ + b’\rfloor$. También asumimos que $k’ \geq 0$ y $b’\geq 0$. En caso contrario se debería sustituir $x’=-t$ y sumar $\lceil|b’|\rceil$ a $b’$.

Discutamos cómo podemos evaluar una suma $\sum\limits_{x=0}^{n - 1} \lfloor k \cdot x + b\rfloor$. Tenemos dos casos:

• $k \geq 1$ o $b \geq 1$.

  Entonces deberíamos empezar sumando puntos por debajo de $y=\lfloor k \rfloor \cdot x + \lfloor b \rfloor$. Su cantidad es igual a

  $$\sum\limits_{x=0}^{n - 1} \lfloor k \rfloor \cdot x + \lfloor b \rfloor=\dfrac{(\lfloor k \rfloor(n-1)+2\lfloor b \rfloor) n}{2}.$$

  Ahora solo nos interesan los puntos $(x;y)$ tales que $\lfloor k \rfloor \cdot x + \lfloor b \rfloor < y \leq k\cdot x + b$. Esta cantidad es la misma que el número de puntos tales que $0 < y \leq (k - \lfloor k \rfloor) \cdot x + (b - \lfloor b \rfloor)$. Así redujimos nuestro problema a $k’= k - \lfloor k \rfloor$, $b’ = b - \lfloor b \rfloor$ y tanto $k’$ como $b’$ son menores que $1$ ahora. Aquí hay una imagen: acabamos de sumar los puntos azules y restamos la función lineal azul de la negra para reducir el problema a valores más pequeños de $k$ y $b$:

• $k < 1$ y $b < 1$.

  Si $\lfloor k \cdot n + b\rfloor$ es igual a $0$, podemos devolver $0$ de forma segura. Si no es el caso, podemos decir que no hay puntos de retículo tales que $x < 0$ y $0 < y \leq k \cdot x + b$. Eso significa que tendremos la misma respuesta si consideramos un sistema de referencia nuevo en el que $O’=(n;\lfloor k\cdot n + b\rfloor)$, el eje $x’$ está dirigido hacia abajo y el eje $y’$ está dirigido hacia la izquierda. Para este sistema de referencia nos interesan los puntos de retículo en el conjunto

  $$\left{(x;y)~\bigg|0 \leq x < \lfloor k \cdot n + b\rfloor, 0 < y \leq \dfrac{x+(k\cdot n+b)-\lfloor k\cdot n + b \rfloor}{k}\right}$$

  lo que nos devuelve al caso $k>1$. Se puede ver el nuevo punto de referencia $O’$ y los ejes $X’$ e $Y’$ en la imagen de abajo:

Análisis de complejidad

Tenemos que contar a lo sumo $\dfrac{(k(n-1)+2b)n}{2}$ puntos. Entre ellos contaremos $\dfrac{\lfloor k \rfloor (n-1)+2\lfloor b \rfloor}{2}$ en el primer paso. Podemos considerar que $b$ es despreciablemente pequeño porque podemos empezar haciéndolo menor que $1$. En ese caso podemos decir que contamos alrededor de $\dfrac{\lfloor k \rfloor}{k} \geq \dfrac 1 2$ de todos los puntos. Así terminaremos en $O(\log n)$ pasos.

Implementación

Aquí hay una función simple que calcula el número de puntos enteros $(x;y)$ tales que $0 \leq x < n$ y $0 < y \leq \lfloor k x+b\rfloor$:

int count_lattices(Fraction k, Fraction b, long long n) {
    auto fk = k.floor();
    auto fb = b.floor();
    auto cnt = 0LL;
    if (k >= 1 || b >= 1) {
        cnt += (fk * (n - 1) + 2 * fb) * n / 2;
        k -= fk;
        b -= fb;
    }
    auto t = k * n + b;
    auto ft = t.floor();
    if (ft >= 1) {
        cnt += count_lattices(1 / k, (t - t.floor()) / k, t.floor());
    }
    return cnt;
}

Aquí Fraction es alguna clase que maneja números racionales. En la práctica parece que si todos los denominadores y numeradores son a lo sumo $C$ en valor absoluto entonces en las llamadas recursivas serán a lo sumo $C^2$ si se siguen dividiendo numeradores y denominadores por su máximo común divisor. Dada esta hipótesis podemos decir que se pueden usar doubles y exigir una precisión de $\varepsilon^2$ donde $\varepsilon$ es la precisión con la que se dan $k$ y $b$. Eso significa que en floor se deberían considerar números como enteros si difieren a lo sumo en $\varepsilon^2$ de un entero.