/* * Group of functions similar in spirit to Brent's algorithm, but with a * more laissez-faire approach to convergence checking in the short term * and with better use of function evaluation history. * * These functions manipulate lists of x,y points. These lists are implemented * as two-dimensional arrays, with the first index selecting the point, and the * second index being "1" for X and "2" for Y. Note that (unlike C) arrays are * "1"-based rather than "0"-based. All point lists are always sorted by * increasing X value. */ /* * Return the index of the point following the best zero-crossing in the list * of points. Return "error" if there is no zero-crossing. The "best" * crossing is the one with the tightest bounds on the X-axis. * * Note that this function assumes that none of the points is exactly * zero Y-value. After all, if a point has a zero Y-value, why didn't you * just accept that point as the root and be done with it??? */ bent_find_zero_cross(b):= block( [i,zcdx:inf,zci:0], for i:2 step 1 thru length(b) do ( if b[i][2]*b[i-1][2] < 0 and b[i][1] - b[i-1][1] < zcdx then ( zci:i, zcdx:b[i][1] - b[i-1][1], 0), 0), if zci = 0 then ( print("bent_find_zero_cross: no zero crossing: ", b), return(error), 0), return(zci), 0)$ /* * Locate monotonic region covering at most four points surrounding the * zero-crossing in point-list "a" preceding the point indexed by "zci". * (In other words, "zci" might have been returned by bent_find_zero_cross()). */ bent_locate_monotonic(b, zci):= block( [i,highi,mult:1,lenb:length(b),lowi], if b[zci][2] < 0 then mult:-1, i:zci-2, lowi:do ( if i < 1 then return(i+1), if b[i][2]*mult >= b[i+1][2]*mult then return(i+1), i:i-1, 0), i:zci, highi:do ( if i = lenb then return(i), if b[i][2]*mult >= b[i+1][2]*mult then return(i), i:i+1, 0), if highi-lowi >= 4 then ( if zci = highi then lowi:lowi+1 else if zci = 2 then highi:highi-1 else if b[highi-1][1]-b[lowi][1] > b[highi][1]-b[lowi+1][1] then ( lowi:lowi+1, 0) else ( highi:highi-1, 0), 0), c:[], for i:lowi step 1 thru highi do ( c:append(c,[b[i]]), 0), return (c), 0)$ /* * Incorporate the specified point "p" into the list "a". Returns the * resulting list, or error if the point is a duplicate either in X or Y * value. */ bent_incorp_pt(a, p):= block( [b,i], for i:1 step 1 thru length(a) do ( if a[i][1] = p[1] then ( error("Warning: bent_incorp_pt: dup. point"), 0), 0), b:sort(append(a,[p])), return (b), 0)$ /* * Perform inverse quadratic interpolation on the specified triple of points. * Return the X-value of the interpolated point, or error out if the point * could not be computed within the stated precision. */ bent_inv_quad_interp(a3):= block( [a:a3[1][1],b:a3[2][1],c:a3[3][1], fa:a3[1][2],fb:a3[2][2],fc:a3[3][2], p,q,r,s,t], r:fb/fc, s:fb/fa, t:fa/fc, p:s*(t*(r-t)*(c-b) - (1.b0-r)*(b-a)), q:(t-1.b0)*(r-1.b0)*(s-1.b0), return (b+p/q), 0)$ /* * Perform bisection on the specified points. Return the interpolated X value, * or error out if the point could not be computed within the stated precision. */ bent_bisect(a, b):=(a[1]+b[1])/2.b0$ /* * Perform either inverse quadratic interpolation or bisection to the specified * triple of points in a1, which must be a subset of the points in a. * The approach is to try the interpolation, and if it fails or gives a value * outside of the smallest current bracketed zero crossing, do bisection * instead. */ bent_3pt(a1, a, zci):= block([xml], errormsg:false, xml:errcatch(bent_inv_quad_interp(a1)), errormsg:true, if xml = [] then ( xm:bent_bisect(a[zci-1], a[zci]), 0) else ( xm: xml[1], if xm <= a[zci-1][1] or xm >= a[zci][1] then ( xm:bent_bisect(a[zci-1], a[zci]), 0), 0), return(xm), 0)$ /* * Locate a root of the specified function given two points x1 and x2 * bracketing a root. The points may be either simple X-values or they may * be an [x,y] list. The latter is useful in cases where function evaluation * is expensive. The "tol" specifies a maximum tolerance for the X-values * bracketing the root. */ bent(bent_func_name, x1, x2, tol):= block( [a, a1, a2, bent_cr:bent_converge_ratio^(bent_hist_size - 1), do_bisect:false, dxhist, dxhistidx, dxhistidxnext:1, iter, len, xm, ym, zci], num_bent_func_evals:0, if numberp(x1) then ( a:[[x1,bent_func_name(x1)]], num_bent_func_evals:num_bent_func_evals + 1, 0) else if length(x1) = 2 then ( a:[x1], 0) else ( error("Warning: bent: x1 bad format"), 0), if numberp(x2) then ( a:append(a, [[x2,bent_func_name(x2)]]), num_bent_func_evals:num_bent_func_evals + 1, 0) else if length(x2) = 2 then ( a:append(a, [x2]), 0) else ( error("Warning: bent: x2 bad format"), 0), if a[1][2] * a[2][2] >= 0 then ( if zerop(a[1][2]) then return(a[1]), if zerop(a[2][2]) then return(a[2]), error("Warning: bent: x1 and x2 do not bracket zero crossing"), 0), return (for iter:1 step 1 thru bent_itmax do ( /* if (bent_brent_verbose) then print(iter, a), */ zci:bent_find_zero_cross(a), /* if (bent_brent_verbose) then print(iter, "zci = ",zci), */ if abs(a[zci][1] - a[zci-1][1]) < tol then ( if abs(a[zci][2]) < abs(a[zci-1][2]) then ( return(a[zci]), 0) else ( return(a[zci-1]), 0), 0), a1:bent_locate_monotonic(a, zci), /* if (bent_brent_verbose) then print (iter, "a1 = ", a1), */ dxhistidx:dxhistidxnext, dxhistidxnext:dxhistidxnext + 1, if dxhistidxnext > bent_hist_size then dxhistidxnext:1, dxhist[dxhistidx]:abs(a[zci][1] - a[zci-1][1]), if iter > bent_hist_size then ( if dxhist[dxhistidxnext] / dxhist[dxhistidx] < bent_cr then ( do_bisect:true, 0) else ( do_bisect:false, 0), 0), len:length(a1), if len = 2 or do_bisect then ( xm:bent_bisect(a[zci - 1],a[zci]), 0) else if len = 3 then ( xm:bent_3pt(a1, a, zci), 0) else if len = 4 then ( if a1[3][1] - a1[1][1] < a1[4][1] - a1[2][1] then ( a2:rest(a1, -1), 0) else ( a2:rest(a1, 1), 0), xm:bent_3pt(a2, a, zci), 0) else ( error("Warning: bent: list grew too long!"), 0), /* if (bent_brent_verbose) then print (iter, "a2 = ", a2), */ ym:bent_func_name(xm), num_bent_func_evals:num_bent_func_evals + 1, if zerop(ym) then return([xm, ym]), a:bent_incorp_pt(a1, [xm, ym]), 0)), 0)$ bent_converge_ratio:1.9b0$ bent_hist_size:8$ bent_itmax:100$ bent_brent_verbose:false$