- sum += term;
- i++;
- } while (i < num_perms);
-
- return sum / num_perms;
-}
-
-ex symmetrize(const ex & e, exvector::const_iterator first, exvector::const_iterator last)
-{
- return symm(e, first, last, false);
-}
-
-ex antisymmetrize(const ex & e, exvector::const_iterator first, exvector::const_iterator last)
-{
- return symm(e, first, last, true);
-}
-
-ex symmetrize(const ex & e, const lst & l)
-{
- exvector v;
- v.reserve(l.nops());
- for (unsigned i=0; i<l.nops(); i++)
- v.push_back(l.op(i));
- return symm(e, v.begin(), v.end(), false);
+ if ((fx[side]<0 && fx[!side]<0) || (fx[side]>0 && fx[!side]>0)) {
+ // Oops, the root isn't bracketed any more.
+ // Restore, and perform a bisection!
+ xx[side] = xxprev;
+ fx[side] = fxprev;
+
+ // Ah, the bisection! Bisections converge linearly. Unfortunately,
+ // they occur pretty often when Newton-Raphson arrives at an x too
+ // close to the result on one side of the interval and
+ // f(x-f(x)/f'(x)) turns out to have the same sign as f(x) due to
+ // precision errors! Recall that this function does not have a
+ // precision goal as one of its arguments but instead relies on
+ // x converging to a fixed point. We speed up the (safe but slow)
+ // bisection method by mixing in a dash of the (unsafer but faster)
+ // secant method: Instead of splitting the interval at the
+ // arithmetic mean (bisection), we split it nearer to the root as
+ // determined by the secant between the values xx[0] and xx[1].
+ // Don't set the secant_weight to one because that could disturb
+ // the convergence in some corner cases!
+ static const double secant_weight = 0.96875; // == 31/32 < 1
+ numeric xxmid = (1-secant_weight)*0.5*(xx[0]+xx[1])
+ + secant_weight*(xx[0]+fx[0]*(xx[0]-xx[1])/(fx[1]-fx[0]));
+ numeric fxmid = ex_to<numeric>(f.subs(x==xxmid).evalf());
+ if (fxmid.is_zero()) {
+ // Luck strikes...
+ return xxmid;
+ }
+ if ((fxmid<0 && fx[side]>0) || (fxmid>0 && fx[side]<0)) {
+ side = !side;
+ }
+ xxprev = xx[side];
+ fxprev = fx[side];
+ xx[side] = xxmid;
+ fx[side] = fxmid;
+ }
+ } while (xxprev!=xx[side]);
+ return xxprev;