/* * Copyright (c) 1993 David I. Bell * Permission is granted to use, distribute, or modify this source, * provided that this copyright notice remains intact. * * 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; } } } global lib_debug; if (!isnum(lib_debug) || lib_debug>0) print "solve(low, high, epsilon) defined";