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>

=={{header|11l}}==

V sqd = Complex(b^2 - 4*a*c) ^ 0.5

R ((-b + sqd) / (2 * a),

V (r1, r2) = quad_roots(a, b, c)

print(r1, end' ‘ ’)

{{out}}

=={{header|Ada}}==

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

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:

=={{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|FreeBASIC}}==

{{libheader|GMP}}

' compile with: fbc -s console

Data 1, -1e100, 1

Data 1, -1e200, 1

{{out}}

<pre>1: is the naieve way

=={{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

(1, -10e5, 1),

(1, -10e9, 1)

{{Out}}

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

Works in both languages.

solve(1.0, -10.0e5, 1.0)

end

else [r1 := (-b-d)/(2.0*a), c/(a*r1)]

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

{{out}}

=={{header|IDL}}==

print, "input a, press enter, input b, press enter, input c, press enter"

e= y/z

=={{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}}==

'''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}}

=={{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)

=={{header|Mathematica}}/{{header|Wolfram Language}}==

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;

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

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

=={{header|Nim}}==

const Epsilon = 1e-15

let roots = quadRoots(a, b, c)

let plural = if roots.kind == solDouble: "" else: "s"

{{out}}

=={{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}}

<span style="color: #008080;">procedure</span> <span style="color: #000000;">solve_quadratic</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">t3</span><span style="color: #0000FF;">)</span>

<span style="color: #004080;">atom</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">a</span><span style="color: #0000FF;">,</span><span style="color: #000000;">b</span><span style="color: #0000FF;">,</span><span style="color: #000000;">c</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">t3</span><span style="color: #0000FF;">,</span>

<span style="color: #7060A8;">papply</span><span style="color: #0000FF;">(</span><span style="color: #000000;">tests</span><span style="color: #0000FF;">,</span><span style="color: #000000;">solve_quadratic</span><span style="color: #0000FF;">)</span>

<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}}==

r <- sqrt(b * b - 4 * a * c + 0i)

if (abs(b - r) > abs(b + r)) {

qroots(1, -1e9, 1)

Using the builtin '''polyroot''' function (note the order of coefficients is reversed):

[1] -1+1i 1-1i

polyroot(c(1, -1e9, 1))

=={{header|Racket}}==

(define (quadratic a b c)

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

;'(0.99999999999995 -1.00000000000005)

;(quadratic 1 0.0000000000001 1)

=={{header|Raku}}==

Raku has complex number handling built in.

[1, 2, 1],

[1, 2, 3],

given

($b ** 2 - 4 * $a * $c ).Complex.sqrt.narrow

{{out}}

<pre>Roots for 1, 2, 1 => (-1, -1)

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

parse arg a b c . /*obtain the specified arguments: A B C*/

call quad a,b,c /*solve quadratic function via the sub.*/

g=g*.5'e'_%2; do j=0 while h>9; m.j=h; h=h%2+1; end /*j*/

do k=j+5 to 0 by -1; numeric digits m.k; g=(g+x/g)*.5; end /*k*/

{{out|output|text=&nbsp; when using the input of: &nbsp; &nbsp; <tt> 1 &nbsp; -10e5 &nbsp; 1 </tt>}}

<pre>

=== Version 2 ===

* 26.07.2913 Walter Pachl

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

Return (r+0)

exit: Say arg(1)

{{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>

{{trans|Go}}

{{libheader|Wren-complex}}

var quadratic = Fn.new { |a, b, c|

]

{{out}}

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>