7
\$\begingroup\$

Write the shortest program to solve a Quartic equation.

A quartic equation is a polynomial equation of the form:

ax⁴ +bx³+cx²+dx+e=0

A solution for x is a number such that the above evaluates to 0.


Rules

  • The program must take in the 5 floating point coefficients of the equation ax^4+bx^3+cx^2+dx+e as "a b c d e" separated by spaces.
  • It must then output a single real solution. If there exist no real solutions then it outputs an "n".
  • No bonus for extra solutions.
  • No use of CAS (Computer algebra systems) This includes symbolic manipulation properties of Mathematica and related software.
  • No connecting to external libraries, web calculators etc...
  • You may use formulas or your own manipulation
  • Native math libraries that do not solve equations may be used, eg. .NET: System.Math
  • You may deduct 50 from your char count if recursion is a main technique used.
  • As long as there is no difference with the answer obtained algebraically, numerical approximation is welcome however adds 25 chars to your solution. This really relies on the accuracy of the floating-point numbers. (Which should be 64 bit).

The winner will be picked in 1 week. (26/02/2014)

\$\endgroup\$
  • 2
    \$\begingroup\$ Hints at math.stackexchange.com/questions/785/… \$\endgroup\$ – ComFreek Feb 19 '14 at 18:45
  • 3
    \$\begingroup\$ I would remove "Program has to return answer in under a second". It may be under a second on my machine but 30s on yours. \$\endgroup\$ – Danny Feb 19 '14 at 18:52
  • 1
    \$\begingroup\$ You may deduct -50 from your char count - isn't this the same as adding +50 to the char count? Who would choose to do that? \$\endgroup\$ – Digital Trauma Feb 19 '14 at 18:59
  • 1
    \$\begingroup\$ As written, the question requires solutions for the general case. \$\endgroup\$ – Peter Taylor Feb 19 '14 at 19:08
  • 1
    \$\begingroup\$ Use Newton Raphson with care, it can get stuck in the wrong place and divide by zero. See my comments on existing answers (unfortunately it is only possible to tag one person per comment, otherwise I would only have made a general comment here.) \$\endgroup\$ – Level River St Feb 19 '14 at 22:19
5
\$\begingroup\$

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.

enter image description here

\$\endgroup\$
  • \$\begingroup\$ Nice. BTW, theoretically my answer always finds a root if it exists since I use the bisection method. Practically it might need too much time if the smallest root is a very-very big number. \$\endgroup\$ – Eelvex Mar 11 '14 at 0:40
3
\$\begingroup\$

J, 222

247 - recursion + numerical

NB. multiple lines for visibility
f=:[:-(]$:~[:}.[-%&{.*],#&0@l)`[@.(0>l=:-&#)
r=:1={.
s=:[:=/+/@(2=&*/\])@(,.]*_1^i.@#)@({.@>@(([;];f;j;f f j=.]f f)([:}:(i._5)*])))@|.
u=:0,0([(]*_1^r@k)>:@]^:(1 1-:k=:[:*(0 p.~[)*(p.(,-)))^:_)~]
h=:{.@((n@],]{~[:r[:*p.*(p.n=:{.+-:@-~/))^:_ u)`('n'"_)@.s@|.

Methods

  • Using Sturm's theorem to find number of roots.

    f is polynomial division (recursive)

    s is Sturm's theorem output: 0 for no real roots, 1 otherwise. It evaluates Sturm's polynomials ([;];f;j;f f j=.]f f) and checks their behaviour at +/- infinity.

  • Using bisection method to avoid divisions by zero, local maxima etc.

    u scans from 0 to +/- infinity to find an appropriate interval

    ((n@],]{~[:r[:*p.*(p.n=:{.+-:@-~/))^:_ u) evaluates the polynomial at the edges and midpoint then halfs the interval accordingly

Verification of results

Using J's buildin polynomial root finder as realRoots:

ti,>([;h;realRoots) each samples
┌──────────────┬────────┬─────────────────┐
│a b c d e     │output  │all real roots   │
├──────────────┼────────┼─────────────────┤
│_1 2 3 4 5    │_1.11029│3.37192 _1.11029 │
├──────────────┼────────┼─────────────────┤
│_1 2 3 4 _5   │0.728727│3.18248 0.728727 │
├──────────────┼────────┼─────────────────┤
│1 2 3 4 5     │n       │                 │
├──────────────┼────────┼─────────────────┤
│1 1 2 _1 _1   │0.728262│0.728262 _0.52236│
├──────────────┼────────┼─────────────────┤
│1 2 3 4 0     │0       │_1.65063 0       │
├──────────────┼────────┼─────────────────┤
│1 2 0 0 1     │_1      │_1.83929 _1      │
├──────────────┼────────┼─────────────────┤
│_1 _1 _1 _1 _1│n       │                 │
└──────────────┴────────┴─────────────────┘

BTW, a simple implementation of Newton's method would be:

f=:3 :'((-(y&p.%y&p.D.1))^:_)0'@|.
\$\endgroup\$
  • \$\begingroup\$ Thanks for the comment on my solution. Yours is the only other practical solution here, though I do not understand your code at all. I'd love to see a clearer explanation of what's going on, especially the "scan for interval" (I assume you are looking for a change of sign) and bisection (I wonder how it handles multiple sign changes within the selected interval.) \$\endgroup\$ – Level River St Mar 11 '14 at 1:34
  • \$\begingroup\$ @steveverrill, the bisection will never fail no matter what, because from an interval (+,-) you choose the midpoint that it's either +,-,or 0. In any case you can either form a new interval (+,-) or you have a solution. The weak point of my algorithm is the scanning which might go on and on. My scanning makes dx steps left and right (x-dx, x+dx) and compares them to x=0. If (x,x+-dx) <= 0, we have an interval to bisect. It would be tremendounsly faster and robust if I used Sturm's polynomial to make an initial guess but I was too tired :) \$\endgroup\$ – Eelvex Mar 11 '14 at 12:44
2
\$\begingroup\$

C — 340 bytes + 25 points for numerical approximation = 365

To compile gcc -o quartic quartic.c -lm, to run ./quartic a b c d e. Uses Newton–Raphson method, no big whoop. Using deterministic methods, I can't see getting anywhere near this compact.

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#define C(x) atof(argv[x])
int main(int argc,char *argv[]){int i;double a,b,c,d,e,x,y;a=C(1);b=C(2);c=C(3);d=C(4);e=C(5);for(x=0,y=1,i=0;i<100;x=y,i++){y=x-((((a*x+b)*x+c)*x+d)*x+e)/(((4*a*x+3*b)*x+2*c)*x+d);if(fabs(x-y)<1e-09){printf("%.9f\n",y);return 0;}}printf("n\n");return 1;}

e.g.

$ ./quartic 1 2 3 4 5
n

$ ./quartic -1 2 3 4 5
-1.110290760

$ ./quartic -1 2 3 4 -5
0.728726879

Verified via Wolfram|Alpha, results check out.

\$\endgroup\$
  • 1
    \$\begingroup\$ Troll FTW - naming it quintic to solve a "quartic" equation ;-) \$\endgroup\$ – Digital Trauma Feb 19 '14 at 20:57
  • \$\begingroup\$ Pretty sure you could work recursion into this fairly easily to grab the -50 bonus \$\endgroup\$ – Digital Trauma Feb 19 '14 at 20:58
  • \$\begingroup\$ @DigitalTrauma - fixed it... \$\endgroup\$ – user15259 Feb 19 '14 at 21:09
  • 2
    \$\begingroup\$ This solution is incomplete as it runs the risk of getting caught in a local minimum/maximum and dividing by zero without showing conclusively that there are no real roots. See my comment on DigitalTrauma's answer. \$\endgroup\$ – Level River St Feb 19 '14 at 22:13
  • \$\begingroup\$ @steveverrill - You're welcome to provide your own solution. \$\endgroup\$ – user15259 Feb 20 '14 at 1:28
1
\$\begingroup\$

bc (driven by bash), score:137 (162 bytes -50 recursion bonus +25 numerical penalty)

bc -ql<<Q
n=0
define r(x){if(n++>99){print "n
";halt}
y=x-($1*x^4+$2*x^3+$3*x*x+$4*x+$5)/(4*$1*x^3+3*$2*x*x+2*$3*x+$4)
if(y==x){print x,"
";halt}else r(y)}r(0)
Q

This is pretty much the same algorithm as the answer, though it uses recursion to spice it up.

Run as follows:

$ ./quartic.sh 1 2 3 4 5
n
$ ./quartic.sh -1 2 3 4 5
-1.11029075952058816917
$ ./quartic.sh 1 2 3 4 -5
.68412431945306971026
$ 

Of course Newton-Raphson has its limitations, and won't necessarily give all real roots in all cases where the exist.

According to WolframAlpha, this one has real roots at -1.9144 and 1.3731:

$ ./quartic.sh 1 1 -2 -1 -1
n
$ 

Also there is a very real risk of divide-by-zero in some cases. But this is code-golf, so who cares ;-)

\$\endgroup\$
  • 1
    \$\begingroup\$ Divide by zero is a major issue. If your code fails to find a solution, it won't be for lack of iterations. It will be because it gets caught in one of the local maxima/minima, where dy/dx falls to zero. Newton-Raphson will try to divide by zero and the program will crash. With a quadratic such as x^2+1 (only one maximum/minimum) this error could be trapped and used as conclusive proof of no real roots. But a quartic has up to 3 local maxima/minima so you don't know if there's really no roots or if you made a bad first guess for x. \$\endgroup\$ – Level River St Feb 19 '14 at 22:02
1
\$\begingroup\$

C, 182 + 25 = 207

This is not the most eleant solution but it is certainly creative. Using iterger pointer conversions, it loops through every possibile value that a double can hold, tries it out, and checks if it equals 0.

int main(int z,char **v){double a[5],f=0,r=1;int i=5;while(i-->0)scanf("%lf",a+i);while(++(*(int*)&f)&&r){i=5,r=0;while(i-->0)r+=a[i]*pow(f,i);}if(r)printf(“%f",r);else printf("n");}
\$\endgroup\$
1
\$\begingroup\$

Ruby — 132 characters + 25 numerical = 157

require'matrix';f=->a,b,c,d,e{a=-a.to_f;Matrix[[0,0,0,e/a],[1,0,0,d/a],[0,1,0,c/a],[0,0,1,b/a]].eigen.eigenvalues.find(&:real?)||?n}

I think this toes the line of using native math libraries, but it's a pretty cool demonstration of what can be done with Ruby, as well as an intriguing way to look at solving univariate polynomial equations through the lens of linear algebra.

Basically, f just constructs the companion matrix of the equivalent monic quartic and computes its eigenvalues, which is the same thing as finding the roots of the original quartic.

\$\endgroup\$
0
\$\begingroup\$

Javascript/ES6 148

Since only a solution is required, I used a Newton solver.

f=s=>{a=s.split(' ');for(x=i=0;i++<1e5;x-=y)z=x*x,y=(a[0]*z*z+a[1]*x*z+a[2]*z+a[3]*x+ +a[4])/(4*a[0]*z*x+3*a[1]*z+2*a[2]*x+ +a[3]);return y|0?'n':x}

If you don't have Firefox, you can test it without fat arrow :

function f(s){a=s.split(' ');for(x=i=0;i++<1e5;x-=y)z=x*x,y=(a[0]*z*z+a[1]*x*z+a[2]*z+a[3]*x+ +a[4])/(4*a[0]*z*x+3*a[1]*z+2*a[2]*x+ +a[3]);return y|0?'n':x}

Test : f('1 -2 5 -12 4') returns 0.3882914410047445

Warning : This can fail for some equations (if d=0 for example).

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.