czt of a short sinusoid
You are not analyzing a sine wave, you are analyzing a windowed sine wave.

$$
y_1(t) = \sin(2\pi f_0t /f_s) w_{\rm rect}(n)
$$
where $f_0 = 510,000$, $f_s = 4,000,000$, and
$$
w_{\rm rect}(n) = \left\{  \begin{array}{cl}
1 &amp; \mbox{where } 0 \le t \le (2^8 - 1)/fs\\
0 &amp; \mbox{otherwise}
 \end{array}\right.
$$

One way to get closer to what you expect is to choose the duration so that your sine wave is "continuous" at each end.  That would mean setting

points = 40000; //2^8


for example.

If I do that, then I get



See this question and its answers for addressing a similar problem.

The reason I chose 40,000 points is to ensure that there is no discontinuity at the ends of the data. This is what happens when you plot [y2 y2] for 40,000 points vs 256 points as in your example.

Notice the discontinuity in the 256 case (second plot).





I agree that choosing a larger number of points is "interpolating the problem away".  There is no real way around that, however.  There will only be 255 zero points in the FFT of a signal of length 256, so if you're expecting to see zeroes everywhere except $f_0$, you're going to be disappointed.



Let's look at the analytic expression for the DFT of a windowed sine wave:
$$
X(\omega) = \sum_{n=0}^{N-1}\sin(2\pi \hat{f}_0 n)e^{-j2\pi nk/N} \\
= \sum_{n=0}^{N-1} \frac{1}{2}\left[e^{+j2\pi \hat{f}_0 n} - e^{-j2\pi \hat{f}_0 n} \right]e^{-j2\pi n\omega/N}\\
= \sum_{n=0}^{N-1} \left [ e^{+j2\pi n ( \hat{f}_0 - \omega )}  - e^{-j2\pi n ( \hat{f}_0 + \omega )}  \right]\\
= \frac{1 - e^{+j2\pi ( \hat{f}_0 - \omega ) N} }{1 - e^{+j2\pi ( \hat{f}_0 - \omega ) } } + \frac{1 - e^{-j2\pi ( \hat{f}_0 + \omega ) N} }{1 - e^{-j2\pi ( \hat{f}_0 + \omega ) } }\\
= \frac{e^{+j\pi ( \hat{f}_0 - \omega ) N}}{e^{+j\pi ( \hat{f}_0 - \omega ) }}
\frac{(e^{-j\pi ( \hat{f}_0 - \omega ) N} - e^{+j\pi ( \hat{f}_0 - \omega ) N})}
{e^{-j\pi ( \hat{f}_0 - \omega )} - e^{+j\pi ( \hat{f}_0 - \omega )} }
+
\frac{e^{-j\pi ( \hat{f}_0 + \omega ) N}}{e^{-j\pi ( \hat{f}_0 + \omega ) }}
\frac{(e^{+j\pi ( \hat{f}_0 + \omega ) N} - e^{-j\pi ( \hat{f}_0 + \omega ) N})}
{e^{+j\pi ( \hat{f}_0 + \omega )} - e^{-j\pi ( \hat{f}_0 + \omega )} }\\
= e^{+j\pi ( \hat{f}_0 - \omega ) (N-1)} \frac{\sin(\pi ( \hat{f}_0 - \omega ) N/2)}{\sin(\pi ( \hat{f}_0 - \omega )/2)}
+
e^{-j\pi ( \hat{f}_0 + \omega ) (N-1)} \frac{\sin(\pi ( \hat{f}_0 + \omega ) N/2)}{\sin(\pi ( \hat{f}_0 + \omega )/2)}
$$
which is where the sinc-like shape comes from.