#C++ A fast and robust solution#
This is by no means the shortest answer, but so far no other answer has confirmed that it will work with curves with the full 3 stationary points (maxima/minima.) These correspond to curves with large values of b and c of opposite sign to a. None of the test data in the other solutions has local maxima/minima apart from the global one. The nearest is Eelvex's 1 2 0 0 1 which has a kink on one side where dy/dx falls to zero.
This program handles all the test data shown in other solutions, plus the following interesting examples, amongst others:
1 0 1 0 0 --> multiple root 0
1 0 -1 0 0 --> multiple root at 0, but the program prefers the simple root at -1
1 0 -1 0 0.25 --> multiple roots at +/- 0.7071 = sqrt(0.5)
It cannot handle the case a=0, but that is not a quartic so it is out of scope. The only problem I have found but not eliminated is b=c=d=0, but that is an easy problem to fix. I will do it when I golf the code (if I get round to that.) I expect it to come in around 700-800 chars.
How it works
The program divides the whole formula by a (which does not affect the roots, but ensures that we know that y= +infinity at x = +/- infinity )
It then calculates the formula for dy/dx (a cubic) and solves this analytically for dy/dx=0. The code to do this is quite short (the explanation in the comments is much longer than the code itself.)
It calculates the values of y at the up to 3 stationary points where dy/dx=0 and looks for the lowest (global minimum.) If this is positive, the whole curve is above the x axis and there are no roots.
If it is zero, then it is a root and this is reported directly. Such a root would be problematic for newton raphson.
If it is less than zero, the program makes a guess for the root on the opposite side of the global minimum from the other stationary points (to avoid getting trapped in the local maximum if there is one.)
Then it switches to newton-raphson (recursive, as you insisted.) This method is fast and simple, and once the right region of the curve has been identified, does not run the risk of dividing by zero.
I did run into other problems with it though. It quickly got very close to a solution (good enough for a human) but kept searching. With the rapid initial convergence, I assumed that the problem was that it couldn't get to f(x)==0 because f'(x) was large and f(x)/f'(x) was being rounded to zero in the calculation. So I put in an extra test condition for f(x)/f'(x) but it still did the same. Finally I limited it to just 10 iterations, which is more than enough for all the tests I have done. When I golf it I may put it up to 99 iterations.
There are certainly other methods which might work better, but I think that might go into another answer (which I don't think I have time for.)
#include "stdafx.h"
#include "math.h"
double a,b,c,d,e,j,k,l,p,q,t,z;
//dy/dx= i* x3 + j*x2 + k*x + l. No variable needed for i as it is normalised to 4.
//depressed form of dy/dx: t3+pt+q. A subsitution is used to eliminate the squared term, see below.
//z=q^2/4 + p^3/27. z is proportional to the "discriminant." The sign tells the number of roots of dy/dx
// +ve z means 1 real root, -ve z means 3 real roots, 0 means one simple root plus a double root.
// recursive newton-raphson to depth i
void s(double x, int i){
double m=(x + b)*x*x*x + c*x*x + d*x + e;
double n=m/(4*x*x*x + j*x*x + k*x + l);
if (m == 0 | n == 0 | i == 0) { printf("solution %f \n", x); }
else {
printf("searching %f \n", x );
s(x-n, i-1);
}
}
void _tmain(){
while (true){
double r[8] = { 0, 0, 0, 0, 1, 1, 1, 1 };
//r[1,2,3] and r[5,6,7] store x and y values of maxima and minima.
//r[0] and r[4] store dummy values for comparison.
scanf_s("%lf %lf %lf %lf %lf", &a, &b, &c, &d, &e);
printf("data for minima and maxima \n");
printf(" z x dy/x y hi/low \n");
//If there is only one root, the others display data for x=0.
//hi/low indicate which y values are highest and lowest: 1st, 2nd, 3rd or dummy.
//Divide by a to give a "monic polynomial" with a=1.
//This does not affect the roots but saves a lot of typing of powers of a.
// differentiate to i*x3 + j*x2 + k*x + l. As a=1, we know i=4.
b /= a; c /= a; d /= a; e /= a;
j = 3*b; k = 2*c; l = d;
// to prepare to solve this cubic, convert to depressed cubic t3+pt+q, eliminating x2. x=t-j/(3*i)
// p=(3*i*k - j2)/(3*i2) q=(2*j3 - 9*i*j*k + 27*i2*l)/27*i3. But we know i=4.
p = (12 * k - j*j) / 48;
q = (2 * j*j*j - 36 * j*k + 432 * l) / 1728;
z = q*q / 4 + p*p*p / 27;
//u and v indicate which of the stationary points are highest and lowest.
//solve dy/dx, store the roots in r[1,2,3] and the y values in r[5,6,7].
int u=0; int v=0;
for (int g = 1; g < 4; g++){
//if z>0 use cardano's method to store the one real root of dy/x in r[3] and dummy values in r[1,2]
//if z>=0 use viete's trigonometric method to find the multiple real roots.
//use of different methods for each case avoids the need to handle complex numbers.
if (z > 0){ r[1] = 0; r[2] = 0; r[3] = cbrt(-q / 2 + sqrt(z)) + cbrt(-q / 2 - sqrt(z)) - j / 12; }
else r[g] = sqrt(-p / .75)*cos(acos(-q / sqrt(-p * p * p * 4 / 27)) / 3 - g*acos(-.5)) - j / 12;
r[g+4] = (r[g]+b)*pow(r[g],3) + c*r[g]*r[g] + d*r[g] + e;
if (r[g + 4] > r[g+3])u=g;
if (r[g + 4] < r[g+3])v=g;
printf("%f %f %f %f %d %d \n \n", z, r[g], 4 * pow(r[g], 3) + j*pow(r[g], 2) + k*r[g] + l, r[g+4], u, v);
}
//because we divided by a, the new a=1 and y tends to infinity at large |x|.
//so if the lowest stationary point has y>0, the whole curve is above the x axis and there is no solution.
//if the lowest stationary point has y=0 it is tangent to the x axis and would cause problems for newton-raphson
//if the lowest stationary point has y<0, use newton raphson, 10 iterations. ensure 1st guess is on the side
//opposite the highest stationary point, just in case that point also has y<0, to avoid getting trapped.
if (r[v+4] >0) printf("n\n");
if (r[v+4]==0) printf("multiple root %f \n",r[v]);
if (r[v+4] <0) s(r[v]+(r[v]-r[u]),10);
}
}