The big trick in polynomial interpolation problems uses the function:
$$g(x)=(x-x_1)(x-x_2)\dots(x-x_n)$$
Note that this function is zero at each $x_n$ (and nowhere else). Also, $g'(x)=\sum_{i=1}^n \frac{g(x)}{x-x_i}$ (where the division is intended to cancel out the extra factor, so that it is defined even at each $x_i$).
If I divide $g(x)$ by $g(x^*)$ for some $x^*$ which is not equal to any $x_n$, I get a function that is $1$ at $x^*$ and $0$ at $x_i$. Adding functions of this form multiplied by scalars allows us to solve the problem of finding a polynomial which goes through $x_i,y_i$ (just take $\sum_{i=1}^ny_ig_i(x)$ where $g_i$ is zero for $x_j,j\ne i$ and one at $x_i$).
To get the derivative right at some specified $x^*$, just add to this expression $a\frac{g^*(x)}{{g^*}'(x^*)}$ where $g^*(x)$ is zero at each $x_i$ and at $x^*$; since $g^*(x)$ has a single root at $x^*$ its derivative there is nonzero, so the derivative of this part of the function is $a$, and we can set $a$ to be the desired derivative minus the derivative of the point-interpolation part of the function.
To give an example, let's find a polynomial such that $f(-2)=1$, $f(-1)=0$, $f'(0)=4$, $f(2)=2$. (I just made these numbers up.) First let's ignore the derivative part and construct functions that are zero at two of the points and one at the remaining point:
$$g_1(x)=\frac{(x+1)(x-2)}{(-2+1)(-2-2)}=\frac{x^2-x-2}{4}$$
$$g_2(x)=\frac{(x+2)(x-2)}{(-1+2)(-1-2)}=\frac{-x^2+4}{3}$$
$$g_3(x)=\frac{(x+2)(x+1)}{(2+2)(2+1)}=\frac{x^2+3x+1}{12}$$
Then $h(x)=g_1(x)+0g_2(x)+2g_3(x)=\frac{5x^2+3x-2}{12}$ satisfies $h(-2)=1$, $h(-1)=0$, $h(2)=2$. But it probably has the wrong derivative: $h'(x)=\frac{10x+3}{12}$, so $h'(0)=\frac14$ is not what we want. So let's construct a function which won't affect our other values, but has a nonzero derivative at $0$:
$$g^*(x)=(x+2)(x+1)(x-2)=x^3+x^2-4x-4\implies {g^*}'(x)=3x^2+2x-4$$
So ${g^*}'(0)=-4$. (It is possible for this to come out zero, in which case you can multiply the formula by $x-x^*$ to guarantee nonzero derivative at the expense of a higher degree.) Now the $h$ function is off from what we want by $4-\frac14=\frac{15}4$, while the adjustment function has derivative $-4$, so $-\frac{15}{16}{g^*}(x)$ has derivative $4-\frac14$ at $0$ and $f(x)=h(x)-\frac{15}{16}{g^*}(x)=-\frac{15}{16}x^3-\frac{25}{48}x^2+4x+\frac{43}{12}$ should do the trick.
Let
- $g(x) = x^5 f\left(\frac1x\right) = x^5+5x^4-2x^3+3x^2-4x+1$.
- $S = \{ a,b,c,d,e \}$ be the roots of $f(x)$.
- $T = \{ \frac1a, \frac1b, \frac1c, \frac1d, \frac1d \}$ be the roots of $g(x)$.
- For $I \subset S$ and $J \subset T$, let $\lambda_I = \prod_{\lambda \in I}\lambda$ and $\mu_J = \prod_{\mu \in J}\mu$.
The polynomial we seek equals to
$\quad\displaystyle\;F(x) \stackrel{def}{=} \prod_{I \subset S,|I| = 3}(x - \lambda_I)$.
Define a similar polynomial for $g$,
$\quad\displaystyle\;G(x) \stackrel{def}{=} \prod_{J \subset T,|J| = 2}(x - \mu_J)$.
By Vieta's formula, we have $abcde = -1$, this implies
$$F(x) = \prod_{I\subset S,|I|=3} \left(x + \frac{1}{\lambda_{S \setminus I}}\right) = \prod_{J\subset T,|J|=2}(x + \mu_J) = G(-x)$$
The problem comes down to given $g(x)$, how to compute $G(x)$ whose roots are
product of distinct pairs of roots of $g(x)$.
It will be hard to relate the coefficients of $g$ and $G$ directly. However, there is a simple relation between the power sums. More precisely, for any
$k \in \mathbb{Z}_{+}$, let
- $P_k(g) \stackrel{def}{=} \sum_{\mu \in T} \mu^k$ be the sum of roots of $f(x)$ raised to power $k$.
- $P_k(G) \stackrel{def}{=} \sum_{J \subset T,|J|=2} \mu_J^k$ be the sum of roots of $G(x)$ raised to power $k$.
We have
$$P_k(G) = \frac12( P_k(g)^2 - P_{2k}(g))\tag{*1}$$
To make following descriptions more generic, let $n = 5$ and $m = \frac{n(n-1)}{2}$.
Define coefficients $\alpha_k, \beta_k$ as follow:
$$g(x) = x^n - \sum\limits_{k=1}^n \alpha_k x^{n-k}
\quad\text{ and }\quad
G(x) = x^m - \sum\limits_{k=1}^m \beta_k x^{m-k}$$
Following are the steps to compute coefficients $\beta_k$ from coefficients $\alpha_k$ manually.
- Compute $P_k(g)$ using
Newton's identity for $1 \le k \le 2m$.
$$P_k(g) = \sum_{j=1}^{\min(n,k-1)} \alpha_j P_{k-j}(g) + \begin{cases}
k \alpha_k, & k \le n\\
0, & \text{otherwise}\end{cases}
$$
Compute $P_k(G)$ from $P_k(g)$ using $(*1)$.
Compute $\beta_k$ from $P_k(G)$ using Newton's identities again:
$$\beta_k = \frac1k\left( P_k(G) - \sum_{j=1}^{k-1} \beta_j P_{k-j}(G) \right)$$
I am lazy, I implement above logic in maxima (the CAS I use) and compute these numbers. The end result is
$$F(x) = x^{10}-2x^9+19x^8-112x^7+82x^6+97x^5-15x^4+58x^3+3x^2+3x+1$$
If one has access to a CAS, there is a quicker way to get the result.
For example, in maxima, one can compute the resultant between $g(t)$ and $g\left(-\frac{x}{t}\right)$ using the command resultant(g(t), g(-x/t), t)).
The resultant of two polynomials is essentially the GCD of them over the polynomial ring. It vanishes when and only when the two polynomials share a root. When the resultant between $g(t)$ and $g\left(-\frac{x}{t}\right)$ vanishes, $x$ either equals to $-\mu^2$ for a root $\mu \in T$ or $-\mu\nu$ for some $\mu, \nu$ in $T$.
If one ask maxima to factor output of above command, the result is
$$-(x^5+29x^4-34x^3+3x^2+10x+1)F(x)^2$$
The first factor is nothing but $\prod\limits_{\mu \in T}(x + \mu^2)$, this confirm the expression we get for $F(x)$ is the product $\prod\limits_{J \subset T,|J| = 2}(x + \mu_J)$ we desired.
Best Answer
You could factorise the polynomial as $p(x-a)^2(x-b)$ and use this to express coefficients in terms of the roots and solve. Or note that $f$ and $f'$ have a common root and use the division algorithm to locate it.