C++ A fast and robust solution, 704 chars + numeric - recursion = 679
#include "stdafx.h"
#include "math.h"
double a,b,c,e,j,k,l,p,q,t,z;void s(double x, int i){if(i==99)printf("%f\n",x);else s(x-((x+b)*x*x*x+c*x*x+l*x+e)/(4*x*x*x+j*x*x + k*x + l),i+1);}
void _tmain(){double r[8];scanf_s("%lf%lf%lf%lf%lf",&a,&b,&c,&l,&e);b/=a;c/=a;l/=a;e/= a;j=3*b;k=2*c;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;int u=0,v=0,g;
for(g=1;g<4;g++){r[g]=z>0|p==0?cbrt(-q/2+sqrt(z))+cbrt(-q/2-sqrt(z))-j/12: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]+l*r[g]+e;if(r[g+4]>r[g+3]|g==1)u=g;if(r[g + 4]<r[g+3]|g==1)v=g;}
if(r[v+4]>0)printf("n\n");if(r[v+4]==0)printf("%f\n",r[v]);if(r[v+4]<0)s(2*r[v]-r[u]+(r[v]-r[u]==0),0);}
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 y is positive here , the whole curve is above the x axis and there are no roots.
If y is zero, then it is a root and this is reported directly. Such a root would be problematic for newton raphson.
If y 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 searches using 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. Then I limited it to just 10 iterations, which is more than enough for all the tests I have done.
The current version goes up to 99 iterations. The ungolfed version has medium verbosity output (with uncommentable code for high verbosity) as seen in the screenshot and has its iterations terminated when f(x)==0. The golfed version has low verbosity (in accordance with specified requirements) and always runs to 99 iterations. The golfed version has other minor differences, such as no while loop and simplified (but less clear) expressions.
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.)
Ungolfed version
#include "stdafx.h"
#include "math.h"
double a,b,c,d,e,j,k,l,p,q,t,z;
//explanation of variables
//dy/dx= i* x3 + j*x2 + k*x + l. No i variable needed 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=q2/4 + p3/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;
if (m == 0 | i == 99) printf("z= %f x= %f y= %f iteration %d\n", z,x,m,i);
else s(x-m/(4*x*x*x + j*x*x + k*x + l), i+1);
}
void _tmain(){
while (true){
double r[8];
//r[1,2,3] and r[5,6,7] store x and y values of the maxima and minima.
//r[0] and r[4] are dummies to handle [subscript-1] references.
scanf_s("%lf %lf %lf %lf %lf", &a, &b, &c, &d, &e);
//uncomment for verbose output: printf("data for minima and maxima \n");
//uncomment for verbose output: printf(" z x dy/x y hi/low \n");
//hi/low indicate which y values are highest and lowest: 1st, 2nd, 3rd.
//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,v=0,g;
for (g = 1; g < 4; g++){
//if z>0 or p==0 use cardano's method to store the one real root of dy/x three times in r[1,2,3]
//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.
r[g] =z > 0 | p==0? cbrt(-q / 2 + sqrt(z)) + cbrt(-q / 2 - sqrt(z)) - j / 12 :
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]|g==1)u=g;
if (r[g + 4]<r[g+3]|g==1)v=g;
// uncomment for verbose output
//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, 99 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.
//special case: if r[v]-r[u] == 0 (implies only one stationary point) then add 1!
if (r[v+4] >0) printf("n\n");
if (r[v+4]==0) printf("z= %f x= %f multiple root", z,r[v]);
if (r[v+4] <0) s(r[v]+(r[v]-r[u])+(r[v]-r[u]==0),0);
}
}
Output
All data tested by other answers is here, as well as many additional examples. The sign of the value z is positive when there is only one stationary point where dy/dx=0, negative when there are three (2 minima and a maximum if a is positive), and zero when there are two (a minimum if a is positive, plus an inflexion.) This is the only answer so far that presents test cases for three stationary points.
The y values are shown for easy checking of the program. y is zero at the root, and as can be seen, an acceptable zero is always found.
The golfed version is low verbosity, displaying only the root per question spec.
