Computing the inverse geodesic length in planar graphs and graphs of bounded treewidth

TL;DR

This paper presents efficient algorithms for computing the inverse geodesic length, also known as the Harary index, in planar graphs and graphs with bounded treewidth, improving computational complexity for these classes.

Contribution

It introduces algorithms with sub-quadratic and near-linear time complexities for calculating the inverse geodesic length in specific graph classes, extending existing distance sum techniques.

Findings

Inverse geodesic length can be computed in roughly O(n^{9/5}) time for planar graphs.

In graphs with bounded treewidth, it can be computed in O(n log^{O(k)} n) time.

Techniques from distance sum computations and rational function evaluations are adapted for the inverse case.

Abstract

The inverse geodesic length of a graph $G$ is the sum of the inverse of the distances between all pairs of distinct vertices of $G$. In some domains it is known as the Harary index or the global efficiency of the graph. We show that, if $G$ is planar and has $n$ vertices, then the inverse geodesic length of $G$ can be computed in roughly $O(n^{9/5})$ time. We also show that, if $G$ has $n$ vertices and treewidth at most $k$, then the inverse geodesic length of $G$ can be computed in $O(n \log^{O(k)}n)$ time. In both cases we use techniques developed for computing the sum of the distances, which does not have "inverse" component, together with batched evaluations of rational functions.

Equations

$$\mathrm{IGL}(G)\text{leavevmode}=\text{leavevmode}\sum _{uv\in \left(\genfrac{}{}{0ex}{}{V}{2}\right)}\frac{1}{d_G(u,v)},$$

$$W(G)\text{leavevmode}=\text{leavevmode}\sum _{uv\in \left(\genfrac{}{}{0ex}{}{V}{2}\right)}{d}_G(u,v),\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\mathrm{diam}(G)\text{leavevmode}=\text{leavevmode}\max \{d_G(u,v)\mid uv\in (\genfrac{}{}{0ex}{}{V}{2})\}.$$

$$\displaystyle\sum_{u\in U}d_{G}(a,u)\leavevmode\nobreak\$$

$$\displaystyle\max_{u\in U}d_{G}(a,u)\leavevmode\nobreak\$$

$$\sum _{u\in U}\frac{1}{d_G(a,u)}\text{leavevmode}=\text{leavevmode}\sum _{u\in U}\frac{1}{d_G(a,s)+d_G(s,u)}$$

$${\rho }_{s,U}(x)\text{leavevmode}=\text{leavevmode}\sum _{u\in U}\frac{1}{x+d_G(s,u)}.$$

$$\operatorname{IGL}(G)\leavevmode\nobreak\ =\leavevmode\nobreak\ \frac{1}{2}\,\sum_{u\in V}\Bigl{(}\deg_{G}(u)+(n-1-\deg_{G}(u))\tfrac{1}{2}\Bigr{)}\leavevmode\nobreak\ =\leavevmode\nobreak\ \frac{n(n-1)}{2}-\frac{3|E(G)|}{4},$$

$$\mathrm{IGL}(A,B)\text{leavevmode}=\text{leavevmode}\sum _{a\in A,b\in B}\frac{1}{d_G(a,b)}.$$

$$\mathrm{IGL}(A,b)\text{leavevmode}=\text{leavevmode}\mathrm{IGL}(A,\{b\})\text{leavevmode}=\text{leavevmode}\sum _{a\in A}\frac{1}{d_G(a,b)}.$$

$$\forall a\in A,b\in B:\text{leavevmode}\text{leavevmode}\text{leavevmode}d(a,b)=\min \{d_G(a,s)+d_G(s,b)\mid s\in S\}.$$

$$\frac{P(x)}{Q(x)}\text{leavevmode}=\text{leavevmode}\sum _{i=1}^n{R}_i(x)\text{leavevmode}=\text{leavevmode}\sum _{i=1}^n\frac{P_i(x)}{Q_i(x)}.$$

$$P(x)\text{leavevmode}=\text{leavevmode}\sum _{i=1}^n\left(P_i(x)\cdot \prod _{j\ne i}{Q}_j(x)\right)\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{and}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}Q(x)\text{leavevmode}=\text{leavevmode}\prod _{i=1}^n{Q}_i(x).$$

$${\tilde{R}}_j(x)\text{leavevmode}:=\text{leavevmode}{R}_{2j-1}(x)+R_{2j}(x)\text{leavevmode}=\text{leavevmode}\frac{P_{2j-1}(x)Q_{2j}(x)+P_{2j}(x)Q_{2j-1}(x)}{Q_{2j-1}(x)Q_{2j}(x)}$$

$$\displaystyle O\left(\sum_{k=1}^{\lceil\log_{2}n\rceil}\left(\frac{n}{2^{k}}\cdot\bigl{(}2^{k}+T_{\rm prod}(2^{k})\bigr{)}\right)\right)\leavevmode\nobreak\$$

$$=\text{leavevmode}O(n)\cdot \sum _{k=1}^{\lceil {\log }_{2}n\rceil }(k\log k)$$

$$=\text{leavevmode}O(n{\log }^{2}n\mathrm{loglog}n).$$

$$\displaystyle O(1+m/n)O\left(T_{\rm prod}(d)\cdot\log d\right)\leavevmode\nobreak\$$

$$=\text{leavevmode}O((n+m){\log }^{2}n\mathrm{loglog}n).$$

$$R_i(x)\text{leavevmode}=\text{leavevmode}\frac{1}{d_G(u_i,v_i)+x}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{ for some ui,vi∈V(G)}$$

$$\sum _v|P(v)|\text{leavevmode}=\text{leavevmode}O(nd\cdot B(n,d)),$$

$$\mathrm{ISW}(R,\delta )\text{leavevmode}=\text{leavevmode}\sum _{p\in P\cap R}\frac{1}{\delta +\omega (p)}.$$

$$\sum _{i\in I}|J_i|\text{leavevmode}=\text{leavevmode}|\Pi |\text{leavevmode}=\text{leavevmode}\sum _{j\in [m]}|I_j|\text{leavevmode}=\text{leavevmode}\sum _{j\in [m]}O(2^{d}B(n,d))\text{leavevmode}=\text{leavevmode}O(2^{d}m\cdot B(n,d)).$$

$${\rho }_i(x)\text{leavevmode}=\text{leavevmode}\sum _{p\in P_i}\frac{1}{x+\omega (p)}.$$

$$O((|P_i|+|J_i|){\log }^2|P_i|\mathrm{loglog}|P_i|)\text{leavevmode}=\text{leavevmode}O((|P_i|+|J_i|){\log }^{2}n\mathrm{loglog}n)$$

$$\displaystyle\sum_{i\in I}O(|P_{i}|+|J_{i}|)\log^{2}n\operatorname{loglog}n)\leavevmode\nobreak\$$

$$=\text{leavevmode}\left(O(nd\cdot B(n,d))+O(2^{d}m\cdot B(n,d))\right)O({\log }^{2}n\mathrm{loglog}n)$$

$$=\text{leavevmode}O((n+m)2^{d}B(n,d){\log }^{2}n\mathrm{loglog}n).$$

$$\mathrm{ISW}(R_j,{\delta }_j)\text{leavevmode}=\text{leavevmode}\sum _{p\in P\cap R_j}\text{leavevmode}\frac{1}{{\delta }_j+\omega (p)}\text{leavevmode}=\text{leavevmode}\sum _{i\in I_j}\sum _{p\in P_i}\text{leavevmode}\frac{1}{{\delta }_j+\omega (p)}\text{leavevmode}=\text{leavevmode}\sum _{i\in I_j}\text{leavevmode}{\rho }_i({\delta }_j).$$

$$A(b,i)\text{leavevmode}:=\text{leavevmode}\{a\in A\mid$$

$$d_G(a,b)=d_G(a,s_i)+d_G(s_i,b),$$

$$d_G(a,b)\le d_G(a,s_j)+d_G(s_j,b)\text{leavevmode}\text{leavevmode}\text{leavevmode}\forall j>i\}.$$

$${\varphi }^{(i)}=({\varphi }_j^{(i)}{)}_{j\in [k]}:A$$

$$a$$

$$I^{(i)}_{j}(b)\leavevmode\nobreak\ =\leavevmode\nobreak\ \begin{cases}\bigl{(}-\infty,\,d_{G}(b,s_{j})-d_{G}(b,s_{i})\bigr{)}&\text{if $j<i$,}\\ \mathbb{R}&\text{if $j=i$,}\\ \bigl{(}-\infty,\,d_{G}(b,s_{j})-d_{G}(b,s_{i})\bigr{]}&\text{if $j>i$.}\end{cases}$$

$$\displaystyle d_{G}(a,b)\leavevmode\nobreak\$$

Full text

Keywords: IGL, distances in graphs, planar graphs, bounded treewidth, algebraic tools

1 Introduction

Let $G=(V,E)$ be an undirected graph with $n$ vertices and (abstract) positive edge-lengths $\lambda\colon E\rightarrow\mathbb{R}_{>0}$. The length of a walk in $G$ is the sum of the edge-lengths along the walk. The distance between two vertices $u$ and $v$ of $G$, denoted by $d_{G}(u,v)$, is the minimum length over all paths in $G$ from $u$ to $v$. A particularly important case is when $\lambda(e)=1$ for all edges $e\in E$. In this case, the distance between two vertices $u$ and $v$ is the minimum number of edges over the $u$-$v$ paths in $G$. (Usually this is called the unweighted case or unit-length case.)

In this paper we are interested in the inverse geodesic length of $G$, defined as

[TABLE]

where $\binom{V}{2}$ denotes all the unordered pairs of vertices of $G$. In this definition we use the convention that $1/d_{G}(u,v)=0$ when there is no path from $u$ to $v$. As a consequence, $\operatorname{IGL}(G)$ is equal to the sum of the inverse geodesic length over its connected components.

The inverse geodesic length has been considered in different contexts. In Chemical Graph Theory it is considered for unit-length edges and called the Harary index of graph $G$, a common topological index [11, 24]; sometimes it is defined as twice $\operatorname{IGL}(G)$. In the context of Network Analysis, the inverse geodesic length goes under the name of efficiency, a concept introduced by Latora and Marchiori [20, 21] that has been heavily used. In this context, $\operatorname{IGL}(G)$ is usually normalized dividing by $\binom{n}{2}$, that is, the average inverse geodesic length is considered. The author learned the concept with the algorithmic work of Gaspers and Lau [14] and refers to the more careful discussion about related work therein.

The problem is related to but different from the computation of the Wiener index $W(G)$ and the diameter $\operatorname{diam}(G)$, defined as

[TABLE]

Our contribution.

We present algorithms to efficiently compute $\operatorname{IGL}(G)$ for two standard graph classes. The challenge to obtain efficient algorithms is to avoid computing the distance between all pairs of vertices. More precisely, we show the following results, where $n$ denotes the number of vertices of the graph:

• •

For each fixed $k\geq 1$, the inverse geodesic length of graphs with treewidth at most $k$ can be computed in $O(n\log^{k+2}n\operatorname{loglog}n)$ time. Thus, the running time is near-linear for graphs of bounded treewidth.

• •

When we consider the running time parameterized by $n$ and the treewidth $k$, the inverse geodesic length can be computed in $n^{1+\varepsilon}2^{O_{\varepsilon}(k)}$ time, for any $\varepsilon>0$.111The constants hidden in the $O_{\varepsilon}$-notation depend on $\varepsilon$.

• •

For planar graphs, the inverse geodesic length can be computed in $\tilde{O}(n^{9/5})$ time222Here and in the rest of the paper we use the notation $\tilde{O}$ to hide logarithmic factors..

To achieve these results we build on the techniques used to compute the Wiener index and the diameter in planar graphs [8, 16], namely the use of additively-weighted Voronoi diagrams in planar graphs, and in graphs with small treewidth [1, 7, 9], namely the use of orthogonal range searching.

Previous papers computing the Wiener index or the diameter rely on the fact that the length of a shortest path is the sum of the lengths of subpaths that compose it. This is not true for the inverse of distances. More precisely, those works decompose the problem into subproblems described by a triple $(a,s,U)$, where $a$ and $s$ are vertices, $U$ is a subset of vertices, and the shortest path from $a$ to each vertex of $U$ goes through vertex $s$. Moreover, a preprocessing step is used to construct a data structure so that $|U|$, $\sum_{u\in U}d_{G}(s,u)$ and $\max_{u\in U}d_{G}(s,u)$ are obtained in sublinear time (in $|U|$). This preprocessing pays off because there are several subproblems with triples of the form $(\cdot,s,U)$ for the same $s$ and $U$. Then we use that

[TABLE]

There is no such simple decomposition of

[TABLE]

that would be useful for $\operatorname{IGL}(G)$. That is the main obstacle to apply the techniques.

Our main new technique is to associate to $(\cdot,s,U)$ a rational function

[TABLE]

Then, for the triple $(a,s,U)$, we are interested in evaluating $\rho_{s,U}(d_{G}(a,s))$. We use computer algebra to make such evaluations at several values in near-liner time, as done by Aronov, Katz and Moroz [3, 23] in completely different settings. The key insight is that $\rho_{s,U}(x)$ can be computed efficiently because, after expansion, it is a rational function whose numerator and denominator have degree $|U|$. Note that all the evaluations of $\rho_{s,U}(x)$ are computed simultaneously, which is also a change in the approach compared to previous works for the Wiener index or the diameter. Thus, we do not use a data structure to handle each triple $(a,s,U)$, but we treat all the triples of the form $(\cdot,s,U)$ batched.

We believe that the use of tools from Computational Geometry and Computer Algebra for solving graph-distance problems, such as computing the $\operatorname{IGL}$, shows an appealing interaction between different areas.

Comparison to related work.

The efficient computation of $\operatorname{IGL}$ for graphs with small treewidth was considered by Gaspers and Lau [14]. They show that $\operatorname{IGL}(G)$ for graphs with unit-length edges and treewidth at most $k$ can be computed in $2^{O(k)}n^{3/2+\varepsilon}$. In the case of unit-length trees the running time is reduced to $O(n\log^{2}n)$. Their approach is based on computing the distance distribution: for any positive integer $P$, one can compute the number of pairs of vertices at distance $i$ for $i=1,\dots,P$ in $2^{O(k)}n^{1+\varepsilon}\sqrt{P}$ time. In the case of trees the distance distribution can be computed $O(n\log n+P\log^{2}n)$ time. From the distance distribution it is easy to compute the $\operatorname{IGL}(G)$ in linear time.

We improve the result of Gaspers and Lau in two ways. First, we reduce the degree of the main polynomial term of the running time from roughly $3/2$ to $1$ for graphs of small treewidth. Second, our results work for arbitrary (positive) weighted graphs, not only unit length (or small integers). For the case of unit-length trees, their result is stronger than ours. Our algorithm does not compute the distance distribution. In fact, for arbitrary weights it would be too large. Their algorithm for computing the distance distribution relies on the fast product of polynomials. They do not use rational functions, which is a key insight in our approach. They use the degree to encode distances, so the edge lengths must be small positive integers, while in our approach the degree of the rational functions depends (implicitly) on the number of vertices.

For graphs with small treewidth, our algorithm to compute the $\operatorname{IGL}$ has a few more logarithmic factors than the algorithms to compute the Wiener index or the diameter [7, 9]. This slight increase in the running time is due to the manipulation and evaluation of rational functions.

We are not aware of any previous algorithmic result for computing the $\operatorname{IGL}$ of planar graphs. Since the distance between all pairs of vertices in a planar graph can be computed in $O(n^{2})$ time [13, 17], the $\operatorname{IGL}$ of planar graphs can be computed in $O(n^{2})$ time. We provide the first algorithm computing the $\operatorname{IGL}$ of planar graphs in subquadratic time, even in the special case of unit-length edges.

For planar graphs, the running time of computing the $\operatorname{IGL}$ is $\tilde{O}(n^{9/5})$, while the Wiener index and the diameter can be computed in $\tilde{O}(n^{5/3})$ time with the algorithm of Gawrychowski et al. [16]. The difference is quite technical and difficult to explain at this level of detail. For readers familiar with that work, we can point out that, to apply our new approach, the dynamic tree used in [16] to encode the bisectors should be combined with the rational functions. The natural way to do this would be to associate a rational function to each node of the dynamic tree, and that leads to an overhead of $O(r)$ when working with a piece of an $r$-division. A more careful treatment of the dynamic tree may lead to an improvement, but we do not see how. The approach we use to manipulate (the inverse of) the distances for each piece is similar to the one used in Cabello [8], where the bisectors are computed and manipulated explicitly. It should be noted that, at this point, the bottleneck in our algorithm is not the computation of the Voronoi diagram, but the time to manipulate the outcome of the Voronoi diagram.

For arbitrary graphs, there is no constant $\delta_{0}>0$ such that the $\operatorname{IGL}$ can be computed in $O(n^{2-\delta_{0}})$ time, unless the strong exponential time hypothesis (SETH) fails. This holds also for sparse graphs. Indeed, Roditty and Vassilevska Williams [25] show that, for arbitrary graphs with $n$ vertices and $O(n)$ edges, one cannot compute the diameter in $O(n^{2-\delta_{0}})$ time, for some constant $\delta_{0}>0$, unless the SETH fails. In fact, their proof shows that for undirected, unweighted graphs we cannot distinguish in $O(n^{2-\delta_{0}})$ time between sparse graphs that have diameter $2$ or larger. Since a unit-length graph $G=(V,E)$ has diameter $2$ if and only if

[TABLE]

we cannot compute the $\operatorname{IGL}$ in $O(n^{2-\delta_{0}})$ time, unless SETH fails.

Computation model.

We assume a model of computation where each arithmetic operation in the ring generated by the edge-lengths takes constant time and can be carried out exactly. Since we are using the Fast Fourier Transform (FFT), the actual running time depends on whether we assume that the ring has primitive roots of the unit available. The assumptions are explained in more detail in Section 2.1.

The output and several intermediate computations are represented as fractions, that is, as pairs $(a,b)$ of numbers representing the fraction $a/b$. The numerators and the denominators of these fractions are in the ring generated by the edge-lengths. Thus, if we assume that the input edge-lengths are integers, the fractions will have integer numerators and denominators, if the input edge-lengths are rationals, the fractions will have rational numerators and denominators, and if the edge-lengths are arbitrary real numbers, then the fractions will have numerators and denominators with real numbers.

The model of computation used by Gaspers and Lau [14] is Word RAM, which is a weaker model and thus a stronger result. They can use this model because they assume graphs with small integer edge-weights and they compute, for each possible distance, the number of pairs of nodes at such distance. After this, one still has to compute $\operatorname{IGL}(G)$ explicitly and it is not obvious how to bound the bit-length of the computation. In our approach, we should also bound the bit-length of intermediate integers appearing through the computation. This is not straightforward. To get a feeling of this endeavor, note for example that we should understand the number of bits needed to write expressions like the Harmonic numbers $H_{n}=\sum_{i=1}^{n}\tfrac{1}{i}$ exactly as a single fraction $a_{n}/b_{n}$.

Roadmap.

We provide some common preliminaries in the next section. In Section 3 we provide the algorithm for graphs with small or bounded treewidth, while in Section 4 we provide the algorithm for planar graphs. Sections 3 and 4 are independent and can be read in any order. The structure of the exposition in both cases is parallel. We suggest to the readers to start with the section dedicated to graphs that are closer to their expertise.

2 Preliminaries

For each positive integer $k$ we use the notation $[k]=\{1,\dots,k\}$.

Let $G$ be a graph with positive edge lengths. Since the graph $G$ under consideration will always be clear from the context, we will often drop the dependency on $G$ from the notation. (We do keep using $d_{G}(\cdot,\cdot)$ for the distance.)

Let $A$ and $B$ be disjoint subsets of vertices in $G$. We define

[TABLE]

For the special case when one of the sets has a single element, for example $B=\{b\}$, we write

[TABLE]

A separation in $G$ is a partition of $V(G)$ into pairwise disjoint sets of vertices $A,B,S$ such that there is no edge of $G$ with one vertex in $A$ and the other vertex in $B$. The set $S$ is the separator. We will systematically use $a$ for vertices in $A$, $b$ for vertices in $B$, and $s$ for vertices in $S$. Note that any path from $a\in A$ to $b\in B$ must necessarily pass through a vertex of $S$. In particular,

[TABLE]

We use $\sqcup$ to denote the disjoint union. Thus, if we write $A=\sqcup_{i\in I}A_{i}$, this means that the sets in $\{A_{i}\mid i\in I\}$ are pairwise disjoint and $A=\cup_{i\in I}A_{i}$.

2.1 Algebraic computations.

A polynomial $P(x)=\sum_{i=0}^{d}a_{i}x^{i}$ given in coefficient representation is the list of coefficients $(a_{0},\dots,a_{d})$ of length $1+\deg(P)$. A rational function $R(x)$ is a function $R(x)=\frac{P(x)}{Q(x)}$ where $P(x)$ and $Q(x)$ are polynomials. The coefficient representation of such rational function is just the coefficient representation of $P(x)$ and $R(x)$. The following result is obtained by using fast multiplication of polynomials followed by fast multipoint evaluation of polynomials. We refer to the book by von zur Gathen and Gerhard [15, Chapters 8 and 10] for a comprehensive treatment of the tools of Computer Algebra that we will use.

Lemma 1** **(Corollary A.4 in Aronov and Katz [3];

see also Lemmas 6 and 7 in Moroz and Aronov [23]).

Given a set of $n$ rational functions $R_{i}(x)=P_{i}(x)/Q_{i}(x)$ of constant degree each, in coefficient representation, and a set of $m$ values $x_{j}$, one can compute the $m$ values $\sum_{i}R_{i}(x_{j})$ in time $O((n+m)\log^{2}n\operatorname{loglog}n)$.

We provide the proof to be able to discuss the model of computation in more detail.

Proof.

As a first step, we compute polynomials $P(x)$ and $Q(x)$ such that

[TABLE]

More precisely, we will compute

[TABLE]

We do this using a recursive algorithm. If $n=1$, we then have $P(x)=P_{1}(x)$ and $Q(x)=Q_{1}(x)$. If $n\geq 2$, we compute, for each $j=1,\dots,\lfloor n/2\rfloor$, the rational function

[TABLE]

in coefficient representation. If $n$ is odd we set $\tilde{R}_{\lceil n/2\rceil}(x)=R_{n}(x)$. We then compute recursively the rational function $\sum_{j=1}^{\lceil n/2\rceil}\tilde{R}_{j}(x)$ and return it.

To analyze the running time note that, after $k$ levels of the recursion, we have $O(n/2^{k})$ rational functions, each with a denominator and a numerator of degree $O(2^{k})$ and coefficients in the ring generated by the coefficients of $\cup_{i}\{P_{i}(x),Q_{i}(x)\}$. Let $T_{\rm prod}(d)$ be the time needed to multiply two polynomials of degree $d$ in coefficient representation, assuming that each arithmetic operation in a ring that contains the coefficients of the polynomials takes constant time. Using the Fast Fourier Transform we have $T_{\rm prod}(d)=O(d\log d\operatorname{loglog}d)$; see [15, Theorem 8.23]. The time needed to compute $\frac{P(x)}{Q(x)}\leavevmode\nobreak\ =\leavevmode\nobreak\ \sum_{i}R_{i}(x)$ is then

[TABLE]

Once we have the polynomials $P(x)$ and $Q(x)$ such that $\frac{P(x)}{Q(x)}\leavevmode\nobreak\ =\leavevmode\nobreak\ \sum_{i}R_{i}(x)$ in coefficient representation, we evaluate $P(x)$ and $Q(x)$ at $x=x_{j}$ for all $j$. Note that $P(x)$ and $Q(x)$ have degree $O(n)$. Evaluating a polynomial of degree $d$ at $d$ points takes $O\left(T_{\rm prod}(d)\cdot\log d\right)$ time; see [15, Corollary 10.8]. To evaluate $P(x)$ and $Q(x)$ at $x=x_{j}$ for all $j$, we make $\lceil m/n\rceil$ groups of $n$ values each, and evaluate at each group. The running time for this is

[TABLE]

Note: the $\operatorname{loglog}n$ term may seem superfluous, but this depends on the assumptions. If we assume that the ring supports the FFT, that is, we can manipulate primitive roots of the unity in constant time, and divisions can be carried out exactly, then $T_{\rm prod}(d)=O(d\log d)$; see [15, Theorem 8.18]. In this case the $\operatorname{loglog}n$ factor in the final running time disappears. If we are assuming that the coefficients and the values $x_{j}$ are arbitrary real numbers, this assumption would mean that arithmetic operations involving the complex roots of the unit $\cos(\pi/2^{k})+i\sin(\pi/2^{k})$ can be manipulated exactly in constant time. In our stated running times we are making the more conservative assumption that only operations in the ring generated by $\{d_{G}(u,v)\mid u,v\in V(G)\}$ are available, and divisions are not performed. ∎

All our use of computer algebra is encoded in Lemma 1. More precisely, we will be using Lemma 1 for rational functions of the form

[TABLE]

and we will be evaluating $\sum_{i}R_{i}(x)$ at values of the form $x=d_{G}(u_{j},v_{j})$ for some $u_{j},v_{j}\in V(G)$. Thus, looking into the proof of Lemma 1 and the proofs in [15, Chapters 8 and 10], we see that we are making arithmetic operations in the ring generated by $\{d_{G}(u,v)\mid u,v\in V(G)\}$, without computing inverses and without ever simplifying fractions. Thus, we carry exact results, assuming that arithmetic operations are performed exactly. The running times are expressed assuming that each such arithmetic operation in the ring takes constant time.

3 Bounded treewidth

In this section we provide the near-linear time algorithm for graphs with bounded treewidth. First we review what we need for the usual orthogonal range searching queries, and then we extend this to provide our new, key data structure for so-called inverse shifted queries. We describe how to compute the interaction across a separator, and finally apply the divide-and-conquer approach of Cabello and Knauer [9] to derive the final result.

3.1 Preliminaries on orthogonal range searching

In our algorithm for small treewidth we will use the range trees used for orthogonal range searching. They are explained in the textbook [10, Chapter 5] or the survey [2], but the usual analysis in Computational Geometry assumes constant dimension. We use the analysis by Bringmann et al. [7], which uses the dimension $d$ as a parameter and is based on the results of Monier [22].

A box in $\mathbb{R}^{d}$ is the product of intervals, some of them possibly infinite and some of them possibly open or closed. (We only consider such axis-parallel boxes in this paper.) We will use $B(n,d)=\binom{d+\lceil\log n\rceil}{d}$ to bound the time and size of the range trees. First we note a bound on $B(n,d)$, and then provide a multivariate bound on range trees. Note that the stated bounds do not assume that $d$ is constant.

Lemma 2** (Lemma 5 in Bringmann et al. [7]).**

$B(n,d)=O(\log^{d}n)$* and $B(n,d)=n^{\varepsilon}2^{O_{\varepsilon}(d)}$ for each $\varepsilon>0$.*

Theorem 3**.**

Given a set $P$ of $n$ points in $\mathbb{R}^{d}$, there is a family of sets $\mathcal{P}=\{P_{i}\mid i\in I\}$ and a data structure with the following properties:

• •

$P_{i}\subseteq P$* for each $P_{i}\in\mathcal{P}$;*

• •

all the sets of $\mathcal{P}$ together have $O(nd\cdot B(n,d))$ points, counting with multiplicity; that is, $\sum_{P_{i}\in\mathcal{P}}|P_{i}|=O(nd\cdot B(n,d))$;

• •

for each box $R\subset\mathbb{R}^{d}$, the data structure finds in $O(2^{d}B(n,d))$ time indices $I_{R}\subset I$ such that $|I_{R}|=O(2^{d}B(n,d))$ and $P\cap R=\bigsqcup_{i\in I_{R}}P_{i}$;

• •

the family $\mathcal{P}$ and the data structure can be computed in $O(nd\cdot B(n,d))$ time.

Proof.

We consider the following variant of range trees. We assume that the reader is familiar with range trees and provide only a sketchy description. We describe a tree $T(d,Q)$ that is defined recursively on the dimension $d$ and size of the point set $Q\subset\mathbb{R}^{d}$ under consideration. When the dimension $d$ is $1$, we make a balanced binary search tree $T(1,Q)$ that stores the elements of $Q$ at the leaves.

When $d>1$ and $|Q|>1$, we split $Q$ into two sets $Q_{\ell}$ and $Q_{r}$ of roughly the same size, such that each element of $Q_{\ell}$ has smaller $d$-coordinate than each element of $Q_{r}$. We recursively construct $T(d,Q_{\ell})$, $T(d,Q_{r})$ and $T(d-1,Q)$. Finally, we make a node, the root of $T(d,Q)$, whose left child is the root $T(d,Q_{\ell})$, its right child is the root of $T(d,Q_{r})$, and its associated data structure is $T(d-1,Q)$. When $d>1$ and $|Q|=1$, we make a node without children and with a pointer to the associated data structure $T(d-1,Q)$, which is built recursively.

The data structure is $T(d,P)$. For each node $v$ of $T(d,P)$, let $P(v)$ denote the set of points stored in $T(d,P)$ under $v$. This is a so-called canonical subset. (In fact it suffices to consider nodes $v$ in trees $T(1,Q)$, but that is not really relevant in our discussion.) For each node $v$ of $T(d,P)$, we add the canonical subset $P(v)$ to $\mathcal{P}$.

The analysis by Bringmann et al. [7, Lemma 6] shows that

[TABLE]

where the sum is over all the nodes $v$ of $T(d,P)$. This gives an upper bound on the size of the data structure and the construction time. (For the construction, instead of computing medians, as Bringmann et al. suggest, it is better to sort the points once in each dimension and pass to the recursive calls the points sorted in $d$ lists, one for each dimension.)

The analysis of Bringmann et al. for queries [7, Lemma 8] identifies for each box $R$ the nodes of the tree such that $P\cap R$ is the disjoint union of $O(2^{d}B(n,d))$ canonical subsets from $\mathcal{P}$ and they can be found in $O(2^{d}B(n,d))$ time. ∎

3.2 Inverse shifted queries

Let $P$ be a set of $n$ points in $\mathbb{R}^{d}$ and assume that each point $p$ of $P$ has a positive weight $\omega(p)$. An inverse shifted-weight query (ISW query) is specified by a pair $(R,\delta)$, where $R$ is an (axis-parallel) box and $\delta$ is a non-negative value, the shift. We want to compute

[TABLE]

Since $\omega(p)$ is positive and $\delta\geq 0$, no fraction has [math] in the denominator and $\operatorname{ISW}(R,\delta)$ is well defined.

We do not know how to answer ISW queries online, but we can solve the batched version.

Theorem 4**.**

Let $P$ be a set of $n$ points in $\mathbb{R}^{d}$ and assume that each point $p\in P$ has a positive weight $\omega(p)$. Consider $\operatorname{ISW}$ queries specified by $(R_{1},\delta_{1}),\dots,(R_{m},\delta_{m})$, where each $R_{j}$ is a box and $\delta_{j}\geq 0$. The values $\operatorname{ISW}(R_{1},\delta_{1}),\dots,\operatorname{ISW}(R_{m},\delta_{m})$ can be computed in $O((n+m)2^{d}\cdot B(n,d)\cdot\log^{2}n\operatorname{loglog}n)$ time.

Proof.

We compute the family $\mathcal{P}=\{P_{i}\mid i\in I\}$ and the data structure given in Theorem 3. This takes $O(nd\cdot B(n,d))=O(n2^{d}\cdot B(n,d))$ time. For each $R_{j}$ we use the data structure to compute the set $I_{j}$ of $O(2^{d}\cdot B(n,d))$ indices such that $P\cap R_{j}=\bigsqcup_{i\in I_{j}}P_{i}$. For all $j\in[m]$ together this takes $O(m2^{d}\cdot B(n,d))$ time.

Define the pairs $\Pi=\{(i,j)\in I\times[m]\mid i\in I_{j}\}$ and, for each $i\in I$, let $J_{i}=\{j\in[m]\mid i\in I_{j}\}$. Note that $I_{j}$ are the fibers of $\Pi$ when we fix the second coordinate, while $J_{i}$ are the fibers of $\Pi$ when we fix the first coordinate. From the data structure we can compute $\Pi$ and thus also the fibers $J_{i}$ for all $i$. We have

[TABLE]

For each $i\in I$ we define the rational function

[TABLE]

For each $i\in I$, we use Lemma 1 to evaluate $\rho_{i}(\delta_{j})$ for each $j\in J_{i}$. This means that, for each $i\in I$, we spend

[TABLE]

time. For all $i\in I$ together we spend a total time of

[TABLE]

For each query $(R_{j},\delta_{j})$, since $P\cap R_{j}=\bigsqcup_{i\in I_{j}}P_{i}$, we have

[TABLE]

Since the values $\rho_{i}(\delta_{j})$ have been already computed for all $(i,j)\in\Pi$, the computation of $\operatorname{ISW}(R_{j},\delta_{j})$ for all $j$ takes additional $O(|\Pi|)=O(2^{d}m\cdot B(n,d))$ time. ∎

3.3 Across a separation

Let $G$ be a graph with $n$ vertices and let $A,B,S$ be a separation of $G$. Let $k$ be the size of the separator; thus $k=|S|$. Our objective here is to compute $\operatorname{IGL}(A,B)$ efficiently. For the rest of the section, we assume that $G$ and the partition $A,B,S$ are fixed.

We fix an order on the vertices of the separator $S$ and denote them by $s_{1},\dots,s_{k}$. For each $b\in B$, we want to partition $A$ into groups depending on which vertices of $S$ belong to a shortest path from $b$. To obtain a partition, we have to be careful to assign each vertex of $A$ to one group. For this we use the minimum index of the vertex in the separator that is contained in some shortest path. More precisely, for each $b\in B$ and each $i\in[k]$ we define the set

[TABLE]

It is easy to see that the sets $A(b,i)$, $i\in[k]$, form a partition of $A$ for each $b\in B$.

Consider some index $i\in[k]$. Define the map $\phi^{(i)}=(\phi^{(i)}_{j})_{j\in[k]}$ by

[TABLE]

Note that the $i$-th coordinate of $\phi^{(i)}(a)$ is always [math], and thus it does not provide any information. We keep it to maintain slightly simpler notation.

For each vertex $b\in B$, we define the box $R^{(i)}(b)=I^{(i)}_{1}(b)\times\dots\times I^{(i)}_{k}(b)\subset\mathbb{R}^{k}$, where $I^{(i)}_{j}(b)$ is

[TABLE]

Lemma 5**.**

A vertex $a\in A$ belongs to $A(b,i)$ if and only if $\phi^{(i)}(a)\in R^{(i)}(b)$.

Proof.

Vertex $a$ belongs to $A(b,i)$ if and only if

[TABLE]

Since each path from $a$ to $b$ passes through some vertex of $S$, this can be rewritten as

[TABLE]

Rearranging terms this is equivalent to

[TABLE]

This last condition is precisely the condition for $\phi^{(i)}(a)\in R^{(i)}(b)$. ∎

Lemma 6**.**

Given a graph $G$ with $n$ vertices and $m$ edges, and a separation $A,B,S$ of $G$, where $S$ is the separator and has size $|S|=k\geq 2$, we can compute $\operatorname{IGL}(A,B)$ in $O(km+n2^{k}k\cdot B(n,k-1)\cdot\log^{2}n\operatorname{loglog}n)$ time.

Proof.

Let $s_{1},\dots,s_{k}$ be an enumeration of the vertices of $S$. We compute shortest path trees from each vertex $s_{i}$ of $S$. This takes $O(k(m+n\log n))$ time. After this, each distance $d_{G}(s_{i},v)$, where $s_{i}\in S$ and $v\in V(G)$, is available in constant time.

Consider a fixed index $i\in[k]$. We compute the set of points $P^{(i)}=\phi^{(i)}(A)$ and assign to point $\phi^{(i)}(a)$ weight $\omega(\phi^{(i)}(a))=d_{G}(s_{i},a)$. For each $b\in B$ we compute the description of the box $R^{(i)}(b)$. This takes $O(kn)$ time, for a single $i\in[k]$.

Since for each $b\in B$ and each $a\in A(b,i)$ the shortest path from $b$ to $a$ goes through $s_{i}$, we have

[TABLE]

Recall that $P^{(i)}$ is effectively a point set in $k-1$ dimensions because the $i$-th coordinate is identically [math]. Using Theorem 4 for $P^{(i)}$ we can compute $\operatorname{IGL}(b,A(b,i))=\operatorname{ISW}(R^{(i)}(b),d_{G}(s_{i},b))$ for all points $b\in B$ together in

[TABLE]

time. We repeat the procedure for each $i\in[k]$, which adds a multiplicative factor of $O(k)$ to the running time. After this, we have the values $\operatorname{IGL}(b,A(b,i))$ for all $i\in[k]$ and $b\in B$.

Finally we use that, for each $b\in B$, the sets $A(b,i)$, $i\in[k]$, form a partition of $A$ to obtain

[TABLE]

Since the $O(k|B|)=O(kn)$ values used in the last expression are already available, the result follows. ∎

Note that in the proof of the previous Lemma we can actually compute the values $\operatorname{IGL}(b,A)$ for all $b\in B$ within the same time bound.

3.4 Final algorithm

We use the recursive approach of Cabello and Knauer [9]. The idea is that graphs of small treewidth have small, balanced separators with the property that adding edges between all pairs of vertices of the separator does not increase the treewidth. Given such a separation $A,B,S$ of the graph we use Lemma 6 to compute $\operatorname{IGL}(A,B)$, we also compute $\operatorname{IGL}(S,V(G)\setminus S)$ using shortest paths from $S$, and then solve recursively the problems for $A$ and $B$. For the recursive problems, we add edges between the vertices of the separator $S$ whose length is equal to the length of the shortest path in $G$. Thus, the distances between vertices in $A$ and $B$ are correct in the recursive calls. Also, adding those edges does not increase the treewidth of the graphs used in the recursion. A bit of care is needed to avoid that distances between vertices of $S$ are counted more than once; it is easy to handle this by subtracting the terms that are counted twice.

Now we have two regimes depending on whether we want to assume that the treewidth is constant, as done in [9], or whether we want to consider the treewidth a parameter, as done in [7]. This difference affects the time to find a tree decomposition and a balanced separator. In both cases we use that an $n$-vertex graph with treewidth $k$ has $O(kn)$ edges [5].

Theorem 7**.**

Let $k\geq 2$ be an integer constant. For graphs $G$ with $n$ vertices and treewidth at most $k$, we can compute $\operatorname{IGL}(G)$ in $O(n\log^{k+2}n\operatorname{loglog}n)$ time.

Proof.

We obtain a tree decomposition of width $k$ in $O(n)$ time using the algorithm of Bodlaender [4], which takes $O(2^{O(k^{3}\log k)}n)=O(n)$ time, when $k$ is a constant. From the tree decomposition we can obtain a separation $A,B,S$, such that $A$ and $B$ contain $\Theta(n)$ vertices each, the separator $S$ has size $k$, and adding edges between the vertices of $S$ does not increase the treewidth; see [9, Lemma 3]. We use Lemma 6 to compute $\operatorname{IGL}(A,B)$. Using the estimate of Lemma 2 and that the number of edges is $m=O(kn)=O(n)$, we spend

[TABLE]

time to compute $\operatorname{IGL}(A,B)$. Because of the divide-and-conquer approach and because each $A$ and $B$ have a constant fraction of the vertices, the recursion increases the running time by another logarithmic factor. ∎

Theorem 8**.**

For graphs $G$ with $n$ vertices and treewidth at most $k$, we can compute $\operatorname{IGL}(G)$ in $n^{1+\varepsilon}2^{O_{\varepsilon}(k)}$ time, for any $\varepsilon>0$.

Proof.

Bodlaender et al. [6] give an algorithm that, for graphs of treewidth at most $k$, finds a tree decomposition of width $3k+4$ in $2^{O(k)}n\log n$ time. Given such a decomposition, we can obtain a separation $A,B,S$, such that $A$ and $B$ contain $\Theta(n-k)$ vertices each, the separator $S$ has size $k^{\prime}=(3k+4)+1=O(k)$, and adding edges between the vertices of $S$ does not increase the treewidth; see for example [5, Theorem 19]. We use Lemma 6 to compute $\operatorname{IGL}(A,B)$. Using the estimate of Lemma 2 and that the number of edges is $m=O(kn)$, we spend

[TABLE]

time to compute $\operatorname{IGL}(A,B)$, for each $\varepsilon>0$. The recursive calls add another logarithmic factor, which is absorbed by the polynomial term $n^{1+\varepsilon}$. (Actually, the logarithmic factor does not appear because the exponent of the polynomial is a constant strictly larger than $1$ and that term dominates the recursive formulation of the running time.) ∎

4 Planar graphs

First we provide some tools for planar graphs, including $r$-divisions, duality and Voronoi diagrams. Then we introduce a data structure to compute inverse shifted queries; this is our essential, new contribution. Finally we explain how to combine the data structure with pieces of an $r$-division, and how the whole algorithm works.

In the following we assume that $G=(V,E)$ is a fixed planar (or plane) graph and often drop the dependency on $G$ in the notation.

4.1 Preliminaries on planar graphs

A planar graph with $n$ vertices has $O(n)$ edges. An embedding of a planar graph in the plane can be computed in linear time. A plane graph is a planar graph together with a fixed embedding. The embedding is usually described combinatorially using rotation systems: it specifies for each vertex the circular ordering (clockwise, say) of the edges incident to the vertex. Adding edges of sufficiently large length to $G$, so that the distances between vertices are not changed, we can assume that $G$ is triangulated, that is, all faces of $G$ are triangles. Such addition of edges, without introducing parallel edges, can be done in linear time. Henceforth, we will assume that the planar graph $G$ is triangulated, and thus $3$-connected. (When $n=3$, the graph is not $3$-connected, but the case $n=O(1)$ is trivial.)

Let $G$ be a plane graph. A piece of $G$ is a subgraph of $G$ without isolated vertices. We assume in $P$ the embedding inherited from $G$. The boundary of $P$, denoted by $\partial P$, is the set of vertices of $P$ that are incident to some edge of $E\setminus E(P)$. The interior of $P$ is $V(P)\setminus\partial P$. Each path from a vertex of $P$ to a vertex of $V\setminus V(P)$ passes through some vertex of the boundary, $\partial P$. Thus, $V\setminus V(P)$, $V(P)\setminus\partial P$, $\partial P$ is a separation of $G$. (It may be very unbalanced.) A hole in a piece $P$ is a face of $P$ that is not a face of $G$. Note that each boundary vertex of a piece $P$ must be incident to some hole.

For any positive integer $r$, an $r$-division with few holes of an $n$-vertex plane graph $G$ is a collection $\mathcal{P}=\{P_{1},\dots,P_{k}\}$ of pieces of $G$ with the following properties:

• •

each edge of $G$ belongs to precisely one piece $P_{i}\in\mathcal{P}$;

• •

each piece $P_{i}\in\mathcal{P}$ has at most $r$ vertices, $O(\sqrt{r})$ boundary vertices, and $O(1)$ holes;

• •

$k=|\mathcal{P}|=O(n/r)$, that is, there are $O(n/r)$ pieces in $\mathcal{P}$.

The first item implies that, if a vertex $v\in\partial P_{i}\in\mathcal{P}$, then $v$ belongs to the boundary of each piece that contains it. Indeed, if such vertex $v$ belongs to two distinct pieces $P_{i},P_{j}\in\mathcal{P}$, then $v$ is incident to an edge of $P_{i}$ and an edge of $P_{j}$ (pieces do not have isolated vertices), and thus belongs to the boundary of $P_{j}$ and the boundary of $P_{i}$. As a consequence, the interiors of distinct pieces are disjoint. In particular, the sets $V(P_{1})\setminus\partial P_{1}$, $\dots$, $V(P_{k})\setminus\partial P_{k}$ and $\bigcup_{i\in[k]}\partial P_{i}$, form a pairwise disjoint partition of $V$.

Klein, Mozes and Sommer [18] give an algorithm to compute in $O(n)$ time an $r$-division with few holes, given a plane, triangulated graph $G$ and a value $r$. A simpler algorithm with running time $O(n\log n)$ was already known before [12, 19].

Let $G$ be a plane graph and denote by $F(G)$ its set of faces. The dual (multi)graph of $G$, denoted as $G^{*}$, has vertex set $V(G^{*})=F(G)$ and edge set

[TABLE]

We use $e^{*}$ to denote the edge of $G^{*}$ dual to the edge $e\in E(G)$. The embedding of $G$ defines naturally an embedding of the dual graph $G^{*}$, as follows. Each face $f\in F(G)$ is represented by a point $p_{f}$ chosen inside the face $f$, and the dual edge $fg=e^{*}$ is drawn as a curve with endpoints $p_{f}$ and $p_{g}$ that crosses the embedding of $G$ exactly once, namely in the interior of the drawing of $e$. See Figure 1. In general plane graphs, the dual graph is a multigraph because two faces can share several edges. To keep the notation simpler, we will talk about $G^{*}$ as a simple graph; adapting it to multigraphs is straightforward.

Let $\gamma$ be a cycle in the dual graph $G^{*}$ and let $s$ be a vertex of $G$. Let $A$ and $B$ be the two connected components of $\mathbb{R}^{2}\setminus\gamma$, and assume without loss of generality that $s$ belongs to $A$. We define the exterior of $\gamma$ with respect to $s$ as $\operatorname{ext}_{G}(\gamma,s)=V(G)\cap B$. See Figure 2. Thus, the exterior of $\gamma$ (with respect to $s$) is the set of vertices of $G$ such that any path connecting them to $s$ must intersect $\gamma$. While in the plane there is a usual meaning of exterior and interior of a closed Jordan curve, our definition is the natural parallel to the concept in the $2$-dimensional sphere and will be more convenient in our application.

4.2 Voronoi diagrams in plane graphs

The following material is taken from Gawrychowski et al. [16], where there is a comprehensive treatment. See also Cabello [8] for the initial treatment, where only a particular case is discussed.

Let $G$ be a plane graph with (abstract) edge lengths and let $P$ be a piece of $G$. (In fact, the graph $G$ is irrelevant for the forthcoming discussion; we use $P$ and $G$ to be consistent with the use later.) Let $S$ be a subset of vertices of $P$, the sites, and let $\delta\colon S\rightarrow\mathbb{R}_{\geq 0}$ assign non-negative weights to the sites. Usually we only talk about $S$ and treat $\delta$ as implicitly given with $S$. The (additively-weighted) Voronoi diagram of $P$ with respect to $S$ (and $\delta$) is the family of subsets of vertices

[TABLE]

Here, $\operatorname{cell}_{P}(s,S)$ is the Voronoi cell for $s$. Note that we are using the distances in $P$; the graph $G$ is irrelevant. A Voronoi cell may be the empty set.

It is convenient that the cells of the Voronoi diagram form a partition of $V(P)$. For this, we can break ties in several consistent ways. For example, we may index the sites and say that, in case of ties, the vertex gets assigned to the closest site with smallest index. (There are other valid options.) Formally, we would enumerate $S$ as $S=\{s_{1},\dots,s_{k}\}$ and define for each site $s_{i}$ the Voronoi cell

[TABLE]

It is easy to see that the Voronoi cell for $s\in S$ is star-shaped from $s$: for each vertex $v\in\operatorname{cell}_{P}(s,S)$, the vertices in any shortest path in $P$ from $s$ to $v$ are also contained in $\operatorname{cell}_{P}(s,S)$. In particular, $\operatorname{cell}_{P}(s,S)$ forms a connected subgraph of $P$. In the forthcoming claims we will assume (for simplicity) that each Voronoi cell is nonempty. Empty Voronoi cells have to be treated as a special, degenerate case.

Since each Voronoi cell is connected, each (non-empty) Voronoi cell can be described by a collection of cycles in the dual graph $P^{*}$, as follows. Here it is important that $P$ is a plane graph, not just planar. For each site $s\in S$, there is a family of cycles $H_{P}(s,S)$ in the dual graph $P^{*}$ such that

[TABLE]

See Figure 3. Recall that $\sqcup$ denotes that the sets in $\{\operatorname{ext}_{P}(\gamma,s)\mid\gamma\in H_{P}(s,S)\}$ are pairwise disjoint (for each single $s$). We call the family of cycles $H_{P}(s,S)$ the dual description of $\operatorname{cell}(s,S)$. (We use $H_{P}(s,S)$ because of their role as “holes”.) The existence of such cycles follows from the following fact for planar graphs $P$: the removal of $E^{\prime}\subset E(P)$ in $P$ leaves exactly two connected components if and only if the dual edges $(E^{\prime})^{*}=\{e^{*}\mid e\in E^{\prime}\}$ define a simple cycle in the dual graph $P^{*}$. More precisely, a dual edge $e^{*}$ belongs to some cycle of $H_{P}(s,S)$ if and only if $e$ has one vertex in $\operatorname{cell}_{P}(s,S)$ and one vertex outside $\operatorname{cell}_{P}(s,S)$. Thus, each connected component of $G-\operatorname{cell}_{P}(s,S)$ defines one cycle of $H_{P}(s,S)$. The behavior is parallel to what occurs with Voronoi diagrams in the plane for some non-Euclidean metrics; think for example of a Voronoi diagram in triangulated terrain, where some Voronoi cells may have holes. (Voronoi cells are not necessarily simply connected in such an scenario.) A Voronoi cell in such a Voronoi diagram is described by a collection of curves, where each curve describes a “hole”, that is, a part that does not belong to the Voronoi cell.

A bisector is the cycle in the dual graph separating the two (non-empty) Voronoi cells of a two-site Voronoi diagram. That is, given two additively weighted sites $s$ and $t$, we observe that $H_{P}(s,\{s,t\})$ consists of a single cycle because both $\operatorname{cell}_{P}(s,\{s,t\})$ and $\operatorname{cell}_{P}(t,\{s,t\})$ are connected. We denote such bisector as $\beta(s,t)$ and note that it consists of the set of edges dual to

[TABLE]

When we have multiple sites, the cycles in the dual representation $H_{P}(s,S)$ describing the boundary of $\operatorname{cell}(s,S)$ are made of portions of bisectors. (Portions in general, but some of them could be a whole bisector.) This is the case because, for each dual edge $e^{*}$ belonging to some cycle in $H_{P}(s,S)$, one vertex of $e$ belongs to $\operatorname{cell}_{P}(s,S)$ while the other vertex of $e$ lies in $\operatorname{cell}_{P}(t,S)$ for some $t\in S\setminus\{s\}$. Each cycle in $H_{P}(s,S)$ can be described by a sequence of patches, where each patch is a maximal portion of a bisector $\beta(s,t)$ for some site $t\in S\setminus\{s\}$. Such a patch is described by a tuple $(s,t,e_{1},e_{2})$, where $t\in S\setminus\{s\}$ and where $e_{1},e_{2}\in\beta(s,t)$ are the first and last edge, respectively, in the portion of the bisector (in some prescribed direction for $\beta(s,t)$). Thus, the description of one patch requires $O(1)$ space.

The combinatorial complexity of $H_{P}(s,S)$ or $\operatorname{cell}_{P}(s,S)$ is the number of patches required for all cycles in $H_{P}(s,S)$. The combinatorial complexity of the Voronoi diagram defined by the sites $S$ is the sum of the combinatorial complexities of $H_{P}(s,S)$ over all $s\in S$. Using properties of planar graphs one can see that the combinatorial complexity of Voronoi diagram of $S$ is $O(|S|)$. Here it is important to note that the size of the piece is not relevant. By the computation of the Voronoi diagram we mean the computation of the cycles $H_{P}(s,S)$ for all $s\in S$, where each cycle of each $H_{P}(s,S)$ is given as a sequence of patches. (In the description of the Voronoi diagram one usually considers also the adjacency relations between cells; we will not use such information explicitly in our algorithm, but it is implicit in the description of the patches because they contain information about which bisectors are being used.)

Cabello provided a randomized algorithm to compute Voronoi diagrams efficiently, after an expensive preprocessing, in a particular case when the sites are cofacial. Gawrychowski et al. improved the result in several aspects, allowing more general families of sites, making the preprocessing faster and the construction deterministic, as follows.

Theorem 9** (Gawrychowski et al. [16]).**

Let $P$ be a plane graph with $r$ vertices, let $S$ be a set of $b$ vertices in $P$, assume that the vertices of $S$ can be covered with $O(1)$ faces of $P$, and such covering with $O(1)$ faces is given. After $\tilde{O}(rb^{2})$ preprocessing time we can handle the following queries: for any given weights $\delta(s)$ assigned to each $s\in S$ at query time, the additively-weighted Voronoi diagram of $P$ with respect to $S$ can be computed in $\tilde{O}(b)$ time.

We will be using the result for pieces $P$ in an $r$-division of a plane graph $G$. In such case we will have $b=O(\sqrt{r})$, the preprocessing will take $\tilde{O}(r^{2})$ time, and the computation of each Voronoi diagram takes $\tilde{O}(\sqrt{r})$ time.

4.3 Inverse shifted queries

Let $P$ be a piece with $r$ vertices in a plane graph $G$, and let $s$ be a fixed vertex of $P$. Assume that each vertex $v$ of $P$ has a prescribed non-negative weight $\omega(v)$. As a warm up for the forthcoming technique, we provide the following result.

Lemma 10**.**

For $t$ given positive values $\delta_{1},\dots,\delta_{t}$, we can compute the values

[TABLE]

in $O((r+t)\log^{2}r\operatorname{loglog}r)$ time.

Proof.

Consider the rational function $\rho(x)=\sum_{v\in V(P)}\,\frac{1}{x+\omega(v)}$. We evaluate $\rho(x)$ at $x=\delta_{i}$ for all $i\in[t]$. Using Lemma 1, we can compute all these values in $O((r+t)\log^{2}r\operatorname{loglog}r)$ time. ∎

In general we will be considering similar sums as in the previous result, but for subsets of vertices defined by cycles in the dual graph. Recall that for a cycle $\gamma$ in the dual graph $P^{*}$ and a vertex $s$ in $P$, $\operatorname{ext}_{P}(\gamma,s)$ is the set of vertices of $P$ that are on the opposite side of $\gamma$ than $s$.

Given a cycle $\gamma$ in $P^{*}$, and a positive weight $\delta$ (the shift), we define the inverse shifted-weight (ISW) values

[TABLE]

Note that we are adding the inverse of $\delta+\omega(v)$, and we are only taking into account the vertices of $P$ outside $\gamma$. Since $\omega(v)$ is non-negative an $\delta$ is positive, the denominator is always positive, and thus the fractions are well-defined.

We are interested in the efficient computation of the values $\operatorname{ISW}_{P,s}$ for several pairs of the form $(\gamma,\delta)$. More precisely, we will handle the batched version: computing $\operatorname{ISW}_{P,s}(\gamma,\delta)$ for several given pairs $(\gamma,\delta)$, where $P$ and $s$ remain constant.

In our application we will consider only cycles $\gamma$ with a particular structure that we explain now. A cycle $\gamma$ in $P^{*}$ is $s$-star-shaped if the following holds: for each $u,v\in V(P)$, if $u$ belongs to the shortest path in $P$ from $s$ to $v$ and $u\in\operatorname{ext}_{P}(\gamma,s)$, then $v\in\operatorname{ext}_{P}(\gamma,s)$. (The intuition is that the interior of $\gamma$ with respect to $s$ has to be star-shaped from $s$.)

Let $\Xi(s)$ be the family of cycles in $P^{*}$ that are $s$-star-shaped. We will be interested in $\operatorname{ISW}_{P,s}$ queries for some cycles in $\Xi(s)$. We will assume that the cycles $\gamma$ in $\Xi(s)$ are oriented in such a way that $s$ lies in the region of $\mathbb{R}^{2}\setminus\gamma$ to the right of $\gamma$, as we walk along $\gamma$. Thus, if $s$ is in the bounded region of $\mathbb{R}^{2}\setminus\gamma$, then $\gamma$ is oriented clockwise, and otherwise $\gamma$ is oriented counterclockwise.

In the following we consider the dual graph $P^{*}$ as an oriented graph that contains both orientations of each edge. The word arc is used for oriented edges. For a primal arc $u\mathord{\shortrightarrow}v$ with face $a$ to the left and face $b$ to the right, we define its dual arc $(u\mathord{\shortrightarrow}v)^{*}=a\mathord{\shortrightarrow}b$. Thus, the dual arc crosses the primal arc from left to right. See Figure 4. For simplicity we assume that each arc has different endpoints. It is easy to enforce this or modify the discussion to the general case. Recall that $r$ is the number of vertices in $P$.

Lemma 11**.**

We can associate to each dual arc $a\mathord{\shortrightarrow}b$ of $P^{*}$ a set $U(a\mathord{\shortrightarrow}b)$ of vertices of $P$ with the following properties:

• •

For each oriented cycle $\gamma\in\Xi(s)$,

[TABLE]

• •

The sets $U(a\mathord{\shortrightarrow}b)$ for all $a\mathord{\shortrightarrow}b\in E(P^{*})$ together can be computed in $O(r^{2})$ time.

Proof.

Fix a shortest path tree $T_{s}$ from $s$ in $P$ and orient the edges as arcs away from the root $s$. For each arc $u\mathord{\shortrightarrow}v$ of $T_{s}$, let $T_{u\mathord{\shortrightarrow}v}$ be the subtree of $T_{s}-uv$ that contains $v$. For each arc $u\mathord{\shortrightarrow}v$ of $T_{s}$, we define $U((u\mathord{\shortrightarrow}v)^{*})=V(T_{u\mathord{\shortrightarrow}v})$. For all other arcs $a\mathord{\shortrightarrow}b$ of $P^{*}$ we define $U(a\mathord{\shortrightarrow}b)=\emptyset$. See Figure 5 for the definition and intuition in the forthcoming argument.

The description of the sets $U(\cdot)$ we gave readily leads to an algorithm to compute them with time complexity $O(r^{2})$ because $P$ has $O(r)$ edges.

Consider a cycle $\gamma\in\Xi(s)$, that is, $\gamma$ is $s$-star-shaped. Then the edges of $T_{s}$ crossed by $\gamma$ are not in any ancestor-descendant relation. That is, in any $s$-to-leaf path in $T_{s}$ there is at most one edge whose dual edge appears in $\gamma$. Indeed, if there would be two edges $e_{1}$ and $e_{2}$ of $T_{s}$ crossed by $\gamma$, then along the path in $T_{s}$ between $e_{1}$ and $e_{2}$ would get an alternation interior-exterior-interior with respect to $\gamma$ and $s$, and we would get a contradiction with the property that $\gamma$ is $s$-star-shaped.

Since the edges of $T_{s}$ crossed by $\gamma\in\Xi(s)$ are not in any ancestor-descendant relation, it follows from the construction that the sets $U(a\mathord{\shortrightarrow}b)$ for all $a\mathord{\shortrightarrow}b\in E(\gamma)$ are pairwise disjoint.

Consider any vertex $v\in V(P)$. If $v\in\operatorname{ext}_{P}(\gamma,s)$, then the path in $T_{s}$ from $s$ to $v$ crosses $\gamma$ exactly once, and thus $v$ belongs to exactly one set $U(a\mathord{\shortrightarrow}b)$ with $a\mathord{\shortrightarrow}b\in E(\gamma)$. If $v\notin\operatorname{ext}_{P}(\gamma,s)$, then the path in $T_{s}$ from $s$ to $v$ does not cross $\gamma$, and thus $v$ does not belong to any $U(a\mathord{\shortrightarrow}b)$ with $a\mathord{\shortrightarrow}b\in E(\gamma)$. ∎

Using ideas similar to those used in segment trees [10, Section 10.3] we obtain the following.

Lemma 12**.**

Assume that the sets $U(a\mathord{\shortrightarrow}b)$ of Lemma 11 are already computed. Let $\gamma$ be an oriented cycle in $\Xi(s)$. Assume that we are given $t$ pairs $(\pi_{1},\delta_{1}),\dots,(\pi_{t},\delta_{t})$, such that each $\pi_{i}$ is a subpath of $\gamma$ specified by the starting and ending arc, and each $\delta_{i}>0$. In $O((r+t)\log^{3}r\operatorname{loglog}r)$ time we can compute the values

[TABLE]

Proof.

Let $e_{1},\dots,e_{m}$ be the arcs along $\gamma$, and note that $m=O(r)$. We make a rooted balanced binary search tree $T$ with $m$ leaves such that the leaf storing $e_{i}$ has search key $i$. A left-to-right traversal of $T$ traverses the leaves containing $e_{1},\dots,e_{m}$ in that order. See Figure 6. Each node $z$ of $T$ represents a subpath of $\gamma$, denoted by $\gamma[z]$, which is the concatenation of arcs stored in leaves under $z$. Alternatively, we can define $\gamma[z]$ recursively: for each leaf $z$ that stores $e_{i}$, we define $\gamma[z]=e_{i}$, and for an internal node $z$ with left child $z^{\prime}$ and right children $z^{\prime\prime}$ the path $\gamma[z]$ is the concatenation of $\gamma[z^{\prime}]$ and $\gamma[z^{\prime\prime}]$. The paths $\gamma[z]$ for all nodes $z$ of $T$ are called canonical subpaths of $\gamma$.

For each node $z$ of $T$, let $r_{z}$ be $\sum_{e_{i}\in\gamma[z]}|U(e_{i})|$. Note that for each leaf $z$ of $T$ that stores $e_{i}$ we have $r_{z}=|U(e_{i})|$, while for each internal node $z$ of $T$ with children $z^{\prime}$ and $z^{\prime\prime}$ we have $r_{z}=r_{z^{\prime}}+r_{z^{\prime\prime}}$. Since the sets $U(e_{1}),\dots,U(e_{m})$ are pairwise disjoint because of Lemma 11 and canonical paths for nodes at the same level of the tree $T$ are arc-disjoint, we have $r_{z}\leq r$ for all $z$ and $\sum_{z\in V(T)}r_{z}=O(r\log m)=O(r\log r)$.

Each subpath $\pi_{i}$ of $\gamma$ is the union of $O(\log m)=O(\log r)$ canonical subpaths. That is, for each subpath $\pi_{i}$ there is a subset $Z_{i}$ of $O(\log r)$ nodes in $T$ such that $\pi_{i}$ is the concatenation of the canonical subpaths $\gamma[z]$, $z\in Z_{i}$. Furthermore, the nodes $Z_{i}$ describing those canonical subpaths can be identified in $O(\log r)$ time with a bottom-up traversal from the first and the last arcs of $\pi_{i}$. (Recall that we assume that $\pi_{i}$ is specified by the starting and ending arc.) This property is parallel to the property often used for segment trees [10, Lemma 10.10]. We compute in $O(t\log r)$ time the sets $Z_{1},\dots,Z_{t}$. For each node $z$ of $T$, let $I_{z}=\{i\in[t]\mid z\in Z_{i}\}$. The non-empty sets $I_{z}$ for all nodes $z$ together can be computed also in $O(t\log r)$ time. Note that

[TABLE]

For each node $z$ of $T$ we define the rational function

[TABLE]

Note that $\rho_{z}$ is the sum of

[TABLE]

rational functions, each of bounded degree.

For each node $z$ of $T$ we compute $\rho_{z}(\delta_{i})$ for each $i\in I_{z}$; that is, we evaluate $\rho_{z}(x)$ at $x=\delta_{i}$ for all $i\in I_{z}$. Using Lemma 1 this takes $O((|I_{z}|+r_{z})\log^{2}r_{z}\operatorname{loglog}r_{z})=O((|I_{z}|+r_{z})\log^{2}r\operatorname{loglog}r)$ time for each node $z$ of $T$. Thus, for all nodes $z$ of $T$ and all $i\in I_{z}$ together, we compute $\rho_{z}(\delta_{i})$ in time

[TABLE]

Note that for each given pair $(\pi_{i},\delta_{i})$, since $\pi_{i}$ is the concatenation of $\gamma[z]$, $z\in Z_{i}$, we have

[TABLE]

Since the values $\rho_{z}(\delta_{i})$ are already available, for each $i\in[t]$ we spend $O(|Z_{i}|)=O(\log r)$ additional time. The result follows. ∎

Here we have the final result that we will use.

Theorem 13**.**

Let $P$ be a piece with $r$ vertices in a plane graph $G$, and let $s$ be a vertex of $P$. Let $\Gamma=\{\gamma_{1},\dots,\gamma_{\ell}\}$ be a family of oriented cycles in $\Xi(s)$. Assume that we are given $t$ pairs $(\pi_{1},\delta_{1}),\dots,(\pi_{t},\delta_{t})$, where each $\delta_{i}>0$, each $\pi_{i}$ is an element of $\Xi(s)$, and each $\pi_{i}$ is given as a concatenation of $k_{i}$ subpaths of cycles from $\Gamma$. Set $k=\sum_{i}k_{i}$. We can compute $\operatorname{ISW}_{P,s}(\pi_{1},\delta_{1}),\dots\operatorname{ISW}_{P,s}(\pi_{t},\delta_{t})$ in $O(r^{2}+(k+\ell r)\log^{3}r\operatorname{loglog}r)$ time.

Proof.

Let $\Pi_{i}$ be the family of subpaths of cycles from $\Gamma$ that are given to describe $\pi_{i}$. Let $\Pi=\cup_{i}\Pi_{i}$ be all the subpaths used in the description; we consider $\Pi$ as a multiset. For each $\varpi\in\Pi$, let $i(\varpi)\in[t]$ be the index such that $\varpi\in\Pi_{i(\varpi)}$ and let $j(\varpi)\in[\ell]$ be the index such that $\varpi$ is a subpath of $\gamma_{j(\varpi)}$. Let $\Gamma_{j}=\{\varpi\in\Pi\mid j(\varpi)=j\}$, that is, the subpaths of $\gamma_{j}$ that appear in the description of some $\pi_{i}$. Note that

[TABLE]

First we compute the sets $U(a\mathord{\shortrightarrow}b)$ described in Lemma 11 for all dual arcs $a\mathord{\shortrightarrow}b$ in $P^{*}$. This takes $O(r^{2})$ time.

For each $j\in[\ell]$, we use Lemma 12 with the cycle $\gamma_{j}$ and the pairs $(\varpi,\delta_{i(\varpi)})$, for all $\varpi\in\Gamma_{j}$. This means that in time $O((r+|\Gamma_{j}|)\log^{3}r\operatorname{loglog}r)$ we compute the sums

[TABLE]

For all $j\in[\ell]$ together we spend

[TABLE]

time and we obtain $\operatorname{partsum}(\varpi,\delta_{i(\varpi)})$ for all $\varpi\in\Pi$. This means that we obtain $\operatorname{partsum}(\varpi,\delta_{i})$ for all $\varpi\in\Pi_{i}$ and all $i\in[t]$.

For each $(\pi_{i},\delta_{i})$, we use the properties of $U(a\mathord{\shortrightarrow}b)$ described in Lemma 11, namely that $\operatorname{ext}_{P}(\pi_{i},s)$ is the disjoint union of the sets $U(a\mathord{\shortrightarrow}b)$ over all arcs $a\mathord{\shortrightarrow}b\in E(\pi_{i})$, to obtain

[TABLE]

Thus, we can compute $\operatorname{ISW}_{P,s}(\pi_{i},\delta_{i})$ using an additional $O(|\Pi_{i}|)=O(k_{i})$ time, for each $i$. ∎

Lemma 10 and Theorem 13 will be combined to compute the sum of inverses over the vertices of a Voronoi cell using simple inclusion-exclusion; see the end of the proof of Lemma 14.

Remark.

In our discussion we have assumed that we want to add the shifted-inverse over vertices of $V(P)$ outside $\gamma$ with respect to a reference $s\in V(P)$. In fact, for any fixed $U\subset V(P)$, we could also be considering the inverse-shifted queries

[TABLE]

For example, if we would have red and blue vertices, we could consider sums over only the blue vertices.

Remark.

Similar results hold for arbitrary cycles, not just cycles from $\Xi(s)$, with worse running time. In Lemma 12 we have used that the sets $U(a\mathord{\shortrightarrow}b)$ are pairwise disjoint when the arcs $a\mathord{\shortrightarrow}b$ come from a cycle of $\Xi(s)$. For arbitrary cycles one can use inclusion-exclusion; see for example Cabello [8, Section 3.1] for the relevant idea. In such case, we can assume only that each $U(a\mathord{\shortrightarrow}b)$ has $O(r)$ vertices and, in the proof of Lemma 12, $\sum_{z\in V(T)}r_{z}=O(r^{2})$. The running time of Lemma 12 for arbitrary cycles becomes $O((r^{2}+t\log r)\log^{2}r\operatorname{loglog}r)$. Similarly, Theorem 13 holds for arbitrary cycles $\gamma_{1},\dots,\gamma_{\ell}$, but the running time becomes $O((k\log r+\ell r^{2})\log^{2}r\operatorname{loglog}r)$.

4.4 Across a separation

The next result handles the interaction between a piece $P$ of a plane graph $G$ and its complement. The boundary $\partial P$ of $P$ plays the role of the separation because each path from $V(P)\setminus\partial P$ to $V(G)\setminus V(P)$ must pass through $\partial P$. We use Voronoi diagrams for planar graphs, as Cabello [8] or Gawrychowski et al. [16], to group the vertices of $V(P)$ depending on which boundary vertex is used in the shortest path to them (for a fixed source). To collect the relevant data from each Voronoi cell, we use the data structure developed in Section 4.3.

Lemma 14**.**

Let $G$ be a plane graph with $n$ vertices and let $P$ be a piece of $G$ with $r$ vertices, $b$ boundary vertices, and $O(1)$ holes. In $\tilde{O}(bn+b^{2}r^{2})$ time we can compute

[TABLE]

Proof.

First we introduce some notation, express some properties, and then look into the algorithmic part. Let $A=V(G)\setminus V(P)$, $S=\partial P$ and $\Omega=A\times S$. Note that $\Omega$ has $O(bn)$ pairs. The set $S$ is a separator between $A$ and $V(P)\setminus S$. We will systematically use $a$ to index vertices of $A$, $s$ for vertices of $S$, and $v$ for vertices of $P$.

Consider one fixed $a\in A$ and assign weight $\delta(s)=d_{G}(s,a)$ to each $s\in S$. Consider the corresponding Voronoi diagram $\{\operatorname{cell}_{P}(s,S)\mid s\in S\}$, and recall that each Voronoi cell has a dual description: a family of cycles $H_{P}(s,S)$ in $P^{*}$ such that

[TABLE]

This Voronoi diagram depends on $a$. To make this dependency on $a$ explicit, because we will be considering several different vertices $a\in A$ together, we define $X(a,s)=\operatorname{cell}_{P}(s,S)$ and $H(a,s)=H_{P}(s,S)$. This notation drops the dependency on $P$ and $S$, which we may assume constant throughout this proof. On the other hand, the notation makes it clear that we consider a different Voronoi diagram for each $a\in A$, and each such Voronoi diagram has a cell for each $s\in S$ whose additive weight $\delta(s)$ was set to $d_{G}(a,s)$. With the new notation we have

[TABLE]

As we explained before, we may and we will assume that the Voronoi cells form a partition of $V(P)$. Therefore, for each $a\in A$, the sets $X(a,s)$, where $s\in S$, form a partition of $V(P)$. Since $S$ separates $A$ from $V(P)\setminus S$, it follows from the definition of additively-weighted Voronoi diagram that

[TABLE]

Note that in the second part we are using distances in $P$, which is consistent with the distance used for defining the Voronoi diagram in $P$. The values we want to compute can be rewritten as

[TABLE]

We will explain how to compute the values \operatorname{IGL}\bigl{(}a,X(a,s)\bigr{)} for all $(a,s)\in\Omega$ together in $\tilde{O}(bn+b^{2}r^{2})$ time. From this, the result follows, as the desired values can be computed in $O(|\Omega|)=O(bn)$ additional time.

For each two distinct boundary vertices $s,s^{\prime}$ of $S$, let $\tilde{\Gamma}(s,s^{\prime})$ be the family of bisectors between $s$ and $s^{\prime}$. Note that $\tilde{\Gamma}(s,s^{\prime})\subset\Xi(s)$ because they are bisectors. Define $\Gamma(s)=\bigcup_{s^{\prime}\in S,s^{\prime}\neq s}\tilde{\Gamma}(s,s^{\prime})$ for each $s\in S$, that is, all the bisectors in $P$ between $s$ and each possible other boundary vertex. Obviously, $\Gamma(s)\subset\Xi(s)$. Each cycle in each $H(a,s)$ is described by patches from $\Gamma(s)$. For each $(a,s)\in\Omega$, let $k(a,s)$ be the number of patches used to describe all the cycles of $H(a,s)$. Since the combinatorial complexity of each Voronoi diagram is linear in $|S|$, we have $\sum_{s\in S}k(a,s)=O(|S|)=O(b)$ for each $a\in A$.

A simple counting shows that the family $\Gamma(s)$ has in total $O(br^{2})$ arcs, counted with multiplicity. Indeed, each family $\tilde{\Gamma}(s,s^{\prime})$ consists of $O(r)$ cycles in the dual graph $P^{*}$, and each such cycle has $O(r)$ arcs. (Since the cycles of $\tilde{\Gamma}(s,s^{\prime})$ are nested and any two consecutive ones in the nested order differ by at least one vertex, the bound follows.) Thus, $\tilde{\Gamma}(s,s^{\prime})$ has a total of $O(r^{2})$ arcs, counted with multiplicity, which implies that $\Gamma(s)$ has in total $O(br^{2})$ arcs because $|S|=b$.

We move now to the computational part. We will consider the cycles in $H(a,s)$ for all pairs $(a,s)\in\Omega$, first from the perspective of fibers with constant $a$, and then from the perspective of fibers with constant $s$.

First we compute the relevant distances. From each $s$ in $S$ we compute two shortest path trees, one in $G$ and one in $P$. This takes $\tilde{O}(bn)$ time because $|S|=b$ and $G$ has $O(n)$ edges. Note that now, we have $d_{G}(s,u)$ for each $(s,u)\in S\times V(G)$ and $d_{P}(s,v)$ for each $(s,v)\in S\times V(P)$. We also compute the cycles in $\tilde{\Gamma}(s,s^{\prime})$ for all $s,s^{\prime}$. This can be done in $\tilde{O}(r^{2})$ for each $s,s^{\prime}\in S$ easily; see Cabello [8] for an explicit computation or Gawrychowski et al. [16] for a faster, implicit representation. The idea is that a vertex $v\in V(P)$ changes sides of the Voronoi cells when the additive weights of the sites satisfy $\delta(s)-\delta(s^{\prime})=d_{P}(s^{\prime},v)-d_{P}(s,v)$. We easily obtain also $\Gamma(s)$ for all $s\in S$. We have spent $\tilde{O}(bn+b^{2}r^{2})$ time.

We use Theorem 9 to preprocess $P$ in $\tilde{O}(rb^{2})$ time. We can do this because the vertices of $S$ are incident to the holes of $P$, and thus $S$ is covered by $O(1)$ faces of $P$, and those faces can be computed from the embedding of $G$ in $O(n)$ time. After this preprocessing, we can handle the following queries: for any given weights $\delta(s)$ assigned to each $s\in S$ at query time, the additively-weighted Voronoi diagram of $S$ in $P$ can be computed in $\tilde{O}(b)$ time. Note that this Voronoi diagram uses distances in $P$

We use the data structure to compute $H(a,s)$ for all $(a,s)\in\Omega$. More precisely, for each $a\in A$, we set the weights $\delta(s)=d_{G}(s,a)$ for all $s\in S$, and query the data structure once to obtain $H(a,s)$ for all $s\in S$. Since we are querying the data structure $|A|<n$ times, we obtain $H(a,s)$ for all $(a,s)\in\Omega$ in $\tilde{O}(bn)$ time.

Now we switch to the fibers with constant $s$. Consider one fixed $s\in S$. We assign to each vertex $v$ of $P$ weight $\omega(v)=d_{P}(s,v)$. Consider the family of cycles $\Gamma(s)$, which has $\ell=O(br)$ cycles of $\Xi(b)$ with a total of $O(br^{2})$ arcs, counted with multiplicity. Using Theorem 13 for all the pairs

[TABLE]

we obtain, for each cycle $\gamma\in H(a,s)$, the value

[TABLE]

in time

[TABLE]

Still for the same fixed $s$ and the same weights $\omega(\cdot)$, we use Lemma 10 with the values $\{\delta_{1},\dots,\delta_{|A|}\}=\{d_{G}(a,s)\mid a\in A\}$ to compute for each $a\in S$ the value

[TABLE]

Note that $d_{G}(a,s)+d_{P}(s,v)$ is not $d_{G}(a,v)$, in general. This takes $\tilde{O}(n+r)=\tilde{O}(n)$ time, for a fixed $s\in S$.

We repeat the procedure for each $s\in S$. In total, we spend

[TABLE]

time and we obtain the values

[TABLE]

for all $\gamma\in H(a,s)$ and all $(a,s)\in\Omega$.

Combining properties (1) and (2), we have for each $(a,s)\in\Omega$

[TABLE]

Thus, we can obtain each single \operatorname{IGL}\bigl{(}a,X(a,s)\bigr{)} in $O(1+|H(a,s)|)$ time, using the data already computed. Since $\sum_{s\in S}|H(a,s)|\leq\sum_{s\in S}k(a,s)=O(b)$, this computation for all $(a,s)\in\Omega$ together takes $O(bn)$ time. With this we have computed \operatorname{IGL}\bigl{(}a,X(a,s)\bigr{)} for all $(a,s)\in\Omega$, and we have finished the proof because of equation (3). ∎

4.5 Final algorithm

Using Lemma 14 for each piece in a $r$-division with a few holes, we obtain the final result.

Theorem 15**.**

Let $G$ be a planar graph with $n$ vertices and positive, abstract edge lengths. In $\tilde{O}(n^{9/5})$ time we can compute

[TABLE]

Proof.

First, we compute an $r$-division of $G$ with few holes. Such division has $k=O(n/r)$ pieces $\mathcal{P}=\{P_{1},\dots,P_{k}\}$, where each piece $P_{i}\in\mathcal{P}$ has $O(r)$ vertices, $O(\sqrt{r})$ boundary vertices, and $O(1)$ holes. Moreover, different pieces only intersect at boundary vertices. Such an $r$-division can be computed in $\tilde{O}(n)$ time [12, 18, 19].

Let $B$ be the set of vertices that are boundary vertices in some piece. Note that $B$ contains $O(n/r)\cdot O(\sqrt{r})=O(n/\sqrt{r})$ vertices. The vertex set $V(G)$ is the disjoint union of $B,V(P_{1})\setminus\partial P_{1},\dots,V(P_{k})\setminus\partial P_{k}$.

We compute a shortest path tree from each vertex of $B$. Each shortest path tree takes linear time [17], and thus we spend $O(n^{2}/\sqrt{r})$ time in this computation. After this, for each vertex $u$ of $G$ and each vertex $b$ of $B$, we can retrieve $d_{G}(u,b)$ in constant time. In particular, any sum of inverse of distances between pair of vertices, where are at least one vertex is in $B$, can be computed in time proportional to the number of terms in the sum. It follows that the desired values $\operatorname{IGL}(a,V(G)\setminus\{a\})$ can be computed now for all $a\in B$ with additional time $O(|B|\cdot|V(G)|)=O(n^{2}/\sqrt{r})$.

Consider one fixed piece $P_{i}\in\mathcal{P}$. We use Lemma 14 to compute in $\tilde{O}(r^{1/2}n+r^{3})$ time the values \operatorname{IGL}\bigl{(}a,V(P_{i})\bigr{)} for all $a\in V(G)\setminus V(P_{i})$. We also compute $\operatorname{IGL}(a,\partial P_{i})$ by adding the distances from $\partial P_{i}\subset B$; this takes $O(r^{1/2}n)$ time. With this we obtain the values

[TABLE]

Finally, we compute the distance (in $G$) between all pairs of vertices in $P_{i}$ in $\tilde{O}(r^{2})$ using usual machinery. For example, we may add edges between the $O((\sqrt{r})^{2})=O(r)$ pairs of boundary vertices $\partial P_{i}$ and compute distances in the resulting graph, which has $O(r)$ vertices and edges. Repeating this for the $O(n/r)$ pieces of $\mathcal{P}$ together we have spent $\tilde{O}(n^{2}/\sqrt{r}+nr^{2})$ time.

Recall that $V(G)$ is the disjoint union of $B,V(P_{1})\setminus\partial P_{1},\dots,V(P_{k})\setminus\partial P_{k}$. Thus, for each $a\in V(G)\setminus B$, there is a unique index $j(a)$ such that $a\in P_{j(a)}\in\mathcal{P}$, and therefore

[TABLE]

The $k-1=O(n/r)$ terms in the first sum have been computed, the $|B|=O(n/\sqrt{r})$ terms in the second sum have also been computed, and the $|V(P_{j(a)})\setminus\partial P_{j(a)}|=O(r)$ terms in the last sum have also been computed. We conclude that, for a single $a\in V(G)\setminus B$, we can recover $\operatorname{IGL}(a,V(G)\setminus\{a\})$ in $O((n/\sqrt{r})+r)$ time. Doing this for each vertex $a\in V(G)\setminus B$ takes $O(n^{2}/\sqrt{r}+nr)$.

The whole algorithm, as explained, takes $\tilde{O}(n^{2}/\sqrt{r}+nr^{2})$. Setting $r=n^{2/5}$ we obtain a running time of $\tilde{O}(n^{9/5})$. ∎

Corollary 16**.**

Let $G$ be a planar graph with $n$ vertices and positive, abstract edge lengths. In $\tilde{O}(n^{9/5})$ time we can compute $\operatorname{IGL}(G)$.

Proof.

We use the Theorem to compute $\operatorname{IGL}(a,V(G)\setminus\{a\})$ for all $a\in V(G)$. Then we note that

[TABLE]

∎

5 Conclusions

Let us remark a few extensions of the results we presented. The extensions are applicable for graphs of bounded treewidth and for planar graphs without affecting the asymptotic running times of each case.

Firstly, the algorithm can be adapted to work for directed planar graphs. In such case it is meaningful to compute

[TABLE]

We provide two clarifying remarks on how to achieve this. The first observation is that distances to a separator vertex can be computed by reversing the arcs and computing the distances from the separator vertex. The second observation is that Voronoi diagrams can also be considered for such "distances" defined in directed graphs; this has been done also in previous works [8, 16].

As a second extension, we note that we can also compute the sum of the distances between some marked vertices. This means that, for any given $U\subset V(G)$, we can compute in the same asymptotic running time the sum

[TABLE]

For the case of planar graphs, where this may be less obvious, see the remark after Theorem 13.

As a final extension, note that we can handle any rational function of constant degree that depends on the distances. For example, the same approach can be used for computing

[TABLE]

Indeed, this only affects the rational (or polynomial) functions that we have to consider and evaluate (at the shifts), but Lemma 1 keeps being applicable. On the other hand, current methods developed here or in previous works do not seem applicable to compute in subquadratic time values like

[TABLE]

even assuming a strong model of computation, like Real RAM, where sums of square roots can be manipulated in constant time.

Acknowledgments

I am grateful to the reviewers for pointing out that the model of computation should be clarified and several additional corrections.