First note that in any optimal solution we can always have $x_1 = a_1$ as we wish to maximize distances, and $x_1 = a_1$ unconditionally does that without sacrificing distance anywhere.

Secondly note that while our 'smallest differences' metric we want to optimize appears a global problem, it really is not if we also sort our $x$. We can have the equivalent definition:

$$\Delta = \min_{1 \leq i < K} \left( x_{i+1} - x_i \right)$$

Now suppose our candidate solution has value $\Delta$. Then given $x_i$ we can find the leftmost $x_{i+1}$:

$$x_{i+1} = \min_y \{ y \mid y \in S \wedge y \geq x_i + \Delta\} $$

So we find $x_1, x_2, \dots$ until we are done. Repeatedly finding which interval in $S$ (if any) will contain $x_{i+1}$ can be done in a combined $O(K)$ time during this process if you just keep track of the highest used interval so far. If there is no valid interval for $x_{i+1}$, it means $\Delta$ is not a valid outcome, but if we find al values up until $x_K$ it is.

So now we have a $O(K + N)$ time test to see whether a certain $\Delta$ is valid. But we can say with certainty that $\Delta \in [0, b_K - a_1]$. So we can do a binary search on this interval using our test to find a range which contains our optimal $\Delta$.

Since each step in this binary search halves the size of the interval that contains the optimal result, we can find optimal $\Delta$ (and corresponding $x_i$) up to a precision of $2^{-p}(b_K - a_1)$ in $O(p(K+ N))$ time.

First note that in any optimal solution we can always have $x_1 = a_1$ as we wish to maximize distances, and $x_1 = a_1$ unconditionally does that without sacrificing distance anywhere.

Secondly note that while our 'smallest differences' metric we want to optimize appears a global problem, it really is not if we also sort our $x$. We can have the equivalent definition:

$$\Delta = \min_{1 \leq i < K} \left( x_{i+1} - x_i \right)$$

Now suppose our candidate solution has value $\Delta$. Then given $x_i$ we can find the leftmost $x_{i+1}$:

$$x_{i+1} = \min_y \{ y \mid y \in S \wedge y \geq x_i + \Delta\} $$

So we find $x_1, x_2, \dots$ until we are done. Repeatedly finding which interval in $S$ (if any) will contain $x_{i+1}$ can be done in a combined $O(K)$ time during this process if you just keep track of the highest used interval so far. If there is no valid interval for $x_{i+1}$, it means $\Delta$ is not a valid outcome, but if we find al values up until $x_K$ it is.

So now we have a $O(K + N)$ time test to see whether a certain $\Delta$ is valid. But we can say with certainty that $\Delta \in [0, b_K - a_1]$. So we can do a binary search on this interval using our test to find a range which contains our optimal $\Delta$.

Since each step in this binary search halves the size of the interval that contains the optimal result, we can find optimal $\Delta$ (and corresponding $x_i$) up to a precision of $2^{-p}(b_K - a_1)$ in $O(p(K+ N))$ time.

If you're a bit smarter you can terminate early and exactly as an optimal solution will always contain a limiting sequence of steps of the minimum distance $\Delta$ starting at the left part of an interval and ending at the right part of an interval. Thus we have $a_i + k\cdot \Delta = b_j$. In other words, if our remaining interval contains only a single value of the form $\frac{b_j - a_i}{k}$ for some $i \leq j$ and $k$, then that is the optimal value and you can terminate the search.