One can use Green's Theorem if the offending point, say $\vec r_0$, is excluded from the region of integration.
To do that, we deform the boundary contour with a "keyhole" contour that encircles the excluded point $\vec r_0$.
This reduced the problem to evaluating the line integral
$$\lim_{\epsilon \to 0}\int_0^{2\pi}\vec F(\vec r_0+\epsilon \hat r)\cdot \hat \phi \epsilon \, d\phi \tag 1$$
For example, suppose $\vec F(\vec r)=\frac{-\hat xy+\hat yx}{x^2+y^2}$.
NOTE:
The vector $\vec F$ here represents the static magnetic field $\vec H$ (modulo the constant $I/2\pi$) from a line current $I$ aligned along the $z$-axis.
Clearly, for all $(x,y)\ne (0,0)$, the first partial derivatives of $\vec F$ satisfy the relationship
$$\frac{\partial F_y}{\partial x}=\frac{\partial F_x}{\partial y} \tag 2$$
NOTE:
We remark that $(2)$ is the $z$ component of the curl of $\vec F$, which by Maxwell's Equation for static magnetic fields is $\hat z\cdot \nabla \times \vec H=0$ for $\vec r \ne 0$.
However, $\vec F$ is singular at the origin. We can still use $(2)$ in Green's Theorem to evaluate the line integral of $\vec F$ on any contour $C$ that encloses the origin by using the contour deformation given in $(1)$. Proceeding we find with $\vec r_0=0$
$$\begin{align}
\oint_C \vec F(\vec r)\cdot \,d\vec \ell&=\lim_{\epsilon \to 0}\int_0^{2\pi}\vec F(\vec r_0+\epsilon \hat r)\cdot \hat \phi \epsilon \, d\phi\\\\
&=\lim_{\epsilon \to 0}\int_0^{2\pi} \left(\frac{\hat \phi}{\epsilon}\right)\cdot \hat \phi \epsilon \,d\phi\\\\
&=2\pi
\end{align}$$
as expected from Ampere's Law (after multiplying by $I/2\pi$).
You are not suppose to have that negative in the second component. You are computing a work integral over the ellipse. Unless stated otherwise, it is always assumed that curves are given the counter clock-wise orientation. Note that orientation really does matter since the force field may be stronger (or weaker) in the opposing direction.
\begin{align*} &\int_{0}^\pi (3 \sin^2t + 2,16 \cos t) \cdot (- \sin t, b\cos t) \ dt \\ & = \int_{0}^\pi -3b^2 \sin^3t- 2 \sin t + 16 b\cos^2t \ dt \\ & = -4b^2+8 \pi b - 4 \end{align*}
Best Answer
Let $E$ denote the interior of the ellipse and $\partial E$ the boundary (i.e. the ellipse itself). You're trying to calculate
$$\color{blue}{\oint_{\partial E}\mathbf{F}\cdot d\mathbf{r}}$$
Your first hope is that Green's theorem applies, so you can instead calculate
$$\int_E(Q_x-P_y)\,dA$$
The problem is that the hypotheses of Green's theorem are not satisfied, for the functions $P$ and $Q$ are not continuously differentiable throughout the region $E$. As you observed, there is no way we may continuously extend the field to be defined at the point $(1,0)$, so that point is a singularity.
The way to proceed in such cases is to put a little closed curve around the singularity but still lying entirely inside the region $E$ of interest; in this case a unit circle will be convenient (the calculation will show you why). So let $D$ be the unit disk centered at $(1,0)$, and let $\partial D$ be the unit circle centered there, which has parametrization $C(t)=(\cos t +1, \sin t)$ for $t\in[0,2\pi)$.
You can easily check by direct calculation that
$$\color{red}{\oint_{\partial D}\mathbf{F}\cdot d\mathbf{r}}=\int_0^{2\pi}(\cos^2t+\sin^2t)\,dt=2\pi$$
Now, consider the region $E-D$, the complement of the disk in the elliptical region. In this region, $P$ and $Q$ are continuously differentiable everywhere, so we may use Green's theorem, which says
$$\int_{E-D}(Q_x-P_y)\,dA=\oint_{\partial(E-D)}\mathbf{F}\cdot d\mathbf{r}\tag{$\star$}$$
The left side vanishes, because $Q_x-P_y=0$ in the region $E-D$. What about the right side? What is the "boundary" $\partial(E-D)$ of the elliptical annular region $E-D$? It is $\partial E - \partial D$, where we interpret the negative sign on $\partial D$ to mean reversing the standard counterclockwise orientation.
So $(\star)$ says
$$0=\color{blue}{\oint_{\partial E}\mathbf{F}\cdot d\mathbf{r}}-\color{red}{\oint_{\partial D}\mathbf{F}\cdot d\mathbf{R}}$$
The blue integral is the one the problem asks for, and the red integral we calculated above.
To summarize, the strategy for calculating line integrals of curl-free fields around closed curves that encircle a single singularity is to replace the original closed curve with an easier one (but still contained within the original one). Green's theorem in the region between the curves guarantees the line integral around the easier curve equals the line integral around the original curve.
What makes a good candidate for an "easier curve"? Well, you're trying to find a curve whose parametrization will make the integrand simple enough to calculate. In your problem, integrating around the ellipse gives a "hard" integrand. But integrating around a circle centered at the singularity gives an "easy" integrand.
Why a circle in this case? Because your field is closely related to the pullback of the length element of a circle of radius $r$ (the length element is $ds=r d\theta$); thus it was, in some sense, tailor-made to be integrated around a circle. It is not always so easy to spot what sort of curve will make the line integral easier, but you should remember the general form $(-\frac{y}{r^2}, \frac{x}{r^2})$, which is $\frac{1}{r}$ times the circle's unit tangent vector. The line integral of the unit tangent vector around the circle gives $2\pi r$, and the factor $\frac{1}{r}$ scales the result to $2\pi$. (In other words, you are just integrating $d\theta=\frac{1}{r}ds$ around the circle.) So the radius of the circle you use around the singularity doesn't matter. (This make sense, because the region between two concentric circles centered at the singularity doesn't contain a singularity and the field is conservative there, so by the above argument, the line integral of the field around them should be the same.)