[m-dev.] for review: math__solve_quadratic/3

[Mark Anthony BROWN]

Hi,

The well-known general formula for solving quadratics is prone to
numerical error when there are roots close to zero.  There is a better
formula, but I can never remember it offhand, so I thought it might be
a useful addition to the library.

Cheers,
Mark.

Estimated hours taken: 1

Add a function to the standard library to solve quadratic equations.

library/math.m:
	Export the new function, as well as the type that it returns.

NEWS:
	Mention the new function.

tests/hard_coded/Mmakefile:
tests/hard_coded/solve_quadratic.m:
tests/hard_coded/solve_quadratic.exp:
	Test case for the new feature.

Index: NEWS
===================================================================
RCS file: /home/mercury1/repository/mercury/NEWS,v
retrieving revision 1.176
diff -u -r1.176 NEWS
--- NEWS	2000/10/03 02:21:37	1.176
+++ NEWS	2000/10/31 04:24:02
@@ -124,6 +124,8 @@
 	  truncated to the specified precision, but for ISO C's
 	  printf(), they are rounded rather than being truncated.
 
+* We've added a new function, math__solve_quadratic/3.
+
 Changes to the Mercury implementation:
 
 * We've implemented a new back-end for the Mercury compiler.
Index: library/math.m
===================================================================
RCS file: /home/mercury1/repository/mercury/library/math.m,v
retrieving revision 1.24
diff -u -r1.24 math.m
--- library/math.m	1999/12/21 10:33:58	1.24
+++ library/math.m	2000/11/01 03:59:11
@@ -77,7 +77,7 @@
 :- mode math__truncate(in) = out is det.
 
 %---------------------------------------------------------------------------%
-% Power/logarithm operations
+% Polynomial roots
 
 	% math__sqrt(X) = Sqrt is true if Sqrt is the positive square
 	% root of X.
@@ -86,6 +86,21 @@
 :- func math__sqrt(float) = float.
 :- mode math__sqrt(in) = out is det.
 
+:- type math__quadratic_roots
+	--->	no_roots
+	;	one_root(float)
+	;	two_roots(float, float).
+
+	% math__solve_quadratic(A, B, C) = Roots is true if Roots are
+	% the solutions to the equation Ax^2 + Bx + C.
+	%
+	% Domain restriction: A \= 0
+:- func math__solve_quadratic(float, float, float) = quadratic_roots.
+:- mode math__solve_quadratic(in, in, in) = out is det.
+
+%---------------------------------------------------------------------------%
+% Power/logarithm operations
+
 	% math__pow(X, Y) = Res is true if Res is X raised to the
 	% power of Y.
 	%
@@ -187,8 +202,9 @@
 %---------------------------------------------------------------------------%
 
 :- implementation.
+:- import_module float.
 
-% These operations are all implemented using the C interface.
+% These operations are mostly implemented using the C interface.
 
 :- pragma c_header_code("
 
@@ -305,6 +321,49 @@
 #endif
 	SquareRoot = sqrt(X);
 ").
+
+%
+% math__solve_quadratic(A, B, C) = Roots is true if Roots are
+% the solutions to the equation Ax^2 + Bx + C.
+%
+% Domain restrictions:
+% 		A \= 0
+%
+math__solve_quadratic(A, B, C) = Roots :-
+	%
+	% This implementation is designed to minimise numerical errors;
+	% it is adapted from "Numerical recipes in C".
+	%
+	DSquared = B * B - 4.0 * A * C,
+	compare(CmpD, DSquared, 0.0),
+	(
+		CmpD = (<),
+		Roots = no_roots
+	;
+		CmpD = (=),
+		Root = -0.5 * B / A,
+		Roots = one_root(Root)
+	;
+		CmpD = (>),
+		D = sqrt(DSquared),
+		compare(CmpB, B, 0.0),
+		(
+			CmpB = (<),
+			Q = -0.5 * (B - D),
+			Root1 = Q / A,
+			Root2 = C / Q
+		;
+			CmpB = (=),
+			Root1 = -0.5 * D / A,
+			Root2 = -Root1
+		;
+			CmpB = (>),
+			Q = -0.5 * (B + D),
+			Root1 = Q / A,
+			Root2 = C / Q
+		),
+		Roots = two_roots(Root1, Root2)
+	).
 
 %
 % math__pow(X, Y) = Res is true if Res is X raised to the
Index: tests/hard_coded/Mmakefile
===================================================================
RCS file: /home/mercury1/repository/tests/hard_coded/Mmakefile,v
retrieving revision 1.97
diff -u -r1.97 Mmakefile
--- tests/hard_coded/Mmakefile	2000/10/31 02:01:19	1.97
+++ tests/hard_coded/Mmakefile	2000/11/01 04:34:10
@@ -100,6 +100,7 @@
 	rtc_bug \
 	rtti_strings \
 	shift_test \
+	solve_quadratic \
 	space \
 	string_alignment \
 	string_alignment_bug \
Index: tests/hard_coded/solve_quadratic.exp
===================================================================
RCS file: solve_quadratic.exp
diff -N solve_quadratic.exp
--- /dev/null	Tue Jul 25 14:12:01 2000
+++ solve_quadratic.exp	Wed Nov  1 15:24:59 2000
@@ -0,0 +1,16 @@
+Two roots: ok, ok.
+Two roots: ok, ok.
+Two roots: ok, ok.
+Two roots: ok, ok.
+Two roots: ok, ok.
+Two roots: ok, ok.
+One root: ok.
+One root: ok.
+One root: ok.
+No roots.
+No roots.
+No roots.
+Two roots: ok, ok.
+Two roots: ok, ok.
+One root: ok.
+Two roots: ok, ok.
Index: tests/hard_coded/solve_quadratic.m
===================================================================
RCS file: solve_quadratic.m
diff -N solve_quadratic.m
--- /dev/null	Tue Jul 25 14:12:01 2000
+++ solve_quadratic.m	Wed Nov  1 15:31:15 2000
@@ -0,0 +1,75 @@
+:- module solve_quadratic.
+:- interface.
+:- import_module io.
+
+:- pred main(io__state::di, io__state::uo) is det.
+
+:- implementation.
+:- import_module float, math.
+
+main -->
+	% Two roots (i.e. B^2 > 4AC)
+	check_quad(1.0, 3.0, -2.0),
+	check_quad(1.0, -5.0, 3.0),
+	check_quad(2.0, -5.0, 3.0),
+	check_quad(2.0, -5.0, -6.0),
+	check_quad(-15.0, 68.0, 34.5),
+	check_quad(-4.0, -17.0, 3.5),
+
+	% One root (i.e. B^2 = 4AC)
+	check_quad(5.0, 10.0, 5.0),
+	check_quad(4.0, -8.0, 4.0),
+	check_quad(-3.0, -18.0, -27.0),
+
+	% No roots (i.e. B^2 < 4AC)
+	check_quad(4.0, 3.0, 2.0),
+	check_quad(1.0, 1.0, 1.0),
+	check_quad(-1.0, -2.0, -2.0),
+
+	% Special cases
+	check_quad(1.0, -2.0, 0.0),
+	check_quad(1.0, 0.0, -2.0),
+	check_quad(2.0, 0.0, 0.0),
+	check_quad(-100.0, 0.0, 0.0001).
+
+
+:- pred check_quad(float, float, float, io__state, io__state).
+:- mode check_quad(in, in, in, di, uo) is det.
+
+check_quad(A, B, C) -->
+	{ Ss = solve_quadratic(A, B, C) },
+	(
+		{ Ss = no_roots },
+		io__write_string("No roots.\n")
+	;
+		{ Ss = one_root(R) },
+		io__write_string("One root: "),
+		check_result(A, B, C, R),
+		io__write_string(".\n")
+	;
+		{ Ss = two_roots(R1, R2) },
+		io__write_string("Two roots: "),
+		check_result(A, B, C, R1),
+		io__write_string(", "),
+		check_result(A, B, C, R2),
+		io__write_string(".\n")
+	).
+
+:- pred check_result(float, float, float, float, io__state, io__state).
+:- mode check_result(in, in, in, in, di, uo) is det.
+
+check_result(A, B, C, R) -->
+	{ Val = ((A * R + B) * R) + C },
+	(
+		%
+		% This test is pretty conservative, since I don't know
+		% how much error should be expected.
+		%
+		{ abs(Val) < 1E-6 }
+	->
+		io__write_string("ok")
+	;
+		io__write_string("problem: Val = "),
+		io__write_float(Val)
+	).
+
--------------------------------------------------------------------------
mercury-developers mailing list
Post messages to:       mercury-developers at cs.mu.oz.au
Administrative Queries: owner-mercury-developers at cs.mu.oz.au
Subscriptions:          mercury-developers-request at cs.mu.oz.au
--------------------------------------------------------------------------