/* * solve - solve f(x) = 0 to within the desired error value for x * * Copyright (C) 1999  David I. Bell * * Calc is open software; you can redistribute it and/or modify it under * the terms of the version 2.1 of the GNU Lesser General Public License * as published by the Free Software Foundation. * * Calc is distributed in the hope that it will be useful, but WITHOUT * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY * or FITNESS FOR A PARTICULAR PURPOSE.	 See the GNU Lesser General * Public License for more details. * * A copy of version 2.1 of the GNU Lesser General Public License is * distributed with calc under the filename COPYING-LGPL.  You should have * received a copy with calc; if not, write to Free Software Foundation, Inc. * 59 Temple Place, Suite 330, Boston, MA  02111-1307, USA. * * @(#) $Revision: 29.2 $ * @(#) $Id: solve.cal,v 29.2 2000/06/07 14:02:25 chongo Exp $ * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/solve.cal,v $ * * Under source code control:	1990/02/15 01:50:37 * File existed as early as:	before 1990 * * Share and enjoy!  :-)	http://www.isthe.com/chongo/tech/comp/calc/ *//* * Solve the equation f(x) = 0 to within the desired error value for x. * The function 'f' must be defined outside of this routine, and the low * and high values are guesses which must produce values with opposite signs. */define solve(low, high, epsilon){	local flow, fhigh, fmid, mid, places;	if (isnull(epsilon))		epsilon = epsilon();	if (epsilon <= 0)		quit "Non-positive epsilon value";	places = highbit(1 + int(1/epsilon)) + 1;	flow = f(low);	if (abs(flow) < epsilon)		return low;	fhigh = f(high);	if (abs(flow) < epsilon)		return high;	if (sgn(flow) == sgn(fhigh))		quit "Non-opposite signs";	while (1) {		mid = bround(high - fhigh * (high - low) / (fhigh - flow), places);		if ((mid == low) || (mid == high))			places++;		fmid = f(mid);		if (abs(fmid) < epsilon)			return mid;		if (sgn(fmid) == sgn(flow)) {			low = mid;			flow = fmid;		} else {			high = mid;			fhigh = fmid;		}	}}