Fast way to calculate a few points in convolution?
I'll write names of time domain sequences in lower case and their frequency-domain counterparts in upper case, for example $X = \operatorname{DTFT}(x)$, where $\operatorname{DTFT}$ stands for discrete-time Fourier transform. I chose this transform, because it is unaffected by zero-padding of the time-domain sequences, allowing equations that mix frequency-domain counterparts of sequences of different lengths. With the upper-case frequency-domain symbols, time-domain convolution is denoted by multiplication, and a time-domain delay of $N$ samples is denoted by $z^{-N},$ where $z = e^{i\omega},$ with $\omega$ as continuous frequency. No frequency-domain thinking is required here. You can mentally translate the equations while reading to time-domain convolutions, delays, and additions. Addition in frequency domain is addition in time domain.

The frequency-domain counterparts $X$ of $Y$ of time-domain sequences $x$ and $x$ each of length $2N,$ with $N=500,$ can be written in partitioned form $X = X_0 + X_1z^{-N}$ and $Y = Y_0 + Y_1z^{-N},$ where $\{X_0, X_1, Y_0, Y_1\}$ are frequency-domain counterparts of time-domain partitions $\{x_0, x_1, y_0, y_1\}$ each of length $N.$ The partitions are non-overlapping sub-sequences of $x$ and $y$, stripped of delay which is then reintroduced in the equations. Convolution of $x$ and $y$ can be written in frequency domain as:

$$\begin{align}&amp;XY\\
= &amp;(X_0 + X_1z^{-N})(Y_0 + Y_1z^{-N})\\
= &amp;(X_0 + X_1z^{-N})Y_0 + (X_0 + X_1z^{-N})Y_1z^{-N}\\
= &amp;X_0Y_0 + X_1Y_0z^{-N} + X_0Y_1z^{-N} + X_1Y_1z^{-N}z^{-N}\\
= &amp;X_0Y_0 + (X_1Y_0 + X_0Y_1)z^{-N} + X_1Y_1z^{-2N}\end{align}\tag{1}$$

Now all convolutions between sequences are convolutions between partitions of length $N$. Each result of such a convolution is of length $2N-1$. Denoting the frequency-domain counterpart of the first $N$ samples of the convolution result by $H()$ as in "head" and the frequency-domain counterpart of the delay-stripped sub-sequence formed by the remaining $N-1$ samples by $T()$ as in "tail":

$$\begin{align}= &amp;H(X_0Y_0) + T(X_0Y_0)z^{-N}\\
+ &amp;\big(H(X_1Y_0) + T(X_1Y_0)z^{-N}\big)z^{-N} + \big(H(X_0Y_1) + T(X_0Y_1)z^{-N}\big)z^{-N}\\
+ &amp;\big(H(X_1Y_1) + T(X_1Y_1)z^{-N}\big)z^{-2N}\\
= &amp;H(X_0Y_0)\\
+ &amp;\big(T(X_0Y_0) + H(X_1Y_0) + H(X_0Y_1)\big)z^{-N}\\
+ &amp;\big(T(X_1Y_0) + T(X_0Y_1) + H(X_1Y_1)\big)z^{-2N}\\
+ &amp;T(X_1Y_1)z^{-3N}\end{align}\tag{2}$$

Now we have neatly split the full convolution result into length $N$ partitions. You were interested in the partition of length $N$ beginning at index $2N,$ or in frequency domain:

$$T(X_1Y_0) + T(X_0Y_1) + H(X_1Y_1)$$

The time-domain counterparts of these three terms can can be calculated by three FFT convolutions of length $2N$ compared to your single FFT convolution of length $4N.$ I don't know if that is of any use. You'd still be throwing away a total of $3N$ samples of convolution results, as originally. 

But wait! The time-domain counterpart of $T(X_0Y_1) + H(X_1Y_1)$ can be calculated by a circular FFT convolution of length $2N$ between $x$ and $y_1.$ First, writing that as linear (not circular) convolution:

$$\begin{align}&amp;XY_1\\
= &amp;X_0Y_1 + X_1Y_1z^{-N}\\
= &amp;H(X_0Y_1) + T(X_0Y_1)z^{-N} + H(X_1Y_1)z^{-N} + T(X_1Y_1)z^{-2N}\\
= &amp;H(X_0Y_1) + \big(T(X_0Y_1) + H(X_1Y_1)\big)z^{-N} + T(X_1Y_1)z^{-2N}\end{align}$$

Wrapping that around to make it circular, $z^{-2N}$ wraps to $z^{0} = 1:$

$$H(X_0Y_1) + T(X_1Y_1) + \big(T(X_0Y_1) + H(X_1Y_1)\big)z^{-N}$$

The second half of that has what we wanted, $T(X_0Y_1) + H(X_1Y_1).$ We still need to calculate the time-domain counterpart of $T(X1, Y0)$ by another length $2N$ FFT convolution so we are down to two FFT convolutions of length $2N,$ with only a total of $2N$ samples of the results wasted.

We can also confirm that user Fat32's suggestion of using a length $3N$ FFT convolution works, by wrapping the latter part of Eq. 2 around whereby $z^{-3N}$ becomes $z^{0} = 1$:

$$\begin{align}&amp;H(X_0Y_0) + T(X_1Y_1)\\
+ &amp;\big(T(X_0Y_0) + H(X_1Y_0) + H(X_0Y_1)\big)z^{-N}\\
+ &amp;\big(T(X_1Y_0) + T(X_0Y_1) + H(X_1Y_1)\big)z^{-2N}\end{align}$$

The middle line has what we wanted. So a single length $3N$ FFT convolution is also an option, with only a total of $2N$ result samples wasted.