Roots of a quadratic function: Difference between revisions

(For double-precision floats, set <math>b = -10^9</math>.)

Consider the following implementation in [[Ada]]:

with Ada.Numerics.Elementary_Functions; use Ada.Numerics.Elementary_Functions;

 

begin

Put_Line ("X1 =" & Float'Image (R (1)) & " X2 =" & Float'Image (R (2)));

{{out}}

<pre>X1 = 1.00000E+06 X2 = 0.00000E+00</pre>

 

'''Task''': do it better. This means that given <math>a = 1</math>, <math>b = -10^9</math>, and <math>c = 1</math>, both of the roots your program returns should be greater than <math>10^{-11}</math>. Or, if your language can't do floating-point arithmetic any more precisely than single precision, your program should be able to handle <math>b = -10^6</math>. Either way, show what your program gives as the roots of the quadratic in question. See page 9 of

 

=={{header|Ada}}==

with Ada.Numerics.Elementary_Functions; use Ada.Numerics.Elementary_Functions;

 

begin

if B < 0.0 then

return (X, C / (A * X));

else

return (C / (A * X), X);

end if;

begin

Put_Line ("X1 =" & Float'Image (R (1)) & " X2 =" & Float'Image (R (2)));

Here precision loss is prevented by checking signs of operands. On errors, Constraint_Error is propagated on numeric errors and when roots are complex.

{{out}}

{{works with|ALGOL 68G|Any - tested with release [http://sourceforge.net/projects/algol68/files/algol68g/algol68g-1.18.0/algol68g-1.18.0-9h.tiny.el5.centos.fc11.i386.rpm/download 1.18.0-9h.tiny]}}

{{wont work with|ELLA ALGOL 68|Any (with appropriate job cards) - tested with release [http://sourceforge.net/projects/algol68/files/algol68toc/algol68toc-1.8.8d/algol68toc-1.8-8d.fc9.i386.rpm/download 1.8-8d] - probably need to "USE" the compl ENVIRON}}

BEGIN

 

ESAC

OD

{{out}}

<pre>

=={{header|AutoHotkey}}==

ahk forum: [http://www.autohotkey.com/forum/viewtopic.php?p=276617#276617 discussion]

MsgBox % quadratic(u,v, 1,3,2) ", " u ", " v

MsgBox % quadratic(u,v, -2,4,-2) ", " u ", " v

x2 := c/a/x1

Return 2

 

=={{header|BBC BASIC}}==

READ a$, b$, c$

PRINT "For a = " ; a$ ", b = " ; b$ ", c = " ; c$ TAB(32) ;

PRINT "the complex roots are " ; -b/2/a " +/- " ; SQR(-d)/2/a "*i"

ENDCASE

{{out}}

<pre>For a = 1, b = -1E9, c = 1 the real roots are 1E9 and 1E-9

=={{header|C}}==

Code that tries to avoid floating point overflow and other unfortunate loss of precissions: (compiled with <code>gcc -std=c99</code> for <code>complex</code>, though easily adapted to just real numbers)

#include <stdlib.h>

#include <complex.h>

 

return 0;

{{out}}<pre>(-1e+12 + 0 i), (-1 + 0 i)

(1.00208e+07 + 0 i), (9.9792e-08 + 0 i)</pre>

 

#include <math.h>

#include <complex.h>

x[0] = (-b + csqrt(delta)) / (2.0*a);

x[1] = (-b - csqrt(delta)) / (2.0*a);

{{trans|C++}}

{

b /= a;

x[1] = c/sol;

}

 

{

complex double x[2];

 

return 0;

 

<pre>x1 = (1.00000000000000000000e+20, 0.00000000000000000000e+00)

 

=={{header|C sharp|C#}}==

using System.Numerics;

 

Console.WriteLine(Solve(1, -1E20, 1));

}

{{out}}

<pre>((1E+20, 0), (1E-20, 0))</pre>

 

=={{header|C++}}==

#include <utility>

#include <complex>

std::pair<complex, complex> result = solve_quadratic_equation(1, -1e20, 1);

std::cout << result.first << ", " << result.second << std::endl;

{{out}}

(1e+20,0), (1e-20,0)

=={{header|Clojure}}==

 

"Compute the roots of a quadratic in the form ax^2 + bx + c = 1.

Returns any of nil, a float, or a vector."

(zero? sq-d) (f +)

(pos? sq-d) [(f +) (f -)]

 

{{out}}

nil

user=> (quadratic 1.0 2.0 1.0)

user=> (quadratic 1.0 3.0 1.0)

[2.618033988749895 0.3819660112501051]

 

=={{header|Common Lisp}}==

(list

(/ (+ (- b) (sqrt (- (expt b 2) (* 4 a c)))) (* 2 a))

 

=={{header|D}}==

 

CommonType!(T1, T2, T3)[] naiveQR(T1, T2, T3)

writefln(" Naive: [%(%g, %)]", naiveQR(1.0L, -10e5L, 1.0L));

writefln("Cautious: [%(%g, %)]", cautiQR(1.0L, -10e5L, 1.0L));

{{out}}

<pre>With 32 bit float type:

Naive: [1e+06, 1e-06]

Cautious: [1e+06, 1e-06]</pre>

 

=={{header|Elixir}}==

def roots(a, b, c) do

IO.puts "Roots of a quadratic function (#{a}, #{b}, #{c})"

Quadratic.roots(1, -1.0e10, 1)

Quadratic.roots(1, 2, 3)

 

{{out}}

 

=={{header|ERRE}}==

 

PROCEDURE SOLVE_QUADRATIC

DATA(1,3,2)

DATA(3,4,5)

{{out}}

<pre>For a= 1 ,b=-1E+09 ,c= 1 the real roots are 1E+09 and 1E-09

=={{header|Factor}}==

{{trans|Ada}}

b sq a c * 4 * - sqrt :> sd

b 0 <

[ b neg sd + a 2 * / ]

[ b neg sd - a 2 * / ] if :> x

 

--- Data stack:

1.0e+20

 

Middlebrook method

a c * sqrt b / :> q

1 4 q sq * - sqrt 0.5 * 0.5 + :> f

b neg a / f * c neg b / f / ;

 

 

--- Data stack:

1.0e+20

 

=={{header|Forth}}==

Without locals:

frot frot

( c a b )

2e f/ fover f/

( c a r1 )

With locals:

b b f* 4e a f* c f* f-

fdup f0< if abort" imaginary roots" then

\ test

1 set-precision

 

=={{header|Fortran}}==

===Fortran 90===

{{works with|Fortran|90 and later}}

 

IMPLICIT NONE

WRITE(*,*) "The roots are complex:"

WRITE(*,"(2(A,2E23.15,A))") "Root1 = ", croot1, "j ", " Root2 = ", croot2, "j"

{{out}}

Coefficients are: a = 0.300000000000000E+01 b = 0.400000000000000E+01 c = 0.133333333330000E+01

===Fortran I===

Source code written in FORTRAN I (october 1956) for the IBM 704.

COMPUTE ROOTS OF A QUADRATIC FUNCTION - 1956

READ 100,A,B,C

332 FORMAT(3HX1=,E12.5,4H,X2=,E12.5)

999 STOP

 

=={{header|GAP}}==

local d;

d := Sqrt(b*b - 4*a*c);

# This works also with floating-point numbers

QuadraticRoots(2.0, 2.0, -1.0);

 

=={{header|Go}}==

 

import (

test(c[0], c[1], c[2])

}

{{out}}

<pre>

 

=={{header|Haskell}}==

 

type CD = Complex Double

 

quadraticRoots :: (CD, CD, CD) -> (CD, CD)

where

d = sqrt $ b ^ 2 - 4 * a * c

mapM_

(print . quadraticRoots)

{{Out}}

<pre>((-0.6666666666666666) :+ 0.0,(-0.6666666666666666) :+ 0.0)

(999999.999999 :+ 0.0,1.000000000001e-6 :+ 0.0)

(1.0e10 :+ 0.0,1.0e-10 :+ 0.0)</pre>

 

=={{header|Icon}} and {{header|Unicon}}==

 

Works in both languages.

solve(1.0, -10.0e5, 1.0)

end

procedure solve(a,b,c)

d := sqrt(b*b - 4.0*a*c)

write(a,"*x^2 + ",b,"*x + ",c," has roots ",roots[1]," and ",roots[2])

 

{{out}}

->

</pre>

 

=={{header|IS-BASIC}}==

100 PROGRAM "Quadratic.bas"

110 PRINT "Enter coefficients a, b and c:":INPUT PROMPT "a= ,b= ,c= ":A,B,C

220 PRINT "The complex roots are ";-B/2/A;"+/- ";STR$(SQR(-D)/2/A);"*i"

230 END SELECT

 

=={{header|J}}==

'''Solution''' use J's built-in polynomial solver:

p.

'''Example''' using inputs from other solutions and the unstable example from the task description:

 

> {:"1 p. coeff

_0.666667 _0.666667

_0.333333j0.471405 _0.333333j_0.471405

1e6 1e_6

 

Of course <code>p.</code> generalizes to polynomials of arbitrary order (which isn't as great as that might sound, because of practical limitations). Given the coefficients <code>p.</code> returns the multiplier and roots of the polynomial. Given the multiplier and roots it returns the coefficients. For example using the cubic <math>0 + 16x - 12x^2 + 2x^3</math>:

+-+-----+

|2|4 2 0|

+-+-----+

p. 2 ; 4 2 0 NB. return coefficients

 

Exploring the limits of precision:

 

100000 1e_5

1 _1e5 1 p. 1{::p. 1 _1e5 1 NB. test roots

1e_5 _1e_15

1 _1e5 1 p. 1e5 1e_5-1e_5 _1e_15 NB. test displayed roots with adjustment

 

When these "roots" are plugged back into the original polynomial, the results are nowhere near zero. However, double precision floating point does not have enough bits to represent the (extremely close) answers that would give a zero result.

Middlebrook formula implemented in J

 

'a b c' =. y

q=. b %~ %: a * c

 

q_r 1 _1e6 1

 

=={{header|Java}}==

private static class Complex {

double re, im;

}

}

{{out}}

<pre>

'''Section 1''': Complex numbers (scrolling window)

<div style="overflow:scroll; height:200px;">

def real(z): if (z|type) == "number" then z else z[0] end;

 

| [ ($r * cos), ($r * sin)]

end

'''Section 2:''' quadratic_roots(a;b;c)

# if none, emit empty;

# otherwise always emit two.

negate(divide(c; multiply(b; $f)))

end

'''Section 3''':

Produce a table showing [i, error, solution] for solutions to x^2 - 10^i + 1 = 0

def pow(i): . as $in | reduce range(0;i) as $i (1; . * $in);

def poly(a;b;c): plus( plus( multiply(a; multiply(.;.)); multiply(b;.)); c);

;

 

{{Out}} (scrolling window)

<div style="overflow:scroll; height:200px;">

$ jq -M -r -n -f Roots_of_a_quadratic_function.jq

0: error = +0 x=[0.5,0.8660254037844386]

11: error = +0 x=[1e-11,-0]

12: error = 1 x=[1000000000000,0]

 

=={{header|Julia}}==

Alternative solutions might make use of Julia's Polynomials or Roots packages.

 

 

function quadroots(x::Real, y::Real, z::Real)

for (x, y, z) in zip(a, b, c)

@printf "The roots of %.2fx² + %.2fx + %.2f\n\tx₀ = (%s)\n" x y z join(round.(quadroots(x, y, z), 2), ", ")

 

{{out}}

The roots of 10.00x² + 1.00x + 1.00

x₀ = (-0.05 + 0.31im, -0.05 - 0.31im)</pre>

 

=={{header|Kotlin}}==

{{trans|Java}}

 

data class Equation(val a: Double, val b: Double, val c: Double) {

 

equations.forEach { println("$it\n" + it.quadraticRoots) }

{{out}}

<pre>Equation(a=1.0, b=22.0, c=-1323.0)

 

=={{header|lambdatalk}}==

1) using lambdas:

 

x1 = 1

x2 = 1

 

=={{header|Liberty BASIC}}==

'assume a<>0

print quad$(a,b,c)

quad$=str$(x/(2*a)+sqr(D)/(2*a));" , ";str$(x/(2*a)-sqr(D)/(2*a))

end if

 

=={{header|Logo}}==

localmake "d sqrt (:b*:b - 4*:a*:c)

if :b < 0 [make "d minus :d]

output list (:d-:b)/(2*:a) (2*:c)/(:d-:b)

 

=={{header|Lua}}==

like that from the Complex Numbers article. However, this should be enough to demonstrate its accuracy:

 

if b < 0 then return qsolve(-a, -b, -c) end

val = b + (b^2 - 4*a*c)^(1/2) --this never exhibits instability if b > 0

for i = 1, 12 do

print(qsolve(1, 0-10^i, 1))

The "trick" lies in avoiding subtracting large values that differ by a small amount, which is the source of instability in the "normal" formula. It is trivial to prove that 2c/(b + sqrt(b^2-4ac)) = (b - sqrt(b^2-4ac))/2a.

 

=={{header|Maple}}==

 

solve(1.0*x^2-10.0^9*x+1.0,x,explicit,allsolutions);

 

{{out}}

<pre> (1/2) (1/2)

1.000000000 10 , 1.000000000 10 </pre>

 

Possible ways to do this are (symbolic and numeric examples):

Solve[x^2 - 10^5 x + 1 == 0, x]

Root[#1^2 - 10^5 #1 + 1 &, 1]

FindInstance[x^2 - 10^5 x + 1 == 0, x, Reals, 2]

FindRoot[x^2 - 10^5 x + 1 == 0, {x, 0}]

gives back:

 

 

=={{header|MATLAB}} / {{header|Octave}}==

 

=={{header|Maxima}}==

 

/* 2 2

bfloat(%);

/* [x = 1.0000000000000000009999920675269450501b-9,

 

=={{header|МК-61/52}}==

ИП2 / - Вx <-> КвКор НОП x>=0 28 ИП3

x<0 24 <-> /-/ + / Вx С/П /-/ КвКор

 

''Input:'' a С/П b С/П c С/П

=={{header|Modula-3}}==

{{trans|Ada}}

 

IMPORT IO, Fmt, Math;

BEGIN

IF b < 0.0D0 THEN

RETURN Roots{x, c / (a * x)};

ELSE

RETURN Roots{c / (a * x), x};

END;

r := Solve(1.0D0, -10.0D5, 1.0D0);

IO.Put("X1 = " & Fmt.LongReal(r[1]) & " X2 = " & Fmt.LongReal(r[2]) & "\n");

 

=={{header|OCaml}}==

 

| RealRoots of float * float

| ComplexRoots of Complex.t * Complex.t ;;

in

RealRoots (r, c /. (r *. a))

 

{{out}}

- : quadroots = RealRoots (-1.4142135623730951, 1.4142135623730949)

 

 

# quadsolve 1.0 (-1.0e5) 1.0 ;;

 

=={{header|Octave}}==

=={{header|PARI/GP}}==

{{works with|PARI/GP|2.8.0+}}

 

{{trans|C}}

Otherwise, coding directly:

b /= a;

c /= a;

[sol,c/sol]

)

 

Either way,

gives one root around 0.000000001000000000000000001 and one root around 999999999.999999999.

 

=={{header|Pascal}}==

some parts translated from Modula2

 

var

end;

end.

{{out}}

<pre>

=={{header|Perl}}==

When using [http://perldoc.perl.org/Math/Complex.html Math::Complex] perl automatically convert numbers when necessary.

 

($x1,$x2) = solveQuad(1,2,3);

my $root = sqrt($b**2 - 4*$a*$c);

return ( -$b + $root )/(2*$a), ( -$b - $root )/(2*$a);

 

=={{header|Phix}}==

{{trans|ERRE}}

<pre>

for a=1,b=-1e+9,c=1 the real roots are 1e+9 and 1e-9

 

=={{header|PicoLisp}}==

 

(de solveQuad (A B C)

(mapcar round

(solveQuad 1.0 -1000000.0 1.0)

{{out}}

<pre>-> ("999,999.999999" "0.000001")</pre>

 

=={{header|PL/I}}==

declare (c1, c2) float complex,

(a, b, c, x1, x2) float;

put data (x1, x2);

end;

 

=={{header|Python}}==

{{libheader|NumPy}}

This solution compares the naïve method with three "better" methods.

 

import math

print('\nNumpy:')

for c in testcases:

{{out}}

<pre>

 

=={{header|R}}==

} else {

}

}

 

 

=={{header|Racket}}==

(define (quadratic a b c)

(let* ((-b (- b))

;'(0.99999999999995 -1.00000000000005)

;(quadratic 1 0.0000000000001 1)

 

=={{header|REXX}}==

<br>(from a default of &nbsp; '''9''') &nbsp; to &nbsp; '''200''' &nbsp; to accommodate when a root is closer to zero than the other root.

 

 

This REXX version supports &nbsp; ''complex numbers'' &nbsp; for the result.

<pre>

a = 1

root2 = 0.000001

</pre>

 

<pre>

a = 1

c = 1

 

The following output is when R4 REXX is used.

 

<pre>

a = 1

root2 = 0.0000000001

</pre>

<pre>

a = 3

c = 1

 

</pre>

<pre>

a = 1

 

=== Version 2 ===

* 26.07.2913 Walter Pachl

**********************************************************************/

Numeric Digits prec

Return (r+0)

{{out}}

<pre>

 

=={{header|Ring}}==

x1 = 0

x2 = 0

x2 = (-b - sqrtDiscriminant) / (2.0*a)

return [x1, x2]

 

=={{header|Ruby}}==

{{works with|Ruby|1.9.3+}}

The CMath#sqrt method will return a Complex instance if necessary.

 

def quadratic(a, b, c)

p quadratic(-2, 7, 15) # [-3/2, 5]

p quadratic(1, -2, 1) # [1]

{{out}}

<pre>

 

=={{header|Run BASIC}}==

print "FOR 4,5,6 => ";quad$(4,5,6)

 

quad$ = str$(x/(2*a)+sqr(d)/(2*a))+" , "+str$(x/(2*a)-sqr(d)/(2*a))

end if

FOR 4,5,6 => -0.625 +i1.05326872 , -0.625 -i1.05326872</pre>

 

=={{header|Scala}}==

Using [[Arithmetic/Complex#Scala|Complex]] class from task Arithmetic/Complex.

object QuadraticRoots {

def solve(a:Double, b:Double, c:Double)={

}

}

Usage:

(1.0, 22.0, -1323.0), // two distinct real roots

(6.0, -23.0, 20.0), // with a != 1.0

if(roots._1 != roots._2) println("x2="+roots._2)

println

{{out}}

<pre>a=1.00000 b=22.0000 c=-1323.00

 

=={{header|Scheme}}==

(if (= a 0)

(if (= b 0) 'fail (- (/ c b)))

 

(quadratic 1 -1e5 1)

 

=={{header|Seed7}}==

{{trans|Ada}}

include "float.s7i";

include "math.s7i";

r := solve(1.0, -10.0E5, 1.0);

writeln("X1 = " <& r.x1 digits 6 <& " X2 = " <& r.x2 digits 6);

 

{{out}}

 

=={{header|Sidef}}==

[1, 2, 1],

[1, 2, 3],

say ("Roots for #{coefficients}",

"=> (#{quadroots(coefficients...).join(', ')})")

{{out}}

<pre>

 

=={{header|Stata}}==

: polyroots((-2,0,1))

1 2

+-------------------------------+

1 | -1.41421356i 1.41421356i |

=={{header|Tcl}}==

{{tcllib|math::complexnumbers}}

namespace import math::complexnumbers::complex math::complexnumbers::tostring

 

report_quad -2 7 15 ;# {5, -3/2}

report_quad 1 -2 1 ;# {1}

{{out}}

<pre>3x**2 + 4x + 1.3333333333333333 = 0

 

TI-89 BASIC has built-in numeric and algebraic solvers.

returns

<pre>x=1.E-9 or x=1.E9</pre>

 

=={{header|zkl}}==

zkl doesn't have a complex number package.

{{trans|Elixir}}

println("Roots of a quadratic function %s, %s, %s".fmt(a,b,c));

d,a2:=(b*b - 4*a*c), a+a;

println(" the complex roots are %s and \U00B1;%si".fmt(-b/a2,sd/a2));

}

quadratic(a,b,c)

{{out}}

<pre>